%The following code is written for Matlab, working with at least version
%R2016a.
%This code generates the figures included in the paper L. Cini and J.I.
%Mackenzie, "Analytical thermal model for end-pumped solid-state lasers",
%Applied Physics B VX IsY pp ZZZ-ZZZ DOI:AAAAAAAAAAAA
%Data and the code (below) is located at DOI:10.5258/SOTON/D0117 and should
%be cited when used according to creative common copyright license CC BY.


% Initial parameters for the pump
P=25; % Incident pump power (W)
eta=0.25; %proportion of the power absorbed converted to heat
a=0.3; %Top hat beam radius (mm)
b=1.25; %Rod radius (mm)
L=5; %Rod length (mm)
alpha=350*10^-3;% Absorption coefficient (1/mm)
k_0=15.09*10^-3; %Thermal conductivity of ~1at.% NdYAG at a temperature T0 (W/mm.K)
T_0=164.17; %Thermal conductivity fitting curve common reference temperature (K) 
m_rt=-0.75; %Room temperature power coefficient for thermal conductivity equation (1) in the paper
m_ct=-1.77;%Cryogenic temperature power coefficient for thermal conductivity equation (1) in the paper
h=2*10^-2; %Thermal conductance at the rod boundary (W/m^2.K)
T_c=300; %Heatsink temperature (K)

%Top hat pump distribution parameters
x_1=0:a/100:a;
x_2=a:a/100:b;
S_1=eta*P*alpha/(pi*a^2)*10^9; %Thermal loading density coefficient (W/m^3)

%Gaussian pump distribution parameters
x=0:a/100:b;
w=0.3;% units mm
S_g=2*eta*P*alpha/(pi*w^2)*exp(-2*(x.^2/w^2))*10^9; %Thermal loading density coefficient(W/m^3)

%Super Gaussian pump distribution parameters
S_r_1=(2/pi)^(3/2)*eta*P*alpha/(w^2)*exp(-2*(x.^4/w^4))*10^9;
S_r_2=4*(2)^(1/4)*eta*P*alpha/(pi*w^2*3.6256)*exp(-2*(x.^8/w^8))*10^9; %Thermal loading density coefficient(W/m^3)

%Donut beam pump distribution parameters
d_1=3/4*a;
d_2=sqrt(a^2+d_1^2);
xd_1=0:a/100:d_1;
xd_2=d_1:a/100:d_2;
xd_3=d_2:a/100:b;
S_d=eta*P*alpha/(pi*(d_2^2-d_1^2))*10^9; %Thermal loading density coefficient(W/m^3)

figure(3)

h1=plot(x_1,ones(size(x_1))*S_1,'k');hold on
h2=plot(x,S_g,'-.r');
h3=plot(xd_2,ones(size(xd_2))*S_d,'--m');
h6=plot(x,S_r_1,'b');
h7=plot(x,S_r_2,'g');
plot([a a], [0 S_1],'k',[d_1 d_1], [0 S_d],'--m',[d_2 d_2], [0 S_d],'--m');
xlim([0 b/2.5]);
set(gca,'xTick',[0 1/10 2/10 3/10 4/10 5/10])
ylabel('$S(r,0)$ [Wm$^{-3}$]','Interpreter','Latex');
xlabel('$r$ [mm]','Interpreter','Latex');
legend([h1 h2 h6 h7 h3],{'TH','G','SG (n=4)','SG (n=8)','D'},'Interpreter','Latex');

%Parameters to create Figure 4 contour map
x_1_cm=0:a/1500:a;
x_2_cm=a:a/1500:b;
z=0:L/200:L;
[X_1,Z_1] = meshgrid(x_1_cm,z);
[X_2,Z_2] = meshgrid(x_2_cm,z);
T_1=(eta*alpha*P.*exp(-alpha.*Z_1)*T_0^m_rt*(m_rt+1)/(4*k_0*pi).*(1-X_1.^2/a^2+log(b^2/a^2))+(T_c+eta*alpha*P.*exp(-alpha.*Z_1)/(2*pi*b*h)).^(m_rt+1)).^(1/(m_rt+1));
T_2=(eta*alpha*P.*exp(-alpha.*Z_2)*T_0^m_rt*(m_rt+1)/(4*k_0*pi).*log(b^2./X_2.^2)+(T_c+eta*alpha*P.*exp(-alpha*Z_2)/(2*pi*b*h)).^(m_rt+1)).^(1/(m_rt+1));

