disp(' '); disp('>>>>>>>>>>>>Start here <<<<<<<<<<<<<'); disp(' ')
clear
disp('Parameters for MPX5 at sea level')
aircraft='MPX-5'
g=32.17      % accel of gravity, ft/sec^2
h=0      % altitude, ft
Rho=rhofun(h)% air density, Slugs/ft^3
%Rho=0.00237;  % air density, Slugs/ft^3
S=4.333;     % ft^2
HalfRhoS=0.5*Rho*S;
W=19
mass=W/g;    % mass
PhiT=0;      % Thrust inclination, rad
Cl0=0.075711;
Clalpha=3.8357;
Cd0=0.018;
K=0.0464;
BhpMax=.1     % maximum horsepower
EtaP=.65       % Propeller efficiency

% Establish guess for the initial conditions (these are not too important
% as they will be reset later).
v0=50;  % Trim on backside of power required curve
%v0=80;  % Trim on frontside of power required curve
gamma0=0;


% Reference about which to trim (these are important).
vR=v0;            % ft/sec
gammaR=gamma0;          % rad

XR=[vR,gammaR]';       %Reference about which to trim
IX=[1,2]';             % Fix x but vary u to acvheive trim
YR=[vR,gammaR,vR*sin(gammaR)]'; %Reference about which to trim
IY=[1,2,3]';           % Fix y but vary u to acvheive trim
%    [X,U,Y,DX]=TRIM('SYS',X0,U0,Y0,IX,IU,IY)
[XR,UR,YR,DX]=trim('PhugoidSubProp',XR,[],YR,IX,[],IY);
disp('Reference conditions about which we will linearize')   
XR        % Reference (trim) state vector about which we will linearize
UR        % Reference (trm) control vector
YR        % This is the reference (trim) output vector 

% Reset the initial conditions for the simulation based on a perturbation from trim
disp('Initial conditions for the simulation')
v0=XR(1)
gamma0=XR(2)
alphaTrim=UR(1)
BhpTrim=UR(2)


if (1)  % Change to 0 to a 1 if you want to check the derivatives at trim
	flag=1;
	t=0;
    disp('Check derivatives of state vector at our trim conditions')
    [xdot] = PhugoidDynProp(t,XR,UR,flag,g,HalfRhoS,mass,PhiT,Cd0,K,Cl0,Clalpha,EtaP,v0,gamma0)
    %[sys,x0,str,ts] =PhugoidDynProp(t,x,u,flag,g,HalfRhoS,mass,PhiT,Cd0,K,Cl0,Clalpha,EtaP,v0,gamma0)
end



% Power Required Curve
% PLOT POWER REQUIRED CURVE
v=40:.1:100;
CL=W./(.5*Rho*v.*v*S);
CD=Cd0+K*CL.*CL;
Power=.5*Rho*v.^3.*S.*CD;
figure(1)
plot(v,Power)
str1=['Power Required Curve for the ',aircraft];
title(str1); xlabel('trim speed (f/s)');ylabel('Power (ft-lbf/s)')
hold on;
CLt=W/(.5*Rho*v0*v0*S);
CDt=Cd0+K*CLt*CLt;
Powert=.5*Rho*v0^3*S*CDt;
plot(v0,Powert,'rx'); hold off
figure(2)
plot(v,CL)
hold on; plot(v0,CLt,'rx'); 
xlabel('Trim speed (ft/sec)')
ylabel('CL')
str2=['Trim lift coefficient for the ',aircraft];
title(str2);
hold off
