function [wf_to_matrix,wf_c_matrix,wf_cr_matrix_A,wf_cr_matrix_B,wf_cr_matrix_C,wf_cr_matrix_D,...
      h_to_matrix,h_c_matrix,h_cr_matrix,v_to_matrix,v_c_matrix,v_cr_matrix,...
   	P_to_in,P_c_in,P_cr_in_A,P_cr_in_B,P_cr_in_C,P_cr_in_D]=prop_inst_funct(P_to_un,W_to)

% P_to_un is the uninstalled power for one engine

% atmospheric data
rho_o=2.377*10^-3;
mu_o=3.737*10^-7;
To=518.69;
gamma=1.405;
g=32.2;
R=1716;
at=0.00357;

% inputs 

%W_to = 60000;			% assume a constant takeoff weight
V_max = 400; 			% true airspeed (kts)
m_flow_eng=rho_o*(173/144)*(V_max*1.6878);	% need in slugs/sec -- Ae=173in^2, assume Vmax=400kts
P_el =0.0007*W_to;	% (shp) for twin engine turbofan Table 6.1 pg 144
P_mech=0.0006*W_to;	% (shp) for twin engine turbofan Table 6.1 pg 144
ninl=.85;				% assumption for ducted spinner -- pg 172
P_to_engine = 10000/4;	% Turboprop uninstalled thrust
D_p = 14.5;		% prop diameter in ft

% input AIAA turboprop data
% takeoff condition
   h_to_matrix = [0 5000 10000; 0 5000 10000; 0 5000 10000; 0 5000 10000];
	v_to_matrix = [0 0 0; 100 100 100; 200 200 200; 300 300 300];  
	wf_to_matrix = [1437 1342 1236; 1440 1339 1247; 1406 1313 1226; 1342 1257 1194];
	p_to_matrix = [2500 2500 2468; 2500 2500 2500; 2500 2500 2500; 2500 2500 2500];
% max climb condition
	h_c_matrix=[0 5000 10000 15000 20000 25000 30000 40000; 0 5000 10000 15000 20000 25000 30000 40000;...
      0 5000 10000 15000 20000 25000 30000 40000;0 5000 10000 15000 20000 25000 30000 40000;...
      0 5000 10000 15000 20000 25000 30000 40000];
	v_c_matrix=[0 0 0 0 0 0 0 0; 100 100 100 100 100 100 100 100; 200 200 200 200 200 200 200 200;...
      300 300 300 300 300 300 300 300; 400 400 400 400 400 400 400 400];
	wf_c_matrix=[1437 1342 1170 1016 873 739 572 276.5; 1440 1339 1191 1035 889 754 588 292.2; ...
      1406 1313 1226 1072 924 789 630 340; 1342 1257 1194 1133 987 842 699 410.7; ...
      1239 1180 1160 1206 1075 907 790 496.5];
	p_c_matrix= [2250 2250 2250 2074 1821 1565 1212 619; 2250 2250 2250 2126 1866 1609 1256 662;...
      2250 2250 2250 2250 1994 1726 1388 793.7; 2250 2250 2250 2250 2223 1926 1621 1020.3; ...
      2250 2250 2250 2250 2250 2202 1941 1667.5];
% cruise condition
	h_cr_matrix = [20000 25000 30000 40000;20000 25000 30000 40000;20000 25000 30000 40000;...
      20000 25000 30000 40000;20000 25000 30000 40000]; 
	v_cr_matrix = [0 0 0 0; 100 100 100 100; 200 200 200 200; 300 300 300 300; 400 400 400 400];
	wf_cr_matrix_A = [778 664 548 318.3; 791 678 561 331.7; 823 704 590 356.2; 868 752 638 407.7; 930 806 707 479.8];
	p_cr_matrix_A = [1585 1382 1147 714.3; 1626 1420 1182 743.3; 1744 1524 1285 829.2; 1926 1697 1459 993.5; ...
      2177 1940 1706 1234.5];
	wf_cr_matrix_B = [683 593 492 302.8; 696 604 503 311.5; 717 627 527 334.7; 752 662 572 392; 788 709 624 461];
	p_cr_matrix_B = [1359 1218 1021 692.3; 1397 1250 1053 757.3; 1497 1347 1150 810.8; 1659 1494 1306 956.8; ...
      1883 1706 1532 1180.5];
	wf_cr_matrix_C = [588 516 434 281.67; 601 527 442 284.83; 617 545 463 310.67; 638 569 498 358.33; 662 595 550 434.33];
	p_cr_matrix_C = [1112 1026 876 650.7; 1147 1056 906 674.8; 1238 1141 994 758.3; 1374 1271 1135 901.5; ...
      1559 1453 1336 1114.8];
	wf_cr_matrix_D = [466 421 354 245.67; 474 429 365 259.17; 490 442 379 270.5; 508 458 402 297; 524 478 429 334.5];
	p_cr_matrix_D = [768 764 671 588.83; 794 791 697 615.2; 879 865 776 685.5; 1012 976 897 789.2; 1205 1131 1072 936.5];
   
