% Analytical trim by balancing lift and weight 
% and setting pitch moment to zero
clear all
format compact
disp(' '); disp('Start Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<')

disp(' Data Input section')
echo on

% Data is selected to duplicate Roskam figures 4.9 <<<<<<<<<<<
disp(' '); disp(' Second we will duplicate figures 4.9, 4.10, and 4.11(right)')
CL0=.17
CLalpha=.1      % Per degree
CLdeltaE=.25/30 % Per degree  <<< Change here
CLiH=3*CLdeltaE % Per degree

CM0=.025
CMalpha=-.015    % Per degree
CMdeltaE=-.0062  %Per degree
CMiH=3*CMdeltaE %Per degree
iH=0             % Degree

mom_ref_pt=.25   % Moment reference point in % chord
forward_cg=.15   % forward cg limit in % chord
aft_cg=.30       % aft cg limit in % chord
alpha_stall=10   % angle of attack (deg) for stall
Cl_PlotMax=1.2   % maximum Cl for plots in figure 2 (the aircraft trim diagram)
alpha_PlotMax=10 % maximum angle of attack (deg) for figure 1

deltaE=[-30 -20 -10 0 10]
% End of data input section


% Data below if for S-211 (cruise 2 modified) page 494-495. Data above is
% unchanged from the data that duplicates figures 4.5 and 4.9.
altitudeFT=35000;  % altitude, ft
W=4000;     % weight, lbf
S=136;      % wing area, ft^2
V=450;      % trim speed, ft/sec
V=284;      % trim speed, ft/sec
% End of data input section
echo off

rho=rhofun(altitudeFT)
CL1=2*W/(rho*V*V*S)

% Plotting information
color=['-bo-gx-r+-c*-md-yv-k^'];
s1=['De=',num2str(deltaE(1)),' deg.'];
s2=['De=',num2str(deltaE(2)),' deg.'];
s3=['De=',num2str(deltaE(3)),' deg.'];
s4=['De=',num2str(deltaE(4)),' deg.'];
s5=['De=',num2str(deltaE(5)),' deg.'];
% End plotting information

% First find CLmax for the different elevator deflections.
CLmax=CL0+CLalpha*alpha_stall+CLiH*iH+CLdeltaE*deltaE
% Now compute stick fixed neutral point (aka the airplane aerodynamic center)
disp('Another name for the airplane aerodynamic center is the')
disp(' stick fixed neutral point.')
XbarNP=-CMalpha/CLalpha+mom_ref_pt
disp('Recall that the neutral point is an aerodynamic property')
disp(' of the aircraft and does not depend on c.g. position.')

% Perform the analytical solution of the two trim equations
% See equations 4.45
disp('Compute a trim solution for the one speed specified above.')
a=[CLalpha CLdeltaE
    CMalpha CMdeltaE];
b=[CL1-CLiH*iH-CL0
    -CM0-CMiH*iH];
x=inv(a)*b;

alpha_trim=x(1)
dE_trim=x(2)
disp('Check for trim (CLtrim=CL1, CM=0)')
CLtrim=CL0+CLalpha*alpha_trim+CLiH*iH+CLdeltaE*dE_trim
CMtrim=CM0+CMalpha*alpha_trim+CMiH*iH+CMdeltaE*dE_trim

% Now do this repeatedly for a range of speeds
Vmat=260:10:540; % a vector
CL1mat=2*W./(rho*Vmat.*Vmat*S); % a vector

alpha_trim_mat=zeros(size(Vmat));
dE_trim_mat=zeros(size(Vmat));

% Compute trim conditions at various speeds
%  when the c.g. is at the moment reference point
for i=1:length(Vmat)
    V=Vmat(i);
    CL1=CL1mat(i);
    a=[CLalpha CLdeltaE;
    CMalpha CMdeltaE];
    b=[CL1-CLiH*iH-CL0
        -CM0-CMiH*iH];
    x=inv(a)*b;
    alpha_trim_mat(i)=x(1);
    dE_trim_mat(i)=x(2);
end

%alpha_trim_mat
%dE_trim_mat


figure(1); clf
subplot(311)
plot(Vmat,-dE_trim_mat,color(1:3))
ylabel('MINUS dE trim (deg)')
xlabel('Trim speed (f/s)')
grid on
hold on
title('Trim elevator as a function of trim speed (Like Figure 4.7)')


subplot(312)
plot(Vmat,alpha_trim_mat,color(1:3))
ylabel('alpha trim (deg)')
xlabel('Trim speed (f/s)')
grid on
hold on
title('Trim angle of attack as a function of trim speed')

CLtrim=CL0+CLalpha*alpha_trim_mat+CLiH*iH+CLdeltaE*dE_trim_mat;
subplot(313)
plot(Vmat,CLtrim,color(1:3))
ylabel('CL trim)')
xlabel('Trim speed (f/s)')
grid on
hold on
title('Trim lift coefficient as a function of trim speed')



% Move cg to forward cg location
% Compute trim conditions at various speeds
%  when the c.g. is at the forward limit

% Move the equivalent force system
d=mom_ref_pt-forward_cg;
CM0newcg=CM0-d*CL0;
CMalphanewcg=CMalpha-d*CLalpha;    % Per degree
CMdeltaEnewcg=CMdeltaE-d*CLdeltaE;  %Per degree
CMiHnewcg=CMiH-d*CLiH; %Per degree

