function thrust = enginethrust(t, z)

%	This function gives an instantaneous one-dimensional thrust value


global runstate;
global index;
global dt;
global enDTo;
global enDEo;
global enERATE;
global eoPOX;
global etPC;
global eoMOX;
global ecMFUEL;

global cpropep;

% INITIALIZE
if runstate == 0
   runstate = 1;
   
   rho_fuel = grainChem;
end



%	MAIN PROGRAM

		
Ae = pi * (enDEo/2)^2;
At = pi * (enDTo/2 + t*enERATE)^2;			% a function of Dto and nozzle throat erosion rate
epsilon = Ae/At;

if (ecMFUEL > 0) & (eoMOX > 0)
   DP = eoPOX - etPC(index);
   Cd = .8;
   Ainj = .0015; %m2
   mdot_ox = Cd * Ainj * sqrt(DP);
   mdot_f = mdot_ox / 2;
   mdot_p = mdot_ox + mdot_f;
else
   mdot_ox = 0;
   mdot_f = 0;
   mdot_p = 0;
end


eoMOX(index+1)    = eoMOX(index) - mdot_ox * dt(index);
ecMFUEL(index+1)  = ecMFUEL(index) - mdot_f * dt(index);


%Pc = 7000000;	% needs to be a function of mdot_oxidizer
Pe = 101325;	% needs to be a function of chamber pressure and nozzle dimensions


% FIGURE OUT HOW MUCH FUEL MASS SHOULD BURN THIS ITTERATION
Mhtpb = 83;
Mipdi = 11;
Mcarbon = 5;
Mcastor = 1;

Rp = .7*.0254;
Ap = pi * Rp^2;


% CREATE input.dat for cpropep
if mdot_p > 0
      
   cd 'D:\Development\Websites\rlpotter\www\ryan\projects\6DOF\cpropep';
   
   %dos('erase input.dat');   
   %dos('erase output.dat');
      
   fid = fopen('input.dat','w');
   fprintf(fid,'#run time: %f\r\n', t);
	fprintf(fid,'Propellant HTPB/KClO4/Al\r\n');
	fprintf(fid,'+788		%f g\r\n', Mhtpb);
	fprintf(fid,'+491		%f g\r\n', Mipdi);
   fprintf(fid,'+212		%f g\r\n', Mcarbon);
   fprintf(fid,'+228    %f g\r\n\r\n', Mcastor);
	fprintf(fid,'FR\r\n');
	fprintf(fid,'+chamber_pressure %f kPa\r\n', etPC(index)/1000);
	fprintf(fid,'+supersonic_area_ratio %f\r\n', epsilon);
	fclose(fid);

	% RUN CPROPEP
	cpropep = dos('cpropep.exe -f input.dat -o output.dat');

	% RUN PERL
   %dos ('extract_data.pl');
   cd 'D:\Development\Websites\rlpotter\www\ryan\projects\6DOF';
   
   
   
   
	% EXTRACT THE DATA
   
	%load data.dat
   k = 1.2;
   a = 1500; % m/s

	% OPERATE ON THE DATA

	Pc = mdot_p * a / Ap;
	etPC(index+1) = Pc;

	thrust = At * Pc * sqrt(((2*k^2)/(k-1))*((2/(k+1))^((k+1)/(k-1)))*(1-(Pe/Pc)^((k-1)/k))) ...
   	+ (Pe-pressure(z))*Ae;

else
   thrust = 0;
   etPC(index+1) = 0;
end


% ===========================================================================================
% INTERNAL FUNCTIONS
% ===========================================================================================


function rho_fuel = grainChem()

global runstate;

if runstate == 0
   a = 1;
end

rho_fuel = 4;


function rdot = regressionRate()

rdot = .3/.0254;	% m/s



function burn_perimeter = grainGeometry()

%load graingeom.dat;
rdot = regressionRate;

burn_perimeter = 2*pi*.03;		% m



function [T, rho] = saturationCurve()

% these are just the critical point values
T = 154.8;		% K
rho = 1100;		% kg/m^3