%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Chapter 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         For plotting the graphs and calculating some complex            %
%               expressions represented in Chapter 4                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Sound velocity, derivative dV/dx
clear all
clc

syms c_A A r x C_1 C_2 omega

chi=sqrt(c_A*((A^2-r^2/c_A^2)*x^2+2*A*x+1));
%pretty(chi)

zak=sqrt(omega^2-r^2); 
phi=A+(A^2-r^2/c_A^2)*x;
nu=atanh(c_A*phi/abs(r))-atanh(c_A*A/abs(r));
%pretty(phi)
%pretty(nu)

g=C_1*cos(zak*nu/abs(r))-C_2*sin(zak*nu/abs(r));
%pretty(g)

%pretty(diff(chi,x))
%pretty(diff(g,x))

v=chi*g;
diff_V=diff(v);
%pretty(diff_V)
pretty(simplify(diff_V))

%M1=subs(diff_V,x,0)

%% Figure 7
% graph 7.1, reflection coefficient

clc
clear all

syms omega;

L=1;
kappa=7/5;pl=287.058;

c_A=345;
r=abs(345*sqrt(3));
A=1;
%
T_A=c_A^2/(kappa*pl);
T_B=T_A*(A^2-r^2/c_A^2+2*A+1)^2;
c_B=sqrt(kappa*pl*T_B);

len=sqrt(omega^2-r^2);
k_A=omega/c_A;
k_B=omega/c_B;

phi=A+(A^2-r^2/c_A^2)*L;
chi=sqrt(c_A*((A^2-r^2/c_A^2)*L^2+2*A*L+1));
nu=atanh(c_A*phi/abs(r))-atanh(c_A*A/abs(r));
zak=(A^2-r^2/c_A^2)*len*c_A*chi/(c_A^2*phi^2-r^2);
kap=c_A*phi/chi;

G1=sqrt(c_A);
H1=chi*cos(len/r*nu);
H2=-chi*sin(len/r*nu);
M1=sqrt(c_A)*A;
M2=len/sqrt(c_A);
N1=zak*sin(len*nu/abs(r))+kap*cos(len*nu/abs(r));
N2=zak*cos(len*nu/abs(r))-kap*sin(len*nu/abs(r));

R=(k_A*k_B*G1*H2-M1*N2+M2*N1+1i*(k_A*G1*N2-k_B*(H1*M2-H2*M1)))/(k_A*k_B*G1*H2+M1*N2-M2*N1+1i*(k_A*G1*N2+k_B*(H1*M2-H2*M1)));

magnitude=abs(R);

x2=[r:10:7000];
RR2=subs(magnitude,'omega',x2);
hold on
p2=plot(x2,RR2);
hold off
set(gca, 'FontSize', 13)
title('Reflection coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|R^{(v)}| (-)','FontWeight','bold') 
set(p2,'LineWidth',2.5)

matlab2tikz('fig7_1.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},',...
                       'legend style={font=\scshape},']);
                   
%% graph 7.2, transmission coefficient

clc
clear all

syms omega;

L=1;
kappa=7/5;pl=287.058;

c_A=345;
r=abs(345*sqrt(3));
A=1;
%
T_A=c_A^2/(kappa*pl);
T_B=T_A*(A^2-r^2/c_A^2+2*A+1)^2;
c_B=sqrt(kappa*pl*T_B);

len=sqrt(omega^2-r^2);
k_A=omega/c_A;
k_B=omega/c_B;

phi=A+(A^2-r^2/c_A^2)*L;
chi=sqrt(c_A*((A^2-r^2/c_A^2)*L^2+2*A*L+1));
nu=atanh(c_A*phi/abs(r))-atanh(c_A*A/abs(r));
zak=(A^2-r^2/c_A^2)*len*c_A*chi/(c_A^2*phi^2-r^2);
kap=c_A*phi/chi;

G1=sqrt(c_A);
H1=chi*cos(len/r*nu);
H2=-chi*sin(len/r*nu);
M1=sqrt(c_A)*A;
M2=len/sqrt(c_A);
N1=zak*sin(len*nu/abs(r))+kap*cos(len*nu/abs(r));
N2=zak*cos(len*nu/abs(r))-kap*sin(len*nu/abs(r));

Tr=(1i*2*k_A*(H1*N2-H2*N1))/(k_A*k_B*G1*H2+M1*N2-M2*N1+1i*(k_A*G1*N2+k_B*(H1*M2-H2*M1)));

magnitude=abs(Tr);

x2=[r:10:7000];
TTr2=subs(magnitude,'omega',x2);
hold on
p2=plot(x2,TTr2);
hold off
set(gca, 'FontSize', 13)
title('Transmission coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|Tr^{(v)}| (-)','FontWeight','bold') 
set(p2,'LineWidth',2.5)

matlab2tikz('fig7_2.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},',...
                       'legend style={font=\scshape},']);                   
%% Figure 8
% graph 8.1, reflection coefficient

clc
clear all

syms omega;

L=1;
kappa=7/5;pl=287.058;

c_A=345;
r=172.5;
A=-0.008;
%
T_A=c_A^2/(kappa*pl);
T_B=T_A*(A^2-r^2/c_A^2+2*A+1)^2;
c_B=sqrt(kappa*pl*T_B);

len=sqrt(omega^2-r^2);
k_A=omega/c_A;
k_B=omega/c_B;

phi=A+(A^2-r^2/c_A^2)*L;
chi=sqrt(c_A*((A^2-r^2/c_A^2)*L^2+2*A*L+1));
nu=atanh(c_A*phi/abs(r))-atanh(c_A*A/abs(r));
zak=(A^2-r^2/c_A^2)*len*c_A*chi/(c_A^2*phi^2-r^2);
kap=c_A*phi/chi;

