Plotting problem

Hi. I've written this code to plot the graph it produces. However, the top three lines have complex parts to them. How can I plot only the real parts? Code follows:

%Matlab code for "MATLAB 
%Question 1
%
%13/02/12

function [lmda, Ct] = problem

%QUESTION 1%

rs=0.05; %rotor solidity
Cl=2*pi; %specific lift coefficient

pa1=0.1745329250; %pitch angle 1
pa2=0.0872664625; %pitch angle 2 
pa3=-0.0017453292; %pitch angle 3
pa4=-0.0349065850; %pitch angle 4
pa5=-0.0872664625; %pitch angle 5

%--------------------%
i=1; 
    for lmda=0:0.1:20;        
        a2(i)=(((((rs*lmda*Cl)/16)+0.5)^2)-((rs*lmda*Cl*(1-(lmda*pa1)))/8));            
            if a2(i) < 0                
                break            
            end
i=i+1;
   end
%--------------------%
i=1;
    for lmda=(0:0.1:20); 
        a1(i)=(((rs*lmda*Cl)/16)+0.5);  

i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);
        at(i)=a1(i)-sqrt(a2(i));      
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);      
        Ct1(i)=((rs*(lmda^2)*Cl)/2)*(((1-at(i))/lmda)-pa1);
i=i+1;
    end
%----------------------------------------------------------------%
i=1; 
    for lmda=0:0.1:20;        
        a2(i)=(((((rs*lmda*Cl)/16)+0.5)^2)-((rs*lmda*Cl*(1-(lmda*pa2)))/8));
            if a2(i) < 0
              break
            end
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);
        a1(i)=(((rs*lmda*Cl)/16)+0.5);
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);
        at(i)=a1(i)-sqrt(a2(i));      
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);        
        Ct2(i)=((rs*(lmda^2)*Cl)/2)*(((1-at(i))/lmda)-pa2);
i=i+1;
    end
%----------------------------------------------------------------%
i=1; 
    for lmda=0:0.1:20;        
        a2(i)=(((((rs*lmda*Cl)/16)+0.5)^2)-((rs*lmda*Cl*(1-(lmda*pa3)))/8));            
            if a2(i) < 0                
                break            
            end           
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);    
        a1(i)=(((rs*lmda*Cl)/16)+0.5);         
i=i+1;
    end
%--------------------%    
i=1;
    for lmda=(0:0.1:20);    
        at(i)=a1(i)-sqrt(a2(i));        
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);        
        Ct3(i)=((rs*(lmda^2)*Cl)/2)*(((1-at(i))/lmda)-pa3);
i=i+1;
    end
%----------------------------------------------------------------%
i=1; 
    for lmda=0:0.1:20;        
        a2(i)=(((((rs*lmda*Cl)/16)+0.5)^2)-((rs*lmda*Cl*(1-(lmda*pa4)))/8));            
            if a2(i) < 0                
                break
            end
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);    
        a1(i)=(((rs*lmda*Cl)/16)+0.5);  
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);    
        at(i)=a1(i)-sqrt(a2(i));        
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);
        Ct4(i)=((rs*(lmda^2)*Cl)/2)*(((1-at(i))/lmda)-pa4);
i=i+1;
    end

%----------------------------------------------------------------%
i=1; 
    for lmda=0:0.1:20;        
        a2(i)=(((((rs*lmda*Cl)/16)+0.5)^2)-((rs*lmda*Cl*(1-(lmda*pa5)))/8));            
            if a2(i) < 0                
                break            
            end            
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);    
        a1(i)=(((rs*lmda*Cl)/16)+0.5);          
i=i+1;
    end
%--------------------%    
i=1;
    for lmda=(0:0.1:20);    
        at(i)=a1(i)-sqrt(a2(i));        
i=i+1;
    end
%--------------------%
i=1;
    for lmda=(0:0.1:20);
        Ct5(i)=((rs*(lmda^2)*Cl)/2)*(((1-at(i))/lmda)-pa5);
i=i+1;
    end
%----------------------------------------------------------------%
lmda=(0:0.1:20);
figure,plot(lmda,real(Ct1),'k-','Linewidth',2),grid on, hold on 
plot(lmda,real(Ct2),'k--','Linewidth',2)
plot(lmda,real(Ct3),'k-.','Linewidth',2)
plot(lmda,real(Ct4),'r-','Linewidth',2)
plot(lmda,real(Ct5),'r--','Linewidth',2)
axis([0 20 0 1])
xlabel('\lambda');ylabel('C_t');title('Thrust coefficient against tip speed ratio.');
legend('\theta_T=10^{o}','\theta_T=5^{o}','\theta_T=-0.1^{o}','\theta_T=-2^{o}','\theta_T=-5^{o}');
end

