%% FINAL TOOL

clear all
clc

%% User Inputs
Wpay = 15000;
rho_pay = 50;
Ncrw = 0;
ICI = 37.7; 
h_targetorbit = 100000;
Mtr = 5; 
htr = 28500; %m
% Engine selection. Please digit 1 for the selected engine and 0 for the others
SABRE = 1; % Database
turboramjet = 0; %Analytical Method for Turboramjet+Rocket

%% Assumptions
% Atmosphere ISA
            h = 0:100:51000; % m 
            R = 287.05; 
            [T, a, p, rho] = atmosisa(h,extended=true);

% Propellant LOX/LH2
rho_f = 70.52/1000; %70.52/1000; %LH2 [ton/m^3]
hpr = 119954; % heat of reaction, H2 [kJ/kg]
rho_ox = 1141/1000; %LOX [ton/m^3]
MR = 6; % mixure ratio %6

% Value from statistics for the input payload to LEO orbit
GTOW_stat = (15.677775549*(Wpay/1000)^1.1229186457)*1000; %Kg
GTOW = GTOW_stat;
Wdry_stat = (0.867751184852*(GTOW_stat/1000)^0.75533773897032)*1000; %Kg
Spln_stat = 1.640205138*(GTOW_stat/1000) -88.1486838; %m^2;
T_to_stat = 1.158721737*10^-6*(GTOW_stat/1000)^3 - 0.001988945288194*(GTOW_stat/1000)^2 + 7.478778432789740*(GTOW_stat/1000) + 56.7545873127545; %kN
%T_vacuum_stat = 23.8854941423234*(GTOW_stat/1000)^0.837718399300137; %kN;
T_vacuum_stat = 10.943743815826*(GTOW_stat/1000) + 37.5226918895443; %kN
T_ramjet_stat = (5.582098252913430*GTOW_stat + 7382.6467029)/1000; %kN
TW0 = T_to_stat*1000/(GTOW_stat*9.81); % Take off Thrust-Weight 
TvW0 = T_vacuum_stat*1000/(GTOW_stat*9.81); % T vacuum / GTOW
TrW0 = T_ramjet_stat*1000/(GTOW_stat*9.81); % T high speed / GTOW
WR_stat = GTOW_stat/Wdry_stat;


%% Mission Analysis

g = 9.81; 
gamma = 1.35; %adiabatic index

% Airbreathing ascent phase
% equations from Billing, "Design and Development of SSTO Vehicles"
V(1) = 500*0.305; %Take-off Velocity [m/s]
V(2:74) = sqrt(h(2:1:74)/(0.305*2.035*10^-2) + (V(1)/0.305)^2)*0.305; % [m/s] 
q = 0.5*rho(htr/100)*(Mtr*a(htr/100))^2; % constant dynamic pressure [Pa] 
V(200) =  sqrt(2*q./rho(200));
x = [h(1) h(74) h(200)]; y = [V(1) V(74) V(200)];
V(75:199) = spline(x,y,h(75:199));
V(200:length(h)) = sqrt(2*q./rho(200:511));
q = 0.5*rho.*V.^2;
M_h = V./a;
if SABRE == 1
M = M_h(M_h<5);
h = h(1:length(M));
V = V(1:length(M)); rho = rho(1:length(M)); hmax_ab = length(M);
Mturb = M(M<2.2); h_turb = length(Mturb); Msub = M(M<1); h_subc = length(Msub);
theta0 = 1+((gamma-1)/2)*M.^2;
end
if turboramjet == 1
    M = M_h;
    theta0 = 1+((gamma-1)/2)*M.^2;
end

%% Aerodynamics from database

CopiadiCD = importfile3("C:\Users\Tommaso\Downloads\Copia di CD.xlsx", "Foglio1", [1, Inf]);
CopiadiCL = importfile4("C:\Users\Tommaso\Downloads\Copia di CL.xlsx", "Foglio1", [1, Inf]);
for i = 1:1:height(CopiadiCD)
    cL(i) = CopiadiCL(i,3);
    cD(i) = CopiadiCD(i,3);
    LD(i) = cL(i)/cD(i);
    M_data(i) = CopiadiCD(i,1);
    alpha_data(i) = CopiadiCD(i,2);
end
j = 0;
for i = 3:8:83
    j = j+1;
    Cl(j) = cL(i);
    Cd(j) = cD(i);
    M_skylon(j) = M_data(i);
end

