% FUNCTION: fclimb_04 = restricted climb profile
% Calculates the rates of change of aircraft state
% variables during an aircraft climb profile

% This function differs from fclimb_06 because it restrictes the aircraft
% climb profile to some predefined cruise altitude hcruise

% programmer: A. Trani
% date: Sept. 17, 2004

% Calls: FUNCTION  isam.m
%                                drag_03.m
%                                thrust_calculation.m

% Called from: various climb segment profiles

function yprime = fclimb_04(t,y)

global  g tsfc macht Cdoct thrust_table mach_table lapse_rate_factor rhos Vclimb altc hcruise

V = interp1 (altc,Vclimb,y(1));  			     % speed schedule - speed is calibrated airspeed (CAS) or indicated airspeed (IAS)

% Computes standard atmosphere values

% Speed of sound in m/s at different cruise altitudes in meters

[mtrue,a,rho,temp] = isam(y(1),V);          % calculates the true mach number of the aircraft at every point in the climb profile

vtas = mtrue * a;                   % true airspeed (m/s)
q = 0.5 * rho * vtas^2;             % dynamic pressure (N)
phi = 0.0;                              % bank angle (radians)
	  
% Compute drag and L/D ratio (call drag_new)
[drag,LD,cl,cd]=drag_03(mtrue,rho,a,y(2)/g,macht,Cdoct,phi);

% Estimates the thrust during climb

[thrust] = thrust_calculation(thrust_table,mach_table,lapse_rate_factor,mtrue,rho);

% Define the rate equations (3 rate equations for 3 state variables)
% y(1) - Aircraft altitude (m)
% y(2) - Aircraft weight (N)
% y(3) - Distance traveled along the path (m)
% y(4) - Distance traveled wrt to ground (m)

if y(1) < hcruise                                                    % Statement to zero ROC when cruise altitude is reached
    yprime(1) = vtas * (thrust - drag) / (y(2)); 		 % Rate of climb (m/s)
else
    yprime(1) = 0;                                                  % zero rate of climb (m/s)
end

yprime(2) = - tsfc * thrust;								% Mass flow rate (N/s) 
yprime(3) = vtas ;											   % Distance traveled along the flight path (m)

% Calculation of flight path angle (gamma)

flt_path_angle = asin(yprime(1)/ vtas);     % radians

% Calculate distance traveled across horizontal distance

yprime(4) = vtas * cos(flt_path_angle);     % distance traveled along horizontal (meters)

yprime = yprime'; % transposes the array going out
															% MATLAB is very picky about this
															% it likes the array in column format
                                                            
                                           