% input efficiency map
J_eff = [.2 .2 .2 .2 .2 .2 .2 .2 .2; .4 .4 .4 .4 .4 .4 .4 .4 .4; .6 .6 .6 .6 .6 .6 .6 .6 .6;...
      .8 .8 .8 .8 .8 .8 .8 .8 .8; 1 1 1 1 1 1 1 1 1; 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2;...
      1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4; 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6;...
      1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8; 2 2 2 2 2 2 2 2 2 ; 2.2 2.2 2.2 2.2 2.2 2.2 2.2 2.2 2.2;...
      2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8];
CP_eff = [.02 .04 .08 .16 .2 .24 .28 .32 .4;.02 .04 .08 .16 .2 .24 .28 .32 .4;.02 .04 .08 .16 .2 .24 .28 .32 .4;...
      .02 .04 .08 .16 .2 .24 .28 .32 .4; .02 .04 .08 .16 .2 .24 .28 .32 .4; .02 .04 .08 .16 .2 .24 .28 .32 .4;...
      .02 .04 .08 .16 .2 .24 .28 .32 .4; .02 .04 .08 .16 .2 .24 .28 .32 .4; .02 .04 .08 .16 .2 .24 .28 .32 .4;...
   	.02 .04 .08 .16 .2 .24 .28 .32 .4; .02 .04 .08 .16 .2 .24 .28 .32 .4; .02 .04 .08 .16 .2 .24 .28 .32 .4];
nu_eff = [0 .46 .39 .25 0 0 0 0 0; 0 .71 .64 .47 .41 .34 .28 0 0; .63 .81 .77 .65 .58 .53 .45 0 0;...
      .6 .82 .84 .76 .71 .65 .61 0 0; .6 .81 .875 .823 .785 .75 .705 .67 0;...
      .6 .75 .884 .86 .835 .805 .775 .74 .67; .58 .73 .875 .88 .864 .845 .82 .795 .74;...
      .45 .675 .857 .891 .881 .869 .851 .83 .78; .45 .6 .83 .894 .891 .882 .871 .858 .81;...
      .4 .6 .8 .89 .895 .81 .885 .875 .85; 0 0 .755 .882 .892 .8945 .891 .885 .87;...
      0 .3 .525 .82 .865 .883 .891 .894 .892];

% input figure II	% only positive blade angles -- correspond to positive thrust
J_II=[0 0 0 0 0 0; .2 .2 .2 .2 .2 .2; .4 .4 .4 .4 .4 .4; .6 .6 .6 .6 .6 .6; .8 .8 .8 .8 .8 .8; 1 1 1 1 1 1];
CP_II=[0 0.04 0.08 0.16 0.24 0.4; 0 0.04 0.08 0.16 0.24 0.4; 0 0.04 0.08 0.16 0.24 0.4;...
   0 0.04 0.08 0.16 0.24 0.4; 0 0.04 0.08 0.16 0.24 0.4; 0 0.04 0.08 0.16 0.24 0.4];
bld_angle_II=[-3 7.5 13 20 25 25; -1 8.5 14.5 21 26 26; -1 9.5 15.5 22.5 28 28;...
   -1 12.5 17.5 23.5 29 29; 8.5 16.5 20 25.5 30.5 30.5; 15.5 19.5 23 27.5 32 32];