%QUESTION 2%

Answer 1

[Matt Tearle]

You are plotting only the real parts. I think the problem is coming with how you're calculating these functions. You have break statements to terminate your for loops when a2 goes negative. However, the last value of a2 (the negative value) is still there.

Furthermore, you're reusing a2, but you stop calculating after a2 goes negative. This means that the rest of a2 is still there from the previous calculation. Do you really want that?

If you zoom out, so the y-axis goes higher, you will see that the function carries on doing something else. My suspicion is that this comes from calculating at based on old a2 values.

Finally, you really don't need these loops. This code is really nasty to read and debug.

lmda=0:0.1:20;
a1=(((rs*lmda*Cl)/16)+0.5);

Nice and easy. Just remember to use .* and ./ when multiplying/dividing arrays.

[J]

Great, thank you very much! I'm new to Matlab so this style was really helpful.

Answer 2

[J]

Ok, I've tried to edit the code and get rid of loops.

rs = 0.05;
Cl = 2*pi;
pa = [10*(pi/180) 5*(pi/180) -0.1*(pi/180) -2*(pi/180) -5*(pi/180)];
lmda=1:0.1:20;    

for v = 1:5
    for i=1:length(lmda)
        a(i) = (((rs*lmda(i)*(Cl))/16)+0.5)-sqrt((((rs*lmda(i)*(Cl))/16)+0.5).^2)-((rs*lmda(i)*(Cl).*(1-lmda(i).*pa(v)))/8);
        Ct(i) = ((rs*lmda(i).^2*Cl)/2).*(((1-a(i))./1-pa(v)));
    end
    Z = imag(Ct);
        for d = 1:1:201
            if Z(v,d) ~= 0
                Ct(v,d) = NaN;
            end
        end
end

plot (lmda,Ct)
legend('10','5','-0.1','-2','-5')
axis([0 20 0 1])
xlabel('\lambda');ylabel('C_t');title('Thrust coefficient against tip speed ratio.');
legend('\theta_T=10^{o}','\theta_T=5^{o}','\theta_T=-0.1^{o}','\theta_T=-2^{o}','\theta_T=-5^{o}')

Answer 3

[J]

I changed it again, thanks so much for your help!

clear all
clc

rs=0.05;

Cl=2*pi;

pa=[(pi/18) (pi/36) -(pi/1800) -(pi/90) -(pi/36)];

for v=1:5

      lmda=0:0.1:20;

      Ct(v,:) = ((rs*lmda.^2*Cl)/2).*(((1-((rs*lmda*Cl/16+1/2)-...
                ((rs*lmda*Cl/16+1/2).^2-(rs*lmda*Cl.*(1-lmda*pa(v)))/8).^(1/2))))./lmda-pa(v));

      Z= imag(Ct);

          for d = 1:1:201

              if Z(v,d) ~= 0

                  Ct(v,d) = NaN;

              end

          end

end

figure,plot(lmda,Ct)
grid on
axis([0 20 0 1]);
xlabel('\lambda');ylabel('C_t');title('Thrust coefficient against tip speed ratio.');
legend('\theta_T=10^{o}','\theta_T=5^{o}','\theta_T=-0.1^{o}','\theta_T=-2^{o}','\theta_T=-5^{o}');

[Matt Tearle]

And that's why I love MATLAB. Actually, one more thing I love about MATLAB that would help you right now: logical indexing.

All that stuff to NaN out anything with a nonzero imaginary part? Yeah, try this: Ct(imag(Ct)~=0) = NaN; (Do it after the loop, so you have to do it only once.)

And speaking of doing stuff outside the loop, move your lmda = ... before the loop -- no need to do the same thing 5 times. And notice the orange warning on the line Ct(v,:) = ...? It won't make a whole heap of difference in this case, but get into the habit of preallocating space for your arrays, when you know how big they will be. Ct = zeros(5,length(lmda)); before the loop.

[J]

<3 haha teach me more! Why couldn't you be my lecturer!

[Matt Tearle]

Because MathWorks pays more :)

Maybe some of these suggestions might help you:

http://www.mathworks.com/matlabcentral/answers/8026-best-way-s-to-master-matlab