figure(4)

contourf(Z_1,X_1,T_1);hold on
contourf(Z_2,X_2,T_2);
xlim([0 L])
ylim([0 b])
ylabel('$r$ [mm]','Interpreter','Latex');
xlabel('$z$ [mm]','Interpreter','Latex');
c=colorbar;
c.Label.String = '$T(r,z)$ [K]';
c.Label.Interpreter = 'latex';
hold off 

%Equations for top hat pumping (Figure 5)
T_1_true=@(x_1)(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*(1-x_1.^2/a^2+log(b^2/a^2))+(T_c+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1));
T_2_true=@(x_2)(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*log(b^2./x_2.^2)+(T_c+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1));
T_1=@(x_1)(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*(1-x_1.^2/a^2+log(b^2/a^2))+(T_c+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c+eta*alpha*P/(2*pi*b*h));
T_2=@(x_2)(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*log(b^2./x_2.^2)+(T_c+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c+eta*alpha*P/(2*pi*b*h));
T_avg=T_c;
k_0_i=k_0*(T_avg/T_0)^m_rt;
T_th_1=eta*P*alpha/(4*pi*k_0_i)*(1-x_1.^2/a^2+log(b^2/a^2));
T_th_2=eta*P*alpha/(4*k_0_i*pi)*log(b^2./x_2.^2);

figure(5)
h1=plot(x_1,T_th_1,'k');hold on
h2=plot(x_1,T_1(x_1),'--b');
plot(x_2,T_th_2,'k');
plot(x_2,T_2(x_2),'--b');
ylim([0 80]);
xlim([0 b]);
ylabel('$\Delta T(r)$ [K]','Interpreter','Latex');
xlabel('$r$ [mm]','Interpreter','Latex');
legend([h1 h2],{'A. K. Cousins','Analytical expression with $k(T)$'},'Interpreter','Latex');
hold off 

%Top hat pumping Different Coolant Temperature
figure(6)
T_c_1=75;

T_1_ct_1=(eta*alpha*P*T_0^m_ct*(m_ct+1)/(4*k_0*pi)*(1-x_1.^2/a^2+log(b^2/a^2))+(T_c_1+eta*alpha*P/(2*pi*b*h))^(m_ct+1)).^(1/(m_ct+1))-(T_c_1+eta*alpha*P/(2*pi*b*h));
T_2_ct_1=(eta*alpha*P*T_0^m_ct*(m_ct+1)/(4*k_0*pi)*log(b^2./x_2.^2)+(T_c_1+eta*alpha*P/(2*pi*b*h))^(m_ct+1)).^(1/(m_ct+1))-(T_c_1+eta*alpha*P/(2*pi*b*h));
h1=plot(x_1,T_1_ct_1,'k');hold on
plot(x_2,T_2_ct_1,'k');

T_c_2=100;

T_1_ct_2=(eta*alpha*P*T_0^m_ct*(m_ct+1)/(4*k_0*pi)*(1-x_1.^2/a^2+log(b^2/a^2))+(T_c_2+eta*alpha*P/(2*pi*b*h))^(m_ct+1)).^(1/(m_ct+1))-(T_c_2+eta*alpha*P/(2*pi*b*h));
T_2_ct_2=(eta*alpha*P*T_0^m_ct*(m_ct+1)/(4*k_0*pi)*log(b^2./x_2.^2)+(T_c_2+eta*alpha*P/(2*pi*b*h))^(m_ct+1)).^(1/(m_ct+1))-(T_c_2+eta*alpha*P/(2*pi*b*h));
h2=plot(x_1,T_1_ct_2,'--k');
plot(x_2,T_2_ct_2,'--k');

T_c_7=125;

T_1_ct_3=(eta*alpha*P*T_0^m_ct*(m_ct+1)/(4*k_0*pi)*(1-x_1.^2/a^2+log(b^2/a^2))+(T_c_7+eta*alpha*P/(2*pi*b*h))^(m_ct+1)).^(1/(m_ct+1))-(T_c_7+eta*alpha*P/(2*pi*b*h));
T_2_ct_3=(eta*alpha*P*T_0^m_ct*(m_ct+1)/(4*k_0*pi)*log(b^2./x_2.^2)+(T_c_7+eta*alpha*P/(2*pi*b*h))^(m_ct+1)).^(1/(m_ct+1))-(T_c_7+eta*alpha*P/(2*pi*b*h));
h3=plot(x_1,T_1_ct_3,'-.k');
plot(x_2,T_2_ct_3,'-.k');

T_c_3=225;

T_1_rt_1=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*(1-x_1.^2/a^2+log(b^2/a^2))+(T_c_3+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_3+eta*alpha*P/(2*pi*b*h));
T_2_rt_1=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*log(b^2./x_2.^2)+(T_c_3+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_3+eta*alpha*P/(2*pi*b*h));
h4=plot(x_1,T_1_rt_1,'b');
plot(x_2,T_2_rt_1,'b');

T_c_4=300;

T_1_rt_2=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*(1-x_1.^2/a^2+log(b^2/a^2))+(T_c_4+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_4+eta*alpha*P/(2*pi*b*h));
T_2_rt_2=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*log(b^2./x_2.^2)+(T_c_4+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_4+eta*alpha*P/(2*pi*b*h));
h5=plot(x_1,T_1_rt_2,'--b');
plot(x_2,T_2_rt_2,'--b');

T_c_5=375;

T_1_rt_3=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*(1-x_1.^2/a^2+log(b^2/a^2))+(T_c_5+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_5+eta*alpha*P/(2*pi*b*h));
T_2_rt_3=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*log(b^2./x_2.^2)+(T_c_5+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_5+eta*alpha*P/(2*pi*b*h));
h6=plot(x_1,T_1_rt_3,'-.b');
plot(x_2,T_2_rt_3,'-.b');

T_c_6=450;

T_1_rt_4=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*(1-x_1.^2/a^2+log(b^2/a^2))+(T_c_6+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_6+eta*alpha*P/(2*pi*b*h));
T_2_rt_4=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*log(b^2./x_2.^2)+(T_c_6+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_6+eta*alpha*P/(2*pi*b*h));
h7=plot(x_1,T_1_rt_4,'m');
plot(x_2,T_2_rt_4,'m');

ylim([0 110]);
xlim([0 b]);
ylabel('$\Delta T(r)$ [K]','Interpreter','Latex');
xlabel('$r$ [mm]','Interpreter','Latex');
legend([h1 h2 h3 h4 h5 h6 h7],{'$T_c=75$ K','$T_c=100$ K','$T_c=125$ K','$T_c=225$ K','$T_c=300$ K','$T_c=375$ K','$T_c=450$ K'},'Interpreter','Latex');
hold off 

figure(7)
T_g=(eta*P*alpha*(m_rt+1)*T_0^m_rt/(4*pi*k_0)*(log(b^2./x.^2)+expint(2*b^2/w^2)-expint(2*x.^2/w^2))+(T_c+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c+eta*P*alpha/(2*pi*b*h));
T_g_true=@(x)(eta*P*alpha*(m_rt+1)*T_0^m_rt/(4*pi*k_0)*(log(b^2./x.^2)+expint(2*b^2/w^2)-expint(2*x.^2/w^2))+(T_c+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1));
k_0_i=k_0*((T_avg)/T_0)^m_rt;
T_i=eta*P*alpha/(4*pi*k_0_i)*(log(b^2./x.^2)+expint(2*b^2/w^2)-expint(2*x.^2/w^2));
h1=plot(x,T_i,'k');hold on
h2=plot(x,T_g,'--b');
ylim([0 90]);
xlim([0 b]);
ylabel('$\Delta T(r)$ [K]','Interpreter','Latex');
xlabel('$r$ [mm]','Interpreter','Latex');
legend([h1 h2],{'Innocenzi et al.','Analytical expression with $k(T)$'},'Interpreter','Latex');
hold off

figure(8)
T_g_ct_1=(eta*P*alpha*(m_ct+1)*T_0^m_ct/(4*pi*k_0)*(log(b^2./x.^2)+expint(2*b^2/w^2)-expint(2*x.^2/w^2))+(T_c_1+eta*P*alpha/(2*pi*b*h))^(m_ct+1)).^(1/(m_ct+1))-(T_c_1+eta*P*alpha/(2*pi*b*h));
h1=plot(x,T_g_ct_1,'k');hold on

T_g_ct_2=(eta*P*alpha*(m_ct+1)*T_0^m_ct/(4*pi*k_0)*(log(b^2./x.^2)+expint(2*b^2/w^2)-expint(2*x.^2/w^2))+(T_c_2+eta*P*alpha/(2*pi*b*h))^(m_ct+1)).^(1/(m_ct+1))-(T_c_2+eta*P*alpha/(2*pi*b*h));
h2=plot(x,T_g_ct_2,'--k');

T_g_ct_3=(eta*P*alpha*(m_ct+1)*T_0^m_ct/(4*pi*k_0)*(log(b^2./x.^2)+expint(2*b^2/w^2)-expint(2*x.^2/w^2))+(T_c_7+eta*P*alpha/(2*pi*b*h))^(m_ct+1)).^(1/(m_ct+1))-(T_c_7+eta*P*alpha/(2*pi*b*h));
h3=plot(x,T_g_ct_3,'-.k');

T_g_rt_1=(eta*P*alpha*(m_rt+1)*T_0^m_rt/(4*pi*k_0)*(log(b^2./x.^2)+expint(2*b^2/w^2)-expint(2*x.^2/w^2))+(T_c_3+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_3+eta*P*alpha/(2*pi*b*h));
h4=plot(x,T_g_rt_1,'b');

T_g_rt_2=(eta*P*alpha*(m_rt+1)*T_0^m_rt/(4*pi*k_0)*(log(b^2./x.^2)+expint(2*b^2/w^2)-expint(2*x.^2/w^2))+(T_c_4+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_4+eta*P*alpha/(2*pi*b*h));
h5=plot(x,T_g_rt_2,'--b');

T_g_rt_3=(eta*P*alpha*(m_rt+1)*T_0^m_rt/(4*pi*k_0)*(log(b^2./x.^2)+expint(2*b^2/w^2)-expint(2*x.^2/w^2))+(T_c_5+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_5+eta*P*alpha/(2*pi*b*h));
h6=plot(x,T_g_rt_3,'-.b');

T_g_rt_4=(eta*P*alpha*(m_rt+1)*T_0^m_rt/(4*pi*k_0)*(log(b^2./x.^2)+expint(2*b^2/w^2)-expint(2*x.^2/w^2))+(T_c_6+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c_6+eta*P*alpha/(2*pi*b*h));
h7=plot(x,T_g_rt_4,'m');

ylim([0 120]);
xlim([0 b]);
ylabel('$\Delta T(r)$ [K]','Interpreter','Latex');
xlabel('$r$ [mm]','Interpreter','Latex');
legend([h1 h2 h3 h4 h5 h6 h7],{'$T_c=75$ K','$T_c=100$ K','$T_c=125$ K','$T_c=225$ K','$T_c=300$ K','$T_c=375$ K','$T_c=450$ K'},'Interpreter','Latex');
hold off

figure(9)
n=2;
S_n_2=10^9*n*eta*P*alpha/(2*pi*4^(-1/n)*w^2*gamma(2/n))*exp(-2*x.^n/w^n);
h1=plot(x,S_n_2,'k');hold on
n=4;
S_n_4=10^9*n*eta*P*alpha/(2*pi*4^(-1/n)*w^2*gamma(2/n))*exp(-2*x.^n/w^n);
h2=plot(x,S_n_4,'--b');
n=8;
S_n_8=10^9*n*eta*P*alpha/(2*pi*4^(-1/n)*w^2*gamma(2/n))*exp(-2*x.^n/w^n);
h3=plot(x,S_n_8,'-.g');
n=16;
S_n_16=10^9*n*eta*P*alpha/(2*pi*4^(-1/n)*w^2*gamma(2/n))*exp(-2*x.^n/w^n);
h4=plot(x,S_n_16,'m');
n=32;
S_n_32=10^9*n*eta*P*alpha/(2*pi*4^(-1/n)*w^2*gamma(2/n))*exp(-2*x.^n/w^n);
h5=plot(x,S_n_32,'--c');
xlim([0 b/3]);
ylabel('$S(r,0)$ [Wm$^{-3}$]','Interpreter','Latex');
xlabel('$r$ [mm]','Interpreter','Latex');
legend([h1 h2 h3 h4 h5],{'SG (n=2)','SG (n=4)','SG (n=8)','SG (n=16)','SG (n=32)'},'Interpreter','Latex');
hold off

figure(10)
n=2;
T_sg_n_2=(2^(2/n-3)*n*eta*P*alpha*T_0^m_rt*(m_rt+1)/(pi*k_0*w^2*gamma(2/n))*(b^2*hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*b^n/w^n)-x.^2.*hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*x.^n/w^n))+(T_c+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-T_c-eta*P*alpha/(2*pi*b*h);
h1=plot(x,T_sg_n_2,'k');hold on
n=4;
T_sg_n_4=(2^(2/n-3)*n*eta*P*alpha*T_0^m_rt*(m_rt+1)/(pi*k_0*w^2*gamma(2/n))*(b^2*hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*b^n/w^n)-x.^2.*hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*x.^n/w^n))+(T_c+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-T_c-eta*P*alpha/(2*pi*b*h);
h2=plot(x,T_sg_n_4,'--b');
n=8;
T_sg_n_8=(2^(2/n-3)*n*eta*P*alpha*T_0^m_rt*(m_rt+1)/(pi*k_0*w^2*gamma(2/n))*(b^2*hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*b^n/w^n)-x.^2.*hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*x.^n/w^n))+(T_c+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-T_c-eta*P*alpha/(2*pi*b*h);
h3=plot(x,T_sg_n_8,'-.g');
n=16;
T_sg_n_16=(2^(2/n-3)*n*eta*P*alpha*T_0^m_rt*(m_rt+1)/(pi*k_0*w^2*gamma(2/n))*(b^2*hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*b^n/w^n)-x.^2.*hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*x.^n/w^n))+(T_c+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-T_c-eta*P*alpha/(2*pi*b*h);
h4=plot(x,T_sg_n_16,'m');
n=32;
T_sg_n_32=(2^(2/n-3)*n*eta*P*alpha*T_0^m_rt*(m_rt+1)/(pi*k_0*w^2*gamma(2/n))*(b^2*hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*b^n/w^n)-x.^2.*hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*x.^n/w^n))+(T_c+eta*P*alpha/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-T_c-eta*P*alpha/(2*pi*b*h);
h5=plot(x,T_sg_n_32,'--c');
xlim([0 b]);
ylim([0 90]);
ylabel('$\Delta T$ [K]','Interpreter','Latex');
xlabel('$r$ [mm]','Interpreter','Latex');
legend([h1 h2 h3 h4 h5],{'SG (n=2)','SG (n=4)','SG (n=8)','SG (n=16)','SG (n=32)'},'Interpreter','Latex');
hold off

figure(11)
h1=plot(x,T_sg_n_2,'k');hold on
h2=plot(x,T_sg_n_4,'--b');
h3=plot(x,T_sg_n_8,'-.g');
h4=plot(x,T_sg_n_16,'m');
h5=plot(x,T_sg_n_32,'--c');
ylim([40 90]);
ylabel('$\Delta T$ [K]','Interpreter','Latex');
xlabel('$r$ [mm]','Interpreter','Latex');
legend([h1 h2 h3 h4 h5],{'SG (n=2)','SG (n=4)','SG (n=8)','SG (n=16)','SG (n=32)'},'Interpreter','Latex');
hold off

figure(12)
%Top hat 
h1=plot(x_1,T_1(x_1),'k');hold on
plot(x_2,T_2(x_2),'k');
%Gaussian
h2=plot(x,T_g,'-.r');
%Super-Gaussian n=4
h3=plot(x,T_sg_n_4,'b');
%Super-Gaussian n=8
h4=plot(x,T_sg_n_8,'g');
%Donut
Td_1=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*((d_2^2-d_1^2+d_2^2*log(b^2/d_2^2)-d_1^2*log(b^2/d_1^2))/(d_2^2-d_1^2))+(T_c+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c+eta*alpha*P./(2*pi*b.*h));
Td_2=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*((d_2^2-xd_2.^2+d_2^2*log(b^2/d_2^2)-d_1^2*log(b^2./xd_2.^2))/(d_2^2-d_1^2))+(T_c+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c+eta*alpha*P./(2*pi*b.*h));
Td_3=(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*log(b^2./xd_3.^2)+(T_c+eta*alpha*P/(2*pi*b*h))^(m_rt+1)).^(1/(m_rt+1))-(T_c+eta*alpha*P./(2*pi*b.*h));
h5=plot(xd_1,ones(size(xd_1))*Td_1,'--m');
plot(xd_2,Td_2,'--m');
plot(xd_3,Td_3,'--m');

xlim([0 b]);
ylabel('$\Delta T(r)$ [K]','Interpreter','Latex');
xlabel('$r$ [mm]','Interpreter','Latex');
legend([h1 h2 h3 h4 h5 ],{'TH','G','SG (n=4)','SG (n=8)','D'},'Interpreter','Latex');
hold off

figure(13)
h_r=0.002:1e-5:0.20; %units WK/mm^-2
T_h=@(h)(eta*alpha*P*T_0^m_rt*(m_rt+1)/(4*k_0*pi)*(1+log(b^2/a^2))+(T_c+eta*alpha*P./(2*pi*b.*h)).^(m_rt+1)).^(1/(m_rt+1));
semilogx(h_r*100,T_h(h_r),'k');
ylabel('$T_1(0,0)$ [K]','Interpreter','Latex');
xlabel('$h$ [WK$^{-1}$cm$^{-2}$]','Interpreter','Latex');
xlim([0.2 20]);
New_XTickLabel = get(gca,'xtick');
set(gca,'XTickLabel',New_XTickLabel);
zoomH = zoom(gcf); 
set(zoomH,'ActionPostCallback',{@zoom_mypostcallback}); 

figure(14)
A=eta*alpha*P/(2*pi*b);
B=A*T_0^m_rt*(m_rt+1)*b/(2*k_0)*(1+log(b^2/a^2));
dT=@(h_r)-A.*(T_c+A./h_r).^m_rt.*(B+(T_c+A./h_r).^(m_rt+1)).^(-m_rt/(m_rt+1))./h_r.^2/100;
semilogx(h_r*100,dT(h_r),'k');hold on
ylabel('$dT_1(0,0)/dh$ [K$^2$cm$^2$W$^{-1}$]','Interpreter','Latex');
xlabel('$h$ [WK$^{-1}$cm$^{-2}$]','Interpreter','Latex');
xlim([0.2 20]);
New_XTickLabel = get(gca,'xtick');
set(gca,'XTickLabel',New_XTickLabel);
zoomH = zoom(gcf); 
set(zoomH,'ActionPostCallback',{@zoom_mypostcallback}); 
hold off

figure(15)
C_1=eta*P*alpha*T_0^m_rt*(m_rt+1)/(4*pi*k_0)*(1+log(b^2/a^2))+(T_c)^(m_rt+1);
C_2=eta*P*alpha*T_0^m_rt*(m_rt+1)/(4*pi*k_0*a^2);
T_1=(C_1-C_2.*x_1.^2).^(1/(m_rt+1));
T_1_T_f=@(x)C_1^(1/(m_rt+1))-C_2*C_1^(-m_rt/(m_rt+1)).*x.^2/(m_rt+1);
T_1_T_s=@(x)C_1^(1/(m_rt+1))-C_2*C_1^(-m_rt/(m_rt+1)).*x.^2/(m_rt+1)-C_2^2*m_rt*C_1^((-1-2*m_rt)/(m_rt+1))/(2*(m_rt+1)^2).*x.^4;
T_2=(eta*P*alpha*T_0^m_rt*(m_rt+1)/(4*pi*k_0)*log(b^2./x_2.^2)+(T_c)^(m_rt+1)).^(1/(m_rt+1));
x_3=0:a/500:a;
h1=plot(x_1,T_1,'--k');hold on
h2=plot(x_3,T_1_T_f(x_3),'r');
plot(x_2,T_2,'--k');
xlabel('r [mm]','Interpreter','Latex');
ylabel('T(r) [K]','Interpreter','Latex');
legend([h1 h2],{'Analytical solution','Quadratic Taylor expansion'},'Interpreter','Latex');
xlim([0 b/3]);
hold off 

figure(16)
l0_ct=-2.3*10^-6;
l1_ct=3.8*10^-8;
l2_ct=0.18*10^-6;
l3_ct=2.1*10^-8;

l0_rt=3.4*10^-6;
l1_rt=1.7*10^-8;
l2_rt=4.2*10^-6;
l3_rt=0.62*10^-8;

chi_0_ct=l0_ct+l2_ct;
chi_1_ct=l1_ct+l3_ct;
chi_0_rt=l0_rt+l2_rt;
chi_1_rt=l1_rt+l3_rt;

T_c_ct=40:0.1:175;
T_c_rt=225:0.1:475;

B_ct=eta*P*T_0^m_ct*(m_ct+1)*alpha/(4*pi*k_0)*(1+log(b^2/a^2));
B_rt=eta*P*T_0^m_rt*(m_rt+1)*alpha/(4*pi*k_0)*(1+log(b^2/a^2));
D_th_lim_ct=@(T_c)2/(a^2*alpha*(1+log(b^2/a^2)))*(chi_0_ct*(B_ct+T_c.^(m_ct+1)).^(1/(m_ct+1))+chi_1_ct*(B_ct+T_c.^(m_ct+1)).^(2/(m_ct+1))-chi_0_ct*(B_ct*exp(-alpha*L)+T_c.^(m_ct+1)).^(1/(m_ct+1))-chi_1_ct*(B_ct*exp(-alpha*L)+T_c.^(m_ct+1)).^(2/(m_ct+1)))*1e3;
D_th_lim_rt=@(T_c)2/(a^2*alpha*(1+log(b^2/a^2)))*(chi_0_rt*(B_rt+T_c.^(m_rt+1)).^(1/(m_rt+1))+chi_1_rt*(B_rt+T_c.^(m_rt+1)).^(2/(m_rt+1))-chi_0_rt*(B_rt*exp(-alpha*L)+T_c.^(m_rt+1)).^(1/(m_rt+1))-chi_1_rt*(B_rt*exp(-alpha*L)+T_c.^(m_rt+1)).^(2/(m_rt+1)))*1e3;

hb = area([125 300], [110 110]);hold on
hb.FaceColor = [0.95 0.95 0.95];

h1=plot(T_c_ct,D_th_lim_ct(T_c_ct),'k');
h2=plot(T_c_rt,D_th_lim_rt(T_c_rt),'--k');

n=2;
G_ct=2^(2/n-3)*n*eta*P*alpha*T_0^m_ct*(m_ct+1)/(pi*k_0*w^2*gamma(2/n));
G_rt=2^(2/n-3)*n*eta*P*alpha*T_0^m_rt*(m_rt+1)/(pi*k_0*w^2*gamma(2/n));
F=hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*b^n/w^n);
D_sg_n_ct=@(T_c)2/(alpha*b^2*F)*(chi_0_ct*(G_ct*b^2*F+T_c.^(m_ct+1)).^(1/(m_ct+1))+chi_1_ct*(G_ct*b^2*F+T_c.^(m_ct+1)).^(2/(m_ct+1))-chi_0_ct*(G_ct*b^2*F*exp(-alpha*L)+T_c.^(m_ct+1)).^(1/(m_ct+1))-chi_1_ct*(G_ct*b^2*F*exp(-alpha*L)+T_c.^(m_ct+1)).^(2/(m_ct+1)))*1e3;
D_sg_n_rt=@(T_c)2/(alpha*b^2*F)*(chi_0_rt*(G_rt*b^2*F+T_c.^(m_rt+1)).^(1/(m_rt+1))+chi_1_rt*(G_rt*b^2*F+T_c.^(m_rt+1)).^(2/(m_rt+1))-chi_0_rt*(G_rt*b^2*F*exp(-alpha*L)+T_c.^(m_rt+1)).^(1/(m_rt+1))-chi_1_rt*(G_rt*b^2*F*exp(-alpha*L)+T_c.^(m_rt+1)).^(2/(m_rt+1)))*1e3;
h3=plot(T_c_ct,D_sg_n_ct(T_c_ct),'r');
h4=plot(T_c_rt,D_sg_n_rt(T_c_rt),'--r');

n=4;
G_ct=2^(2/n-3)*n*eta*P*alpha*T_0^m_ct*(m_ct+1)/(pi*k_0*w^2*gamma(2/n));
G_rt=2^(2/n-3)*n*eta*P*alpha*T_0^m_rt*(m_rt+1)/(pi*k_0*w^2*gamma(2/n));
F=hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*b^n/w^n);
D_sg_n_ct=@(T_c)2/(alpha*b^2*F)*(chi_0_ct*(G_ct*b^2*F+T_c.^(m_ct+1)).^(1/(m_ct+1))+chi_1_ct*(G_ct*b^2*F+T_c.^(m_ct+1)).^(2/(m_ct+1))-chi_0_ct*(G_ct*b^2*F*exp(-alpha*L)+T_c.^(m_ct+1)).^(1/(m_ct+1))-chi_1_ct*(G_ct*b^2*F*exp(-alpha*L)+T_c.^(m_ct+1)).^(2/(m_ct+1)))*1e3;
D_sg_n_rt=@(T_c)2/(alpha*b^2*F)*(chi_0_rt*(G_rt*b^2*F+T_c.^(m_rt+1)).^(1/(m_rt+1))+chi_1_rt*(G_rt*b^2*F+T_c.^(m_rt+1)).^(2/(m_rt+1))-chi_0_rt*(G_rt*b^2*F*exp(-alpha*L)+T_c.^(m_rt+1)).^(1/(m_rt+1))-chi_1_rt*(G_rt*b^2*F*exp(-alpha*L)+T_c.^(m_rt+1)).^(2/(m_rt+1)))*1e3;

h5=plot(T_c_ct,D_sg_n_ct(T_c_ct),'b');
h6=plot(T_c_rt,D_sg_n_rt(T_c_rt),'--b');

n=8;
G_ct=2^(2/n-3)*n*eta*P*alpha*T_0^m_ct*(m_ct+1)/(pi*k_0*w^2*gamma(2/n));
G_rt=2^(2/n-3)*n*eta*P*alpha*T_0^m_rt*(m_rt+1)/(pi*k_0*w^2*gamma(2/n));
F=hypergeom([2/n,2/n],[2/n+1,2/n+1],-2*b^n/w^n);
D_sg_n_ct=@(T_c)2/(alpha*b^2*F)*(chi_0_ct*(G_ct*b^2*F+T_c.^(m_ct+1)).^(1/(m_ct+1))+chi_1_ct*(G_ct*b^2*F+T_c.^(m_ct+1)).^(2/(m_ct+1))-chi_0_ct*(G_ct*b^2*F*exp(-alpha*L)+T_c.^(m_ct+1)).^(1/(m_ct+1))-chi_1_ct*(G_ct*b^2*F*exp(-alpha*L)+T_c.^(m_ct+1)).^(2/(m_ct+1)))*1e3;
D_sg_n_rt=@(T_c)2/(alpha*b^2*F)*(chi_0_rt*(G_rt*b^2*F+T_c.^(m_rt+1)).^(1/(m_rt+1))+chi_1_rt*(G_rt*b^2*F+T_c.^(m_rt+1)).^(2/(m_rt+1))-chi_0_rt*(G_rt*b^2*F*exp(-alpha*L)+T_c.^(m_rt+1)).^(1/(m_rt+1))-chi_1_rt*(G_rt*b^2*F*exp(-alpha*L)+T_c.^(m_rt+1)).^(2/(m_rt+1)))*1e3;

h7=plot(T_c_ct,D_sg_n_ct(T_c_ct),'g');
h8=plot(T_c_rt,D_sg_n_rt(T_c_rt),'--g');
xlim([40 470]); 
ylim([0 110]);
%filled area
ha = area([150 250], [110 110]);
ha.FaceColor = [0.9 0.9 0.9];

xlabel('$T_c$ [K]','Interpreter','Latex');
ylabel('$D_{th}$ [m$^{-1}$]','Interpreter','Latex');
legend([h1 h2 h3 h4 h5 h6 h7 h8],{'TH CT','TH RT','G CT','G RT','SG (n=4) CT','SG (n=4) RT','SG (n=8) CT','SG (n=8) RT'},'Interpreter','Latex');