% input figure III % only positive blade angles
J_III=[0 0 0 0; .2 .2 .2 .2; .4 .4 .4 .4; .6 .6 .6 .6; .8 .8 .8 .8; 1 1 1 1];
bld_angle_III=[0 10 20 30; 0 10 20 30; 0 10 20 30; 0 10 20 30; 0 10 20 30; 0 10 20 30];
CT_III=[3.5 .145 .21 .21; 1.5 .115 .18 .21; 0 .07 .135 .21; -.04 -.025 .035 .21; -.14 -.1 .09 .19;...
      -.24 -.105 .035 .16];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% change uninstalled takeoff data to installed data
N = 1200; % (rpm) propeller speed at takeoff
for i = 1:4,  % loop through airspeeds 
   v = (i-1) * 100; % true airspeed in kts
   for j = 1:3,  % loop through altitudes
      h = (j-1) * 5000; % altitude in ft
      % calculate new atmospheric properties
      rho = rho_o * (1-at*h/To)^(g/(at*R)-1);
      T = To - at*h;
      a = sqrt(gamma*R*T);
      M = v * 1.6878 / a;
      % find available installed power -- Roskam 'Airplane Design Part VI' 6.4.1.2 pg 195
      m_flow_bleed=0.03*m_flow_eng;					%(slugs/sec)
      P_pneum=(m_flow_bleed/m_flow_eng)*p_to_matrix(i,j);	% (shp) for twin engine turbofan Table 6.1 pg 144
      P_extr = P_el + P_mech + P_pneum;			% electrical, mechanical and pneumatic extraction 6.1.2 pg 145
      J = 101.4 * (v/(N*D_p)); 						% 1200 is propeller speed at takeoff
      CP = p_to_matrix(i,j) * (rho_o/rho) / (2000 * (N^3 * 10^-9) * D_p^5 * 10^-5);
      if v > 50
         np=INTERP2(CP_eff,J_eff,nu_eff,CP,J);		% propeller efficiency
         P_to_in(i,j) = (ninl*p_to_matrix(i,j)-P_extr)*np;
      else
         bld_angle=interp2(CP_II,J_II,bld_angle_II,CP,J);
         CT=interp2(bld_angle_III,J_III,CT_III,bld_angle,J);
         T=29500*(N^2 * 10^-6)*CT/ (rho_o/rho);
         P_to_in(i,j) = T*v*1.6878/550;
      end	% ifelse
   end	% inner loop
end	% outer loop

% change uninstalled climb data to installed data
h_climb=[0 5000 10000 15000 20000 25000 30000 40000];
N = 1100; % (rpm) propeller speed at climb
for i = 1:5, % loop through airspeeds 
   v = (i-1) * 100; % true airspeed in kts
   for j = 1:8,  % loop through altitudes
      h = h_climb(j); % altitude in ft
      % calculate new atmospheric properties
      rho = rho_o * (1-at*h/To)^(g/(at*R)-1);
      T = To - at*h;
      a = sqrt(gamma*R*T);
      M = v * 1.6878 / a;
      % find available installed power -- Roskam 'Airplane Design Part VI' 6.4.1.2 pg 195
      m_flow_bleed=0.03*m_flow_eng;	%(slugs/sec)
      P_pneum=(m_flow_bleed/m_flow_eng)*p_c_matrix(i,j);	% (shp) for twin engine turbofan Table 6.1 pg 144
      P_extr = P_el + P_mech + P_pneum; 			% electrical, mechanical and pneumatic extraction 6.1.2 pg 145
      J = 101.4 * (v/(N*D_p));  						% 1200 is propeller speed at takeoff
      CP = p_c_matrix(i,j) * (rho_o/rho)/ (2000 * (N^3 * 10^-9) * D_p^5 * 10^-5);
      if v > 50
         np=INTERP2(CP_eff,J_eff,nu_eff,CP,J);		% propeller efficiency
         P_c_in(i,j) = (ninl*p_c_matrix(i,j)-P_extr)*np;
      else
         bld_angle=interp2(CP_II,J_II,bld_angle_II,CP,J);
         CT=interp2(bld_angle_III,J_III,CT_III,bld_angle,J);
         T=29500*(N^2 * 10^-6)*CT/ (rho_o/rho);
         P_c_in(i,j) = T*v*1.6878/550;
      end	% ifelse
   end