G1=sqrt(c_A);
H1=chi*cos(len/r*nu);
H2=-chi*sin(len/r*nu);
M1=sqrt(c_A)*A;
M2=len/sqrt(c_A);
N1=zak*sin(len*nu/abs(r))+kap*cos(len*nu/abs(r));
N2=zak*cos(len*nu/abs(r))-kap*sin(len*nu/abs(r));

R=(k_A*k_B*G1*H2-M1*N2+M2*N1+1i*(k_A*G1*N2-k_B*(H1*M2-H2*M1)))/(k_A*k_B*G1*H2+M1*N2-M2*N1+1i*(k_A*G1*N2+k_B*(H1*M2-H2*M1)));

magnitude=abs(R);

x2=[r:10:7000];
RR2=subs(magnitude,'omega',x2);
hold on
p2=plot(x2,RR2);
hold off
set(gca, 'FontSize', 13)
title('Reflection coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|R^{(v)}| (-)','FontWeight','bold') 
set(p2,'LineWidth',2.5)

matlab2tikz('fig8_1.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},',...
                       'legend style={font=\scshape},']);
                   
%% graph 8.2, transmission coefficient

clc
clear all

syms omega;

L=1;
kappa=7/5;pl=287.058;

c_A=345;
r=172.5;
A=-0.008;
%
T_A=c_A^2/(kappa*pl);
T_B=T_A*(A^2-r^2/c_A^2+2*A+1)^2;
c_B=sqrt(kappa*pl*T_B);

len=sqrt(omega^2-r^2);
k_A=omega/c_A;
k_B=omega/c_B;

phi=A+(A^2-r^2/c_A^2)*L;
chi=sqrt(c_A*((A^2-r^2/c_A^2)*L^2+2*A*L+1));
nu=atanh(c_A*phi/abs(r))-atanh(c_A*A/abs(r));
zak=(A^2-r^2/c_A^2)*len*c_A*chi/(c_A^2*phi^2-r^2);
kap=c_A*phi/chi;

G1=sqrt(c_A);
H1=chi*cos(len/r*nu);
H2=-chi*sin(len/r*nu);
M1=sqrt(c_A)*A;
M2=len/sqrt(c_A);
N1=zak*sin(len*nu/abs(r))+kap*cos(len*nu/abs(r));
N2=zak*cos(len*nu/abs(r))-kap*sin(len*nu/abs(r));

Tr=(1i*2*k_A*(H1*N2-H2*N1))/(k_A*k_B*G1*H2+M1*N2-M2*N1+1i*(k_A*G1*N2+k_B*(H1*M2-H2*M1)));

magnitude=abs(Tr);

x2=[r:10:7000];
TTr2=subs(magnitude,'omega',x2);
hold on
p2=plot(x2,TTr2);
hold off
set(gca, 'FontSize', 13)
title('Transmission coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|Tr^{(v)}| (-)','FontWeight','bold') 
set(p2,'LineWidth',2.5)

matlab2tikz('fig8_2.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},',...
                       'legend style={font=\scshape},']);                   
%% Figure 9
% graph 9.1, reflection coefficient

clc
clear all

syms omega;

L=1;
kappa=7/5;pl=287.058;

c_A=345;
r=abs(345*sqrt(10.25-3*sqrt(10)));
A=-2.5+sqrt(10);
%
T_A=c_A^2/(kappa*pl);
T_B=T_A*(A^2-r^2/c_A^2+2*A+1)^2;
c_B=sqrt(kappa*pl*T_B);

len=sqrt(omega^2-r^2);
k_A=omega/c_A;
k_B=omega/c_B;

phi=A+(A^2-r^2/c_A^2)*L;
chi=sqrt(c_A*((A^2-r^2/c_A^2)*L^2+2*A*L+1));
nu=atanh(c_A*phi/abs(r))-atanh(c_A*A/abs(r));
zak=(A^2-r^2/c_A^2)*len*c_A*chi/(c_A^2*phi^2-r^2);
kap=c_A*phi/chi;

G1=sqrt(c_A);
H1=chi*cos(len/r*nu);
H2=-chi*sin(len/r*nu);
M1=sqrt(c_A)*A;
M2=len/sqrt(c_A);
N1=zak*sin(len*nu/abs(r))+kap*cos(len*nu/abs(r));
N2=zak*cos(len*nu/abs(r))-kap*sin(len*nu/abs(r));

R=(k_A*k_B*G1*H2-M1*N2+M2*N1+1i*(k_A*G1*N2-k_B*(H1*M2-H2*M1)))/(k_A*k_B*G1*H2+M1*N2-M2*N1+1i*(k_A*G1*N2+k_B*(H1*M2-H2*M1)));

magnitude=abs(R);

x2=[r:10:7000];
RR2=subs(magnitude,'omega',x2);
hold on
p2=plot(x2,RR2);
hold off
set(gca, 'FontSize', 13)
title('Reflection coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|R^{(v)}| (-)','FontWeight','bold') 
set(p2,'LineWidth',2.5)

matlab2tikz('fig9_1.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},',...
                       'legend style={font=\scshape},']);
                   
%% graph 9.2, transmission coefficient

clc
clear all

syms omega;

L=1;
kappa=7/5;pl=287.058;

c_A=345;
r=abs(345*sqrt(10.25-3*sqrt(10)));
A=-2.5+sqrt(10);
%
T_A=c_A^2/(kappa*pl);
T_B=T_A*(A^2-r^2/c_A^2+2*A+1)^2;
c_B=sqrt(kappa*pl*T_B);

