%start simulation
out = sim('Case1_Model2019b');
%% RESULTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%Liquid cooling
fp=npoints+1;
Q_evap=repelem(out.simoutevap,fp)'; %Heat load in evaporator (W)
W_pump=out.simoutpump; %Pump work (W)

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

%Exterior flow
Taw= out.simoutTaw;
Pr_ext=out.simoutPrext;
mi_ext=out.simoutmiext;
cp_ext=out.simoutcpext;
r=out.simoutr;
V_ext=out.simoutvext;
Tenv=out.simouttenv;
rhoenv=out.simoutrhoenv;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%VCS
Q_cond=repelem(out.simoutcond,fp)'; %Heat load in condenser (W)
W_comp=repelem(out.simoutcomp,fp)'; %Compression work (W)

psat_cond=repelem(out.simoutpsatcond,fp)'; 
psat_evap=repelem(out.simoutpsatevap,fp)'; 
m_dotref=repelem(out.simoutmdotref,fp)'; 
T1_VCS=repelem(out.simoutT1VCS,fp)'; %Outlet evaporator temperature (K)/Inlet compressor temperature (K)
T2s_VCS=repelem(out.simoutT2sVCS,fp)';%Outlet compressor temperature (K)/Inlet condenser temperature (K) -isentropic
T2_VCS=repelem(out.simoutT2VCS,fp)'; %Outlet compressor temperature (K)/Inlet condenser temperature (K) -real
T3_VCS=repelem(out.simoutT3VCS,fp)'; %Outlet condenser temperature (K)/Inlet expansion valve temperature (K)
T4_VCS=repelem(out.simoutT4VCS,fp)';%Outlet expansion valve temperature (K)/Inlet evaporator temperature (K) 
s1_VCS=repelem(out.simouts1VCS,fp)';
s2s_VCS=repelem(out.simouts2sVCS,fp)';
s2_VCS=repelem(out.simouts2VCS,fp)';
s3_VCS=repelem(out.simouts3VCS,fp)';
s4_VCS=repelem(out.simouts4VCS,fp)';
h1_VCS=repelem(out.simouth1VCS,fp)';
h2s_VCS=repelem(out.simouth2sVCS,fp)';
h2_VCS=repelem(out.simouth2VCS,fp)';
h3_VCS=repelem(out.simouth3VCS,fp)';
h4_VCS=repelem(out.simouth4VCS,fp)';
v1_VCS=repelem(out.simoutv1VCS,fp);
v2_VCS=repelem(out.simoutv2VCS,fp);
v3_VCS=repelem(out.simoutv3VCS,fp);
v4_VCS=repelem(out.simoutv4VCS,fp);
%s_VCS=[s1_VCS,s2s_VCS,s2_VCS,s3_VCS,s4_VCS];
%T_VCS=[T1_VCS,T2s_VCS,T2_VCS,T3_VCS,T4_VCS];

