# Plot temperature profile from monitor with analytic result

# simulation with a single material only
clear;
if (getnamed('solid_2','enabled')==0) {

T_left = 350;
T_right = 280;
k = 148;  # thermal conductivity of Silion in W/mK

temp = getresult('HEAT::monitor','temperature');
T_sim = pinch(temp.T);
x_sim = temp.x;

x = linspace(min(x_sim),max(x_sim),101);
T = linspace(T_left,T_right,101);  # linear change in temperature assuming constant k
 
plot(x,T);
holdon;
plot(x_sim,T_sim);  
legend('analytical','simulation');
setplot('x label','x (m)');
setplot('y label','Temp (K)');

flux = getresult('HEAT','boundaries');
Q_sim = flux.P_left/getnamed('HEAT simulation region','z span')/getnamed('HEAT','norm length');  # flux in W/m^2

Q = k * (T_left - T_right) / getnamed('HEAT simulation region','x span');

?'Analytic heat flux = '+num2str(Q)+' W/m^2';
?'Heat flux from simulation = '+num2str(Q_sim)+' W/m^2';
}

# simulation with two materials

else {
clear;
T_left = 350;
T_right = 280;
k1 = 148;  # thermal conductivity of Silion in W/mK
k2 = 1.38;  # thermal conductivity of SiO2 in W/mK

temp = getresult('HEAT::monitor','temperature');
T_sim = pinch(temp.T);
x_sim = temp.x;

x = linspace(min(x_sim),max(x_sim),101);
k_net = 2*k1*k2/(k1+k2);
Q = k_net * (T_left - T_right) / getnamed('HEAT simulation region','x span');

T_mid = T_left - Q / k1 *getnamed('HEAT simulation region','x span')/2;
T1 = linspace(T_left,T_mid,51);
T2 = linspace(T_mid,T_right,51);
T = [T1;T2(2:51)];
 
plot(x,T);
holdon;
plot(x_sim,T_sim);  
legend('analytical','simulation');
setplot('x label','x (m)');
setplot('y label','Temp (K)');


flux = getresult('HEAT','boundaries');
Q_sim = flux.P_left/getnamed('HEAT simulation region','z span')/getnamed('HEAT','norm length');  # flux in W/m^2

?'Analytic heat flux = '+num2str(Q)+' W/m^2';
?'Heat flux from simulation = '+num2str(Q_sim)+' W/m^2';
}