[fitresult_cd, gof_cd] = createFit(M_skylon, Cd);
[fitresult_cl, gof_cl] = createFit(M_skylon, Cl);
Cd = fitresult_cd(M)';
Cl = fitresult_cl(M)';
% [fitresult_a, gof_a] = createFit(M_skylon, alpha);
% alpha = fitresult_a(M(M<=17))';
L_D = Cl./Cd;

if turboramjet == 1
%% Propulsion
% Turbojet
cp = gamma*R/(gamma-1)/1000;
Tt4_turb = 1550; % T total inlet turbine, hypothetical value
thetaT = Tt4_turb./T;
pi_c = 10; % compression ratio compressor, hypothetical value
tau_c = pi_c^((gamma-1)/gamma); 
f_turb = 0.0291; % stoichiometric;
%f_turb = cp*T.*(thetaT-tau_c*theta0)/hpr;

Specific_Thrust_turb = sqrt((2*theta0/(gamma-1)).*(thetaT./(theta0*tau_c)-1)*(tau_c-1) + thetaT.*M.^2./(theta0*tau_c)) - M;
Isp_turb = Specific_Thrust_turb.*a./(g*f_turb);
SFC_turb = 1./(g*Isp_turb);

% Ramjet
Tt4_ram = 1650; % T total exit combustor, hypotetical value
tau_max = Tt4_ram./T;
tau_b = tau_max./theta0;
f = 0.0291; % stoichiometric; 

Specific_Thrust_ram = M.*(sqrt(tau_b)-1);
Isp_ram = Specific_Thrust_ram.*a./(g*f);
SFC_ram = 1./(g*Isp_ram);

% Scramjet
% energetic formulation from Heiser & Pratt, "Hypersonic airbreathing propulsion" 
nke = 0.75; 
eta_b = 0.9;
Isp_scram = (1./(g*f)).*a.*M.*(sqrt(nke*(1+f).*(1+(eta_b*f*hpr./(cp*T.*theta0))))-1);
SFC_scram = 1./(g*Isp_ram);

% Rocket
eps = 40; %area ratio
Tc = 3550; pc = 70; 
Mmol = 12; yr = 1.21; %LOX/LH2 
Rr = 8314.5/Mmol;
RHO = sqrt(yr)*(2/(yr+1))^((yr+1)/(2*(yr-1)));
c_star = sqrt(Rr*Tc)/RHO;
pe_pc_test = 0.3; tollpe_pc = 0.01;
while 1 
pe_pc = (((RHO/eps)^2*(yr-1)/(2*yr))/(1-pe_pc_test^((yr-1)/yr)))^(yr/2);
if abs(pe_pc - pe_pc_test) < tollpe_pc
    break
else
    pe_pc_test = pe_pc_test-0.01;
end
if pe_pc_test < 0
    break
end
end
p_bar = p/(1*10^5);
pe = pe_pc*pc;
cf = RHO*sqrt((2*yr/(yr-1))*(1 - (pe/pc)^((yr-1)/yr))) + eps*(pe/pc - p_bar/pc);
Isp_rocket = c_star*cf/g;

% OVERALL
% NB: Hypotetical Mach and altitude for transitions, i'll insert analytical methods
for i = 1:1:length(h)
if Isp_turb(i) > Isp_ram(i)
   Isp(i) = Isp_turb(i);
   Mtr_turboram = M(i);
   h_turb = i;
end
if turboramjet == 1
    if Isp_ram(i) > Isp_turb(i) && M(i) <= 5
    Isp(i) = Isp_ram(i);
    hmax_ab = i;
    end
    if M(i) > 5
    Isp(i) = Isp_rocket(i);
    end
end
end

if turboramjet == 1
M = M(M<=5);
end
V = V(1:length(M)); rho = rho(1:length(M)); T = T(1:length(M)); a = a(1:length(M)); h_subc = length(M(M<1));
end

%% ITERATION 
N_engine = 1; 
if turboramjet == 1
    N_engine = 2;
end
while 1
z = 1; % loop control parameter
zz = 0;
tau = 0.25; % First iterative variable
while 1
    Kw = - 93.831*tau.^3 + 59.920*tau.^2 - 5.648*tau + 2.821; %wing body

    %L_D estimation from Curran Model
    F = (tau^0.333)*(Kw^0.75);
    A = 3; B = 2; %empirical coefficient
    L_D_est = A*((M+B)./M)*(1.11238 - 0.1866*F);