Delta_TLMcond=out.simoutDELTATLMcond;
Delta_T1cond=repelem(out.simoutDELTAT1cond,fp)';
Delta_T2cond=out.simoutDELTAT2cond;
UA_cond=out.simoutUAcond;
Delta_TLMevap=out.simoutDTMevap;
Delta_T1evap=out.simoutDT1evap;
Delta_T2evap=repelem(out.simoutDT2evap,fp)';
UA_evap=out.simoutUAevap;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Ram air inlet
W_fan=[out.simoutfan(end,1),repelem(0,59),out.simoutfan(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_comp_elect=W_comp./eta_compressor;
wce=W_comp_elect(1);
COP=Q_evap./W_comp_elect; %VCS coefficient of performance - COP (W)
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=wce+wfe+wpe;

%Energy balance of liquid cooling verification (tolerance due to due to floating points)
B1_LIQUID=Q_equip*10^3;
B2_LIQUID=Q_evap;
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


%Energy balance of VCS verification 
B1_VCS=W_comp+Q_evap;
B2_VCS=Q_cond;
if abs(B1_VCS-B2_VCS)<5*10^(-10)
        disp ('  Energy Balance of VCS verified')
else
    disp('  Energy Balance of VCS NOT verified')
end


%% Heat Results

Qequip=repelem(Q_equip*10^3,fp)';
% figure
% subplot(2,2,1);
% plot(Qequip,'-s',...
%     'LineWidth',2,...
%     'MarkerSize',7,...
%     'Color','b',...
%     'MarkerEdgeColor','b',...
%     'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel ('Heat Load (W)');
% 
% subplot(2,2,2);
% plot(Q_evap,'-s',...
%     'LineWidth',2,...
%     'MarkerSize',7,...
%     'Color','b',...
%     'MarkerEdgeColor','b',...
%     'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel (' Evaporator Heat Flux (W)');
% 
% subplot(2,2,3);
% plot(Q_cond,'-s',...
%     'LineWidth',2,...
%     'MarkerSize',7,...
%     'Color','b',...
%     'MarkerEdgeColor','b',...
%     'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel ('Condenser Heat Flux (W)');

%% Work Results

% figure
% subplot(2,2,1);
% plot(W_comp_elect,'-s',...
%     'LineWidth',2,...
%     'MarkerSize',7,...
%     'Color','b',...
%     'MarkerEdgeColor','b',...
%     'MarkerFaceColor','b');
% xlabel ('Flight Duration (min)'); 
% ylabel ('Compressor electric power consumption (W)');
% 
% subplot(2,2,2);
% 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,3);
% 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,4);
% 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_condenser=max(0.0008*Q_cond - 0.2165);
m_evaporator=max(0.0005*Q_evap + 2.629);
m_compressor= max(0.001*W_comp + 8.6225);
al_density = 2780; %Alluminum density - SHX material [kg/m³]
cu_density=8960; %Copper density [kg/m^3]


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_ref_ducts=max(cu_density*pi*D_vcs*L_duct_vcs*t_duct_vcs); 
rho_R314amax=1/(min([v1_VCS,v2_VCS,v3_VCS,v4_VCS])); %Max density is on stage 3 (liquid)
m_ref=rho_R314amax.*pi.*D_vcs.^2.*L_duct_vcs/4; %use of max density to calculate the refrigerant masS

m_fan=10.542*flowrate_ramair_max + 0.9366;

m_TMS=2*(m_pump+m_condenser+m_evaporator+m_compressor+m_liquid_ducts+m_liquid+m_ref_ducts+m_ref+m_fan);

%% Save to Excel
%folder = 'C:\Users\Utente\Desktop\Laurea Magistrale\Thesis\Matlab_model\Technical & Architectures\Case 1';
baseFileName='Case1TMSResults.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(Delta_T1evap,Delta_T2evap,Delta_TLMevap,UA_evap,Q_evap,W_pump,Wpump_elect,Q_dot_pump);
writetable(RESULTS_Liquid,fileName,'Sheet',2)

RESULTS_VCS=table(psat_evap,m_dotref,T1_VCS,T2_VCS,T2s_VCS,T3_VCS,T4_VCS,s1_VCS,s2_VCS,s2s_VCS,s3_VCS,s4_VCS,h1_VCS,h2_VCS,h2s_VCS,h3_VCS,h4_VCS,Delta_TLMcond,Delta_T1cond,Delta_T2cond,UA_cond,Q_cond,W_comp,W_comp_elect,COP,psat_cond);
writetable(RESULTS_VCS,fileName,'Sheet',3)

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

RESULTS_Mass=table(m_pump,m_evaporator,m_liquid_ducts,m_liquid,m_condenser,m_compressor,m_ref,m_ref_ducts,m_fan,m_TMS);
writetable(RESULTS_Mass,fileName,'Sheet',5)


%winopen(fileName)

