% Aircraft trim diagram
clear all
%close all
disp(' '); disp('Start Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<')
disp(' Data Input section')

% Load script for particular aircraft
disp('Load aircraft: Units expected are degree or per degree')
Trim_Data_Cessna182

% 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


disp(' First we will plot the CL and Cm versus alpha for nominal cg position')
alpha=0:1:alpha_PlotMax;
dCMdCL=CMalpha/CLalpha;
CM0bar=CM0-dCMdCL*CL0;
CMiHbar=CMiH-dCMdCL*CLiH;
CMdeltaEbar=CMdeltaE-dCMdCL*CLdeltaE;

StaticMarginNominal=-CMalpha/CLalpha
Xneutralpoint=StaticMarginNominal+mom_ref_pt

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,s5,'Location',LOC)
LOC='SouthEast';
subplot(212); legend(s1,s2,s3,s4,s5,'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
strw=['CM about ',num2str(mom_ref_pt),'c'];
ylabel(strw)
strq=['Cm and CL vs alpha: ',aircraft];
title(strq)
text2(.02,.92,['iH= ',num2str(iH),' deg.'])


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


% Data is selected to duplicate Roskam figures 4.9 <<<<<<<<<<<
disp(' '); disp(' Second we will plot the aircraft trim diagram and include other c.g. data')
alpha=0:1:alpha_PlotMax;
dCMdCL=CMalpha/CLalpha;
CM0bar=CM0-dCMdCL*CL0;
CMiHbar=CMiH-dCMdCL*CLiH;
CMdeltaEbar=CMdeltaE-dCMdCL*CLdeltaE;

% Plot aircraft trim diagram
figure(2) % <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
clf

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,s5,'Location',LOC)
subplot(121); legend(s1,s2,s3,s4,s5,'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

str1=['CM about ',num2str(mom_ref_pt),' c'];
xlabel(str1)
ylabel('CL')
%title(strq)
grid on
text(.18,.75,['iH= ',num2str(iH),' deg.'])
strXF=['fwd cg xbar=',num2str(forward_cg)];
strXR=['aft cg xbar=',num2str(aft_cg)];
text(.16,1.18,strXF)
text(-.02,1.18,strXR)
stralpha=['alpha stall=',num2str(alpha_stall), ' deg.'];
text(.09,.95,stralpha)

subplot(121);z2=axis; axis([z2(1) z2(2) 0 Cl_PlotMax]);
ylabel('CL')
xlabel('alpha (deg)')
strq=['Aircraft Trim Diagram: ',aircraft];
title(strq)
grid on
text2(.1,.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');
axis([-.2 .2 0 Cl_PlotMax])
set(gca, 'XDir', 'reverse'); % reverse the plotting direction on the x axis

% End alpha_stall line plotting