%% SABRE
if SABRE == 1
M_isp = [0.00778267917923067
0.0719042052302153
0.168092952961613
0.360463989769488
0.428163610646444
0.535002680341789
0.737965911219327
0.830531353540309
0.905303201555242
1.14017864639510
1.23266658485703
1.33227196104139
1.54575634078447
1.64540046889835
1.73079680425754
1.96920513333893
2.12944436191719
2.36443606254561
2.42853821263184
2.52823400998508
2.73832113723996
2.81311881987457
2.91279524126305
3.13713661992753
3.22615626069715
3.31162364126047
3.53236109047930
3.61784138835246
3.69973713274473
3.92405913544445
4.00956526793730
4.07365450071368
4.29080739645161
4.37980120260155
4.47591244647390
4.71086539517280
4.78917012742926
4.88528782995653
4.99919266813493]';

Isp = [26333.9382940109
25970.9618874773
25390.1996370236
24264.9727767695
23793.1034482758
23357.5317604356
22704.1742286751
22486.3883847550
22268.6025408348
22268.6025408348
22486.3883847550
22704.1742286751
22921.9600725953
22921.9600725953
22994.5553539020
23139.7459165154
22595.2813067150
21941.9237749546
21687.8402903811
21397.4591651542
20707.8039927405
20344.8275862069
20163.3393829401
19364.7912885662
19074.4101633394
18747.7313974591
18203.2667876588
17803.9927404718
17549.9092558983
16860.2540834846
16315.7894736842
16134.3012704174
15735.0272232305
15589.8366606170
15444.6460980036
15009.0744101633
14936.4791288566
14754.9909255898
14609.8003629764]'/9.81;

[fitresult_isp, gof_isp] = createFit3(M_isp, Isp);
Isp = fitresult_isp(M)';

M_thrust =[0.00342416434331318
0.109956599623499
0.184513763482972
0.397570863625016
0.500549627589753
0.564539022522818
0.784944348263935
0.880955637127682
0.980466204378580
1.20425425223193
1.29659531350530
1.37472945993521
1.60207117910400
1.68023640720722
1.77615445105103
2.01064236494236
2.08178572501629
2.16712822951537
2.40181040386490
2.47295376393884
2.56184216933503
2.78233296967776
2.86056036112760
2.95301797867590
3.14507940849504
3.23399112514622
3.34067896879297
3.52920226813331
3.63232089962795
3.73189363022548
3.92753204261505
4.01642821842957
4.11245504812997
4.34003765026694
4.42536461392936
4.51781446105933
4.74539706319630
4.80940199896602
4.90542105824810
5.00138572460187]';

Fu = N_engine*8*[153207.007052432
155384.278268005
157345.687262420
161918.889265323
164096.471297629
163434.897881163
161229.031526141
159471.362900293
159244.000460009
161848.333866902
163370.869634138
164675.989096549
167280.011686709
167710.492862094
168576.739097326
170086.842195237
168112.689714764
166574.613110872
162618.226914864
160644.074434392
159324.346585523
154713.224941178
153394.429542509
151637.071733394
147028.436622914
145052.730058776
142856.809839215
137811.166401955
136052.876142640
134076.237128302
128811.312485819
127272.925065194
125077.937295833
120903.513161535
119802.755701156
118264.057463798
114089.633329500
112990.740769520
111014.412571916
110568.701376608]';

[fitresult_thrust, gof_thrust] = createFit3(M_thrust, Fu);
Fu = fitresult_thrust(M)';


%% Rocket Mode

eps = 40; %area ratio
Tc = 3550; pc = 70; 
Mmol = 12; yr = 1.21; %LOX/LH2 
Rr = 8314.5/Mmol;
RHO = sqrt(yr)*(2/(yr+1))^((yr+1)/(2*(yr-1)));
c_star = sqrt(Rr*Tc)/RHO;
pe_pc_test = 0.3; tollpe_pc = 0.01;
while 1 
pe_pc = (((RHO/eps)^2*(yr-1)/(2*yr))/(1-pe_pc_test^((yr-1)/yr)))^(yr/2);
if abs(pe_pc - pe_pc_test) < tollpe_pc
    break
else
    pe_pc_test = pe_pc_test-0.01;
end
if pe_pc_test < 0
    break
end
end
p_bar = p/(1*10^5);
pe = pe_pc*pc;
cf = RHO*sqrt((2*yr/(yr-1))*(1 - (pe/pc)^((yr-1)/yr))) + eps*(pe/pc - p_bar/pc);
Isp_rocket = c_star*cf(length(cf))/g;
At = 0.154;
Thrust_r = cf*At*pc*10^5;
Fu_rocket = Thrust_r(hmax_ab)*N_engine;
end

