disp(' '); disp('>>>>>>>>>>>>Start here <<<<<<<<<<<<<'); disp(' ')
clear
% Parameters for Predator at sea level
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=112        % wing area, ft^3
HalfRhoS=0.5*Rho*S;
W=2250      % weight, lbf
mass=W/g;    % mass
PhiT=0;      % Thrust inclination, rad
Cd0=0.0209    % parasite drag (non-dim)
K=0.03     % square ther in drag polar (non-dim)
Cl0=0.738     % lift coef at zero alpha (non-dim)
Clalpha=8.1581 % lift curve slope, per rad
BhpMax=81     % maximum horsepower
EtaP=.7       % Propeller efficiency

% Establish guess for the initial conditions (these are not too important
% as they will be reset later).
v0=180;
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)
[XR2,UR2,YR2,DX]=trim('PhugoidSubProp',XR,[],YR,IX,[],IY);
disp('Reference conditions about which we will linearize')   
XR2        % Reference (trim) state vector about which we will linearize
UR2        % Reference (trm) control vector
YR2        % 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=XR2(1)
gamma0=XR2(2)
alphaTrim=UR2(1)
BhpTrim=UR2(2)

% Check derivatives of state vector at our trim conditions
if (0)  % Change to 0 to a 1 if you want to check the derivatives at trim
	flag=1;
	t=0;
	[sys,x0,str,ts] = PhugoidDynProp(t,XR,UR,flag,g,HalfRhoS,mass,PhiT,Cd0,K,Cl0,Clalpha,EtaP,v0,gamma0)
	% Executing this code will generate a harmless warning.
end
