% FUNCTION: fclimb_06
%
% Calculates the rates of change of aircraft state
% variables during the aircraft climb trajectory

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

% Calls: FUNCTION    densityAltitudeoffISA - estimates off ISA conditions
 %                               drag_03.m                - calculates total drag    
%                                thrust_calculation.m  - calculates available thrust

% Outputs:

% An aircraft state vector (y) with the following elements (each column is a variable):

% y(1) - aircraft altitude (m)
% y(2) - aircraft weight (N)
% y(3) - distance traveled along the flight path (m)
% y(4) - distance traveled in the horizontal plane (m)

% a time vector (t) containinig discrete times of computation of states (y)

function yprime = fclimb_06(t,y)

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

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,dAlt] = densityAltitudeoffISA(deltaTemp,y(1),V);
%[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);
[drag,LD,cl,cd]=drag_03(mtrue,rho,a,y(2)/g,macht,Cdoct,phi);

% Estimates the thrust during climb

[thrust] = thrust_calculationNoLoss(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)

yprime(1) = vtas * (thrust - drag) / (y(2)) ;		 % Rate of climb (m/s)
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
                           
% print the results

fid = fopen ('climb_output.txt','a');
fprintf(fid, '%g %4.0f %4.0f %4.0f %4.0f %4.0f %4.0f %4.0f %4.0f %4.0f %4.0f %4.5f %4.5f\n',...
    t,drag,thrust,yprime(1), yprime(2), yprime(3), yprime(4), y(1), y(2), y(3), y(4), mtrue, flt_path_angle);
status = fclose(fid);
                                                            
                                           