%% WR weight ratio

WRprova = exp((V(hmax_ab))/(g*(Isp(1)+Isp(hmax_ab))/2));
WR_to = 1.01; % statistical value 
% Turbo climb
deltah = 100; % segment of climb [m]
flightpath_angle_sub = 8*pi/180; %hypotetical value
j = 1;
for i = 2:1:h_turb
    Vm = (V(i)+V(i-1))/2;
    Ispm = (Isp(i)+Isp(i-1))/2;
    L_Dm = (L_D(i) + L_D(i-1))/2;
    WR_climb_sub(j) = exp(deltah/(Vm*sin(flightpath_angle_sub)*Ispm*L_Dm));
    j = j+1;
end
WR_turbo = prod(WR_climb_sub);

% Subsonic cruise, buffer
R_sub = 0; % [m]
V_cruise_sub = V(h_subc); Isp_cruise_sub = Isp(h_subc); L_D_cruise_sub = L_D(h_subc);
WR_cruise_sub = exp(R_sub/(V_cruise_sub*Isp_cruise_sub*L_D_cruise_sub));

% Turboram climb
flightpath_angle_super = 3*pi/180; % Hypothetical value
j = 1;
for i = h_turb+1:1:hmax_ab
    Vm = (V(i)+V(i-1))/2;
    Ispm = (Isp(i)+Isp(i-1))/2;
    L_Dm = (L_D(i) + L_D(i-1))/2;
    WR_climb_super(j) = exp(deltah/(Vm*sin(flightpath_angle_super)*Ispm*L_Dm));
    j = j+1;
end
WR_super = prod(WR_climb_super);

% Total Airbreathing Weight ratio 
WR_ab = WR_to*prod(WR_climb_sub)*prod(WR_climb_super)*WR_cruise_sub;
% Rocket climb 
mu = 3.986*10^5;
r0 = 6371+hmax_ab/10;
r1 = h_targetorbit/1000+r0-hmax_ab/10; 
r2 = 400+r0-hmax_ab/10;
a_hohmann1 = (r1+r0)/2;
E_hohmann1 = - mu/(2*a_hohmann1);
V0 = sqrt(2*(E_hohmann1 + mu/r0));
V1_1 = sqrt(2*(E_hohmann1 + mu/r1));
a_hohmann2 = (r1+r2)/2;
E_hohmann2 = - mu/(2*a_hohmann2);
V1_2 = sqrt(2*(E_hohmann2 + mu/r1));
deltaVtot = V0 - V(hmax_ab)/1000; % + (V1_2 - V1_1);
deltaVp = 0.05*deltaVtot;
deltaVtot_rocket = deltaVtot + deltaVp;
if SABRE == 1
Isp_rocket_m = Isp_rocket;
end
if turboramjet == 1
    Isp_rocket_m = (Isp(hmax_ab+1)+Isp(length(h)))/2;
end
WR_r = exp(deltaVtot_rocket*1000/(g*Isp_rocket_m));

% Total Weight Ratio
WR = WR_ab*WR_r;

o_f = (6/7)*(1/WR_ab - (1/WR_ab)/WR_r)/((1-1/WR_ab) + (1/WR_ab - (1/WR_ab)/WR_r)/7);
rho_ppl = rho_f*(1+o_f)/(1+(rho_f/rho_ox)*o_f);
Istr = 10*(rho_ppl*1000/(WR-1))/ICI;

k = 0;
Spln = 1; %[m^2]
while 1
while 1

%% Weight Budget
% typical values from Czysz, "Chapter 16, Transatmospheric Launcher Sizing"
fcprv = 0.45; %0.45:0.01:0.5; %ton/person
Wcprv = fcprv*Ncrw; %ton
Cun = 1.9; % 1.9:0.1:2.1; %ton
fmnd = 1.05; %1.05:0.05:1.45; %ton/person
Csys = Cun + fmnd*Ncrw;
Wcrw = 0.14*Ncrw;
X = 0.697; %0.63:0.01:0.71; % (1/(1+mu) - fsys)
fsys = 0.2;
mu = 1/(X+fsys) -1;
ETW = (-36.25*X + 30.537)*TW0*WR; % suggested relation with X from Czysz, "Chapter 16, Transatmospheric Launcher Sizing"
%ETW = circa 19.8; % SABRE, data from REL
% NB: the ETW value is for the general airbreathing propulsion system and depends on thrust-weight at take off, i haven't found yet pratical sizing equation for specific configutations like ATR engine or dual mode ram-scram engine.
Wdry1 = (Istr*Kw*Spln/1000 + Wcprv + Csys + (TW0*WR/ETW)*(Wpay/1000+Wcrw))/(X-TW0*WR/ETW);

