% Step2_Trim.m       Version 9.2.1 1/8/07
%
% OBJECTIVE: Trim the aircraft and develop initial conditions for the
%            SIMULINK nonlinear simulation.
%
% INPUTS: the array called constant
% OUTPUTS: In the Matlab workspace will be defined initial conditions
% on the state and controls (xIC and uIC). The array constant, xIC and uIC
% are all required by the SIMULINK nonlinear simulation.
%
%

% Quick trim the Predator
disp(' '); disp(' ')
disp(' >>>>>>>>>>>Start here, Step2_Trim <<<<<<<<<<<<<<<')
disp(' ')

disp(aircraft); % Display aircraft identification
% Establish a trim condition (initial condition).
Vt=constant(58); % Trim airspeed (ft/sec)
Hp=constant(59); % Trim altitude (ft)
% Alternately you can make changes to the trim airspeed and 
% altitude here if the nondimensional stability and control
% derivatives in the array called constant are not strongly
% dependent on trim airspeed ar alitude, as is often the case.

% At this point in the script if you want you can change the c.g. location as follows.
% constant(57)=.25;  		
% XbarCG,  nondimensional, measured aft from leading edge of wing mean aerodynamic chord.

s1=['Trim airspeed=',num2str(Vt),' ft/sec, '];
s2=['Trim altitude=',num2str(Hp),' ft, '];
s3=['c.g.=',num2str(constant(57)),' cbar.'];
disp(' '); disp([s1 s2 s3]); disp(' ');

% Define the steady winds
disp('NED omponents of the steady wind')
VwindKnots=0;
Vwx=-VwindKnots*1.687808 %ft/sec, - sign means headwind
Vwy=0
Vwz=0 

% Location of accelerometer package
Xa(1)=0; % x-location of accelerometer wrt c.g. (ft)
Xa(2)=0; % y-location of accelerometer wrt c.g. (ft)
Xa(3)=0; % z-location of accelerometer wrt c.g. (ft)

% Trim the aircraft
disp('Trim conditions')
echo on
% States x
%   x=[U,V,W,P,Q,R,Phi,Theta,Psi,Xprime,Yprime,Hprime]
%   Units are [ft/sec, ft/sec, ft/sec, rad/sec, rad/sec, rad/sec, rad, rad, rad, ft, ft, ft]
% Inputs u
%   u=[deltaE,deltaA,deltaR,bhp, Vwx, Vwy, Vwz]]
%   Units are [rad,rad,rad,bhp, ft/sec, ft/sec, ft/sec]
echo off
[xIC,uIC,CL,CD,CM,alphadeg]=QuickTrim2(Vt,Hp,constant,Vwx,Vwy,Vwz)

% Check to insure that the derivatives of the state vector (sys) are small 
% for this trim condition.
flag=1;;t=0;
[Xdot] = aircraft9(t,xIC,uIC,flag,constant,xIC,Xa);
disp('Trim Check: The first 9 derivatives (Xdot) should be small at trim')
Xdot
disp('  Now you are ready to run the nonlinear simulation using')
disp('  SIMULINK with the model called FlatEarth_MATLAB7.mdl from MATLAB 7')
disp('  or to execute Step3_Simulink_SimMATLAB7.m.')


% PLOT POWER REQUIRED CURVE
w=constant(1);
rho=rhofun(constant(59));
%rho=.002
S=constant(20);
%S=3
Cdo=constant(25);
%Cdo=.025
k=constant(26);
%k=.1
v=20:.1:120;
CL=w./(.5*rho*v.*v*S);
CD=Cdo+k*CL.*CL;
Power=.5*rho*v.^3.*S.*CD;
plot(v,Power)
title('Power Required Curve'); xlabel('trim speed (f/s)');ylabel('Power (ft-lbf/s)')
hold on;
CLt=w/(.5*rho*Vt*Vt*S);
CDt=Cdo+k*CLt*CLt;
Powert=.5*rho*Vt^3*S*CDt;
plot(Vt,Powert,'rx'); hold off
figure(2)
plot(v,CL)
hold on; plot(Vt,CLt,'rx'); 
xlabel('Trim speed (ft/sec)')
ylabel('CL')
hold off