len=sqrt(omega^2-r^2);
k_A=omega/c_A;
k_B=omega/c_B;

phi=A+(A^2-r^2/c_A^2)*L;
chi=sqrt(c_A*((A^2-r^2/c_A^2)*L^2+2*A*L+1));
nu=atanh(c_A*phi/abs(r))-atanh(c_A*A/abs(r));
zak=(A^2-r^2/c_A^2)*len*c_A*chi/(c_A^2*phi^2-r^2);
kap=c_A*phi/chi;

G1=sqrt(c_A);
H1=chi*cos(len/r*nu);
H2=-chi*sin(len/r*nu);
M1=sqrt(c_A)*A;
M2=len/sqrt(c_A);
N1=zak*sin(len*nu/abs(r))+kap*cos(len*nu/abs(r));
N2=zak*cos(len*nu/abs(r))-kap*sin(len*nu/abs(r));

Tr=(1i*2*k_A*(H1*N2-H2*N1))/(k_A*k_B*G1*H2+M1*N2-M2*N1+1i*(k_A*G1*N2+k_B*(H1*M2-H2*M1)));

magnitude=abs(Tr);

x2=[r:10:7000];
TTr2=subs(magnitude,'omega',x2);
hold on
p2=plot(x2,TTr2);
hold off
set(gca, 'FontSize', 13)
title('Transmission coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|Tr^{(v)}| (-)','FontWeight','bold') 
set(p2,'LineWidth',2.5)

matlab2tikz('fig9_2.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},',...
                       'legend style={font=\scshape},']);                   
%% Figure 10
% graph 10.1, reflection coefficient

clc
clear all

syms omega;

L=1;
kappa=7/5;pl=287.058;

c_A=345;
r=abs(sqrt(33-24*sqrt(2)));
A=-2.5+2*sqrt(2);
%
T_A=c_A^2/(kappa*pl);
T_B=T_A*(A^2-r^2/c_A^2+2*A+1)^2;
c_B=sqrt(kappa*pl*T_B);

len=sqrt(omega^2-r^2);
k_A=omega/c_A;
k_B=omega/c_B;

phi=A+(A^2-r^2/c_A^2)*L;
chi=sqrt(c_A*((A^2-r^2/c_A^2)*L^2+2*A*L+1));
nu=atanh(c_A*phi/abs(r))-atanh(c_A*A/abs(r));
zak=(A^2-r^2/c_A^2)*len*c_A*chi/(c_A^2*phi^2-r^2);
kap=c_A*phi/chi;

G1=sqrt(c_A);
H1=chi*cos(len/r*nu);
H2=-chi*sin(len/r*nu);
M1=sqrt(c_A)*A;
M2=len/sqrt(c_A);
N1=zak*sin(len*nu/abs(r))+kap*cos(len*nu/abs(r));
N2=zak*cos(len*nu/abs(r))-kap*sin(len*nu/abs(r));

R=(k_A*k_B*G1*H2-M1*N2+M2*N1+1i*(k_A*G1*N2-k_B*(H1*M2-H2*M1)))/(k_A*k_B*G1*H2+M1*N2-M2*N1+1i*(k_A*G1*N2+k_B*(H1*M2-H2*M1)));

magnitude=abs(R);

x2=[r:10:7000];
RR2=subs(magnitude,'omega',x2);
hold on
p2=plot(x2,RR2);
hold off
set(gca, 'FontSize', 13)
title('Reflection coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|R^{(v)}| (-)','FontWeight','bold') 
set(p2,'LineWidth',2.5)

matlab2tikz('fig10_1.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},',...
                       'legend style={font=\scshape},']);
                   
%% graph 9.2, transmission coefficient

clc
clear all

syms omega;

L=1;
kappa=7/5;pl=287.058;

c_A=345;
r=abs(sqrt(33-24*sqrt(2)));
A=-2.5+2*sqrt(2);
%
T_A=c_A^2/(kappa*pl);
T_B=T_A*(A^2-r^2/c_A^2+2*A+1)^2;
c_B=sqrt(kappa*pl*T_B);

len=sqrt(omega^2-r^2);
k_A=omega/c_A;
k_B=omega/c_B;

phi=A+(A^2-r^2/c_A^2)*L;
chi=sqrt(c_A*((A^2-r^2/c_A^2)*L^2+2*A*L+1));
nu=atanh(c_A*phi/abs(r))-atanh(c_A*A/abs(r));
zak=(A^2-r^2/c_A^2)*len*c_A*chi/(c_A^2*phi^2-r^2);
kap=c_A*phi/chi;

G1=sqrt(c_A);
H1=chi*cos(len/r*nu);
H2=-chi*sin(len/r*nu);
M1=sqrt(c_A)*A;
M2=len/sqrt(c_A);
N1=zak*sin(len*nu/abs(r))+kap*cos(len*nu/abs(r));
N2=zak*cos(len*nu/abs(r))-kap*sin(len*nu/abs(r));

Tr=(1i*2*k_A*(H1*N2-H2*N1))/(k_A*k_B*G1*H2+M1*N2-M2*N1+1i*(k_A*G1*N2+k_B*(H1*M2-H2*M1)));

magnitude=abs(Tr);

x2=[r:10:7000];
TTr2=subs(magnitude,'omega',x2);
hold on
p2=plot(x2,TTr2);
hold off
set(gca, 'FontSize', 13)
title('Transmission coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|Tr^{(v)}| (-)','FontWeight','bold') 
set(p2,'LineWidth',2.5)

matlab2tikz('fig10_2.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},',...
                       'legend style={font=\scshape},']);                   

%% Sound pressure                   
% second exact analytical solution, derivative dP/dx 
clear all
clc