for i=1:length(Vmat)
    V=Vmat(i);
    CL1=CL1mat(i);
    a=[CLalpha CLdeltaE;
    CMalphanewcg CMdeltaEnewcg];
    b=[CL1-CLiH*iH-CL0
        -CM0newcg-CMiHnewcg*iH];
    x=inv(a)*b;
    alpha_trim_mat(i)=x(1);
    dE_trim_mat(i)=x(2);
end

%alpha_trim_mat
%dE_trim_mat
subplot(311)
plot(Vmat,-dE_trim_mat,color(4:6))


subplot(312)
plot(Vmat,alpha_trim_mat,color(4:6))

CLtrim=CL0+CLalpha*alpha_trim_mat+CLiH*iH+CLdeltaE*dE_trim_mat;
subplot(313)
plot(Vmat,CLtrim,color(4:6));


% Move cg to aft cg location
% Compute trim conditions at various speeds
%  when the c.g. is at the aft limit

% Move the equivalent force system
d=mom_ref_pt-aft_cg;
CM0newcg=CM0-d*CL0;
CMalphanewcg=CMalpha-d*CLalpha;    % Per degree
CMdeltaEnewcg=CMdeltaE-d*CLdeltaE;  %Per degree
CMiHnewcg=CMiH-d*CLiH; %Per degree

for i=1:length(Vmat)
    V=Vmat(i);
    CL1=CL1mat(i);
    a=[CLalpha CLdeltaE;
    CMalphanewcg CMdeltaEnewcg];
    b=[CL1-CLiH*iH-CL0
        -CM0newcg-CMiHnewcg*iH];
    x=inv(a)*b;
    alpha_trim_mat(i)=x(1);
    dE_trim_mat(i)=x(2);
end

%alpha_trim_mat
%dE_trim_mat
subplot(311)
plot(Vmat,-dE_trim_mat,color(7:9))
str1=['CG at MRP, xbarcg= ',num2str(mom_ref_pt)];
str2=['CG Forward, xbarcg= ',num2str(forward_cg)];
str3=['CG Aft, xbarcg= ',num2str(aft_cg)];
legend(str1,str2,str3)
text2(.01,.07,'Note: trim elevator sets trim speed and trim angle of attack')
text2(.01,.9,'ELEVATOR IS REALLY THE SPEEDERATOR')
% Plot max and min elevator lines
z=axis;plot([z(1) z(2)],[max(deltaE) max(deltaE)],'-k')
       plot([z(1) z(2)],[min(deltaE) min(deltaE)],'-k')

subplot(312)
plot(Vmat,alpha_trim_mat,color(7:9))
legend(str1,str2,str3)
text2(.01,.09,'Note: trim elevator sets trim speed and trim algle of attack')

CLtrim=CL0+CLalpha*alpha_trim_mat+CLiH*iH+CLdeltaE*dE_trim_mat;
subplot(313)
plot(Vmat,CLtrim,color(7:9))
legend(str1,str2,str3)
hold off
% Done with figure 1

disp('Compute dDeltaE/dCL vs cg position')
% move cg to forward cg location
d=mom_ref_pt-forward_cg;
CM0newcg=CM0-d*CL0;
CMalphanewcg=CMalpha-d*CLalpha;    % Per degree
CMdeltaEnewcg=CMdeltaE-d*CLdeltaE;  %Per degree
CMiHnewcg=CMiH-d*CLiH; %Per degree
% Use equation 4.53 on page 200
dDeltaEperdCL(1)=-CMalphanewcg/(CLalpha*CMdeltaEnewcg-CMalphanewcg*CLdeltaE);
% move cg to ref point
d=0;
CM0newcg=CM0-d*CL0;
CMalphanewcg=CMalpha-d*CLalpha;    % Per degree
CMdeltaEnewcg=CMdeltaE-d*CLdeltaE;  %Per degree
CMiHnewcg=CMiH-d*CLiH; %Per degree
% Use equation 4.53 on page 200
dDeltaEperdCL(2)=-CMalphanewcg/(CLalpha*CMdeltaEnewcg-CMalphanewcg*CLdeltaE);
% move cg to aft cg location
d=mom_ref_pt-aft_cg;
CM0newcg=CM0-d*CL0;
CMalphanewcg=CMalpha-d*CLalpha;    % Per degree
CMdeltaEnewcg=CMdeltaE-d*CLdeltaE;  %Per degree
CMiHnewcg=CMiH-d*CLiH; %Per degree
% Use equation 4.53 on page 200
dDeltaEperdCL(3)=-CMalphanewcg/(CLalpha*CMdeltaEnewcg-CMalphanewcg*CLdeltaE);


% assemble data for plotting
Xcg=[forward_cg mom_ref_pt aft_cg]
dDeltaEperdCL
% Start figure 2
figure(2); clf
plot(Xcg,dDeltaEperdCL,'-ko'); hold on
xlabel('Xbar cg (fraction of chord)')
ylabel('dDeltaE/dCLtrim')
z=axis;
axis([0 .7 z(3) 0])
%[m,b]=straightline(x,y)
[m1,b1]=straightline(Xcg,dDeltaEperdCL);
disp('Compute neutral point by finding the c.g location on')
disp(' figure 2 where dDeltaEperdCL=0')
xNP=-b1/m1
plot([aft_cg xNP],[dDeltaEperdCL(3) 0],'--ko')
grid on
title('Determination of neutral point from flight test data. Figure 4.6 page 203')
text2(.1,.95,['Stick fixed neutral point is at ',num2str(xNP),' cbar'])
text2(.1,.85,'See equation 4.53 on page 200.')
hold off
% Done with figure 2
