% Delivered to the Government in Accordance with Contract FA8650-10-C-7080
% 
%  (c) 2011 United Technologies Corp.
% Developed with the sponsorship of the Defense Advanced Research Projects Agency (DARPA).
%  
% Permission is hereby granted, free of charge, to any person obtaining a copy of this data, including any software or models in source or binary form, as well as any drawings, specifications, and documentation (collectively "the Data"), to deal in the Data without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Data, and to permit persons to whom the Data is furnished to do so, subject to the following conditions:
%  
% The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Data.
%  
% THE DATA IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS, SPONSORS, DEVELOPERS, CONTRIBUTORS, OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE DATA OR THE USE OR OTHER DEALINGS IN THE DATA.
% 

% This script uses multi-dimensional Newton-Raphson along with the first order
% backward Euler method to compute the transient performance of a 1 spool
% turbojet

% Set Domain Wide Constants (ISA)
global cp_cold cp_hot R gam_cold gam_hot

cp_cold = 1045;           % Specific Heat at Constant Pressure Cold Section
cp_hot = 1244;            % Specific Heat at Constant Pressure Hot Section
gam_cold = 1.4;           % Ratio of Specific Heats Cold Section
gam_hot = 1.3;            % Ratio of Specific Heats Hot Section
R = 287;                  % Gas Constant

G_fuel_design = 0.4174;   % Design fuel flow
fuel_percentage = 5;      % Fuel step increase and decrease percentage

I = 50;                   % Spool Inertia ISA units

w(1) = 630;               % Initial Spool Speed (design)

fprintf('Initial Spool Speed (Rad/s):             %9.2f \n', w(1));

time_step = 0.25;
max_time = 150;
time = 0:time_step:max_time;
n_timesteps = floor(max_time/time_step+1);

d_constant = 1e-2;              % Beta perturbation value for Jacobian calculation
tolerance = 1e-9;

beta_new(1:2,1) = 0.5;          % First guess of auxiliary beta coordinates
Error(1:2,1) = 10*tolerance;    % Initialize Error
torque(1) = 0;                  % Initialize Torque

T_amb = 216.15;                 % Ambient conditions
P_amb = 22729;
Mach = 0.8;

% Fuel flow profile has up step and down step at 50 and 100 seconds
G_fuel(1:floor(n_timesteps/3)) = G_fuel_design;
G_fuel(floor(n_timesteps/3):floor(2*n_timesteps/3))= G_fuel_design*(1+fuel_percentage/100);
G_fuel(ceil(2*n_timesteps/3):n_timesteps)= G_fuel_design*(1-fuel_percentage/100);

t_NL = 0;
t_DF = 0;


for timer = 1:n_timesteps
    
    disp(time(timer))
    
    % Reset beta map variable to middle of map if fuel flow changes
    if timer>1 &&(abs(G_fuel(timer)/G_fuel(timer-1))-1)>tolerance
        beta_new(1:2,1) = 0.5;
    end
    
    % Initiate iteration number
    k = 1;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Initial Guess
    
    
    % Diffuser
    [T1 P1] = Diffuser(T_amb,P_amb, Mach);
    
    % Compressor equations
    [T2, P2, G1, Speed_comp] = HPC(T1, P1, w(timer), beta_new(1));
    
    
    % Burner equations
    [T3 P3] =  Burner(T2, P2, G1, G_fuel(timer));
    
    % Turbine Equations
    [T4, P4, G2, Speed_turb] = HPT(T3, P3, w(timer), beta_new(2));
    
    % Nozzle
    [T5, G3] = Nozzle(T4,P4,P_amb);
    

    % Initialize Error
    Error(1,1) = G3-G1;
    Error(2,1) = G3-G2;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    tic;
    
%     waitforbuttonpress
    
    while max(abs(Error(:,k)))>tolerance
        
        % Count Iterations
        k = k+1;
        
        % Memorize number of iterations
        mem_k(timer) = k;
        
        %Compressor equations
        [T2, P2, G1, Speed_comp] = HPC(T1, P1, w(timer), beta_new(1));
        
        % Burner equations
        [T3 P3] =  Burner(T2, P2, G1, G_fuel(timer));
        
        % Turbine Equations
        [T4, P4, G2, Speed_turb] = HPT(T3, P3, w(timer), beta_new(2));
        
        % Nozzle
        [T5, G3] = Nozzle(T4,P4,P_amb);
        

        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Error calculations
        Error(1,k) = G3-G1;
        Error(2,k) = G3-G2;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Calculate derivitives for to construct Jacobian
        
        
        for perturb = 1:2
            
            % Perturb Apropriate Beta Coordinate and Calculate All Changes
            
            beta_new(perturb) = beta_new(perturb)+d_constant;
            
            %Compressor
            [T2_p, P2_p, G1_p, Speed_p] = HPC(T1, P1, w(timer), beta_new(1));
            
            % Burner
            [T3_p P3_p] =  Burner(T2_p, P2_p, G1_p, G_fuel(timer));
            
            % Turbine
            [T4_p, P4_p, G2_p, Speed_p] = HPT(T3_p, P3_p, w(timer), beta_new(2));
            
            % Nozzle
            [T5_p, G3_p] = Nozzle(T4_p,P4_p,P_amb);
            
            % New Error
            Error_change(1,perturb) = G3_p-G1_p;
            Error_change(2,perturb) = G3_p-G2_p;
            
            % Reset Beta Coordinate Back to Original Value
            beta_new(perturb) = beta_new(perturb)-d_constant;
        end
        
        
        
        for c1 = 1:2
            for c2 = 1:2
                
                % Jacobian calculation
                J(c1,c2) = (Error_change(c1,c2)-Error(c1,k))/(d_constant);
                
            end
        end
        
        
        
        J_mem(1,:,:) = J;                           % Memorize Jacobian
        
        
        Current_Error = Error(:,k);                 % Memorize Current Error
        beta_new = beta_new-J^(-1)*Current_Error;   % Calculate new beta map coordinate
        
        
    end
    
    
    t_DF = t_DF + toc;
    
    % Memorize Relative Corrected Speed for each component
    mem_speed(1,timer) = Speed_comp;
    mem_speed(2,timer) = Speed_turb;
    
    % Memorize the auxiliary beta coordinate
    beta_mem(1:2,timer) = beta_new;
    
    % Torque Calculations
    
    % Compressor Torque
    torque_a_c(timer) =cp_cold*(T1-T2)*G1/w(timer);
    
    % Turbine Torque
    torque_a_t(timer) = cp_hot*(T3-T4)*G1/w(timer);
    
    % Total torque on spool
    torque(timer) = torque_a_t(timer)+torque_a_c(timer);
    
    % Memorize temperature profile within engine for every time step
    T(1,timer) = T1;
    T(2,timer) = T2;
    T(3,timer) = T3;
    T(4,timer) = T4;
    T(5,timer) = T5;
    
    % Memorize pressure profile within engine for every time step
    P(1,timer) = P1;
    P(2,timer) = P2;
    P(3,timer) = P3;
    P(4,timer) = P4;
    P(5,timer) = P_amb;
    
    % Memorize mass flow for each time step
    G(timer) = G1;
    
    % Calculate acceleration
    w_der(timer) = torque(timer)/I;
    
    % Calculate new speed
    if timer < n_timesteps
        w(timer+1) = w(timer)+w_der(timer)*(time(timer+1)-time(timer));
    end
   
    
end

w = w/630;  % Normalize angular velocity with the design point value








