clear
tic
cd T129b

% FLIGHT CONDITIONS

speed=57.075;           %True airspeed [ft/s]
aoa=0;              %Alpha [deg]
density=1.2250;     %Desity of the air [kg/m^3]
beta=0;             %Beta [deg]
P=0;                %Roll angular velocity [deg/s]
Q=0;                %Pitch angular velocity [deg/s]
R=0;                %Yaw angular velocity [deg/s]


flight_con.AS=speed*0.3048;                    
flight_con.alpha=aoa*pi/180;    
flight_con.betha=beta*pi/180;                  
flight_con.P=P*pi/180;                      
flight_con.Q=Q*pi/180;                      
flight_con.R=R*pi/180;                      
flight_con.rho=density;

% WING GEOMETRY
S = 13.5;           %Area of wing Planform
span_par=1;         %Number of semispanwise partitions for this wing geometry
chord=24/12;            %Root chord length [ft]
NACA_r=5400;        %Base chord airfoil NR (4 digits): NACA
NACA_o=5400;        %Outboard airfoil NR (4 digits): NACA
chord_pan=6;        %Number of panels chord wise
bs_twist=0;        %Base chord twist [deg]
bs_twist=0;        %Outboard twist [deg]
dihedral=0;         %Partition dihedral [deg]
span_pan=12;         %Number of panels semi-span wise
taper=1 %linspace(.2,1,7);            %Taper ratio
span=2*S./(chord*(1+taper))./2;             %Span of partition [ft]
sweep=0;            %Quarter chord line sweep [deg]
flap=1;             %Is partition flapped [1 0]
flap_chord=0.2;     %Flap chord in fraction of local chord (0..1)
flap_chord_pan=1;   %Number of chord wise panels on flap
flap_sym=1;         %Does control surfaces deflect symmetrically [1 0]

% COPUTATIONS AND RESULTS
% [1]. Simple solution computation. Forces/Coefficients only
%  
% [2]. Alpha sweep computation 
% [3]. Beta sweep computation 
% [4]. Delta sweep computation 
%  
% [5]. Roll rate sweep computation 
% [6]. Pitch rate sweep computation 
% [7]. Yaw rate sweep computation 
%  
% [8]. Central difference expansion around current state
computation=2;
results=0;

for wing=1:length(taper)
    geometry(wing).nwing= 1;
    geometry(wing).nelem= span_par;
    geometry(wing).ref_point= [0 0 0];
    geometry(wing).symetric= 1;
    geometry(wing).startx= 0;
    geometry(wing).starty= 0;
    geometry(wing).startz= 0;
    geometry(wing).c= chord*0.3048;
    geometry(wing).foil(1,1,1)= NACA_r;
    geometry(wing).foil(1,1,2)= NACA_o;
    geometry(wing).nx= chord_pan;
    geometry(wing).TW(1,1,1)= bs_twist*pi/180;
    geometry(wing).TW(1,1,2)= bs_twist*pi/180;
    geometry(wing).dihed= dihedral;
    geometry(wing).ny= span_pan;
    geometry(wing).b= span(wing)*0.3048;
    geometry(wing).T= taper(wing);
    geometry(wing).SW= sweep*pi/180;
    geometry(wing).flapped= flap;
    geometry(wing).fc= flap_chord;
    geometry(wing).fnx= flap_chord_pan;
    geometry(wing).fsym= flap_sym;
    geometry(wing).flap_vector=zeros(size(geometry(wing).flapped));
    [latt_w(wing),ref_w(wing)]=fLattice_setup(geometry(wing),flight_con);
    Run_ID(wing).name=['wing' num2str(wing)];
    solverloop5b(results,computation,Run_ID(wing).name,latt_w(wing),flight_con,geometry(wing),ref_w(wing),-10,.5,15);
end

% geometryplot(latt_w(wing),geometry(wing),ref_w(wing))
% solverloop5(results,computation,Run_ID(wing).name,latt_w(wing),flight_con,geometry(wing),ref_w(wing));
% resultplot2(1,Run_ID(wing).name)

for wing=1:length(taper)
    cd output
    fname=strcat(Run_ID(wing).name,'-Cx_alpha');
    load(fname)
    cd ..
    
    figure(wing)
    clf
%     set(gcf,'Position',[10 10 1250 750])
    
    subplot(3,2,3), plot(results.alpha_sweep*180/pi,results.CD_f_a);
   xlabel('Alpha [deg]')
   ylabel('CD')
   
   subplot(3,2,5), plot(results.alpha_sweep*180/pi,results.CY_f_a);
   xlabel('Alpha [deg]')
   ylabel('CY')
   
   subplot(3,2,2), plot(results.alpha_sweep*180/pi,results.Cl_f_a);
   xlabel('Alpha [deg]')
   ylabel('Cl')
   
   subplot(3,2,4), plot(results.alpha_sweep*180/pi,results.Cm_f_a);
   xlabel('Alpha [deg]')
   ylabel('Cm')
   
   subplot(3,2,6), plot(results.alpha_sweep*180/pi,results.Cn_f_a);
   xlabel('Alpha [deg]')
   ylabel('Cn')
   
   subplot(3,2,1), plot(results.alpha_sweep*180/pi,results.CL_f_a);
   xlabel('Alpha [deg]')
   ylabel('CL')
  
  axes
  text(0.0,1.05,'Coefficient dependency on alpha')
  axis off  
end
cd ..
toc