%% Volume Budget
% typical values from Czysz, "Chapter 16, Transatmospheric Launcher Sizing"
kvs = 0.02; % 0.02:0.01:0.04;
Vun = 5; %5:1:7; %m^3
fcrw = 11; %11:1:12; %m^3/person
kve = 0.25; %0.25:0.05:0.75; %m^3/ton thrust
kvv = 0.1; %0.1:0.1:0.2; 
kcprv = 5; %5:0.5:6; %m^3/person
kcrew = 0.9; %0.9:0.1:2; %m^3/person

Wdry2 = (tau*Spln^1.5*(1-kvv-kvs) - (kcprv+kcrew)*Ncrw - Wpay/rho_pay - Vun - fcrw*Ncrw)/((WR-1)/rho_ppl + kve*TW0*WR) - Wpay/1000 - 0.14*Ncrw;

%% Iteration
toll = 0.2;
if abs(Wdry1 - Wdry2) < toll
    GTOW = (Wdry1+Wpay/1000+0.14*Ncrw)*WR*1000;
    break
end
    Spln = Spln +1;
    k = k+1;
if k > 2000 % loop control parameter
    break
end
end
if SABRE == 1
if abs(TW0 - Fu(1)/(GTOW*9.81)) < 0.001
    break
end
TW0 = Fu(1)/(GTOW*9.81);
end
if turboramjet == 1
    break
end
end

%% DATA

W_TOT(z) = GTOW; 
Wdry_TOT(z) = Wdry1;
GTOW_Spln = GTOW/Spln;
Swet = Kw*Spln;
Vol_tot = tau*Spln^1.5; 
Wdry = Wdry1*1000;
WLH2 = ((1-1/WR_ab) + (1/WR_ab - (1/WR_ab)/WR_r)/7)*GTOW;
WLOX = (6/7)*(1/WR_ab - (1/WR_ab)/WR_r)*GTOW;
Wprop = WLH2+WLOX;
Wprop_TOT(z) =  Wprop;
Wstr_TOT = Istr*Kw*Spln;
Wengine = TW0(1)*WR*(Wdry+Wpay)/ETW;
Wsubsys = fsys*Wdry;

Vpay = Wpay/rho_pay;
Vempty = kvv*Vol_tot;
Vsubsys = Vun + kvs*Vol_tot;
Vengines = kve*TW0(1)*WR*(Wdry+Wpay)/1000;
Vprop = Wprop/(rho_ppl*1000);
VLOX = WLOX/(rho_ox*1000);
VLH2 = WLH2/(rho_f*1000);

%% Matching Chart
W_S = 100:1:1.5*10^3; % range of values
T_W = 0:0.1:5;

% Take off
l_to = 4000;
s_to = 3.2808*l_to;
TOP_25 = s_to/37.5;
TOP_25 = TOP_25*4.8824;
rho_to = 1.225;
Cl_to = 0.687;
sigma_to = rho_to./rho(1);
T_W_to = (W_S)/(sigma_to*TOP_25*Cl_to);

% Second segment requirement
E_2nd = Cl(1)/Cd(1);
if N_engine == 1
G_2nd = 0.024;
end
if N_engine ==2
G_2nd = 0.024;
end
if N_engine == 3
G_2nd = 0.027;
end
if N_engine == 4
G_2nd = 0.03;
end
T_W_2snd = N_engine/(N_engine-1)*(1/E_2nd + G_2nd)*(1/sigma_to)*ones(1,length(W_S));

% Subsonic climb
th_subclimb = 1; G_subclimb = tan(flightpath_angle_sub);
sigma_subclimb = (sigma_to + rho(h_turb)/rho(1))/2;
Cd_subclimb = (Cd(1)+Cd(h_turb))/2;
q_subclimb = (0.5*rho(1)*V(1)^2 + 0.5*rho(h_turb)*V(h_turb)^2)/2;
T_W_subclimb = (q_subclimb*Cd_subclimb./(g*W_S) + G_subclimb)*(1/(th_subclimb*sigma_subclimb));

% Subsonic cruise (best range)
th_subcruise = 1;
sigma_subcruise = rho(h_subc)/rho(1);
q_subcruise = 0.5*rho(h_subc)*V(h_subc)^2;
Cd_subcruise = Cd(h_subc);
T_W_subcruise = (q_subcruise*Cd(h_subc)./(g*W_S))*(1/(th_subcruise*sigma_subcruise)); 

