function [Q, qw, v_vc0, Tw, h_re, delta_t_tot, Mach, x, Istr_TPS] = function_TPS(W_dry, Spln,L_D_re, cD_re, R, phi)

 
x = 0:1:84; 
beta = W_dry/(cD_re*Spln);
sigma = 5.670367*10^-8;
eps = 0.85;

g0 = 9.81; gamma = 1.35;
h_re = 84000:-1000:1;
[T, a, p, rho, nu] = atmosisa(h_re,extended=true);
nu_e = nu;
r0 = 6371000;
V = sqrt(g0*(r0+h_re)./(1+(rho.*(r0+h_re)*L_D_re)/(2*beta)));
Mach = V./a;
u_e = V;
T0 = T.*(1+(gamma-1)*0.5*Mach.^2);

z = 0; b = 0;
for i = 1:1:length(h_re)
    for j = 1:1:length(x)
% transition above Rex = 3*10^6, fully turbolent achieved Rex = 7*10^6;
Rex(i,j) = u_e(i)*x(j)/nu_e(i); % edge boundary layer = free stream? 
ReT(i) = 10^(6.421*exp(1.209*10^-4*Mach(i)^2.641));
k = 0;
Tw(i,j) = T(i);
while 1

xT = 0;
if Rex(i,j) == 0
%stagnation point
Ms = 3; Ns = 0.5; Cs(i) = (1.83*10^-8)*(R^-0.5)*(1-Tw(i,j)/T0(i));
qw(i,j) = rho(i)^Ns*V(i)^Ms*Cs(i)*10^4; %approx of eq 6.106 
end

if Rex(i,j) < ReT(i) && Rex(i,j) ~= 0
% Laminar flat plate
M = 3.2; N = 0.5; C(i,j) = (2.53*10^-9)*cos(phi)^0.5*sin(phi)*(x(j)^-0.5)*(1-Tw(i,j)/T0(i));
qw(i,j) = rho(i)^N*V(i)^M*C(i,j)*10^4;
xT = Rex(i,j)*nu_e(i)/(u_e(i));
end

if Rex(i,j) >= ReT(i)
    if xT <= 1
        xT = 0;
    end
% Turbolent flat plate 
if V < 3962
    Mt = 3.37;
    Ct(i,j) = (3.89*10^-8)*cos(phi)^1.78*sin(phi)^1.6*(x(j)-xT)^(-1/5)*(Tw(i,j)/556)^(-1/4)*(1-1.11*Tw(i,j)/T0(i)); 
else
    Mt = 3.7;
    Ct(i,j) = (2.2*10^-9)*cos(phi)^2.08*sin(phi)^1.6*(x(j)-xT)^(-1/5)*(1-1.11*Tw(i,j)/T0(i));
end
qw(i,j) = rho(i)^N*V(i)^Mt*Ct(i,j)*10^4;
end

% thermal energy radiated
qw_out(i,j) = sigma*eps*Tw(i,j)^4;

if abs(qw(i,j)-qw_out(i,j)) < 305
    break
end
Tw(i,j) = Tw(i,j) + 0.25;
k = k+1;
if k > 8001
    b = b+1;
    n(b) = i;
    break
end
end
    end
% flight time
v_vc0(i) = sqrt(1/(1+(rho(i)*r0*L_D_re)/(2*beta)));
t_tot(i) = 0.5*sqrt(r0/g0)*L_D_re*log((1+v_vc0(i)^2)/(1-(v_vc0(i))^2));
delta_t_tot(i) = t_tot(1) - t_tot(i);
if i > 1
    z = z+1;
    delta_t(z) = t_tot(i-1) - t_tot(i);
end
end

% heat load
qw_m = zeros(length(h_re)-1,length(x));
for j = 1:1:length(x)
    for i = 1:1:length(h_re)-1
       qw_m(i,j) = (qw(i,j) + qw(i+1,j))/2;
    end
end

