cd T129b

% FLIGHT CONDITIONS

speed=57.075;           %True airspeed [ft/s]
aoa=8: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;         %Number of semispanwise partitions for this wing geometry
chord=20/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=linspace(.4,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

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));
    Run_ID(wing).name=['wing' num2str(wing)];
    stall_CL=0;swp=1;
    computation = 1;
    results=0;
    while stall_CL<1.6 && 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 ..
        stall_CL=iszero(results.CL)
        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)
    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 ..
    
    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
    x_span=results.ystation(:,1)/0.3048;
    y_force=results.ForcePerMeter(:,1)/4.44822162*0.3048;
    plot(x_span,y_force);
    axis([min(min(xwing)) max(max(xwing)) 0 max(y_force)*1.05]);
    Title('Spanload on main wing');
    ylabel('lbf per foot');
    xlabel('Spanstation')
    hold off
    
    subplot(3,2,5)
    hold on, grid on
    y_CL=results.CL_local(:,1);
    plot(x_span,y_CL,'r');
    axis([min(min(xwing)) max(max(xwing)) min(y_CL) max(y_CL)*1.05]);
    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 ..