Make a 4D contour plot.

I'm trying to visualize the efficiency map (ITEn) of an engine from the sweep of three (3) inputs. See code below.The results in my workspace is good. However, plotting this using isosurface has me in doubt. The color on the plot, likewise the colorbar and their corresponding ITEn do not match the ITEn in the workspace. For example, the indicated point in the figure shows ITEn between 37% and 40%, but actual simulation results in 34.65%.

How can I plot this properly?

IVC_sweep = round(linspace(-180,-121,50));
EVO_sweep = round(linspace(105,165,50));
FM_sweep = round(linspace(70,200,50));

for i = 1:length(IVC_sweep)
    for ii = 1:length(FM_sweep)
        for iii = 1:length(EVO_sweep)
        IVC = IVC_sweep(i);
        EVO = EVO_sweep(iii);
        IVO = -11;
        EGR = 30;
        Qf  = FM_sweep(ii); %200;
        Tin = 313;
        pin = 200;
        RPM = 1200*2*pi/60;
        SOI = -40.3;
        u_in = [IVC IVO EVO EGR Qf Tin pin RPM SOI];
        [CA10(i,ii,iii), CA50(i,ii,iii), IMEPg(i,ii,iii), ITEn(i,ii,iii),...
            Pmax(i,ii,iii), MPRR(i,ii,iii), TEVO(i,ii,iii)]  = MX13_3(u_in);
        end
    end
end
figure;
levellist = linspace(30,50,10); % Range of ITEn
for k = 1:length(levellist)
level = levellist(k);
p = patch(isosurface(IVC_sweep,FM_sweep,EVO_sweep,ITEn,level));
p.FaceVertexCData = level;
p.FaceColor = 'flat';
p.EdgeColor = 'none';
p.FaceAlpha = 0.3;

xlabel('IVC [{\circ}CA aTDC]');
ylabel('m_{fuel} [mg/cyc]');
zlabel('EVO [{\circ}CA aTDC]');
title('ITEn [%]')
grid on

end
view(3)
hc=colorbar;
rotate3d on

[Walter Roberson]

Unfortunately MX13_3 is not defined, so I cannot test my ammended version.

[Prince Essemiah]

Hey Walter, MX13_3 is a function that consists of a collection of executables so I cannot just put the code here. Is it possible to send it to you directly? Or I can test your version instead.

[Walter Roberson]

I threw together an example function, and the code seems to work.

Nothing about your calculations have changed. The plotting was made slightly more efficient by moving some code.

The important change for this code is that it creates a custom data cursor to display the coordinates including the ITEn level associated with the point. It takes the current point X Y Z, and finds the closest point to that in the volume, and displays the coordinates it found together with the ITEn stored at the corresponding index. So it might perhaps not have "fixed" anything, but it will at least permit you to see the exact ITEn values at the location being pointed at, and you can cross-reference that to the colorbar for reasonability check.

IVC_sweep = round(linspace(-180,-121,50));
EVO_sweep = round(linspace(105,165,50));
FM_sweep = round(linspace(70,200,50));

for i = 1:length(IVC_sweep)
    for ii = 1:length(FM_sweep)
        for iii = 1:length(EVO_sweep)
            IVC = IVC_sweep(i);
            EVO = EVO_sweep(iii);
            IVO = -11;
            EGR = 30;
            Qf  = FM_sweep(ii); %200;
            Tin = 313;
            pin = 200;
            RPM = 1200*2*pi/60;
            SOI = -40.3;
            u_in = [IVC IVO EVO EGR Qf Tin pin RPM SOI];
            [CA10(i,ii,iii), CA50(i,ii,iii), IMEPg(i,ii,iii), ITEn(i,ii,iii),...
                Pmax(i,ii,iii), MPRR(i,ii,iii), TEVO(i,ii,iii)]  = MX13_3(u_in);
        end
    end
end
fig = figure;
ax = gca(fig);

levellist = linspace(30,50,10); % Range of ITEn
for k = 1:length(levellist)
    level = levellist(k);
    p = patch(ax, isosurface(IVC_sweep,FM_sweep,EVO_sweep,ITEn,level));
    p.FaceVertexCData = level;
    p.FaceColor = 'flat';
    p.EdgeColor = 'none';
    p.FaceAlpha = 0.3;
end

xlabel(ax, 'IVC [{\circ}CA aTDC]');
ylabel(ax, 'm_{fuel} [mg/cyc]');
zlabel(ax, 'EVO [{\circ}CA aTDC]');
title(ax, 'ITEn [%]')
grid(ax, 'on');
view(ax, 3);
hc = colorbar(ax);
rotate3d(ax, 'on');

dcm = datacursormode(fig);
dcm.UpdateFcn = @(obj, info) displayinfo(obj, info, IVC_sweep,FM_sweep,EVO_sweep, ITEn);
dcm.Enable = 'on';

function txt = displayinfo(~, info, IVC_sweep, FM_sweep, EVO_sweep, ITEn);
   x = info.Position(1);
   y = info.Position(2);
   z = info.Position(3);
   d = (IVC_sweep - x).^2 + (FM_sweep - y).^2 + (EVO_sweep - z).^2;
   [~, idx] = min(d);
   truex = IVC_sweep(idx);
   truey = FM_sweep(idx);
   truez = EVO_sweep(idx);
   truew = ITEn(idx);
   txt = ["IVC = " + truex, "FM = " + truey, "EVO = " + truez, "ITEn = " + truew];
end

[Prince Essemiah]

Thank you for the data cursor code. I had to do some modifcations because not only the swept inputs define the ITEn, though that's how the isosurface interprets it. So it'd require patching the surface with all the input parameters in the isosurface command, which with my matlab knowledge is impossible to achieve.

[Walter Roberson]

isosurface() converts all of the input coordinates into appropriate faces and vertex structure. All of the information that you pass to isosurface() is present in the structure that is returned to be used by patch()

isosurface() does not care how the values are calculated. It would be happy, for example, with a 3D projection down from a 4 dimensional space, such as a Maximum Intensity Projection. It does, however, assume that values are continuous, that if it sees a pixel with value 1 and a pixel beside it with value 2, assumes that there must be a location between the two with value between the two.