Program for design point cycle analysis and parametric analysis of turboprop engines and plotting of parametric results, gas path temperatures and pressures.

	%% Requirement:
	% Pratap wants design point cycle analysis for 2500 HP, 2750 HP class engine and
	% 3000HP class engine; Asked on 28 June 2022.
	%% Constant:
	% * TIT = Kept at the same level as 2750 class engine.
	% * Ratio of LPC and HPC kept same
	% * Ideally it should be done on the basis of work function distribution.
	%% Parameters
	% * Overall pressure ratio
	% * Mass flow rate of air
	clc;clear all;close all; commandwindow;
	doPlots = 0;
	ratioOfLPCtoHPCpr = 1.0272;% (18 = 4.3 * 4.18)
	%TIT_k = 1477+273;%Vlaue of PW150A engine.
	TIT_k = 1470;%Value Pratap asked to use.
	whichEngineToRun = 'parametric';
	whichEngineToRun = '2500hp';
	whichEngineToRun = '2750hp';
	whichEngineToRun = '3000hp';
	whichEngineToRun = '5068.4hp';%PW150 value
	switch lower(whichEngineToRun)
	case 'parametric'
	% % Use for parametric study
	targetHp = 3000;%HP max take off power
	rangeOfTotalPR = 14.5:0.01:16;
	rangeOfMassFlow = 10:0.01:12;
	case '2500hp'
	% Picked values for 2500 HP, 0.43982 SFC = 0.95 * 0.463 SFC
	targetHp = 2500;%HP max take off power
	rangeOfTotalPR = 14.4;
	rangeOfMassFlow = 9.72;
	case '2750hp'
	% % % Picked values for 2750 HP, 0.436 SFC = 0.95 * 0.459 SFC
	targetHp = 2750;%HP max take off power
	rangeOfTotalPR = 14.7;
	rangeOfMassFlow = 10.66;
	case '3000hp'
	% % % Picked values for 3000 HP, 0.432 SFC = 0.95 * 0.455 SFC
	targetHp = 3000;%HP max take off power
	rangeOfTotalPR = 15.0;
	rangeOfMassFlow = 11.58;
	case '5068.4hp'
	% % % Picked values for 5000 HP, 0.263 SFC kg/(kW hr)
	targetHp = 5068.4;%HP max take off power
	rangeOfTotalPR = 18.0;
	rangeOfMassFlow = 11.50;
	end
	[prGrid,mDotGrid] = meshgrid(rangeOfTotalPR,rangeOfMassFlow);
	fpt_power_out = prGrid*0;
	sfc_out = prGrid*0;
	mfDot_out = prGrid*0;
	Thspec_out = prGrid*0;
	Powspec_out = prGrid*0;
	for rowIndex = 1:size(prGrid,1)
	for colIndex = 1:size(prGrid,2)
	% Assign parameters
	pAmb = 101325;
	tAmb = 288.15;
	propeller.powerTakeOff_kW = targetHp*0.746;
	%SFC_eshp_TO = 0.263;%kg/(kW*hr);
	machFreestream = 0;
	inlet.pr = 0.99;
	overallPR = prGrid(rowIndex,colIndex);
	lpc.pr = sqrt(overallPR*ratioOfLPCtoHPCpr);
	lpc.eta = 0.78;
	lpc.mDot = mDotGrid(rowIndex,colIndex);
	lpc.exitArea = 0;
	lpc.exitPr = 0.98;
	hpc.pr = sqrt(overallPR/ratioOfLPCtoHPCpr);
	hpc.eta = 0.77;
	hpc.exitArea = 0;
	hpc.bleedOutFractionOfInletAir = 0.01 ;
	hpc.exitPr = 0.98;
	burner.pr = 0.96;
	burner.fhv = 43e6;
	burner.eta = 0.99;
	%burner.mfDot = propeller.powerTakeOff_kW* SFC_eshp_TO/3600;
	hpt.fractionOfBleedAirReturningInFrontOfHptAndDoingWork = 1/2;
	hpt.eta = 0.87;
	hpt.exitPr = 0.99;
	hpt.excessWork = 30*.746*1000;%30HP on HP spool
	lpt.fractionOfBleedAirReturningInFrontOfLptAndDoingWork = 1/4;
	lpt.eta = 0.87;
	lpt.exitPr = 0.99;
	lpt.excessWork = 30*0.746*1000;
	lpt.exitArea = 0;
	fpt.fractionOfBleedAirReturningInFrontOfFptAndDoingWork = 1/4;
	fpt.eta = 0.87;
	fpt.exitPr = 0.99;
	lpSpool.eta = 0.99;
	lpSpool.rpm = 27000;
	hpSpool.eta = 0.99;
	hpSpool.rpm = 31150;
	fpSpool.eta = 0.99;
	nozzle.cd = 0.99;
	nozzle.ct = 0.99;
	gearBox.excessWork = 30*0.746*1000;
	%% Description of the fucntion:
	% This function is part of project for parameter estimation and simulation
	% of RTA
	%
	% Purpose of this function is to run the design point simulation model
	% using the input arguments given and compare the output of the simulation
	% model to the measurements, that are read from the mat file, and calculate
	% the errors.
	%% Input parameters
	% The input parameters are the component parameters of the simulation model
	% at a specific operating point. It could be design or off design rating.
	%
	%%
	%
	% # pAmb : Ambient pressure in Pascals
	% # tAmb : Ambient temperature in K
	% # mach : Free stream mach number
	% # inlet.pr : Inlet total pressure ratio
	% # lpc.pr : Low pressure compressor pressure ratio
	% # lpc.eta : Low pressure compressor efficiency
	% # lpc.mDot : Low pressure compressor actual mass flow rate
	% # splitter.prColdChannel : Splitter pressure ratio in the cold channel,
	% that leads to the bypass duct. This value should be less than 1 always.
	% # splitter.prHotChannel : Splitter pressure ratio in the hot channel,
	% that leads to the core engine. This value should be less than 1 always.
	% # splitter.bpr : Splitter splitting ratio, division of flow into the core
	% and the bypass streams, $\frac{M_{cold}}{M_{hot}}=bypass\ ratio$
	% # bypassDuct.pr : Duct pressure ratio in the bypass dcuct connecting the
	% lpc exit to the aferburner chamber. This value should be less than 1 always.
	% # hpc.pr : High pressure compressor pressure ratio
	% # hpc.eta : High pressure compressor efficiency
	% # hpc.bleedOutFractionOfInletAir : High pressure compressor 7th stage air
	% beed out, as a fraction of air incoming in the engine inlet.
	% # burner.pr : Main combustion chamber pressure ratio
	% # burner.eta : Main combustion chamber efficiency
	% # burner.mfDot : Main combustion chamber fuel mass flow rate
	% # burner.fhv : Main combustion chamber fuel heating value
	% # hpt.excessWork : Additional work produced by the turbine in excess of
	% the work consumedby the compressor.
	% # hpt.eta : Isentropic efficiency of the turbine
	% # hpt.exitPr : Pressure loss down stream of the HP turbine
	% # lpt.excessWork : Additional work produced by the turbine in excess of
	% the work consumedby the compressor.
	% # lpt.eta : Isentropic efficiency of the turbine
	% # lpt.exitPr : Pressure loss down stream of the LP turbine before static
	% pressure measurement port.
	% # lpSpool.eta : Mechanical transmission efficiency of the low pressure
	% spool.
	% # hpSpool.eta : Mechanical transmission efficiency of the high pressure
	% spool.
	% # Afterburner.pr: Afterburner pressure ratio, total pressure loss happens after
	% before nozzle entry and just after lpt. Any loss of total pressure from
	% the shock waves have to be absorbed in this total pressure loss factor.
	% Nozzle itself does not have any total pressure loss factor. It has to be
	% absorbed in after burner only.
	% # Afterburner.mfDot : Fuel flow rate pumped into the afterburner. Heating
	% value is the same as that of the main burner.
	% # Afterburner.eta : Afterburner combustion chamber efficiency of
	% combustion.
	% # nozzle.cd: The ratio between the aerodynamic exit area and geometric
	% exit area of the nozzle .
	% # nozzle.ct: The ratio between the actual thrust generated and the theoretical momentum thrust
	%% Input sanitation
	commandwindow;
	% fprintf('\n Balaji Sankar: ');
	if (inlet.pr >1 || inlet.pr<0.8)
	fprintf('\n In design point code, inlet pr is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	if (lpc.pr >6 || lpc.pr<3)
	fprintf('\n In design point code, LPC pr is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	if (lpc.eta >1 || lpc.eta <0.6)
	fprintf('\n In design point code, LPC eta is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	if (hpc.pr >10 || hpc.pr<1)
	fprintf('\n In design point code, HPC pr is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	if (hpc.eta >1 || hpc.eta <0.6)
	fprintf('\n In design point code, HPC eta is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	% if (splitter.bpr>1 || splitter.bpr <0.41)
	% fprintf('\n In design point code, bypass ratio is unrealistic; pausing executiona and going into debug mode');
	% keyboard;
	% end
	if (burner.pr>1 || burner.pr<0.4)
	fprintf('\n In design point code, burner PR is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	if (hpt.eta >1 || hpt.eta <0.6)
	fprintf('\n In design point code, HPT eta is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	if (lpt.eta >1 || lpt.eta <0.6)
	fprintf('\n In design point code, LPT eta is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	if (lpt.exitPr >1 || lpt.exitPr <0.6)
	fprintf('\n In design point code, LPT exit Pr is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	if (nozzle.cd >1 || nozzle.cd <0.8)
	fprintf('\n In design point code, Coeffecient Discharge of nozzle is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	if (nozzle.ct >1 || nozzle.ct <0.8)
	fprintf('\n In design point code, Thrust coeffecient nozzle is unrealistic; pausing executiona and going into debug mode');
	keyboard;
	end
	%% Inlet
	mach.ambient = machFreestream;
	tP.inletEntry = pAmb*(1+0.2* mach.ambient* mach.ambient)^3.5;
	tK.inletEntry = tAmb*(1+0.2* mach.ambient* mach.ambient);
	m.inletEntry = lpc.mDot;
	theta.inletEntry = tK.inletEntry /288.15;
	delta.inletEntry = tP.inletEntry /101325;
	mc.inletEntry = m.inletEntry*sqrt(theta.inletEntry)/delta.inletEntry;
	tP.inletExit = tP.inletEntry*inlet.pr;
	tK.inletExit = tK.inletEntry ;
	m.inletExit = lpc.mDot;
	theta.inletExit = tK.inletExit /288.15;
	delta.inletExit = tP.inletExit /101325;
	mc.inletExit = m.inletExit*sqrt(theta.inletExit)/delta.inletExit;
	[g.inletExit,~,cp.inletExit] = gordonMcBrideAir_withG(tK.inletExit);
	% LPC
	m.lpcExit = m.inletExit;
	tP.lpcExit = tP.inletExit*lpc.pr;
	tKisen.lpcExit = tK.inletExit*(lpc.pr)^((g.inletExit-1)/g.inletExit);
	tK.lpcExit = tK.inletExit + tK.inletExit*((((lpc.pr)^((g.inletExit-1)/g.inletExit))-1)/lpc.eta);
	work.lpc = m.lpcExit*cp.inletExit*(tK.lpcExit-tK.inletExit);
	theta.lpcExit = tK.lpcExit/288.15;
	delta.lpcExit = tP.lpcExit/101325;
	mc.lpcExit = m.lpcExit*sqrt(theta.lpcExit)/delta.lpcExit;
	[g.lpcExit,~,cp.lpcExit] = gordonMcBrideAir_withG(tK.lpcExit);
	tP.lpcExit = tP.lpcExit*lpc.exitPr;
	% sP.lpcExit = staticFromTotal(m.lpcExit,tP.lpcExit,tK.lpcExit, lpc.exitArea,g.lpcExit,'LPC exit');
	% HPC
	m.hpcExit = m.lpcExit - hpc.bleedOutFractionOfInletAir*m.inletExit;
	tP.hpcExit = tP.lpcExit*hpc.pr;
	tKisen.hpcExit = tK.lpcExit*(hpc.pr)^((g.lpcExit-1)/g.lpcExit);
	tK.hpcExit = tK.lpcExit+ tK.lpcExit*((((hpc.pr)^((g.lpcExit-1)/g.lpcExit))-1)/hpc.eta);
	work.hpc = m.lpcExit*cp.lpcExit*(tK.hpcExit-tK.lpcExit);
	theta.hpcExit = tK.hpcExit/288.15;
	delta.hpcExit = tP.hpcExit/101325;
	mc.hpcExit = m.hpcExit*sqrt(theta.hpcExit)/delta.hpcExit;
	[g.hpcExit,~,cp.hpcExit] = gordonMcBrideAir_withG(tK.hpcExit);
	% sP.hpcExit = staticFromTotal(m.hpcExit,tP.hpcExit,tK.hpcExit, hpc.exitArea,g.hpcExit,'HPC exit');
	tP.hpcExit = tP.hpcExit*hpc.exitPr;
	% Main combustion chamber
	tK.burnerExit = TIT_k;
	[g.burnerExit,~,cp.burnerExit]= gordonMcBrideAir_withG(tK.burnerExit);
	burner.mfDot = ((m.hpcExit * cp.burnerExit * tK.burnerExit) - (m.hpcExit * cp.hpcExit * tK.hpcExit))/(burner.eta*burner.fhv- cp.burnerExit * tK.burnerExit);
	m.burnerExit = m.hpcExit + burner.mfDot;
	tP.burnerExit = tP.hpcExit*burner.pr;
	tK.burnerExit = ((m.hpcExit*cp.hpcExit*tK.hpcExit) + (burner.eta*burner.fhv*burner.mfDot))/(m.burnerExit*cp.burnerExit);
	work.burner = (burner.eta*burner.fhv*burner.mfDot);
	theta.burnerExit = tK.burnerExit/288.15;
	delta.burnerExit = tP.burnerExit/101325;
	mc.burnerExit = m.burnerExit*sqrt(theta.burnerExit)/delta.burnerExit;
	% HPT
	m.hptExit = m.burnerExit + hpt.fractionOfBleedAirReturningInFrontOfHptAndDoingWork*hpc.bleedOutFractionOfInletAir*m.inletExit;
	tK.hptExit = tK.burnerExit - (work.hpc/hpSpool.eta + hpt.excessWork)/(m.hptExit*cp.burnerExit);
	tKisen.hptExit = tK.burnerExit - (tK.burnerExit - tK.hptExit)/hpt.eta;
	tP.hptExit = tP.burnerExit*(tKisen.hptExit/tK.burnerExit)^(g.burnerExit/(g.burnerExit-1));
	tP.hptExit = tP.hptExit*hpt.exitPr;
	work.hpt = (tK.burnerExit - tK.hptExit)*m.hptExit*cp.burnerExit;
	theta.hptExit = tK.hptExit/288.15;
	delta.hptExit = tP.hptExit/101325;
	mc.hptExit = m.hptExit*sqrt(theta.hptExit)/delta.hptExit;
	[g.hptExit,~,cp.hptExit] = gordonMcBrideAir_withG(tK.hptExit);
	hpt.pr = tP.burnerExit/tP.hptExit;
	% LPT
	m.lptExit = m.hptExit + lpt.fractionOfBleedAirReturningInFrontOfLptAndDoingWork*hpc.bleedOutFractionOfInletAir*m.inletExit;
	tK.lptExit = tK.hptExit - (work.lpc/lpSpool.eta + lpt.excessWork)/(m.lptExit*cp.hptExit);
	tKisen.lptExit = tK.hptExit - (tK.hptExit - tK.lptExit)/lpt.eta;
	tP.lptExit = tP.hptExit*(tKisen.lptExit/tK.hptExit)^(g.hptExit/(g.hptExit-1));
	tP.lptExit = tP.lptExit*lpt.exitPr;
	work.lpt = (tK.hptExit - tK.lptExit)*m.lptExit*cp.hptExit;
	theta.lptExit = tK.lptExit/288.15;
	delta.lptExit = tP.lptExit/101325;
	mc.lptExit = m.lptExit*sqrt(theta.lptExit)/delta.lptExit;
	[g.lptExit,~,cp.lptExit] = gordonMcBrideAir_withG(tK.lptExit);
	lpt.pr = tP.hptExit/tP.lptExit;
	% [sP.lptExit,flag,dev] = staticFromTotal(m.lptExit,tP.lptExit,tK.lptExit, lpt.exitArea, g.lptExit,'LPT exit');
	if ((strcmp(flag,'mDot_not_reached_even4_mach0p99_at_LPT_exit')) || ...
	(strcmp(flag,'mDot_not_reached_even4_mach0p05_at_LPT_exit')))
	fprintf('\n *************************************');
	fprintf('\n ILLOGICAL VALUE OF Ps AT LPT EXIT');
	fprintf('\n *************************************');
	fprintf('\n Could not find sP at LPT exit because %s',flag);
	fprintf('\n Could not find sP at LPT exit because %s',flag);
	fprintf('\n Thrust and static pressure values will be');
	fprintf('\n altered so that the error is very high at this guess.');
	fprintf('\n The alteration will be a function of deviation of sP ');
	fprintf('\n from max or min possible value');
	sP.afterburnerEntry = dev*1e7;
	th.momentum = dev*1e7;
	th.pressure = dev*1e7;
	th.pressure = dev*1e7;
	area.nozzleThroat = 777;
	area.nozzleExit = 777;
	area.geometricNozzleExit = 777;
	return;
	end
	%% Power turbine
	m.fptExit = m.lptExit + fpt.fractionOfBleedAirReturningInFrontOfFptAndDoingWork*hpc.bleedOutFractionOfInletAir*m.inletExit;
	tK.fptExit = tK.lptExit - (propeller.powerTakeOff_kW*1000/fpSpool.eta + gearBox.excessWork)/(m.fptExit*cp.lptExit);
	tKisen.fptExit = tK.lptExit - (tK.lptExit - tK.fptExit)/fpt.eta;
	tP.fptExit = tP.lptExit*(tKisen.fptExit/tK.lptExit)^(g.lptExit/(g.lptExit-1));
	tP.fptExit = tP.fptExit*fpt.exitPr;
	work.lpt = (tK.lptExit - tK.fptExit)*m.fptExit*cp.lptExit;
	theta.fptExit = tK.fptExit/288.15;
	delta.fptExit = tP.fptExit/101325;
	mc.fptExit = m.fptExit*sqrt(theta.fptExit)/delta.fptExit;
	[g.fptExit,~,cp.fptExit,~,R.fptExit] = gordonMcBrideAir_withG(tK.fptExit);
	fpt.pr = tP.lptExit/tP.fptExit;
	% [sP.fptExit,flag,dev] = staticFromTotal(m.fptExit,tP.fptExit,tK.fptExit, lpt.exitArea, g.fptExit,'LPT exit');
	%% Nozzle
	m.nozzleThroat = m.fptExit;
	m.nozzleExit = m.nozzleThroat;
	tP.nozzleEntry = tP.fptExit;
	tK.nozzleEntry = tK.fptExit;
	nozzleEntryThroatPressureRatio = tP.nozzleEntry/pAmb;
	mach.nozzleThroat = sqrt((((nozzleEntryThroatPressureRatio )^((g.fptExit -1)/g.fptExit ))-1)*(2/(g.fptExit -1)));
	if(mach.nozzleThroat >=1)
	mach.nozzleThroat = 1;
	end
	% Get exit temperature
	tK.nozzleThroat = tK.nozzleEntry;
	tP.nozzleThroat = tP.nozzleEntry;
	sK.nozzleThroat = tK.nozzleThroat / (1+(((g.fptExit -1)/2)*(mach.nozzleThroat*mach.nozzleThroat)));
	sP.nozzleThroat = tP.nozzleThroat*(sK.nozzleThroat/tK.nozzleThroat)^(g.fptExit /(g.fptExit - 1));
	a.nozzleThroat = sqrt(g.fptExit*R.fptExit*sK.nozzleThroat);
	v.nozzleThroat = a.nozzleThroat*mach.nozzleThroat;
	rho.nozzleThroat = sP.nozzleThroat/(R.fptExit*sK.nozzleThroat);
	area.nozzleThroat = m.nozzleExit/(rho.nozzleThroat*v.nozzleThroat);
	% Now assume expansion to optimum mach number at the exit of the divergent
	% section
	nozzleThroatExitPressureRatio = tP.nozzleThroat/pAmb;
	mach.nozzleExit = sqrt((((nozzleThroatExitPressureRatio )^((g.fptExit -1)/g.fptExit ))-1)*(2/(g.fptExit -1)));
	tK.nozzleExit = tK.nozzleEntry;
	tP.nozzleExit = tP.nozzleEntry;
	sK.nozzleExit = tK.nozzleExit / (1+(((g.fptExit -1)/2)*(mach.nozzleExit*mach.nozzleExit)));
	sP.nozzleExit = tP.nozzleExit*(sK.nozzleExit/tK.nozzleExit)^(g.fptExit /(g.fptExit - 1));
	a.nozzleExit = sqrt(g.fptExit*R.fptExit*sK.nozzleExit);
	v.nozzleExit = a.nozzleExit*mach.nozzleExit;
	rho.nozzleExit = sP.nozzleExit/(R.fptExit*sK.nozzleExit);
	area.nozzleExit = m.nozzleExit/(rho.nozzleExit*v.nozzleExit);
	area.geometricNozzleExit = area.nozzleExit/nozzle.cd;
	th.pressure = (sP.nozzleExit-pAmb)*area.geometricNozzleExit;
	th.momentum = m.nozzleExit * nozzle.ct*(v.nozzleExit - mach.ambient*sqrt(g.inletExit*287*tK.inletEntry));
	th.total = th.pressure + th.momentum;
	theta.nozzleThroat = tK.nozzleThroat/288.15;
	delta.nozzleThroat = tP.nozzleThroat/101325;
	mc.nozzleThroat = m.nozzleThroat*sqrt(theta.nozzleThroat)/delta.nozzleThroat;
	theta.nozzleExit = tK.nozzleExit/288.15;
	delta.nozzleExit = tP.nozzleExit/101325;
	mc.nozzleExit = m.nozzleExit*sqrt(theta.nozzleExit)/delta.nozzleExit;
	[g.nozzleExit,~,cp.nozzleExit] = gordonMcBrideAir_withG(tK.nozzleExit);
	%% Powers
	lpc.power = m.inletExit*cp.inletExit*(tK.inletExit - tK.lpcExit);
	hpc.power = m.lpcExit*cp.lpcExit*(tK.lpcExit - tK.hpcExit);
	burner.power = work.burner;
	hpt.power = m.burnerExit*cp.burnerExit*(tK.burnerExit - tK.hptExit);
	lpt.power = m.hptExit*cp.hptExit*(tK.hptExit - tK.lptExit);
	fpt.power = m.fptExit*cp.lptExit*(tK.lptExit - tK.fptExit);
	%% Output matrices populate
	fpt_power_out(rowIndex,colIndex) = fpt.power;
	sfc_out(rowIndex,colIndex) = burner.mfDot / fpt_power_out(rowIndex,colIndex) * 3600 * 1000; % Approx value = 250 kg/(kw hr)
	Thspec_out(rowIndex,colIndex) = th.total / lpc.mDot; % Approx value = 250 kg/(kw hr)
	Powspec_out(rowIndex,colIndex) = propeller.powerTakeOff_kW / lpc.mDot; % Approx value = 250 kg/(kw hr)
	mfDot_out(rowIndex,colIndex) = burner.mfDot;
	end
	end
	%% Print key results
	if ( size(prGrid,1) == 1)
	list_keyPerformance = [sfc_out;
	tK.lptExit;
	th.total;
	lpc.power;
	hpc.power;
	burner.mfDot];
	display('SFC ITT Thrust LPCpower HPCpower FuelFlowRate');
	fprintf('\n SFC = %f kg/(kW*hr)',list_keyPerformance(1));
	fprintf('\n SFC = %f lb/(hp*hr)',list_keyPerformance(1)/0.746*0.4535);
	fprintf('\n ITT = %f K',list_keyPerformance(2));
	fprintf('\n Th = %f N',list_keyPerformance(3));
	fprintf('\n LPCpower = %f MW',list_keyPerformance(4)/1e6);
	fprintf('\n HPCpower = %f MW',list_keyPerformance(5)/1e6);
	fprintf('\n Mf = %f kg/s',list_keyPerformance(6));
	%% Plot results of parametric search:
	close all;
	% Station wise Plot
	stationNamesLabels={'Inlet Entry','Inlet Exit','LPC Exit','HPC Exit','Combustor Exit','HPT Exit','LPT Exit','FPT Exit','Nozzle Exit'};
	yQuantities={'Mass flow rate $$\frac{kg}{sec}$$','Corrected mass flow rate $$\frac{kg}{sec}$$','Temperature total,( Kelvin)','Pressure Total (Pascals)'};
	yQs = [m.inletEntry m.inletExit m.lpcExit m.hpcExit m.burnerExit m.hptExit m.lptExit m.fptExit m.nozzleExit;
	mc.inletEntry mc.inletExit mc.lpcExit mc.hpcExit mc.burnerExit mc.hptExit mc.lptExit mc.fptExit mc.nozzleExit;
	tK.inletEntry tK.inletExit tK.lpcExit tK.hpcExit tK.burnerExit tK.hptExit tK.lptExit tK.fptExit tK.nozzleExit;
	tP.inletEntry tP.inletExit tP.lpcExit tP.hpcExit tP.burnerExit tP.hptExit tP.lptExit tP.fptExit tP.nozzleExit;];
	titles=char('Mass flow rate','Corrected mass flow rate','Total Temperature ','Total Pressure ');
	fn = char('Mass','CorrMass','Tt','Pt');
	for i=[1 3 4]
	figure('Name',fn(i,:),'units','normalized','position',[0 0 1 1]);
	plot( [1 2 3 4 5 6 7 8 9 ],yQs(i,:),'-ok','linewidth',3);
	set(gca, 'XTick', 1:9, 'XTickLabel', stationNamesLabels);
	set(gca,'fontsize',14);
	xlabel('Station Name','fontsize',14,'fontweight','b','color','k');
	xticklabel_rotate([],45,[],'fontsize',14,'fontweight','b','color','k');
	ylabel(yQuantities(i),'interpreter','latex','fontsize',14,'fontweight','b','color','k');
	fileName=strcat(titles(i,:) ,' variation along the engine');
	title(fileName,'fontsize',16,'fontweight','b','color','k');
	fileName=strcat('.\..\output\pngFiles\',fn(i,:) ,'-var');
	export_fig(strcat(fileName ,'.png'));
	fileName=strcat('.\..\output\figFiles\',fn(i,:) ,'-var');
	saveas(gcf,strcat(fileName ,'.fig'));
	%saveas(gcf,strcat('E:\01Balajis_dd_Work\01_PR\01_genratedDocs\03_sgt\sgtAnalysisReport_pd\', fn(i,:), '-variationStationWise.png'));
	end
	end
	%% Use only if you are doing parametric plot
	if ( size(prGrid,1)>1)
	figure('units','normalized','outerposition',[0 0 1 1]);
	% subplot(131)
	contourf(prGrid,mDotGrid,sfc_out);
	xlabel('Pressure ratio');
	ylabel('Mass flow rate kg/s');
	zlabel('SFC kg/(kW*hr)');
	title (sprintf(' SFC $$\\frac{lb}{shp \\times hr}$$ = f($$ \\dot{m}$$ , $$\\pi$$) , %d HP power',targetHp),'interpreter','latex','fontsize',20,'fontweight','n','color','k');
	colorbar
	set(findall(gcf,'-Property','fontsize'),'fontsize',36)
	end
	%%
	% subplot(132)
	% contourf(prGrid,mDotGrid,Thspec_out);
	% xlabel('Pressure ratio');
	% ylabel('Mass flow rate kg/s');
	% zlabel('Specifice power kW/(kg/s)');
	% title (sprintf('Specific thrust with $$ \\dot{m}$$ and $$\\pi$$ \n %d HP power',targetHp),'interpreter','latex','fontsize',20,'fontweight','n','color','k');
	% colorbar
	% set(findall(gcf,'-Property','fontsize'),'fontsize',36)
	%
	% subplot(133)
	% contourf(prGrid,mDotGrid,Thspec_out./max(max(Thspec_out))+sfc_out./max(max(sfc_out)));
	% xlabel('Pressure ratio');
	% ylabel('Mass flow rate kg/s');
	% zlabel('Specifice power kW/(kg/s)');
	% title (sprintf('Specific thrust with $$ \\dot{m}$$ and $$\\pi$$ \n %d HP power',targetHp),'interpreter','latex','fontsize',20,'fontweight','n','color','k');
	% colorbar
	% set(findall(gcf,'-Property','fontsize'),'fontsize',36)
	% figure('units','normalized','outerposition',[0 0 1 1]);
	% contourf(prGrid,mDotGrid,mfDot_out);
	% xlabel('Pressure ratio');
	% ylabel('Mass flow rate kg/s');
	% zlabel('Fuel flow rate, kg/s');
	% title (sprintf('Variation of fuel flow rate with Mass flow rate and pressure ratio\n 2500 HP propeller useful power'),'fontsize',20,'fontweight','n','color','k');
	% colorbar
	% set(findall(gcf,'-Property','fontsize'),'fontsize',36)
	%
	%
	% figure('units','normalized','outerposition',[0 0 1 1]);
	% contourf(prGrid,mDotGrid,fpt_power_out/746);
	% xlabel('Pressure ratio');
	% ylabel('Mass flow rate kg/s');
	% zlabel('Propeller power, HP');
	% title (sprintf('Variation of fuel flow rate with Mass flow rate and pressure ratio\n 2500 HP propeller useful power'),'fontsize',20,'fontweight','n','color','k');
	% colorbar
	% set(findall(gcf,'-Property','fontsize'),'fontsize',16)
	%

view raw main_dpSim_2500HP_2750HP_3000HP_parametric.m hosted with ❤ by GitHub