syms k_A q x C_1 C_2

f=1/sqrt(q*x+1);
%pretty(f)

g=C_1*cos((sqrt(1/4-k_A^2/q^2))*(log(q*x+1)))+C_2*sin((sqrt(1/4-k_A^2/q^2))*(log(q*x+1)));
%pretty(g)

%pretty(diff(f,x))
%pretty(diff(g,x))

p=f*g;
diff_P=diff(p);
%pretty(diff_P)
pretty(simplify(diff_P))

%% Figure 12
% graph 12.1, reflection coefficient delta<0

clc
clear all

syms q omega

L=1;
q1=1;
q2=2;
q3=3;
kappa=7/5;pl=287.058;
T_A=296;T_B=T_A*(q+1)^2;

c_A=sqrt(kappa*pl*T_A)
c_B=sqrt(kappa*pl*T_B);
k_A=omega/c_A;
k_B=omega/c_B;

%<0
chi=sqrt(1/4-k_A^2/q^2)*log(q*L+1);
H1=1/sqrt(q*L+1)*cos(chi);
H2=1/sqrt(q*L+1)*sin(chi);
M1=-q/2;
M2=q*sqrt(1/4-k_A^2/q^2);
N1=-(q*(q*L+1)^(-3/2))*(0.5*cos(chi)+sqrt(1/4-k_A^2/q^2)*sin(chi));
N2=-(q*(q*L+1)^(-3/2))*(0.5*sin(chi)-sqrt(1/4-k_A^2/q^2)*cos(chi));

R=(k_A*k_B*H2-M1*N2+M2*N1+1i*(k_A*N2+k_B*(-H1*M2+H2*M1)))/(k_A*k_B*H2+M1*N2-M2*N1+1i*(k_A*N2-k_B*(-H1*M2+H2*M1)));
magnitude=abs(R);

sub1=subs(magnitude,'q',q1);
sub2=subs(magnitude,'q',q2);
sub3=subs(magnitude,'q',q3);

x1=[c_A*q1/2+0.0001:10:5001];
x2=[c_A*q2/2+0.0001:10:5001];
x3=[c_A*q3/2+0.0001:10:5001];

RR1=subs(sub1,'omega',x1)
RR2=subs(sub2,'omega',x2)
RR3=subs(sub3,'omega',x3)

hold on
p1=plot(x1,RR1);
p2=plot(x2,RR2,'g');
p3=plot(x3,RR3,'r');
hold off
legend({'$q$ = 1 m$^{-1}$','$q$ = 2 m$^{-1}$','$q$ = 3 m$^{-1}$'},'interpreter','latex','fontsize',17,'FontWeight','bold')

set(gca, 'FontSize', 13)
title('Reflection coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|R^{(p)}| (-)','FontWeight','bold') 
set(p1,'LineWidth',2.5)
set(p2,'LineWidth',2.5)
set(p3,'LineWidth',2.5)

matlab2tikz('fig12_1.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},']);
           
%% graph 12.2, transmission coefficient delta<0

clc
clear all

syms q omega

L=1;
q1=1;
q2=2;
q3=3;
kappa=7/5;pl=287.058;
T_A=296;T_B=T_A*(q+1)^2;

c_A=sqrt(kappa*pl*T_A)
c_B=sqrt(kappa*pl*T_B);
k_A=omega/c_A;
k_B=omega/c_B;

%<0
chi=sqrt(1/4-k_A^2/q^2)*log(q*L+1);
H1=1/sqrt(q*L+1)*cos(chi);
H2=1/sqrt(q*L+1)*sin(chi);
M1=-q/2;
M2=q*sqrt(1/4-k_A^2/q^2);
N1=-(q*(q*L+1)^(-3/2))*(0.5*cos(chi)+sqrt(1/4-k_A^2/q^2)*sin(chi));
N2=-(q*(q*L+1)^(-3/2))*(0.5*sin(chi)-sqrt(1/4-k_A^2/q^2)*cos(chi));

Tr=(1i*2*k_A*(H1*N2-H2*N1))/(k_A*k_B*H2+M1*N2-M2*N1+1i*(k_A*N2-k_B*(-H1*M2+H2*M1)));
magnitude=abs(Tr);

sub1=subs(magnitude,'q',q1);
sub2=subs(magnitude,'q',q2);
sub3=subs(magnitude,'q',q3);

x1=[c_A*q1/2+0.0001:10:5001];
x2=[c_A*q2/2+0.0001:10:5001];
x3=[c_A*q3/2+0.0001:10:5001];

RR1=subs(sub1,'omega',x1)
RR2=subs(sub2,'omega',x2)
RR3=subs(sub3,'omega',x3)

hold on
p1=plot(x1,RR1);
p2=plot(x2,RR2,'g');
p3=plot(x3,RR3,'r');
hold off
legend({'$q$ = 1 m$^{-1}$','$q$ = 2 m$^{-1}$','$q$ = 3 m$^{-1}$'},'interpreter','latex','fontsize',17,'FontWeight','bold')

set(gca, 'FontSize', 13)
title('Transmission coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|Tr^{(p)}| (-)','FontWeight','bold') 
set(p1,'LineWidth',2.5)
set(p2,'LineWidth',2.5)
set(p3,'LineWidth',2.5)

matlab2tikz('fig12_2.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},']);
                   
%% Linear temperature distribution, derivative dP/dx   

clear all
clc

syms T_A A x omega C_1 C_2 b

f=omega*sqrt(T_A*(A*x+1))/b;
%pretty(f)

g=C_1*besselj(0,f)+C_2*bessely(0,f);
%pretty(g)

%pretty(diff(f,x))
%pretty(diff(g,x))

