clear
clc
format short g
cd T129b
results=0;

% FLIGHT CONDITIONS

speed=58;           %True airspeed [ft/s]
aoa=8;              %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 1 1];         %Number of semispanwise partitions for this wing geometry
chord=[20/12 8/12 8/12];            %Root chord length [ft]
NACA_r=[5400;0;0012];        %Base chord airfoil NR (4 digits): NACA
NACA_o=[5400;0;0012];        %Outboard airfoil NR (4 digits): NACA
chord_pan=[5;3;3];        %Number of panels chord wise
bs_twist=[0;0;0];        %Base chord twist [deg]
out_twist=[0;0;0];        %Outboard twist [deg]
dihedral=[0;180;90];         %Partition dihedral [deg]
span_pan=[10;4;4];         %Number of panels semi-span wise
tpr=linspace(.6,1,7);
taper=[tpr;ones(1,length(tpr))*.6;ones(1,length(tpr))*.6];  %Taper ratio
m_spn=S./(chord(1,1).*(1+taper(1,:)));
span=[m_spn;m_spn*.4;ones(1,length(tpr))];             %Span of partition [ft]
sweep=[0;0;0];            %Quarter chord line sweep [deg]
flap=[0;0;0] ;             %Is partition flapped [1 0]
flap_chord=[0;0;0];     %Flap chord in fraction of local chord (0..1)
flap_chord_pan=[0;0;0];   %Number of chord wise panels on flap
flap_sym=[0;0;0];         %Does control surfaces deflect symmetrically [1 0]
main_wing=23.5/12; h_tail=81.66/12; v_tail=67.9/12; %Location of the wings (x-direction)
wing_loc=[main_wing h_tail v_tail]-main_wing; %Location of the wings relative to the main wing

% 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

for wing=1:length(taper(1,:))
    geometry(wing).nwing= length(chord);
    geometry(wing).nelem= span_par;
    geometry(wing).ref_point= [-chord(1)/4 0 0]*0.3048;
    geometry(wing).symetric= [1 1 0];
    geometry(wing).startx= wing_loc*0.3048;
    geometry(wing).starty= [0 0 0]*0.3048;
    geometry(wing).startz= [0 -5/12 0]*0.3048;
    geometry(wing).c= chord*0.3048;
    geometry(wing).foil(:,1,1)= NACA_r;
    geometry(wing).foil(:,1,2)= NACA_o;
    geometry(wing).nx= chord_pan;
    geometry(wing).TW(:,:,1)= bs_twist*pi/180;
    geometry(wing).TW(:,:,2)= out_twist*pi/180;
    geometry(wing).dihed= dihedral*pi/180;
    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));
    Run_ID(wing).name=['wing' num2str(wing)];
    [latt_w(wing),ref_w(wing)]=fLattice_setup(geometry(wing),flight_con);
   
    for comp=[1 8]
    solverloop5(results,comp,Run_ID(wing).name,latt_w(wing),flight_con,geometry(wing),ref_w(wing));
    end
end

% 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

% 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)
    figure(wing)
    clf
    set(gcf,'Position',[10 10 1250 750])
    
    subplot(3,2,1)
    plot(latt_w(wing).XYZ(:,:,2)'/0.3048,latt_w(wing).XYZ(:,:,1)'/0.3048,'k')
    axis equal,hold on
    xlabel('Wing x-coordinate')
    ylabel('Wing y-coordinate')
    zlabel('Wing z-coordinate')
    title('3-D Wing configuration')
    axis tight
    grid on
    hold off
    
  cd output
    fname=strcat(Run_ID(wing).name,'-Cx');
    load(fname)
    cd ..
    
    subplot(3,2,3)
    hold on, grid on
    plot(results.ystation(:,1)/0.3048,results.ForcePerMeter(:,1)/4.44822162*0.3048);
    Title('Spanload on main wing');
    ylabel('lbf per foot');
    xlabel('Spanstation')
    hold off
    
    subplot(3,2,5)
    hold on, grid on
    plot(results.ystation(:,1)/0.3048,results.CL_local(:,1),'r');
    Title('Spanload on main wing');
    ylabel('CL Local');
    xlabel('Spanstation')
    hold off
    
    cd output
    fname=strcat(Run_ID(wing).name,'-Cxx_diff');
    load(fname)
    cd ..     
    
    subplot(3,2,[2 4 6])
    axis off
    text(0,1,'TORNADO CALCULATION RESULTS, Central difference')
