clear
clc
results=0;
% FLIGHT CONDITIONS

speed=50;           %True airspeed [ft/s]
aoa=6.682;              %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 = 17.7778;           %Area of wing Planform
span_par=[1 1 1];         %Number of semispanwise partitions for this wing geometry
chord=[20 10.224 13.064]./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=[8;5;5];        %Number of panels chord wise
bs_twist=[0;-.9;0];        %Base chord twist [deg]
out_twist=[0;-.9;0];        %Outboard twist [deg]
dihedral=[.8;180;90];         %Partition dihedral [deg]
span_pan=[10;4;4];         %Number of panels semi-span wise
tpr=0.6;
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*.419;m_spn*.2279];             %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=88.444/12; v_tail=81.484/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

cd T129b

% 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