%start simulation
out = sim('Case4_Model2019b');
%% RESULTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%Liquid cooling
fp=npoints+1;
W_pump=out.simoutpump; %Pump work (W)
T_OutEquip=out.simouttoutequip; %Equipment EGW outlet temperature (K)
T_InRHX=out.simouttoutequip; %RHX EGW intlet temperature (K)
T_InEquip=out.simouttinequip; %Equipment inlet/RHX outlet EGW temperature (K)
Q_RHX=out.simoutqrhx; %Heat load in RHX (W)
UA_RHX=out.simoutUArhx; %UA values for RHX (W/K)
Delta_T1RHX=repelem(out.simoutDT1RHX,fp)';
Delta_T2RHX=out.simoutDT2RHX;
Delta_TLMRHX=out.simoutDTMLRHX;

%Interior flow
Dh=out.simoutDh;
Re_int=out.simoutRint;
f_int=out.simoutfint;
k_int=out.simoutkint;
mi_int=out.simoutmiint;
cp_int=out.simoutcpint';
rho_int=out.simoutrhoint;
V_int=out.simoutvint;
flowrate_int=out.simoutflowrateint;

%Exterior flow
Taw= out.simoutTaw;
r=out.simoutr;
Tenv=out.simouttenv;
rhoenv=out.simoutrhoenv;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Ram air inlet
W_fan=[out.simoutwfan(end,1),repelem(0,59),out.simoutwfan(end,1)]'; %Hydraulic power of fan (W)
m_dotramair=out.simoutmdotramair; %Ram air requires (kg/s)
mdra=m_dotramair(1);
T_outramair=repelem(out.simouttoutramair,fp)';
T_inramair=out.simouttintramair;


%Effiencies and Performance coeffients
W_fan_elect=W_fan./eta_fan; %Electric power consumption of fan (W)
wfe=W_fan_elect(1);
Q_dot_fan=(1-eta_fan).*W_fan_elect; %Heat load fan (W)
Wpump_elect=W_pump./eta_pump;
wpe=Wpump_elect(1);
Q_dot_pump=(1-eta_pump).*Wpump_elect;
w_tot=wfe+wpe;

%Energy balance of liquid cooling verification (tolerance due to due to floating points)
B1_LIQUID=Q_equip*10^3;
B2_LIQUID=Q_RHX;
if abs(B1_LIQUID-B2_LIQUID)<5*10^(-10)
        disp ('  Energy Balance of Liquid cooling verified')
else
    disp('  Energy Balance of Liquid cooling NOT verified')
end


%% Heat Results
% 
Qequip=repelem(Q_equip*10^3,fp)';
% figure
% subplot(1,2,1);
% plot(Qequip,'-s',...
%     'LineWidth',2,...
%     'MarkerSize',7,...
%     'Color','b',...
%     'MarkerEdgeColor','b',...
%     'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel ('Heat Load (W)');
% 
% 
% subplot(1,2,2);
% plot(Q_RHX,'-s',...
%     'LineWidth',2,...
%     'MarkerSize',7,...
%     'Color','b',...
%     'MarkerEdgeColor','b',...
%     'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel ('RHX Heat Flux (W)');
% 
%% Work Results
% 
% figure
% subplot(2,2,1);
% plot(W_fan_elect,'-s',...
%     'LineWidth',2,...
%     'MarkerSize',7,...
%     'Color','b',...
%     'MarkerEdgeColor','b',...
%     'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel ('Fan electric power consumption (W)');
% 
% subplot(2,2,2);
% plot(Wpump_elect,'-s',...
%     'LineWidth',2,...
%     'MarkerSize',7,...
%     'Color','b',...
%     'MarkerEdgeColor','b',...
%     'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel ('Pump electric power consumption (W)');
% 
% subplot(2,2,3);
% plot(m_dotramair,'-s',...
%     'LineWidth',2,...
%    'MarkerSize',7,...
%    'Color','b',...
%    'MarkerEdgeColor','b',...
%    'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel ('Ram air mass flow (kg/s)');
% 
%% Mission profiles
H_60=out.simoutheight;
Mach_60=out.simoutmach;
% figure
% subplot(1,2,1);
% plot(H_60,'-s',...
%     'LineWidth',2,...
%     'MarkerSize',7,...
%     'Color','b',...
%     'MarkerEdgeColor','b',...
%     'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel ('Height (m)');
% 
% subplot(1,2,2);
% plot(Mach_60,'-s',...
%     'LineWidth',2,...
%     'MarkerSize',7,...
%     'Color','b',...
%     'MarkerEdgeColor','b',...
%     'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel ('Mach number');


%% VCS Temperature vs Entropy
%figure
%plot(s_VCS,T_VCS,'s',...
   % 'LineWidth',2,...
    %'MarkerSize',7,...
    %'Color','b',...
    %'MarkerEdgeColor','b',...
   % 'MarkerFaceColor','b');
%xlabel ('Flight Points'); 
%ylabel ('W Comp (W)');



%% Equipment mass
%Data from Embraer 

flowrate_ramair=zeros(1,61);
for i=1:61
flowrate_ramair(i)=m_dotramair(i)/rhoenv(i);
end

flowrate_ramair_max=max(flowrate_ramair);


m_pump=max(2.*(1942.2.*flowrate_int + 1.619)); %different expression
m_rhx=max(0.0008*Q_RHX - 0.2165);
%m_rhx2=functs(max(UA_RHX),mdot_liquid(1),rho_int(1),V_int(1));
al_density = 2780; %Alluminum density - SHX material [kg/m³]
m_liquid_ducts=max(al_density*pi*Dh*L_duct_liquid*t_duct_liquid);
m_liquid=max(rho_int.*pi.*Dh.^2.*L_duct_liquid)/4;

m_fan=10.542*flowrate_ramair_max + 0.9366;

m_TMS=2*(m_pump+m_rhx+m_liquid_ducts+m_liquid+m_fan);

%% Save to Excel
%folder = 'C:\Users\ASUS\Desktop\MODELAÇÃO\Case 4';
baseFileName='Case4TMSResults.xlsx';
fileName = fullfile(baseFileName);


T_in_Equi=repelem(T_InEquip,fp)';


INPUTS=table(Mach_60,H_60,T_in_Equi,Qequip);
writetable(INPUTS,fileName,'Sheet',1)


RESULTS_Liquid=table(T_OutEquip,T_InRHX,Q_RHX,UA_RHX,Delta_T1RHX,Delta_T2RHX,Delta_TLMRHX,W_pump,Wpump_elect,Q_dot_pump);
writetable(RESULTS_Liquid,fileName,'Sheet',2)

RESULTS_Ramair=table(m_dotramair,T_outramair,T_inramair,W_fan,W_fan_elect,Q_dot_fan);
writetable(RESULTS_Ramair,fileName,'Sheet',3)

RESULTS_Mass=table(m_pump,m_liquid_ducts,m_liquid,m_rhx,m_fan,m_TMS);
%writetable(RESULTS_Mass,fileName,'Sheet',4)


%winopen(fileName)