% Supersonic climb
th_superclimb = 1;
sigma_superclimb = (rho(h_turb) + rho(hmax_ab))/(2*rho(1));
q_superclimb = (0.5*rho(h_turb)*V(h_turb)^2 + 0.5*rho(hmax_ab)*V(hmax_ab)^2)/2;
Cd_superclimb = (Cd(h_turb)+Cd(hmax_ab))/2;
G_superclimb = tan(flightpath_angle_super);
T_W_superclimb = (q_superclimb*Cd_superclimb./(g*W_S) + G_superclimb)*(1/(th_superclimb*sigma_superclimb));

% Landing
s_LFL = 1900;
s_ALD = s_LFL/1.6; 
kapp = 1.7;
Vapp = sqrt(s_ALD*kapp^2);
rho_l = 1.225;
sigma_l = rho_l/rho(1);
Cl_a = 0.864;
W_S_land = rho(1)*sigma_l*Cl_a*Vapp^2/(2*g);
W_S_land_to = W_S_land*WR*ones(1,length(T_W));

%% Supersonic Matching Chart

% Supersonic climb
sigma_superclimb_SUP = (rho(h_turb)+rho(hmax_ab))/(2*rho(h_turb));
T_W_superclimb_SUP = (q_superclimb*Cd_superclimb./(g*W_S) + G_superclimb)*(1/(th_superclimb*sigma_superclimb_SUP));

% Supersonic cruise
sigma_supercruise_SUP = rho(hmax_ab)/rho(h_turb);
th_supercruise = 1;
q_supercruise = 0.5*rho(hmax_ab)*V(hmax_ab)^2;
T_W_supercruise_SUP = (q_supercruise*Cd(hmax_ab)./(g*W_S))*(1/(th_supercruise*sigma_supercruise_SUP));

WSUP = GTOW/(WR_to*prod(WR_climb_sub)*WR_cruise_sub);
ROC_ram = (V(h_turb)+V(hmax_ab))/2*sin(flightpath_angle_super); 
deltaV_ram = V(hmax_ab) - V(h_turb);
h_ram = hmax_ab-h_turb;
Mmax_ram = M(hmax_ab);
% Access to space
[T_W_access_ram] = orbitReachingReq4(tau,WSUP,Wpay,Mmax_ram,h_ram*100,ROC_ram,Istr,rho_pay,deltaV_ram);
T_W_access_ram = T_W_access_ram(W_S);

% Identification of design point
WSUP = GTOW/(WR_to*prod(WR_climb_sub)*WR_cruise_sub);
W_S_Design_SUP = (W_S_land_to(1))*WSUP/GTOW;

% Max design point considering Subsonic W/S

% Supersonic Climb 
T_W_superclimb_Design_SUPGEN = (q_superclimb*Cd_superclimb./(g*W_S_land_to(1)*WSUP/GTOW) + G_superclimb)*(1/(th_superclimb*sigma_superclimb_SUP)); 

% Supersonic Cruise
T_W_supercruise_Design_SUPGEN = (q_supercruise*Cd(hmax_ab)./(g*W_S_land_to(1)*WSUP/GTOW))*(1/(th_supercruise*sigma_supercruise_SUP));

T_W_Design_MAX_SUPGEN(1) = T_W_superclimb_Design_SUPGEN;
T_W_Design_MAX_SUPGEN(2) = T_W_supercruise_Design_SUPGEN;
T_W_Design_MAX_TOTAL_SUPGEN = max(T_W_Design_MAX_SUPGEN);


%% Subsonic Matching Chart

ROC_turb = V(h_turb)/2*sin(flightpath_angle_sub); 
deltaV_turb = V(h_turb);
Mmax_turb = M(h_turb);
% Access to space
[T_W_access_turb] = orbitReachingReq4(tau,GTOW,Wpay,Mmax_turb,h_turb*100,ROC_turb,Istr,rho_pay,deltaV_turb);
T_W_access_turb = T_W_access_turb(W_S);

W_S_Design = W_S_land_to(1);

% Take-off
T_W_to_Design = (W_S_Design)/(sigma_to*TOP_25*Cl_to);

% Second Segment
T_W_2snd_Design = N_engine/(N_engine-1)*(1/E_2nd + G_2nd)*(1/sigma_to);

