% speed to fly in wave, and up / downwind. Uses a + bv + cv^2 which is
% better fit as does not force quadratic at min sink. 

clear all
close all

wet = 0; 
dry = ~wet;

% find S(V) = a + bV + c V^2. Input three values

if dry
    V = [ 75 100 125 150 175 ]'; % km/h
    S = [ (21/32)/(1+3/8)  (3/4)/(1+3/8) (1+1/16)/(1+3/8) (1+11/16)/(1+3/8) (2+1/2)/(1+3/8)]'; % m/s
end;

if wet
    V = [ 100 125 150 175 200]'; % km/h
    S = [ (7/8)/(1+3/8)  (15/16)/(1+3/8) (1+7/32)/(1+3/8) (1+11/16)/(1+3/8) (2+1/4)/(1+3/8)]'; % m/s
end;

V = V*0.539957; % to knots
S = S*1.94384; 

x = [ones(size(V,1),1) V V.^2];
abc = x\S;

a = abc(1);
b = abc(2);
c = abc(3);

Vplot = (30:0.01:120)'; 
Splot = a + b*Vplot + c*Vplot.^2;

figure; 
plot(V,-S,'v');
hold on
plot(Vplot,-Splot,'linewidth',2)
xlabel('V, knots'); 
ylabel('S, knots'); 
title('Dry ASG29 polar'); 

figure
plot(V,V./S,'v'); 
hold on
plot(Vplot,Vplot./Splot,'linewidth',2);
xlabel('V, knots');
ylabel('L/D')
title('Dry ASG29 Polar')

% ***************
% achieved speed up/down and x wind
% ******************

Vw = [0 20 40 60]'; 
M = [2 4 6]';
V = (61:1:140)';
Vmat = V*ones(1,size(Vw,1));
Vwmat = ones(size(V,1),1)*Vw';

SofVmat = (a+b*Vmat+c*Vmat.^2); 

% achieved speed graphs
for Mi = 1:size(M,1); 
    Vau = (Vmat-Vwmat)./(1+SofVmat/M(Mi)); % Upwind achieved average speed
    Vax = ((Vmat.^2-Vwmat.^2).^0.5)./(1+SofVmat/M(Mi)); % xwind achieved average speed
    
    % optimum point 
    Vupmax = Vw + (Vw.^2 +(M(Mi)+a+b*Vw)/c).^0.5; 
    Vxmax = 0*Vupmax;
    for in = 1:size(Vw,1);
        Vxmax(in) = fsolve(@(x) c*x^3 - (2*c*Vw(in)^2 + M(Mi) + a)*x - b*Vw(in)^2,100);
    end
    Supmax = a + b*Vupmax + c*Vupmax.^2;
    Sxmax = a + b*Vxmax + c*Vxmax.^2; 
    Vaupmax = (Vupmax-Vw)./(1+Supmax/M(Mi)); % Upwind achieved average speed
    Vaxmax = ((Vxmax.^2-Vw.^2).^0.5)./(1+Sxmax/M(Mi)); % xwind achieved average speed
    
    
    
    
    figure
    plot(V,[Vau],'linewidth',2)
    hold on
    plot(Vupmax,Vaupmax,'v')
    title(['Achieved speeds, still air to ' num2str(M(Mi)) ' knot wave lift -- upwind '])
    legend('0','20kt','40kt','60kt')  
    xlabel('Airspeed, knots');
    ylabel('Achieved speed, knots');
    
    if M(Mi) == 2; 
    if dry; 
        print -depsc2 wave_up_achieved_dry.eps
    end
    if wet; 
        print -depsc2 wave_up_achieved_wet.eps
    end
    end;
    
    if M(Mi) == 6; 
    if dry; 
        print -depsc2 wave_up_achieved_dry_6.eps
    end
    if wet; 
        print -depsc2 wave_up_achieved_wet_6.eps
    end
    end;
    
    
    figure
    plot(V,[Vax],'linewidth',2)
    hold on
    plot(Vxmax,Vaxmax,'v')
    title(['Achieved speeds, still air to ' num2str(M(Mi)) ' knot wave lift -- crosswind '])
    legend('0','20kt','40kt','60kt') 
    xlabel('Airspeed, knots');
    ylabel('Achieved speed, knots');
  
   if dry; 
        print -depsc2 wave_x_achieved_dry.eps
    end
    if wet; 
        print -depsc2 wave_x_achieved_wet.eps
    end
  
    
    
end



 % Now, maximum calculations
 
Vw = [0 20 40 60]'; 
M = [2 4 6]';
Mmat = M*ones(1,size(Vw,1));
Vwmat = ones(size(M,1),1)*Vw';

Vup = Vwmat + (Vwmat.^2 +(Mmat+a+b*Vwmat)/c).^0.5; 
Mup = c*Vup.^2 - a; 

for in = 1:size(Mmat,1);
    for jn = 1:size(Mmat,2);
        Vx(in,jn) = fsolve(@(x) c*x^3 - (2*c*Vwmat(in,jn)^2 + Mmat(in,jn) + a)*x - b*Vwmat(in,jn)^2,130);
    end;
end;
Mx = c*Vx.^2 - a; 

disp(' ');
disp('Upwind/downwind Speed')
fprintf(' %6.8s   %6.8s   %6.8s   %6.8s  \n', ' ', ' ', 'Wind', 'speed'); 
fprintf(' %6.8s ',' Lift ')
fprintf(' %6.0f ',Vw);
fprintf('\n'); 
for ind = 1:size(M,1);
    fprintf(' %6.1f ', M(ind))
    fprintf(' %6.0f ', Vup(ind,:))
    fprintf('\n');
