Graph not Plotting Correctly??

Hi,

When I run my code I expect to plot a graph which decreases close to 0 as the x axis approaches 10. I ave used the debugger which suggests that the code is correct and giving appropriate values. However when I plot the graph, the trendline becomes erratic at an x value of about 9. I have attached my code for clarity. Does anybody know how I could fix this as Matlab is clearly not plotting correctly at these values?

% Inputs
R=0.4;                                          % Radius of Rotor
B=3;                                            % Number of blades
U=1;                                            % Fluid velocity
Rho=998;                                        % Fluid Density
N=9;                                            % Number of Blade Elements
Cp_estimate=0.5;                                % Estimate power coefficient
Alpha_design=4;                                 % Design alpha
Cl_design=[1;1;1;1;1;1;1;1;1];                  % Design lift coefficient

% Variables
Q=1;
jj=1;                                           % Counter
tolerance=0.1;                                  % Tolerance Value
axial_induction=[0;0;0;0;0;0;0;0;0];
Check=1;                                        % Initial check value
axial_induction_old=[0;0;0;0;0;0;0;0;0];        % Initial value for old axial induction factor
relative_wind=[1;1;1;1;1;1;1;1;1]; 

Cl=[1.3; 1.1; 1; 0.9; 0.86; 0.83; 0.8; 0.75; 0.5];                                          % Lift Coefficients
Cd=[0.027; 0.024; 0.02; 0.019; 0.018; 0.016; 0.013; 0.012; 0.01];                           % Drag Coefficients

r_local=R/9*(1:9)';
r_over_R=r_local / R;

ii=0;                                                   % Counting Interval for integral values

for TSR=8:0.1:10                                      % TSR from 1 to 10
    disp(TSR)                                                                    
    Check=1;
ii=ii+1;

TSR_local=r_over_R .* TSR;    
relative_wind=(2/3)*atan(1./TSR_local);
C=((8*pi.*r_local) ./ (B.*Cl_design)).*(1-cos(relative_wind));
sigma=(B*C) ./ (pi.*R.*2);
sigma_local=(B*C) ./ (pi.*r_local.*2); 

axial_induction= 1 ./ (((4.*(sin(relative_wind).^2)) ./ (sigma_local.*Cl_design.*cos(relative_wind)))+1);
angular_induction= (1-(3.*axial_induction)) ./ ((4.*axial_induction)-1);
angular_velocity=22;

TSR_local = TSR .* (r_local./R);                                                                 % Local Tip Speed Ratio
relative_wind = (2/3) .* atan(1./TSR_local);                                                               % Angle of Relative Fluid

     while abs(Check)>=tolerance

axial_induction_old = axial_induction;        

relative_wind = atan(U.*(1-axial_induction))./(angular_velocity*R).*(1+angular_induction);

F=zeros(size(relative_wind));
for Q=1:length(relative_wind)
    if relative_wind(Q)<0;
    F(Q)=1;
    else
    F(Q)=(2/pi) .* acos(exp(-(((B/2) .* (1-(r_over_R(Q)))) ./ ((r_over_R(Q)) .* sin(relative_wind(Q))))));      % Tip Loss Factor
    end;
end;

alpha=relative_wind-S;

if axial_induction<0.4;
C_T=(sigma_local .* ((1-axial_induction).^2) .* ((Cl.*cos(relative_wind))+(Cd.*sin(relative_wind)))) ./ ((sin(relative_wind)).^2); 
else
    C_T=(8/9)+(4.*F-(40/9)).*axial_induction+((50/9)-4.*F).*axial_induction.^2;
end

for jj=1:length(C_T)
  angular_induction(jj) = 1 ./ ((4.*F(jj).*cos(relative_wind(jj))./(sigma_local(jj).*Cl(jj)))-1);

     if  C_T(jj) < 0.96;
        axial_induction(jj) = 1 / 1+(4*F(jj)*(sin(relative_wind(jj))^2)/(sigma_local(jj)*Cl(jj)*cos(relative_wind(jj))));
     else
         axial_induction(jj) = (1/F(jj)) .* abs(0.142+(0.0203-0.6427.*(0.889-C_T(jj))).^0.5);
     end;
  end;

D=abs((8./(TSR.*9)).*(F.*(sin(relative_wind).^2).*(cos(relative_wind)-((TSR_local).*(sin(relative_wind)))).*(sin(relative_wind)+((TSR_local).*(cos(relative_wind)))).*(1-(Cd./Cl).*cot(relative_wind)).*(TSR_local.^2)));
Cp=sum(D(1:8));

Diff=axial_induction-axial_induction_old;
Check=max(Diff(:));
Check

     end

store_C(:,ii)=C;
store_sigma_local(:,ii)=sigma_local;
store_sigma(:,ii)=sigma;
store_Phi(:,ii)=relative_wind;
store_TSR_local(:,ii)=TSR_local;
store_axial_induction(:,ii)=axial_induction;
store_angular_induction(:,ii)=angular_induction;
store_axial_induction_old(:,ii)=axial_induction_old;
store_relative_wind(:,ii)=relative_wind;
store_Check(:,ii)=Check;
store_Diff(:,ii)=Diff;
store_Cp(:,ii)=Cp;
store_TSR(:,ii)=TSR;
store_F(:,ii)=F;

end

figure(1)
plot(store_TSR,store_Cp)
axis([0 10 0 0.6])
hold all
title('Cp vs Tip Speed Ratio')
xlabel('TSR')
ylabel('Cp')

[Kevin]

[Image Analyst]

I don't know - I can't run it:

Undefined function or variable 'S'.
Error in test3 (line 61)
alpha=relative_wind-S;

What is S? Please give us code that actually runs without errors if we copy and paste.

Answer 1

[Iain]

It is very unlikely that matlab is plotting incorrectly.

It is more likely that there is something going funny in the calculations.

Get matlab to print out TSR and CP when TSR is close to 9. See what you spot.

[Kevin]

Thanks for your comment. I have checked the code using the debugger and focusing around TSR = 8 and 9. The values appear to be correct but when they are plotted the trend line goes crazy.