diff_P=diff(g);
%pretty(diff_P)
pretty(simplify(diff_P))

%% Figure 13
% graph 13.1, reflection coefficient for A>0

clc
clear all

syms A omega;

kappa=7/5;pl=287.058;
L=1;
T_A=296;
A1=456/T_A;
A2=238/T_A;
A3=110/T_A;

T_B=T_A*(A+1);
c_A=sqrt(kappa*pl*T_A);c_B=sqrt(kappa*pl*T_B);
b=abs(T_A*abs(A))/2*sqrt(kappa*pl);
k_A=omega/c_A;
k_B=omega/c_B;

G1=besselj(0,omega/b*sqrt(T_A));
G2=bessely(0,omega/b*sqrt(T_A));
H1=besselj(0,omega/b*sqrt(T_A*(A*L+1)));
H2=bessely(0,omega/b*sqrt(T_A*(A*L+1)));
M1=-omega*A/(2*b)*sqrt(T_A)*besselj(1,omega/b*sqrt(T_A));
M2=-omega*A/(2*b)*sqrt(T_A)*bessely(1,omega/b*sqrt(T_A));
N1=-omega*A/(2*b)*sqrt(T_A/(A*L+1))*besselj(1,omega/b*sqrt(T_A*(A*L+1)));
N2=-omega*A/(2*b)*sqrt(T_A/(A*L+1))*bessely(1,omega/b*sqrt(T_A*(A*L+1)));


R=(k_A*k_B*(G1*H2-G2*H1)-M1*N2+M2*N1+1i*(k_A*(G1*N2-G2*N1)-k_B*(H1*M2-H2*M1)))/(k_A*k_B*(G1*H2-G2*H1)+M1*N2-M2*N1+1i*(k_A*(G1*N2-G2*N1)+k_B*(H1*M2-H2*M1)));

magnitude=abs(R);

sub1=subs(magnitude,'A',A1);
sub2=subs(magnitude,'A',A2);
sub3=subs(magnitude,'A',A3);

x=0:10:7000;

RR1=subs(sub1,'omega',x)
RR2=subs(sub2,'omega',x)
RR3=subs(sub3,'omega',x)

hold on
p1=plot(x,RR1);
p2=plot(x,RR2,'g');
p3=plot(x,RR3,'r');
hold off
legend({'$A$ = $456/T_A$ m$^{-1}$','$A$ = $238/T_A$ m$^{-1}$','$A$ = $110/T_A$ m$^{-1}$'},'interpreter','latex','fontsize',15,'FontWeight','bold')

set(gca, 'FontSize', 13)
title('Reflection coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|R^{(p)}| (-)','FontWeight','bold') 
set(p1,'LineWidth',2.5)
set(p2,'LineWidth',2.5)
set(p3,'LineWidth',2.5)

matlab2tikz('fig13_1.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},']);

%% graph 13.2, transmission coefficient for A>0

clc
clear all

syms A omega

kappa=7/5;pl=287.058;
L=1;
T_A=296;
A1=456/T_A;
A2=238/T_A;
A3=110/T_A;

T_B=T_A*(A+1);
c_A=sqrt(kappa*pl*T_A);c_B=sqrt(kappa*pl*T_B);
b=abs(T_A*abs(A))/2*sqrt(kappa*pl);
k_A=omega/c_A;
k_B=omega/c_B;

G1=besselj(0,omega/b*sqrt(T_A));
G2=bessely(0,omega/b*sqrt(T_A));
H1=besselj(0,omega/b*sqrt(T_A*(A*L+1)));
H2=bessely(0,omega/b*sqrt(T_A*(A*L+1)));
M1=-omega*A/(2*b)*sqrt(T_A)*besselj(1,omega/b*sqrt(T_A));
M2=-omega*A/(2*b)*sqrt(T_A)*bessely(1,omega/b*sqrt(T_A));
N1=-omega*A/(2*b)*sqrt(T_A/(A*L+1))*besselj(1,omega/b*sqrt(T_A*(A*L+1)));
N2=-omega*A/(2*b)*sqrt(T_A/(A*L+1))*bessely(1,omega/b*sqrt(T_A*(A*L+1)));


Tr=(1i*2*k_A*(H1*N2-H2*N1))/(k_A*k_B*(G1*H2-G2*H1)+M1*N2-M2*N1+1i*(k_A*(G1*N2-G2*N1)+k_B*(H1*M2-H2*M1)));

magnitude=abs(Tr);

sub1=subs(magnitude,'A',A1);
sub2=subs(magnitude,'A',A2);
sub3=subs(magnitude,'A',A3);

x=0:10:7000;

RR1=subs(sub1,'omega',x);
RR2=subs(sub2,'omega',x);
RR3=subs(sub3,'omega',x);

hold on
p1=plot(x,RR1);
p2=plot(x,RR2,'g');
p3=plot(x,RR3,'r');
hold off 
legend({'$A$ = $456/T_A$ m$^{-1}$','$A$ = $238/T_A$ m$^{-1}$','$A$ = $110/T_A$ m$^{-1}$'},'interpreter','latex','Position',[0.742878337505775 0.657587425351634 0.154465592972182 0.139025725232622],'fontsize',15,'FontWeight','bold')

set(gca, 'FontSize', 13)
title('Transmission coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|Tr^{(p)}| (-)','FontWeight','bold') 
set(p1,'LineWidth',2.5)
set(p2,'LineWidth',2.5)
set(p3,'LineWidth',2.5)


matlab2tikz('fig13_2.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},']);

%% Figure 14   
% graph 14.1, reflection coefficient for A<0

clc
clear all

syms A omega