% Subsonic Climb
T_W_subclimb_Design = (q_subclimb*Cd_subclimb./(g*W_S_Design) + G_subclimb)*(1/(th_subclimb*sigma_subclimb));
 
% Subsonic Cruise
T_W_subcruise_Design = (q_subcruise*Cd(h_subc)./(g*W_S_Design))*(1/(th_subcruise*sigma_subcruise));

% Max design point (SUB)
T_W_Design_MAX_SUB(1) = T_W_to_Design;
T_W_Design_MAX_SUB(2) = T_W_2snd_Design;
T_W_Design_MAX_SUB(3) = T_W_subclimb_Design;
T_W_Design_MAX_SUB(4) = T_W_subcruise_Design;
T_W_Design_MAX_TOTAL_SUB = max(T_W_Design_MAX_SUB);

%% Rocket Mode Matching Chart
Wmax_ab = GTOW/WR_ab;
ROC_r = 125;
h_orbit = 1000;
h_rocket = h_orbit - hmax_ab;
[T_W_access_rock] = orbitReachingReq5(tau,Wmax_ab,Wpay,WR_r,h_rocket*100,ROC_r,Istr,rho_pay,deltaVtot_rocket*1000);

W_S_Design_ROCK = (W_S_land_to(1))*Wmax_ab/GTOW;
T_W_Design_ROCK = T_W_access_rock(W_S_Design_ROCK);
T_W_access_rock = T_W_access_rock(W_S);

%% Iteration
if GTOW_Spln < W_S_Design && Spln ~= 2002
    break
end
    tau = tau - 0.001;
    z = z+1;
    zz = zz+1;
if z > 150 % loop control parameter
    break
end
end
Wmax_ab = GTOW/WR_ab;
if SABRE == 1
TWsup = (Fu(h_turb)+Fu(hmax_ab))/(2*9.81*WSUP);
TWrock = Fu_rocket/(9.81*Wmax_ab);
end
if turboramjet == 1
    TWsup = TrW0*GTOW/WSUP;
    TWrock = TvW0*GTOW/Wmax_ab;
end
if TW0 >= T_W_Design_MAX_TOTAL_SUB
    if TWsup >= T_W_Design_MAX_TOTAL_SUPGEN
        if TWrock >= T_W_Design_ROCK
    break
        end
    end
end
if SABRE == 1
 N_engine = N_engine +1;
end
 if turboramjet == 1
 TrW0 = TrW0*(TW0+0.01)/TW0;
 TvW0 = TvW0+0.01;
 TW0 = TW0+0.01;
 end
 zz = 1;
end
W_TOT = W_TOT(1:zz); Wprop_TOT = Wprop_TOT(1:zz); Wdry_TOT = Wdry_TOT(1:zz);
TW0 = TW0*ones(1,length(W_S));
TWsup = TWsup*ones(1,length(W_S));
TWrock = TWrock*ones(1,length(W_S));

%% TPS sizing

L_D_re = 1.5; cD_re = 0.5; R = 0.1; phi = 70*pi/180;
[Q, qw, v_vc0, Tw, h_re, delta_t_tot, Mach, x, Istr_TPS] = function_TPS(Wdry1*1000, Spln,L_D_re, cD_re, R, phi);
WTPS = Istr_TPS*Swet;
Wstr = Wstr_TOT-WTPS;

%% PLOT

figure(2)
hold on
axis([0 1500 0 5])
xlabel('W/S [kg/m^2]')
ylabel('T/W')
title('Matching Chart - Turbo-ramjet mode')
plot(W_S,T_W_superclimb_SUP, 'blue','LineWidth',1.5)
plot(W_S,T_W_supercruise_SUP,'LineWidth',1.5)
plot(W_S_land_to*WSUP/GTOW,T_W, 'red','LineWidth',1.5)
plot(W_S,TWsup,'--','LineWidth',1.5)
plot(W_S,T_W_access_ram,'LineWidth',1.5)
plot(W_S_Design_SUP, T_W_Design_MAX_TOTAL_SUPGEN, 'ro','MarkerSize', 10, 'MarkerFaceColor','r')
legend('Climb req.','Cruise req.','Ref Surface','High-speed engine mode','Orbit reaching req.','Design point')

