function [sys,x0,str,ts] = PhugoidDynProp2(t,x,u,flag,g,HalfRhoS,mass,PhiT,Cd0,K,Cl0,Clalpha,EtaP,v0,gamma0)
%   S-file for studying simplified phugoid mode dynamics.
%Units x(1)=y(1)=v [ft/sec]
%      x(2)=y(2)=gamma [rad]
%           y(3)=hdot [ft/sec]
%  g is the acceleration of gravity. 
%  HalfRhoS is one half rho times wing area (S).
%  mass is the aircraft in slugs.
%  PhiT is the thrust inclination angle in rad.
%  Cdo and k are parameters of the drag polar Cd=Cd0+kCl*Cl
%  Cl=Cl0+Clalpha*alpha(rad)
%  v0, gamma0 are the initial condition on v and gamma.
switch flag,      
  case 0,         % Initialization
    [sys,x0,str,ts]=PhugoidDynPropInitializeSizes(v0,gamma0);
  case 1,         % Compute derivatives of continuous states
    sys=PhugoidDynPropDerivatives(t,x,u,g,HalfRhoS,mass,PhiT,Cd0,K,Cl0,Clalpha,EtaP);
  case 2,
    sys=PhugoidDynPropUpdate(t,x,u);
  case 3,
    sys=PhugoidDynPropOutputs(t,x,u);  % Compute output vector
  case 4,                   % Compute time of next sample
    sys=PhugoidDynPropGetTimeOfNextVarHit(t,x,u);
  case 9,                   % Finished. Do any needed 
    sys=PhugoidDynPropTerminate(t,x,u);
  otherwise                 % Invalid input
    error(['Unhandled flag = ',num2str(flag)]);
end

%***********************************************************
%*                    PhugoidDynPropInitializeSizes                   *
%***********************************************************
function [sys,x0,str,ts]=PhugoidDynPropInitializeSizes(v0,gamma0)
% Return the sizes of the system vectors, initial 
% conditions, and the sample times and offets.
sizes = simsizes;   % Create the sizes structure
sizes.NumContStates  = 2;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 3;
sizes.NumInputs      = 2;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;   % at least one sample time is needed
sys = simsizes(sizes);      % load sys with the sizes structure
x0  = [v0,gamma0];    % Specify initial conditions for all states 
str = [];    % str is always an empty matrix
ts  = [0 0]; %initialize the array of sample times


%***********************************************************
%*                    PhugoidDynPropDerivatives                       *
%***********************************************************
function sys=PhugoidDynPropDerivatives(t,x,u,g,HalfRhoS,mass,PhiT,Cd0,K,Cl0,Clalpha,EtaP)
% Compute derivatives of continuous states
%Units x(1)=y(1)=v [ft/sec]
%      x(2)=y(2)=gamma [rad]
%           y(3)=hdot [ft/sec]
%  g is the acceleration of gravity. 
%  HalfRhoS is one half rho times wing area (S).
%  mass is the aircraft in slugs.
%  PhiT is the thrust inclination angle in rad.
%  Cdo and k are parameters of the drag polar Cd=Cd0+kCl*Cl
%  Cl=Cl0+Clalpha*alpha(rad)
%  v0, gamma0 are the initial condition on v and gamma.
%  EtaP is the propeller efficiency
v=x(1);
gam=x(2);
alpha=u(1);
Bhp=u(2);
Thrust=550*Bhp*EtaP/v; % Thrust model for prop
Cl=Cl0+Clalpha*alpha;
Cd=Cd0+K*Cl*Cl;
Drag=HalfRhoS*v*v*Cd;
Lift=HalfRhoS*v*v*Cl;
v_dot = (Thrust/mass)*cos(PhiT+alpha) - Drag/mass - g*sin(gam);
gamma_dot = (Thrust/(mass*v))*sin(PhiT+alpha) + Lift/(mass*v) - (g/v)*cos(gam);
sys = [v_dot,gamma_dot] ;

%***********************************************************
%*                    PhugoidDynPropUpdate                            *
%***********************************************************
function sys=PhugoidDynPropUpdate(t,x,u)
% Compute update for discrete states. If necessary, check for
% sample time hits.
sys = [];    % Empty since this model has no discrete states.

%***********************************************************
%*                    PhugoidDynPropOutputs                           *
%***********************************************************
function sys=PhugoidDynPropOutputs(t,x,u)
% Compute output vector given current state, time, and input
% sys=y=[v(ft/sec), gamma(rad), hdot(ft/sec)]
sys = [x(1),x(2),x(1)*sin(x(2))];

%***********************************************************
%*                    PhugoidDynPropGetTimeOfNextVarHit               *
%***********************************************************
function sys=PhugoidDynPropGetTimeOfNextVarHit(t,x,u)
% Return the time of the next hit for this block.  Note that 
% the result is absolute time.  Note that this function is 
% only used when you specify a variable discrete-time sample
% time [-2 0] in the sample time array in sampleTime = 1;    
sys = [] ;

%************************************************************
%*                    PhugoidDynPropTerminate                          *
%************************************************************
function sys=PhugoidDynPropTerminate(t,x,u)
% Perform any necessary tasks at the end of the simulation
sys = [];