kappa=7/5;pl=287.058;
L=1;
T_A=752;
A1=-218/T_A;
A2=346/T_A;
A3=-456/T_A;

T_B=T_A*(A+1);
c_A=sqrt(kappa*pl*T_A);c_B=sqrt(kappa*pl*T_B);
b=abs(T_A*abs(A))/2*sqrt(kappa*pl);
k_A=omega/c_A;
k_B=omega/c_B;

G1=besselj(0,omega/b*sqrt(T_A));
G2=bessely(0,omega/b*sqrt(T_A));
H1=besselj(0,omega/b*sqrt(T_A*(A*L+1)));
H2=bessely(0,omega/b*sqrt(T_A*(A*L+1)));
M1=-omega*A/(2*b)*sqrt(T_A)*besselj(1,omega/b*sqrt(T_A));
M2=-omega*A/(2*b)*sqrt(T_A)*bessely(1,omega/b*sqrt(T_A));
N1=-omega*A/(2*b)*sqrt(T_A/(A*L+1))*besselj(1,omega/b*sqrt(T_A*(A*L+1)));
N2=-omega*A/(2*b)*sqrt(T_A/(A*L+1))*bessely(1,omega/b*sqrt(T_A*(A*L+1)));


R=(k_A*k_B*(G1*H2-G2*H1)-M1*N2+M2*N1+1i*(k_A*(G1*N2-G2*N1)-k_B*(H1*M2-H2*M1)))/(k_A*k_B*(G1*H2-G2*H1)+M1*N2-M2*N1+1i*(k_A*(G1*N2-G2*N1)+k_B*(H1*M2-H2*M1)));

magnitude=abs(R);

sub1=subs(magnitude,'A',A1);
sub2=subs(magnitude,'A',A2);
sub3=subs(magnitude,'A',A3);

x=0:10:7000;

RR1=subs(sub1,'omega',x);
RR2=subs(sub2,'omega',x);
RR3=subs(sub3,'omega',x);

hold on
p1=plot(x,RR1);
p2=plot(x,RR2,'g');
p3=plot(x,RR3,'r');
hold off
legend({'$A$ = $-218/T_A$ m$^{-1}$','$A$ = $-346/T_A$ m$^{-1}$','$A$ = $-456/T_A$ m$^{-1}$'},'interpreter','latex','fontsize',15,'FontWeight','bold')

set(gca, 'FontSize', 13)
title('Reflection coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|R^{(p)}| (-)','FontWeight','bold') 
set(p1,'LineWidth',2.5)
set(p2,'LineWidth',2.5)
set(p3,'LineWidth',2.5)


matlab2tikz('fig14_1.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},']);

%% graph 14.2, transmission coefficient for A<0

clc
clear all

syms A omega

kappa=7/5;pl=287.058;
L=1;
T_A=752;
A1=-218/T_A;
A2=346/T_A;
A3=-456/T_A;

T_B=T_A*(A+1);
c_A=sqrt(kappa*pl*T_A);c_B=sqrt(kappa*pl*T_B);
b=abs(T_A*abs(A))/2*sqrt(kappa*pl);
k_A=omega/c_A;
k_B=omega/c_B;

G1=besselj(0,omega/b*sqrt(T_A));
G2=bessely(0,omega/b*sqrt(T_A));
H1=besselj(0,omega/b*sqrt(T_A*(A*L+1)));
H2=bessely(0,omega/b*sqrt(T_A*(A*L+1)));
M1=-omega*A/(2*b)*sqrt(T_A)*besselj(1,omega/b*sqrt(T_A));
M2=-omega*A/(2*b)*sqrt(T_A)*bessely(1,omega/b*sqrt(T_A));
N1=-omega*A/(2*b)*sqrt(T_A/(A*L+1))*besselj(1,omega/b*sqrt(T_A*(A*L+1)));
N2=-omega*A/(2*b)*sqrt(T_A/(A*L+1))*bessely(1,omega/b*sqrt(T_A*(A*L+1)));


Tr=(1i*2*k_A*(H1*N2-H2*N1))/(k_A*k_B*(G1*H2-G2*H1)+M1*N2-M2*N1+1i*(k_A*(G1*N2-G2*N1)+k_B*(H1*M2-H2*M1)));

magnitude=abs(Tr);

sub1=subs(magnitude,'A',A1);
sub2=subs(magnitude,'A',A2);
sub3=subs(magnitude,'A',A3);

x=0:10:7000;

RR1=subs(sub1,'omega',x);
RR2=subs(sub2,'omega',x);
RR3=subs(sub3,'omega',x);

hold on
p1=plot(x,RR1);
p2=plot(x,RR2,'g');
p3=plot(x,RR3,'r');
hold off
legend({'$A$ = $-218/T_A$ m$^{-1}$','$A$ = $-346/T_A$ m$^{-1}$','$A$ = $-456/T_A$ m$^{-1}$'},'Position',[0.734042648772236 0.599354158841448 0.166178623718887 0.129856850715746],'interpreter','latex','fontsize',15,'FontWeight','bold')

set(gca, 'FontSize', 13)
title('Transmission coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|Tr^{(p)}| (-)','FontWeight','bold') 
set(p1,'LineWidth',2.5)
set(p2,'LineWidth',2.5)
set(p3,'LineWidth',2.5)


matlab2tikz('fig14_2.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},']);
                     
%% Exponential temperature distribution, derivative dP/dx 
clear all
clc

syms T_A gamma x omega delta C_1 C_2;

f=1/sqrt(T_A*exp(-gamma*x));
%pretty(f)

g=C_1*besselj(1,omega*delta*f)+C_2*bessely(1,omega*delta*f);
%pretty(g)