for j = 1:1:length(x)
Q(j) = sum(qw_m(:,j)'.*delta_t);
end
Q_TPS = Q/(3.281^2*1055);

% Conversion Unit
% 1 [m] = 39.3701 [in]; 1 [J] = [Btu]/1055.056; 1 [m] = 3.2808 [ft]; 1 [kg]=  2.2046 [lb]; 
% T(°F) = ((T(K) - 273.15) * 1.8) + 32; T(K) = (T(°F) + 459,67) * 5/9;

Q_prova_HT = 3900:1:32000;
Q_prova_LW = 0:1:5500;

Q_prova_HT = 3900:1:32000;
Q_prova_LW = 0:1:5500;

TABI_highT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\TABI_highT.mat");
[fitresult1, gof1] = createFit3(TABI_highT(:,1)',TABI_highT(:,2)');
TABI_hT = fitresult1(Q_prova_HT)';
AETB12TUFI_highT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\AETB12TUFI_highT.mat");
[fitresult2, gof2] = createFit3(AETB12TUFI_highT(:,1)',AETB12TUFI_highT(:,2)');
AETB12TUFI_hT = fitresult2(Q_prova_HT)';
AMHC_highT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\AMHC_highT.mat");
[fitresult3, gof3] = createFit3(AMHC_highT(:,1)',AMHC_highT(:,2)');
AMHC_hT = fitresult3(Q_prova_HT)';
SAHC_highT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\SAHC_highT.mat");
[fitresult4, gof4] = createFit3(SAHC_highT(:,1)',SAHC_highT(:,2)');
SAHC_hT = fitresult4(Q_prova_HT)';
SAHC2_highT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\SAHC2_highT.mat");
[fitresult5, gof5] = createFit3(SAHC2_highT(:,1)',SAHC2_highT(:,2)');
SAHC2_hT = fitresult5(Q_prova_HT)';
LI900_highT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\LI900_highT.mat");
[fitresult6, gof6] = createFit3(LI900_highT(:,1)',LI900_highT(:,2)');
LI900_hT = fitresult6(Q_prova_HT)';
AETB8TUFI_highT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\AETB8TUFI_highT.mat");
[fitresult7, gof7] = createFit3(AETB8TUFI_highT(:,1)',AETB8TUFI_highT(:,2)');
AETB8TUFI_hT = fitresult7(Q_prova_HT)';


AFRSI_lowT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\AFRSI_lowT.mat");
[fitresult8, gof8] = createFit3(AFRSI_lowT(:,1)',AFRSI_lowT(:,2)');
AFRSI_lT = fitresult8(Q_prova_LW)';
TABI_lowT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\TABI_lowT.mat");
[fitresult9, gof9] = createFit3(TABI_lowT(:,1)',TABI_lowT(:,2)');
TABI_lT = fitresult9(Q_prova_LW)';
TIHC_lowT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\TIHC_lowT.mat");
[fitresult10, gof10] = createFit3(TIHC_lowT(:,1)',TIHC_lowT(:,2)');
TIHC_lT = fitresult10(Q_prova_LW)';
TIMW_lowT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\TIMW_lowT.mat");
[fitresult11, gof11] = createFit3(TIMW_lowT(:,1)',TIMW_lowT(:,2)');
TIMW_lT = fitresult11(Q_prova_LW)';
LI900_lowT = importdata("C:\Users\Tommaso\OneDrive\Documenti\MATLAB\LI900_lowT.mat");
[fitresult12, gof12] = createFit3(LI900_lowT(:,1)',LI900_lowT(:,2)');
LI900_lT = fitresult12(Q_prova_LW)';

Istr_tps = zeros*length(Q_TPS);
for j = 1:1:length(Q_TPS)
    i = 1; 
    while abs(Q_prova_LW(i) - Q_TPS(j)) > 1
        Istr_tps(j) = min([LI900_lT(i), TIMW_lT(i), TIHC_lT(i), TABI_lT(i), AFRSI_lT(i)]);
        i = i+1;
        if i == length(Q_prova_LW)
            Istr_tps(j) = 0;
            break
        end
    end
end

for j = 1:1:length(Q_TPS)
    i = 1; 
    if Istr_tps(j) == 0
    while abs(Q_prova_HT(i) - Q_TPS(j)) > 1
        Istr_tps(j) = min([TABI_hT(i), AETB12TUFI_hT(i), AETB8TUFI_hT(i), AMHC_hT(i), SAHC_hT(i), SAHC2_hT(i), LI900_hT(i)]);
        i = i+1;
        if i == length(Q_prova_HT)
            Istr_tps(j) = 5;
            break
        end
    end
    end
end

Istr = Istr_tps*3.2808^2/2.2046; Istr_TPS = sum(Istr(1:length(Istr)))/(length(Istr));


   