end

% change uninstalled cruise data A to installed data
h_cruise=[20000 25000 30000 40000];
N = 1000; % (rpm) propeller speed at cruise
for i = 1:5, % loop through airspeeds 
   v = (i-1) * 100; % true airspeed in kts
   for j = 1:4,  % loop through altitudes
      h = h_cruise(j); % altitude in ft
      % calculate new atmospheric properties
      rho = rho_o * (1-at*h/To)^(g/(at*R)-1);
      T = To - at*h;
      a = sqrt(gamma*R*T);
      M = v * 1.6878 / a;
      % find available installed power -- Roskam 'Airplane Design Part VI' 6.4.1.2 pg 195
      m_flow_bleed=0.03*m_flow_eng;	%(slugs/sec)
      P_pneum=(m_flow_bleed/m_flow_eng)*p_cr_matrix_A(i,j);	% (shp) for twin engine turbofan Table 6.1 pg 144
      P_extr = P_el + P_mech + P_pneum; % electrical, mechanical and pneumatic extraction 6.1.2 pg 145
      J = 101.4 * (v/(N*D_p));  % 1200 is propeller speed at takeoff
      CP = p_cr_matrix_A(i,j) * (rho_o/rho) / (2000 * (N^3 * 10^-9) * D_p^5 * 10^-5);   
      if v > 50
         np=INTERP2(CP_eff,J_eff,nu_eff,CP,J);		% propeller efficiency
         P_cr_in_A(i,j) = (ninl*p_cr_matrix_A(i,j)-P_extr)*np;
      else
         bld_angle=interp2(CP_II,J_II,bld_angle_II,CP,J);
         CT=interp2(bld_angle_III,J_III,CT_III,bld_angle,J);
         T=29500*(N^2 * 10^-6)*CT/ (rho_o/rho);
         P_cr_in_A(i,j) = T*v*1.6878/550;
      end	% ifelse
   end
end

% change uninstalled cruise data B to installed data
N = 1000; % (rpm) propeller speed at cruise
for i = 1:5, % loop through airspeeds 
   v = (i-1) * 100; % true airspeed in kts
   for j = 1:4,  % loop through altitudes
      h = h_cruise(j); % altitude in ft
      % calculate new atmospheric properties
      rho = rho_o * (1-at*h/To)^(g/(at*R)-1);
      T = To - at*h;
      a = sqrt(gamma*R*T);
      M = v * 1.6878 / a;
      % find available installed power -- Roskam 'Airplane Design Part VI' 6.4.1.2 pg 195
      m_flow_bleed=0.03*m_flow_eng;	%(slugs/sec)
      P_pneum=(m_flow_bleed/m_flow_eng)*p_cr_matrix_B(i,j);	% (shp) for twin engine turbofan Table 6.1 pg 144
      P_extr = P_el + P_mech + P_pneum; % electrical, mechanical and pneumatic extraction 6.1.2 pg 145
      J = 101.4 * (v/(N*D_p));  % 1200 is propeller speed at takeoff
      CP = p_cr_matrix_B(i,j) * (rho_o/rho) / (2000 * (N^3 * 10e-9) * D_p^5 * 10^-5);
      if v > 50
         np=INTERP2(CP_eff,J_eff,nu_eff,CP,J);		% propeller efficiency
         P_cr_in_B(i,j) = (ninl*p_cr_matrix_B(i,j)-P_extr)*np;
      else
         bld_angle=interp2(CP_II,J_II,bld_angle_II,CP,J);
         CT=interp2(bld_angle_III,J_III,CT_III,bld_angle,J);
         T=29500*(N^2 * 10^-6)*CT/ (rho_o/rho);
         P_cr_in_B(i,j) = T*v*1.6878/550;
      end	% ifelse
   end
end

