Matlab don't show negative values

Question

I am making a matlab file to analyse the performance of an aircraft. What I want to do is see where the excess power is max and where the excess power becomes minimum. I was able to make that matlab doesn't take into consideration values more than the Clmax (by using NaN).
Now I want to do the same for the excess power: I don't want my contour plot to display the negative values of PS. I tried doing the same as I did with the Cl, but without success.

This the code so far:

TA= 667233.2;
W=228000;
S=360;
CDO=0,02;
b=60;
AR=10;
h=0;


v= [60:10:350]
hoogte= [0:1000:17000]

%TEST

for j=1:length(hoogte)
hoogte_temp = hoogte(j)
[rho c] = Atmos(hoogte_temp)
sigma= rho/1.225

for i= 1:length(v)
    v_temp = v(i);
    q= rho*v_temp*v_temp/2
    cl(i,j) = W*9.81/S/q

  if (cl<1.91)
   cd(i,j) = 0.02 + cl(i,j)^2/(AR*3.14*0.8)/beta(i)
    D(i,j)= cd(i,j)*q*S
    PS(i,j)=(TA*sigma-D(i,j))*v_temp/(W*9.81)


  else %als Cl groter is dan 1,91: niet weergeven, ook niet Cd, D en ook niet Ps
    cl(cl>1.91) = NaN;
    cd(i,j)= NaN;
    D(i,j)=NaN;
    PS(i,j)=NaN;
  end

    if(v_temp/c<1)
        beta(i)=sqrt(1-(v_temp/c)^2)

    else
        beta(i) = NaN
    end

    cd(i,j) = 0.02 + cl(i,j)^2/(AR*3.14*0.8)/beta(i)
    D(i,j)= cd(i,j)*q*S
    PS(i,j)=(TA*sigma-D(i,j))*v_temp/(W*9.81)

end

I thought of just adding another if (PS<0) or even if PS<=D PS=NaN, but that does not seem to be working if I rerun the program and look at the graph.

The atmos.m file is a file I made where the atmosphere and its properties (like temperature, altitude and density) vary.
What I want to do now:
Cancel the negative PS values, and find the sigma and v_temp values for which PS becomes 0. This would lead me then to the max cruising altitude.

I'm new to matlab, and everything in this code was done by looking up on the internet and looking youtube movies. Any feedback will greatly be appreciated! I found this:

    data( indices_to_data_not_to_plot )  = NaN;
    surf(x,y,data);  %# or whatever you're using

in another topic, but it's the thing I tried I guess, which does not work, except I use contour instead of surf.

Answer 1

A simple piece of code to clarify what I said in my comment above. This creates a wavy 2D pattern (2D sinc function), then adds contours at different levels - sometimes by forcing the contour level, sometimes by forcing the data.

% contour example
x = linspace(-10,10,20); y=x;
[xx yy]=meshgrid(x,y);
r = sqrt(xx.^2+yy.^2);
z = sin(r)./r;
figure
subplot 121
imagesc(x,y,z); axis xy; axis image
colormap gray
hold all
contour(x,y,z,0:0.1:1,'r')
contour(x,y,z,[0 0],'b') % blue line at "good" zero
zclip = z;
zclip(z<0)=0;
contour(x,y,zclip,[0 0],'g') % green line at "clipped" zero

subplot 122
imagesc(x,y,z); axis xy; axis image
colormap gray
hold all
contour(x,y,z,'w')

This produces the following output:

As you can see, using the vector 0:0.1:1 in the first contour command causes only positive contours to be drawn (in red, from 'r'); using [0 0] as the contour specification in the second call causes a blue ('b') contour to be drawn just at z==0; but when I actually set the negative data to zero, I get the ragged third contour (green). This is because it's difficult to "estimate the zero crossing" when there is no real "crossing"...

On the right you see the result of "conventional" contour plotting, with white contours at "default" spacing.

Note - I plotted the intensity map to show how it could be done, but obviously you can have the contour map plotted without first showing the intensity map...

Please use the above example as a starting point to improving your question (complete code you can run, and that shows the issue and an approach to a solution) - then we'll be able to give you a more relevant answer.

Happy coding.