clear
cd T129b

% FLIGHT CONDITIONS

speed=57.075;           %True airspeed [ft/s]
aoa=0:30;              %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=1%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]

% 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= [0 3 3]*0.3048;
    geometry(wing).starty= [0 0 0]*0.3048;
    geometry(wing).startz= [0 -5/12 -5/12]*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)];
    lift=0;swp=1;
    computation = 1;
    results=0;
%     while lift<55 && swp<=length(aoa) 
        flight_con.alpha=aoa(swp)*pi/180;
        [latt_w(wing),ref_w(wing)]=fLattice_setup(geometry(wing),flight_con);
        solverloop5(results,computation,Run_ID(wing).name,latt_w(wing),flight_con,geometry(wing),ref_w(wing));
%         cd output
%         fname=strcat(Run_ID(wing).name,'-Cx');
%         load(fname)
%         cd ..
%         lift=iszero(results.L/4.44822162)
%         swp=swp+1;
%     end
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(1,:))
    figure(wing)
    clf
    set(gcf,'Position',[0 0 1250 750])
    
    subplot(3,2,1)
    plot3(latt_w(wing).XYZ(:,:,1)'/0.3048,latt_w(wing).XYZ(:,:,2)'/0.3048,latt_w(wing).XYZ(:,:,3)'/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 ..
    
    cdi=iszero(results.CL)^2/(pi*(ref_w(wing).b_ref/0.3048)^2/(ref_w(wing).S_ref/0.3048^2));
    Oswald_e=iszero(results.CD)/cdi;
    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
       
    subplot(3,2,[2 4 6])
    axis off
text(0,1,'Tornado Computation Results ')
text(0,.95,'JID: '); text(0.25,0.95,Run_ID(wing).name)
text(0.4,.95,'Downwash matrix condition: '); text(0.8,0.95,num2str(results.dwcond))

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,.9,'Taper: ');text(0.8,.9,num2str(geometry(wing).T'));
text(0.4,.85,'Reference point pos: ');text(0.8,.85,num2str(geometry(wing).ref_point));
text(0.4,.8,'Oswalds efficiency: ');text(0.8,.8,num2str(Oswald_e));
   
text(0,.7,'Net Wind Forces: (lbf)');			
	text(0.0,.65,'Drag: ');		text(0.1,.65,num2str(iszero(results.D/4.44822162)));
   text(0.0,.6,'Side: ');			text(0.1,.6,num2str(iszero(results.C/4.44822162)));
   text(0.0,.55,'Lift: ');		    text(0.1,.55,num2str(iszero(results.L/4.44822162)));
   
text(0.35,.7,'Net Body Forces: (lbf)');			
   text(0.35,.65,'X: ');	    	text(0.4,.65,num2str(iszero(results.FORCE(1)/4.44822162)));
   text(0.35,.6,'Y: ');			text(0.4,.6,num2str(iszero(results.FORCE(2)/4.44822162)));
   text(0.35,.55,'Z: ');		    text(0.4,.55,num2str(iszero(results.FORCE(3)/4.44822162)));

text(0.7,.7,'Net Body Moments: (lbf-ft)');
   text(0.7,.65,'Roll: ');			text(0.8,.65,num2str(iszero(results.MOMENTS(1)/(4.44822162*0.3048))));
   text(0.7,.6,'Pitch: ');		    text(0.8,.6,num2str(iszero(results.MOMENTS(2)/(4.44822162*0.3048))));   
   text(0.7,.55,'Yaw: ');		    text(0.8,.55,num2str(iszero(results.MOMENTS(3)/(4.44822162*0.3048))));
   
text(0,.4,'CL ');	    text(0.1,.4,num2str(iszero(results.CL)))   
text(0,.35,'CD ');	    text(0.1,.35,num2str(iszero(results.CD))) 
text(0,.3,'CY ');	    text(0.1,.3,num2str(iszero(results.CY)))

text(0.35,.4,'CZ ');    text(0.45,.4,num2str(iszero(results.CZ)))
text(0.35,.35,'CX ');	text(0.45,.35,num2str(iszero(results.CX)))
text(0.35,.3,'CC ');	text(0.45,.3,num2str(iszero(results.CC)))

text(0.7,.4,'Cm ');		text(0.8,.4,num2str(iszero(results.Cm)))
text(0.7,.35,'Cn ');	text(0.8,.35,num2str(iszero(results.Cn)))
text(0.7,.3,'Cl ');		text(0.8,.3,num2str(iszero(results.Cl)))

text(0,.2,'STATE: ');
text(0,.15,'alpha: '); text(.15,.15,num2str(iszero(state.alpha*180/pi)));
text(0,.1,'beta: '); text(.15,.1,num2str(iszero(state.betha*180/pi)));
text(0,.05,'Airspeed: ');text(.15,.05,num2str(iszero(state.AS/0.3048)));
text(0,.0,'Density: ');text(.15,.0,num2str(iszero(state.rho)));

text(0.3,.15,'P: '); text(.45,.15,num2str(iszero(state.P)));
text(0.3,.10,'Q: '); text(.45,.10,num2str(iszero(state.Q)));
text(0.3,.05,'R: '); text(.45,.05,num2str(iszero(state.R)));


text(0.6,.1,'Rudder setting [deg]:'); text(.9,.1,num2str(geometry(wing).flap_vector*180/pi));
end
cd ..