close all
clear all

global k_global;

 [rho S Cp C ] = Constants();
 

%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:1;
Vel = 0.4;


SpringConst = [.001 .025:.05:.2 .4];

SpringConst = 0.05;


for k = 1:length(Vel)
    V = Vel(k);
    qinf = .5*rho*V.^2;
    M0 = S*qinf*C^2;
    F0 = S*qinf*C;
    fh1 = figure('Position',[    1000          84        1168        1414]);
    title_str = ['Vel = ' num2str(V) 'm/s'   '    C/4-Cp = ' num2str((C/4-Cp)/.0254) ' in'];
    for i = 1:length(SpringConst)
        
        k_global = SpringConst(i);
        alpha = NaN(size(theta));
        F = NaN(3,length(alpha));
        
        leg_str{i} = ['k = ' num2str(k_global) ' N-m/deg/wing'];
        
        
        v_w(1,:) = -V*cos(theta*pi/180);
        v_w(2,:) = -V*sin(theta*pi/180);
        
        %for j = 1:length(theta);
        %  [alpha(j) F(:,j)] = OptimumWingAngle(theta(j),V);
        %end
        
        idx = find(theta < -120);
        alpha(idx) = -30;
        for j = idx
            F(:,j) = HorzForce(alpha(j),theta(j));
            F(1:2,j) = F0*F(1:2,j);
        end
        
        idx = find(theta > 120);
        alpha(idx) = 30;
        for j = idx
            F(:,j) = HorzForce(alpha(j),theta(j));  %These are force coefficients!
            F(1:2,j) = F0*F(1:2,j);
        end
        
        idx = find((theta <= 120) & (theta >= -120) )
        
        
        for j = idx
            [alpha(j) F(:,j)] = ActualWingAngle(theta(j),V);
        end
        
        AOA = alpha-theta;
        
        figure(fh1)
        ah(1) = subplot(4,1,1);
        plot(theta,alpha);
        hold all
        
        ah(2) = subplot(4,1,2);
        plot(theta,AOA)
        hold all
        
        ah(3) = subplot(4,1,3);
        plot(theta,F(1,:))
        hold all
        
        ah(4) = subplot(4,1,4);
        plot(theta,F(2,:))
        hold all
        
        
    end


subplot(4,1,1);
grid on;
xlabel('Theta (water angle)');
ylabel('Alpha (wing angle)');
legend(leg_str);
title(title_str);

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

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

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


linkaxes(ah,'x');
filenm = ['Pivot' num2str((C/4-Cp)/.0254,2) '.ps'];
if(k==1)
    print('-fillpage','-dpsc',filenm);
else
    print('-fillpage','-dpsc',filenm,'-append');
end

end

%
%
% 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
%