end

disp(' ');
disp('Crosswind Speed')
fprintf(' %6.8s   %6.8s   %6.8s   %6.8s  \n', ' ', ' ', 'Wind', 'speed'); 
fprintf(' %6.8s ',' Lift ')
fprintf(' %6.0f ',Vw);
fprintf('\n'); 
for ind = 1:size(M,1);
    fprintf(' %6.1f ', M(ind))
    fprintf(' %6.0f ', Vx(ind,:))
    fprintf('\n');
end

Vx

disp(' ');
disp('Upwind/downwind Mc setting')
fprintf(' %6.8s   %6.8s   %6.8s   %6.8s  \n', ' ', ' ', 'Wind', 'speed'); 
fprintf(' %6.8s ',' Lift ')
fprintf(' %6.0f ',Vw);
fprintf('\n'); 
for ind = 1:size(M,1);
    fprintf(' %6.1f ', M(ind))
    fprintf(' %6.1f ', Mup(ind,:))
    fprintf('\n');
end


disp(' ');
disp('Crosswind Mc setting')
fprintf(' %6.8s   %6.8s   %6.8s   %6.8s  \n', ' ', ' ', 'Wind', 'speed'); 
fprintf(' %6.8s ',' Lift ')
fprintf(' %6.0f ',Vw);
fprintf('\n'); 
for ind = 1:size(M,1);
    fprintf(' %6.1f ', M(ind))
    fprintf(' %6.1f ', Mx(ind,:))
    fprintf('\n');
end

% effective polar graphs

if 1;
    Vw = [0 50]';
    V = (30:1:130)';
    Vmat = V*ones(1,size(Vw,1));
    Vwmat = ones(size(V,1),1)*Vw';
    Smat = a + b*(Vmat) + c*(Vmat).^2;
    Vupmat = Vmat - Vwmat; 
    Vxmat = abs((Vmat.^2-Vwmat.^2)).^0.5;
    Vxmat(Vmat<Vwmat)=NaN;
    
    % find mc optimums for plot
    Mplot = 2; 
    Vup = Vw + (Vw.^2 +(Mplot+a+b*Vw)/c).^0.5; 
    Vx = 0*Vup;
    for in = 1:size(Vx,1);
        Vx(in) = fsolve(@(x) c*x^3 - (2*c*Vw(in)^2 + Mplot + a)*x - b*Vw(in)^2,130);
    end;
    Sup = a + b*Vup + c*Vup.^2;
    Sx = a + b*Vx + c*Vx.^2; 
    Mup = c*Vup(2).^2 - a; 
    Mx = c*Vx(2).^2 - a; 

    
    figure;
    plot(Vupmat(:,1),-Smat(:,1),'-b','linewidth',2)
    hold on
    plot(Vupmat(:,2),-Smat(:,2),'--r','linewidth',2)
    plot( [ 0 1.5*(Vup(1)-Vw(1))]', [Mplot Mplot-1.5*(Sup(1)+Mplot)]', '-b','linewidth',2);
    plot( [ 0 1.5*(Vup(2)-Vw(2))]', [Mplot Mplot-1.5*(Sup(2)+Mplot)]', '-r','linewidth',2);   
    plot( [ 0 1.5*(Vup(2))]', [Mup Mup-1.5*(Sup(2)+Mup)]', '-r','linewidth',2);   
    
    plot( [ (Vup(1)-Vw(1))]', [-Sup(1)]', 'vb','linewidth',2);
    plot( [ (Vup(2)-Vw(2))]', [-Sup(2)]', 'vr','linewidth',2); 
    plot( [ Vup(2)]', [-Sup(2)]', 'vr','linewidth',2); 
    
    xlabel('Groundspeed V, knots');
    ylabel('Sink rate, knots');
    title('Dry ASG29 polar upwind');
    axis([-10 130 -10 5]);
    legend('0 wind','50 kt')
    plot([0 0]',[-10 5],'-k');
    plot([-10 120]',[0 0],'-k');
   
    if dry; 
        print -depsc2 wave_up_polar.eps
    end
    
    
    figure
    plot(Vxmat(:,1),-Smat(:,1),'-b','linewidth',2)
    hold on
    plot(Vxmat(:,2),-Smat(:,2),'--r','linewidth',2)
    plot( [ 0 1.5*(Vx(1)^2-Vw(1)^2)^0.5]', [Mplot Mplot-1.5*(Sx(1)+Mplot)]', '-b','linewidth',2);
    plot( [ 0 1.5*(Vx(2)^2-Vw(2)^2)^0.5]', [Mplot Mplot-1.5*(Sx(2)+Mplot)]', '-r','linewidth',2); 
    plot( [ 0 1.5*Vx(2)]', [Mx Mx-1.5*(Sx(2)+Mx)]', '-r','linewidth',2); 
    
    plot( [ (Vx(1)^2-Vw(1)^2)^0.5]', [-Sx(1)]', 'vb','linewidth',2);
    plot( [ (Vx(2)^2-Vw(2)^2)^0.5]', [-Sx(2)]', 'vr','linewidth',2);
    plot( [ Vx(2)]', [-Sx(2)]', 'vr','linewidth',2);   
    xlabel('Groundspeed V, knots');
    ylabel('Sink rate, knots');
    title('Dry ASG29 polar cross wind');
    axis([-10 130 -10 5]);
    legend('0 wind','50 kt')
    plot([0 0]',[-10 5],'-k');
    plot([-10 120]',[0 0],'-k');
    
    if dry; 
        print -depsc2 wave_x_polar.eps
    end
    
end;