text(0,.95,'JID: '); text(0.25,0.95,Run_ID(wing).name)
text(0,.90,'Reference area: ');	text(0.25,0.90,num2str(ref_w(wing).S_ref/0.3048^2 ));	
text(0,.85,'Root chord: ');text(0.25,.85,num2str(geometry(wing).c/0.3048));
text(0,.8,'Reference span: ');text(0.25,.8,num2str(ref_w(wing).b_ref/0.3048));
  
   
text(0.4,.90,'Alpha: '); 	text(.55,.9,num2str(flight_con.alpha*180/pi));
text(0,.75,'Taper: ');text(0.55,.95,num2str(geometry(wing).T));
text(0.4,.85,'Beta: '); 		text(.55,.85,num2str(flight_con.betha*180/pi));
text(0.4,.8,'Airspeed: ');text(.55,.8,num2str(flight_con.AS/0.3048));

text(0.65,.9,'P: ');  		text(.8,.9,num2str(flight_con.P));
text(0.65,.85,'Q: '); 		text(.8,.85,num2str(flight_con.Q));
text(0.65,.8,'R: ');  		text(.8,.8,num2str(flight_con.R));
% 
% subplot(2,2,4)
% axis off
text(0,.7,'CL derivatives : ');			

text(0,.65,'CL-alfa');		text(0.15,.65,num2str(iszero(results.CL_a)));
text(0,.6,'CL-beta');		text(0.15,.6,num2str(iszero(results.CL_b)));
text(0,.55,'CL-P');			text(0.15,.55,num2str(iszero(results.CL_P)));
text(0,.5,'CL-Q');			text(0.15,.5,num2str(iszero(results.CL_Q)));
text(0,.45,'CL-R');			text(0.15,.45,num2str(iszero(results.CL_R)));
   
text(0,.35,'Roll derivatives : ');			
   
text(0,.3,'Cl-alfa');		text(0.15,.3,num2str(iszero(results.Cl_a)));
text(0,.25,'Cl-beta');		text(0.15,.25,num2str(iszero(results.Cl_b)));
text(0,.2,'Cl-P');			text(0.15,.2,num2str(iszero(results.Cl_P)));
text(0,.15,'Cl-Q');			text(0.15,.15,num2str(iszero(results.Cl_Q)));
text(0,.1,'Cl-R');			text(0.15,.1,num2str(iszero(results.Cl_R)));
   
text(0.35,.7,'CD derivatives : ');			

text(0.35,.65,'CD-alfa');	text(0.5,.65,num2str(iszero(results.CD_a)));
text(0.35,.6,'CD-beta');	text(0.5,.6,num2str(iszero(results.CD_b)));
text(0.35,.55,'CD-P');		text(0.5,.55,num2str(iszero(results.CD_P)));
text(0.35,.5,'CD-Q');		text(0.5,.5,num2str(iszero(results.CD_Q)));
text(0.35,.45,'CD-R');		text(0.5,.45,num2str(iszero(results.CD_R)));
   
text(0.35,.35,'Pitch derivatives : ');			

text(0.35,.3,'Cm-alfa');	text(.5,.3,num2str(iszero(results.Cm_a)));
text(0.35,.25,'Cm-beta');	text(0.5,.25,num2str(iszero(results.Cm_b)));
text(0.35,.2,'Cm-P');		text(0.5,.2,num2str(iszero(results.Cm_P)));
text(0.35,.15,'Cm-Q');		text(0.5,.15,num2str(iszero(results.Cm_Q)));
text(0.35,.1,'Cm-R');		text(0.5,.1,num2str(iszero(results.Cm_R)));
   
text(0.7,.7,'CY derivatives : ');			

text(0.7,.65,'CY-alfa');	text(0.85,.65,num2str(iszero(results.CY_a)));
text(0.7,.6,'CY-beta');		text(0.85,.6,num2str(iszero(results.CY_b)));
text(0.7,.55,'CY-P');		text(0.85,.55,num2str(iszero(results.CY_P)));
text(0.7,.5,'CY-Q');			text(0.85,.5,num2str(iszero(results.CY_Q)));
text(0.7,.45,'CY-R');		text(0.85,.45,num2str(iszero(results.CY_R)));
   
text(0.7,.35,'Yaw derivatives : ');			

text(0.7,.3,'Cn-alfa');	text(.85,.3,num2str(iszero(results.Cn_a)));
text(0.7,.25,'Cn-beta');	text(0.85,.25,num2str(iszero(results.Cn_b)));
text(0.7,.2,'Cn-P');			text(0.85,.2,num2str(iszero(results.Cn_P)));
text(0.7,.15,'Cn-Q');		text(0.85,.15,num2str(iszero(results.Cn_Q)));
text(0.7,.1,'Cn-R');			text(0.85,.1,num2str(iszero(results.Cn_R)));

end
cd ..