%pretty(diff(f,x))
%pretty(diff(g,x))

p=f*g;
diff_P=diff(p);
%pretty(diff_P)
pretty(simplify(diff_P))

%% Figure 15
% graph 15.1, reflection coefficient for gam>0

clc
clear all

syms gam omega

kappa=7/5;pl=287.058;
L=1;
gam1=-log(1.1);
gam2=-log(2);
gam3=-log(3);
T_A=296;T_B=T_A*exp(-gam);

c_A=sqrt(kappa*pl*T_A);c_B=sqrt(kappa*pl*T_B);
k_A=omega/c_A;
k_B=omega/c_B;
delta=2/sqrt(kappa*pl*gam^2);

G1=besselj(1,(omega*delta)/sqrt(T_A))/sqrt(T_A);
G2=bessely(1,(omega*delta)/sqrt(T_A))/sqrt(T_A);
H1=besselj(1,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(T_A*exp(-gam*L));
H2=bessely(1,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(T_A*exp(-gam*L));
M1=(gam*omega*delta)*besselj(0,(omega*delta)/sqrt(T_A))/(2*T_A);
M2=(gam*omega*delta)*bessely(0,(omega*delta)/sqrt(T_A))/(2*T_A);
N1=(gam*omega*delta)*besselj(0,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(2*T_A*exp(-gam*L));
N2=(gam*omega*delta)*bessely(0,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(2*T_A*exp(-gam*L));


R=(k_A*k_B*(G1*H2-G2*H1)-M1*N2+M2*N1+1i*(k_A*(G1*N2-G2*N1)-k_B*(H1*M2-H2*M1)))/(k_A*k_B*(G1*H2-G2*H1)+M1*N2-M2*N1+1i*(k_A*(G1*N2-G2*N1)+k_B*(H1*M2-H2*M1)));

magnitude=abs(R);

sub1=subs(magnitude,'gam',gam1);
sub2=subs(magnitude,'gam',gam2);
sub3=subs(magnitude,'gam',gam3);

x=0:10:7000;

RR1=subs(sub1,'omega',x);
RR2=subs(sub2,'omega',x);
RR3=subs(sub3,'omega',x);

hold on
p1=plot(x,RR1);
p2=plot(x,RR2,'g');
p3=plot(x,RR3,'r');
hold off


set(gca, 'FontSize', 13)
title('Reflection coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|R^{(p)}| (-)','FontWeight','bold') 
legend({'$\gamma$ = $-\ln{1.1}$ m$^{-1}$','$\gamma$ = $-\ln{2.0}$ m$^{-1}$','$\gamma$ = $-\ln{3.0}$ m$^{-1}$'},'interpreter','latex','fontsize',17,'FontWeight','bold')
set(p1,'LineWidth',2.5)
set(p2,'LineWidth',2.5)
set(p3,'LineWidth',2.5)


matlab2tikz('fig15_1.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},']);

%% graph 15.2, transmission coefficient for gam>0

clc
clear all

syms gam omega

kappa=7/5;pl=287.058;
L=1;
gam1=-log(1.1);
gam2=-log(2);
gam3=-log(3);
T_A=296;T_B=T_A*exp(-gam);

c_A=sqrt(kappa*pl*T_A);c_B=sqrt(kappa*pl*T_B);
k_A=omega/c_A;
k_B=omega/c_B;
delta=2/sqrt(kappa*pl*gam^2);

G1=besselj(1,(omega*delta)/sqrt(T_A))/sqrt(T_A);
G2=bessely(1,(omega*delta)/sqrt(T_A))/sqrt(T_A);
H1=besselj(1,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(T_A*exp(-gam*L));
H2=bessely(1,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(T_A*exp(-gam*L));
M1=(gam*omega*delta)*besselj(0,(omega*delta)/sqrt(T_A))/(2*T_A);
M2=(gam*omega*delta)*bessely(0,(omega*delta)/sqrt(T_A))/(2*T_A);
N1=(gam*omega*delta)*besselj(0,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(2*T_A*exp(-gam*L));
N2=(gam*omega*delta)*bessely(0,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(2*T_A*exp(-gam*L));


Tr=(1i*2*k_A*(H1*N2-H2*N1))/(k_A*k_B*(G1*H2-G2*H1)+M1*N2-M2*N1+1i*(k_A*(G1*N2-G2*N1)+k_B*(H1*M2-H2*M1)));

magnitude=abs(Tr);

sub1=subs(magnitude,'gam',gam1);
sub2=subs(magnitude,'gam',gam2);
sub3=subs(magnitude,'gam',gam3);

x=0:10:7000;

RR1=subs(sub1,'omega',x);
RR2=subs(sub2,'omega',x);
RR3=subs(sub3,'omega',x);

hold on
p1=plot(x,RR1);
p2=plot(x,RR2,'g');
p3=plot(x,RR3,'r');
hold off


set(gca, 'FontSize', 13)
title('Transmission coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|Tr^{(p)}| (-)','FontWeight','bold') 
legend({'$\gamma$ = $-\ln{1.1}$ m$^{-1}$','$\gamma$ = $-\ln{2.0}$ m$^{-1}$','$\gamma$ = $-\ln{3.0}$ m$^{-1}$'},'interpreter','latex','Position',[0.735067628808477 0.55343266140812 0.167642752562225 0.139059304703476],'fontsize',17,'FontWeight','bold')
set(p1,'LineWidth',2.5)
set(p2,'LineWidth',2.5)
set(p3,'LineWidth',2.5)


matlab2tikz('fig15_2.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},']);
                   
%% Figure 16
% graph 16.1, reflection coefficient for gam<0

clc
clear all

syms gam omega

kappa=7/5;pl=287.058;
L=1;
%gam=-log(0.5);
gam1=-log(0.8);
gam2=-log(0.5);
gam3=-log(0.1);
T_A=296;T_B=T_A*exp(-gam);

c_A=sqrt(kappa*pl*T_A);c_B=sqrt(kappa*pl*T_B);
k_A=omega/c_A;
k_B=omega/c_B;
delta=2/sqrt(kappa*pl*gam^2);

G1=besselj(1,(omega*delta)/sqrt(T_A))/sqrt(T_A);
G2=bessely(1,(omega*delta)/sqrt(T_A))/sqrt(T_A);
H1=besselj(1,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(T_A*exp(-gam*L));
H2=bessely(1,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(T_A*exp(-gam*L));
M1=(gam*omega*delta)*besselj(0,(omega*delta)/sqrt(T_A))/(2*T_A);
M2=(gam*omega*delta)*bessely(0,(omega*delta)/sqrt(T_A))/(2*T_A);
N1=(gam*omega*delta)*besselj(0,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(2*T_A*exp(-gam*L));
N2=(gam*omega*delta)*bessely(0,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(2*T_A*exp(-gam*L));


R=(k_A*k_B*(G1*H2-G2*H1)-M1*N2+M2*N1+1i*(k_A*(G1*N2-G2*N1)-k_B*(H1*M2-H2*M1)))/(k_A*k_B*(G1*H2-G2*H1)+M1*N2-M2*N1+1i*(k_A*(G1*N2-G2*N1)+k_B*(H1*M2-H2*M1)));

magnitude=abs(R);

sub1=subs(magnitude,'gam',gam1);
sub2=subs(magnitude,'gam',gam2);
sub3=subs(magnitude,'gam',gam3);

x=0:10:7000;

RR1=subs(sub1,'omega',x);
RR2=subs(sub2,'omega',x);
RR3=subs(sub3,'omega',x);

hold on
p1=plot(x,RR1);
p2=plot(x,RR2,'g');
p3=plot(x,RR3,'r');
hold off


set(gca, 'FontSize', 13)
title('Reflection coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|R^{(p)}| (-)','FontWeight','bold') 
legend({'$\gamma$ = $-\ln{0.8}$ m$^{-1}$','$\gamma$ = $-\ln{0.5}$ m$^{-1}$','$\gamma$ = $-\ln{0.1}$ m$^{-1}$'},'interpreter','latex','Position',[0.732926340375092 0.123239848086471 0.167642752562225 0.139059304703476],'fontsize',17,'FontWeight','bold')
set(p1,'LineWidth',2.5)
set(p2,'LineWidth',2.5)
set(p3,'LineWidth',2.5)


matlab2tikz('fig16_1.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},']);

%% graph 16.2, transmission coefficient for gam<0

clc
clear all

syms gam omega

kappa=7/5;pl=287.058;
L=1;
gam1=-log(0.8);
gam2=-log(0.5);
gam3=-log(0.1);
T_A=296;T_B=T_A*exp(-gam);

c_A=sqrt(kappa*pl*T_A);c_B=sqrt(kappa*pl*T_B);
k_A=omega/c_A;
k_B=omega/c_B;
delta=2/sqrt(kappa*pl*gam^2);

G1=besselj(1,(omega*delta)/sqrt(T_A))/sqrt(T_A);
G2=bessely(1,(omega*delta)/sqrt(T_A))/sqrt(T_A);
H1=besselj(1,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(T_A*exp(-gam*L));
H2=bessely(1,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(T_A*exp(-gam*L));
M1=(gam*omega*delta)*besselj(0,(omega*delta)/sqrt(T_A))/(2*T_A);
M2=(gam*omega*delta)*bessely(0,(omega*delta)/sqrt(T_A))/(2*T_A);
N1=(gam*omega*delta)*besselj(0,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(2*T_A*exp(-gam*L));
N2=(gam*omega*delta)*bessely(0,(omega*delta)/sqrt(T_A*exp(-gam*L)))/sqrt(2*T_A*exp(-gam*L));


Tr=(1i*2*k_A*(H1*N2-H2*N1))/(k_A*k_B*(G1*H2-G2*H1)+M1*N2-M2*N1+1i*(k_A*(G1*N2-G2*N1)+k_B*(H1*M2-H2*M1)));

magnitude=abs(Tr);

sub1=subs(magnitude,'gam',gam1);
sub2=subs(magnitude,'gam',gam2);
sub3=subs(magnitude,'gam',gam3);

x=0:10:7000;

RR1=subs(sub1,'omega',x);
RR2=subs(sub2,'omega',x);
RR3=subs(sub3,'omega',x);

hold on
p1=plot(x,RR1);
p2=plot(x,RR2,'g');
p3=plot(x,RR3,'r');
hold off


set(gca, 'FontSize', 13)
title('Transmission coefficient')
xlabel('\omega (rad/s)','FontWeight','bold')
ylabel('|Tr^{(p)}| (-)','FontWeight','bold') 
legend({'$\gamma$ = $-\ln{0.8}$ m$^{-1}$','$\gamma$ = $-\ln{0.5}$ m$^{-1}$','$\gamma$ = $-\ln{0.1}$ m$^{-1}$'},'interpreter','latex','fontsize',17,'FontWeight','bold')
set(p1,'LineWidth',2.5)
set(p2,'LineWidth',2.5)
set(p3,'LineWidth',2.5)


matlab2tikz('fig16_2.tikz', ...
   'extraaxisoptions',['title style={font=\Large},'...
                       'xlabel style={font=\Large},'...
                       'ylabel style={font=\Large},']);