% change uninstalled cruise data C to installed data
N = 1000; % (rpm) propeller speed at cruise
for i = 1:5, % loop through airspeeds 
   v = (i-1) * 100; % true airspeed in kts
   for j = 1:4,  % loop through altitudes
      h = h_cruise(j); % altitude in ft
      % calculate new atmospheric properties
      rho = rho_o * (1-at*h/To)^(g/(at*R)-1);
      T = To - at*h;
      a = sqrt(gamma*R*T);
      M = v * 1.6878 / a;
      % find available installed power -- Roskam 'Airplane Design Part VI' 6.4.1.2 pg 195
      m_flow_bleed=0.03*m_flow_eng;	%(slugs/sec)
      P_pneum=(m_flow_bleed/m_flow_eng)*p_cr_matrix_C(i,j);	% (shp) for twin engine turbofan Table 6.1 pg 144
      P_extr = P_el + P_mech + P_pneum; % electrical, mechanical and pneumatic extraction 6.1.2 pg 145
      J = 101.4 * (v/(N*D_p));  % 1200 is propeller speed at takeoff
      CP = p_cr_matrix_C(i,j) * (rho_o/rho) / (2000 * (N^3 * 10^-9) * D_p^5 * 10^-5);
      if v > 50
         np=INTERP2(CP_eff,J_eff,nu_eff,CP,J);		% propeller efficiency
         P_cr_in_C(i,j) = (ninl*p_cr_matrix_C(i,j)-P_extr)*np;
      else
         bld_angle=interp2(CP_II,J_II,bld_angle_II,CP,J);
         CT=interp2(bld_angle_III,J_III,CT_III,bld_angle,J);
         T=29500*(N^2 * 10^-6)*CT/ (rho_o/rho);
         P_cr_in_C(i,j) = T*v*1.6878/550;
      end	% ifelse
   end
end

% change uninstalled cruise data D to installed data
N = 1000; % (rpm) propeller speed at cruise
for i = 1:5, % loop through airspeeds 
   v = (i-1) * 100; % true airspeed in kts
   for j = 1:4,  % loop through altitudes
      h = h_cruise(j); % altitude in ft
      % calculate new atmospheric properties
      rho = rho_o * (1-at*h/To)^(g/(at*R)-1);
      T = To - at*h;
      a = sqrt(gamma*R*T);
      M = v * 1.6878 / a;
      % find available installed power -- Roskam 'Airplane Design Part VI' 6.4.1.2 pg 195
      m_flow_bleed=0.03*m_flow_eng;	%(slugs/sec)
      P_pneum=(m_flow_bleed/m_flow_eng)*p_cr_matrix_D(i,j);	% (shp) for twin engine turbofan Table 6.1 pg 144
      P_extr = P_el + P_mech + P_pneum; % electrical, mechanical and pneumatic extraction 6.1.2 pg 145
      J = 101.4 * (v/(N*D_p));  % 1200 is propeller speed at takeoff
      CP = p_cr_matrix_D(i,j) * (rho_o/rho) / (2000 * (N^3 * 10^-9) * D_p^5 * 10^-5);
      if v > 50
         np=INTERP2(CP_eff,J_eff,nu_eff,CP,J);		% propeller efficiency
         P_cr_in_D(i,j) = (ninl*p_cr_matrix_D(i,j)-P_extr)*np;
      else
         bld_angle=interp2(CP_II,J_II,bld_angle_II,CP,J);
         CT=interp2(bld_angle_III,J_III,CT_III,bld_angle,J);
         T=29500*(N^2 * 10^-6)*CT/ (rho_o/rho);
         P_cr_in_D(i,j) = T*v*1.6878/550;
      end	% ifelse
   end
end

% scale 
ESF = P_to_un / P_to_engine;

%%%%%% ONCE ESF is found, mass flow and power arrays must be scaled %%%%%%%
% takeoff condition 
wf_to_matrix = wf_to_matrix * ESF;
P_to_in = P_to_in * ESF;
% max climb condition
wf_c_matrix = wf_c_matrix * ESF;
P_c_in = P_c_in * ESF;
% cruise condition
wf_cr_matrix_A = wf_cr_matrix_A * ESF;
P_cr_in_A = P_cr_in_A * ESF;
wf_cr_matrix_B = wf_cr_matrix_B * ESF;
P_cr_in_B = P_cr_in_B * ESF;
wf_cr_matrix_C = wf_cr_matrix_C * ESF;
P_cr_in_C = P_cr_in_C * ESF;
wf_cr_matrix_D = wf_cr_matrix_D * ESF;
P_cr_in_D = P_cr_in_D * ESF;
