% Aircraft trim diagram
clear all
disp(' '); disp('Start Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<')

disp(' Data Input section')
% Data is selected to duplicate Roskam figures 4.5 and 4.9
CL0=.2
CLalpha=.1      % Per degree
CLdeltaE=.01    % Per degree
CLiH=3*CLdeltaE % Per degree

CM0=.1
CMalpha=-.02    % Per degree
CMdeltaE=-.005  %Per degree
CMiH=3*CMdeltaE %Per degree

deltaE=[-20 -10 0 10]
iH=0             % Degree
%iH=-6           % 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
% End of data input section


% 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.']
% End plotting information

alpha=0:1:alpha_PlotMax;
dCMdCL=CMalpha/CLalpha
CM0bar=CM0-dCMdCL*CL0
CMiHbar=CMiH-dCMdCL*CLiH
CMdeltaEbar=CMdeltaE-dCMdCL*CLdeltaE



figure(1)
clf
for i=1:length(deltaE)
dE=deltaE(i);
CM=CM0+CMalpha*alpha+CMiH*iH+CMdeltaE*dE;
CL=CL0+CLalpha*alpha+CLiH*iH+CLdeltaE*dE;
subplot(211);plot(alpha,CM,color(3*(i-1)+1:3*(i-1)+3)); hold on
subplot(212);plot(alpha,CL,color(3*(i-1)+1:3*(i-1)+3)); hold on
end
LOC='NorthEast'
subplot(211); legend(s1,s2,s3,s4,'Location',LOC)
LOC='SouthEast'
subplot(212); legend(s1,s2,s3,s4,'Location',LOC)

subplot(211);plot(alpha,0*alpha,'g') %plot zero line on CM plot
hold off
subplot(211);z1=axis;axis([z1(1) z1(2) -.4 .6]); grid on
ylabel('CM c/4')
title('Figure 4.5, page 201 of Roskam')
text2(.02,.92,['iH= ',num2str(iH),' deg.'])


subplot(212);z2=axis;axis([z2(1) z2(2) 0 1]); grid on
ylabel('CL')
xlabel('alpha (deg)')
text2(.02,.92,['iH= ',num2str(iH),' deg.'])
hold off


% Plot aircraft trim diagram
figure(2)
clf

% Temporary plot for scaling purposes only
for i=1:length(deltaE)
dE=deltaE(i);
CL=CL0+CLalpha*alpha+CLiH*iH+CLdeltaE*dE;
CM=CM0bar+dCMdCL*CL+CMiHbar*iH+CMdeltaEbar*dE;
subplot(122);plot(-CM,CL); hold on
subplot(121);plot(alpha,CL); hold on
end
% End temporary plot

% Extend plots
subplot(122); z1=axis
clf
CLexpended=z1(3):.1:z1(4);
for i=1:length(deltaE)
dE=deltaE(i);
CL=CL0+CLalpha*alpha+CLiH*iH+CLdeltaE*dE;
subplot(121);plot(alpha,CL,color(3*(i-1)+1:3*(i-1)+3)); hold on
CMexpanded=CM0bar+dCMdCL*CLexpended+CMiHbar*iH+CMdeltaEbar*dE;
subplot(122);  hold on; plot(-CMexpanded,CLexpended,color(3*(i-1)+1:3*(i-1)+3));
end

LOC='SouthEast'
subplot(122); legend(s1,s2,s3,s4,'Location',LOC)
subplot(121); legend(s1,s2,s3,s4,'Location',LOC)
% end of expanded plotting

subplot(122);plot([0 0],[0 z1(4)],'--k') %plot zero line on CM plot
hold off
subplot(122);z1=axis;axis([z1(1) z1(2) 0 Cl_PlotMax])
% plot the forward cg line of the trim triangle
delta_cg1=mom_ref_pt-forward_cg
Cm_forward=+Cl_PlotMax*delta_cg1
subplot(122);hold on; plot([0 -Cm_forward],[0 Cl_PlotMax],'k') %plot zero line on CM plot
% plot the aft cg line of the trim triangle
delta_cg2=mom_ref_pt-aft_cg
Cm_aft=+Cl_PlotMax*delta_cg2
subplot(122);hold on; plot([0 -Cm_aft],[0 Cl_PlotMax],'k') %plot zero line on CM plot


%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




str1=['MINUS CM about ',num2str(mom_ref_pt),' c']
xlabel(str1)
ylabel('CL')
title('A/C Trim Diag: Figure 4.9, page 207 of Roskam')
grid on
text2(.02,.95,['iH= ',num2str(iH),' deg.'])

subplot(121);z2=axis;axis([z2(1) z2(2) 0 Cl_PlotMax])
ylabel('CL')
xlabel('alpha (deg)')
title('A/C Trim Diag')
grid on
text2(.02,.95,['iH= ',num2str(iH),' deg.'])
hold off





% Plot alpha_stall line in figure 2
CL=CL0+CLalpha*alpha_stall+CLiH*iH+CLdeltaE*deltaE;
CM=CM0+CMalpha*alpha_stall+CMiH*iH+CMdeltaE*deltaE;
subplot(122); hold on; plot(-CM,CL,'k');
% End alpha_stall line plotting




