close all
clear all

global k_global;

 [rho S Cp C ] = Constants();
 

%SpringRate = 0.06;  %N-m/deg/wing
SpringRate = 0.09;  %N-m/deg/wing


%aoa = 0;
%[cl cd cm]=Coef0015(aoa);

spine_pts = [-38.55       35.45         35.45          -38.55        -38.55 ;
    11.28        11.28         5.28           5.28         11.28]; %inches

wing_ctr = 0;%[78.45:-14:-76];
wing_ctr(2,:) = 8.28*ones(1,length(wing_ctr(1,:)));

spine_pts = .0254*spine_pts;  %meters
wing_ctr = 0.0254*wing_ctr;  %meters

Chord = .0254*12.34;
ctr = .0254*2.875;


theta = [-180:1:180];

Vel = .2:.2:.6;

fh1 = figure('Position',[    1000          84        1168        1414]);


for k = 1:length(Vel)
    V = Vel(k);
   
        k_global = 0.0;
        alpha = NaN(size(theta));
        F = NaN(3,length(alpha));
        
        leg_str{k} = ['Vel = ' num2str(Vel(k)) ' m/s'];
        
        v_w(1,:) = -V*cos(theta*pi/180);
        v_w(2,:) = -V*sin(theta*pi/180);
        

        qinf = .5*rho*V.^2;

        F0 = S*qinf*C;
        M0 = S*qinf*C*C;
        
        idx = find(theta < -120);
        alpha(idx) = -30;
        for j = idx
            F(:,j) = HorzForce(alpha(j),theta(j));
        end
        
        idx = find(theta > 120);
        alpha(idx) = 30;
        for j = idx
            F(:,j) = HorzForce(alpha(j),theta(j));
        end
        
        idx = find((theta <= 120) & (theta >= -120) )
        
        for j = idx
          [alpha(j) F(:,j)] = OptimumWingAngle(theta(j),V);
        end
             
        
        AOA = alpha-theta;
        
        figure(fh1)
        ah(1) = subplot(6,1,1);
        plot(theta,alpha);
        hold all
        
        ah(2) = subplot(6,1,2);
        plot(theta,AOA)
        hold all
        
        ah(3) = subplot(6,1,3);
        plot(theta,F(1,:))
        hold all
        
        ah(4) = subplot(6,1,4);
        plot(theta,F(2,:))
        hold all
        
        ah(5) = subplot(6,1,5);
        plot(theta,F(3,:))
        hold all
        
        ah(6) = subplot(6,1,6);
        plot(theta,F(3,:)/M0)
        hold all
        
end
    



subplot(6,1,1);
grid on;
xlabel('Theta (water angle)');
ylabel('Alpha (wing angle)');
legend(leg_str);
title_str =  ['Optimum Wing Angle:'];
title(title_str);

subplot(6,1,2);
grid on;
xlabel('Theta (water angle)');
ylabel('Wing AOA (degrees)');
legend(leg_str);

ah(3) = subplot(6,1,3);
grid on;
xlabel('Theta (water angle)');
ylabel('Horizontal Force (N)');
legend(leg_str);

ah(4) = subplot(6,1,4);
grid on;
xlabel('Theta (water angle)');
ylabel('Vertical Force (N)');
legend(leg_str);



ah(5) = subplot(6,1,5);
title_str =  ['Optimum Wing Moment:  C/4-Cp = ' num2str((C/4-Cp)/.0254) ' in    SpringRate = ' num2str(SpringRate) ' N-m/deg/wing'];
title(title_str);
grid on;
xlabel('Theta (water angle)');
ylabel('Moment (N-m)');
legend(leg_str);

plot(theta,SpringRate*theta);


ah(6) = subplot(6,1,6);
title_str =  ['Optimum Wing Moment Coeff:  C/4-Cp = ' num2str((C/4-Cp)/.0254) ' in    SpringRate = ' num2str(SpringRate) ' N-m/deg/wing'];
title(title_str);
grid on;
xlabel('Theta (water angle)');
ylabel('Moment Coeff M/qinf*C^2*V^2');
legend(leg_str);


linkaxes(ah,'x');

filenm = ['OptimumAngle' num2str((C/4-Cp)/.0254,2) '.ps'];
    print('-fillpage','-dpsc',filenm);


%
%
% figure
% plot(alpha(idx),F(3,(idx)))
% xlabel('alpha(wing angle to spine)');
% ylabel('Required Moment (N-m)');
%
%
% fh = figure('Position',[       100         481        2336        1017]);
% ah1 = axes('OuterPosition',[0 0 1 1]);
%
% set(fh, 'renderer', 'zbuffer');
% set(fh,'DoubleBuffer','on');
% set(ah1,'nextplot','replace');
% set(ah1,'FontSize',18);
%
% hold(ah1,'on');
% axis(ah1,'equal');
% h_spine = fill(spine_pts(1,:),spine_pts(2,:),[.9 .9 .9]);
% h_wingctr = plot(wing_ctr(1),wing_ctr(2),'ro');
% h_Vw = quiver(wing_ctr(1),wing_ctr(2),v_w(1,1),v_w(2,1),'b');
% [pts] = WingPoints([wing_ctr(1) wing_ctr(2)],-alpha(1)*pi/180,Chord,ctr);
% h_wings = plot(pts(1,:),pts(2,:),'k');
%
%
% xlim([-1 1]);
% ylim([-.4 .4]);
% xlabel('m');
% ylabel('m');
%
% for k = 2:length(theta)
%     disp(['theta = ' num2str(theta(k))]);
%     set(h_Vw,'Udata',v_w(1,k));
%     set(h_Vw,'Vdata',v_w(2,k));
%     [pts] = WingPoints([wing_ctr(1) wing_ctr(2)],-alpha(k)*pi/180,Chord,ctr);
%     set(h_wings,'Xdata',pts(1,:));
%     set(h_wings,'Ydata',pts(2,:));
%     pause(.005);
%     drawnow;
% end
%
%
% %end
%