figure(1)
hold on
axis([0 1500 0 5])
xlabel('W/S [kg/m^2]')
ylabel('T/W')
title('Matching Chart - Turbojet mode')
plot(W_S,T_W_subclimb, 'blue','LineWidth',1.5)
plot(W_S,T_W_subcruise,'LineWidth',1.5)
plot(W_S_land_to,T_W, 'red','LineWidth',1.5)
plot(W_S,TW0,'--','LineWidth',1.5)
plot(W_S,T_W_access_turb,'LineWidth',1.5)
plot(W_S,T_W_2snd,'LineWidth',1.5)
plot(W_S,T_W_to,'LineWidth',1.5)
plot(W_S_Design, T_W_Design_MAX_TOTAL_SUB, 'ro', 'MarkerSize', 10, 'MarkerFaceColor','r')
legend('Climb req.','Cruise req.','Ref. Surface','Low-speed engine mode','Orbit reaching req.','2-segment req.','Take-off req.','Design point')

figure(3) 
hold on
axis([0 1500 T_W_Design_ROCK-0.02 T_W_Design_ROCK+0.02])
xlabel('W/S [kg/m^2]')
ylabel('T/W')
title('Matching Chart - Rocket mode')
plot(W_S,T_W_access_rock,'LineWidth',1.5,Color= '#77AC30')
plot(W_S_land_to*Wmax_ab/GTOW,T_W, 'red','LineWidth',1.5)
plot(W_S,TWrock,'--','LineWidth',1.5,Color='#7E2F8E')
plot(W_S_Design_ROCK,T_W_Design_ROCK,'ro', 'MarkerSize', 10, 'MarkerFaceColor','r')
legend('Orbit reaching req.','Ref. Surface','Engine rocket mode','Design point')

if SABRE == 1
figure(4)
hold on
xlabel('Mach')
ylabel('Altitude h [m]')
plot(M,h,'LineWidth',1.5)
end

if turboramjet == 1
figure(4)
hold on
xlabel('Mach')
ylabel('Altitude h [m]')
plot(M,h(1:hmax_ab),'LineWidth',1.5)  
end

figure(5)
hold on
xlabel('x [m]')
ylabel('Q [J/m^2]')
plot(x(2:length(x)),Q(2:length(Q)),'-o','LineWidth',1.5)

figure(6)
for i = 1:12:length(h_re)   
    hold on
    ylim([0 110000])
    xlim([0 85])
    xlabel('x [m]')
    ylabel('qw [W/m^2]')
    plot(x(2:length(x)),qw(i,(2:length(x))),'LineWidth',1.5)
end
legend('Mach 27','Mach 19.5','Mach 10.5','Mach 5','Mach 2.5','Mach 1', 'Mach 0.4')

figure(7)
for i = 1:12:length(h_re)   
    hold on
    ylim([200 1300])
    xlim([0 85])
    xlabel('x [m]')
    ylabel('Tw [K]')
    plot(x(2:length(x)),Tw(i,(2:length(x))),'LineWidth',1.5)
end
legend('Mach 27','Mach 19.5','Mach 10.5','Mach 5','Mach 2.5','Mach 1', 'Mach 0.4')

figure(8)
hold on
xlabel('Mach')
ylabel('Tw stagnation point [K]')
plot(Mach,Tw(:,(1))','-o','LineWidth',1.5)

figure(9)
hold on
xlabel('Mach')
ylabel('qw stagnation point [W/m^2]')
plot(Mach,qw(:,(1))','-o','LineWidth',1.5)

figure(10)
hold on
xlabel('Velocity Ratio')
ylabel('Altitude [m]')
plot(v_vc0,h_re,'LineWidth',1.5)

figure(11)
hold on
xlabel('Time from Entry [s]')
ylabel('Velocity Ratio')
plot(delta_t_tot,v_vc0,'LineWidth',1.5)

if SABRE == 1
figure(12)
hold on
xlabel('Flight Mach')
ylabel('L/D')
plot(M,L_D,'LineWidth',1.5)
plot(M,L_D_est,'LineWidth',1.5)
legend('L/D from CFD data', 'L/D from empirical method')
end
if turboramjet == 1
figure(12)
hold on
xlabel('Flight Mach')
ylabel('L/D')
plot(M,L_D(1:hmax_ab),'LineWidth',1.5)
plot(M,L_D_est,'LineWidth',1.5)
legend('L/D from CFD data', 'L/D from empirical method')
end

zz = 1:1:zz;
figure(13)
hold on
xlabel('N° of iterations for the W/S consistence')
ylabel('GTOW')
plot(zz,W_TOT,'LineWidth',1.5)

figure(14)
hold on
xlabel('N° of iterations for the W/S consistence')
ylabel('W propellant')
plot(zz,Wprop_TOT,'LineWidth',1.5,Color='#D95319')







