%Coating thermal noise of LIGO 4masses (2 ETMs, 2 ITMs) with M4 of 40kg, a thickness of
%0.2m, a diameter of 0.35m and a beam radius of 0.065m

% modified: 150515 after comparison w/ John Miller's Mathematica code
% using "square-sum" method, not "averaging" method
% using 10 Bessel function terms, not 1
% for all 4 masses (2 ETMs, 2 ITMs)

function [total_coating_therm_final] = total_coating_therm_final(freq)

list_constants;


f=freq;  % frequency (Hz)
Kb= constants.kB; % Boltzmann constant (J/K)
T = T_mir; % temperature (K)
omega=2*pi*f;

m = M4;              % mass (kg)

asp_ratio = aspect_ratio; %35/20
rho = specific_weight;     % specific weight
sigma = mirror_sigma;      % w/2a=6.5/35

h = (4*m/(pi*asp_ratio^2*rho))^(1/3);   % mirror thickness (m)
a = asp_ratio*h/2;                      % mirror radius (m)
w = 2*a*sigma;                          % beam radius (m)

phi1 = phi1_SiO2;           % coating loss angle SiO2
phi2 = phi2_Ta205;          % coating loss angle Ta2O5
Y1 = Y1_SiO2;               % coating Young modulus SiO2 (Pa)
Y2 = Y2_Ta205;              % coating Young modulus Ta2O5 (Pa)
ni1 = ni1_SiO2;             % coating Poisson ratio SiO2
ni2 = ni2_Ta205;            % coating Poisson ratio Ta2O5


% zeros bessel function
bessj1 = inline('besselj(1,x)');
for n = 1:10
z1(:,n) = fzero(bessj1,n*pi);
end
J0 = besselj(0,z1);
J2 = besselj(2,z1);
r=sym('r');

k1=z1/a;
Q1=exp(-2*k1*h);
J0r = besselj(0,k1*r);
J2r = besselj(2,k1*r);

%Lame coefficient
lambda=Y1*ni1/((1+ni1)*(1-2*ni1));
mu=Y1/(2*(1+ni1));
lambdac=Y2*ni2/((1+ni2)*(1-2*ni2));
muc=Y2/(2*(1+ni2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   ITM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sigma_itm = 0.053/0.34; 
w_itm = 2*a*sigma_itm;
%w=0.053; % beam radius (m)
d1=182e-9*9; % coating thickness SiO2
d2=131e-9*8; % coating thickness Ta2O5

for ii=1:length(k1)
%constants
p(ii)=exp(-k1(ii)^2*w_itm^2/8)/(pi*a^2*J0(ii)^2);
alpha1=p*(lambda+2*mu)/(k1(ii)*mu*(lambda+mu))*(1-Q1(ii)+2*k1(ii)*h.*Q1(ii))/((1-Q1(ii)).^2-4*k1(ii)^2.*h.^2.*Q1(ii));
beta1=p*(lambda+2*mu)*Q1(ii)/(k1(ii)*mu*(lambda+mu)).*(1-Q1(ii)+2*k1(ii)*h)/((1-Q1(ii)).^2-4*k1(ii)^2.*h.^2.*Q1(ii));
gamma1=-p/(2*k1(ii)*mu*(lambda+mu))*((2*k1(ii)^2*h.^2*(lambda+mu)+2*mu*k1(ii)*h).*Q1(ii)+mu*(1-Q1(ii)))/((1-Q1(ii)).^2-4*k1(ii)^2*h.^2.*Q1(ii));
delta1=-p*Q1(ii)/(2*k1(ii)*mu*(lambda+mu)).*(2*k1(ii)^2*h.^2*(lambda+mu)-2*mu*k1(ii)*h-mu*(1-Q1(ii)))/((1-Q1(ii)).^2-4*k1(ii)^2*h.^2.*Q1(ii));
D(ii)=J0(ii)*p(ii)/z1(ii)^2;
p0=1/(pi*a^2);

%First term of Strain tensors SiO2
A1(ii)=k1(ii)*(gamma1+delta1)/2*(J0r(ii)-J2r(ii));
B1(ii)=k1(ii)*(gamma1+delta1)/2*(J0r(ii)+J2r(ii));
C1(ii)=-1/(lambda+2*mu)*k1(ii)*J0r(ii)*(mu*(alpha1-beta1)+(lambda+2*mu)*(gamma1+delta1));

%First term of Strain tensors Ta
A2(ii)=k1(ii)*(gamma1+delta1)/2*(J0r(ii)-J2r(ii));
B2(ii)=k1(ii)*(gamma1+delta1)/2*(J0r(ii)+J2r(ii));
C2(ii)=-1/(lambdac+2*muc)*k1(ii)*J0r(ii)*(mu*(alpha1-beta1)+(lambdac+2*mu)*(gamma1+delta1));

end
c0=6*a^2/h^2*sum(D);
% Strain tensors SiO2
Err1=sum(A1)+((lambda+2*mu)*c0+lambda*p0)/(2*mu*(3*lambda+2*mu));
Ephiphi1=sum(B1)+((lambda+2*mu)*c0+lambda*p0)/(2*mu*(3*lambda+2*mu));
Ezz1=sum(C1)-(lambda*(lambda+2*mu)*c0+(lambda*lambda+3*lambda*mu+2*mu^2)*p0)/(mu*(3*lambda+2*mu)*(lambda+2*mu));

%Stress tensor SiO2
Trr1=(lambda+2*mu)*Err1+lambda*(Ephiphi1+Ezz1);
Tphiphi1=(lambda+2*mu)*Ephiphi1+lambda*(Ezz1+Err1);
Tzz1=(lambda+2*mu)*Ezz1+lambda*(Err1+Ephiphi1);

%Strain tensors Ta
Err2=sum(A2)+((lambda+2*mu)*c0+lambda*p0)/(2*mu*(3*lambda+2*mu));
Ephiphi2=sum(B2)+((lambda+2*mu)*c0+lambda*p0)/(2*mu*(3*lambda+2*mu));
Ezz2=sum(C2)-(lambdac*(lambda+2*mu)*c0+(lambda*lambdac+3*lambda*mu+2*mu^2)*p0)/(mu*(3*lambda+2*mu)*(lambdac+2*muc));

%Stress tensor Ta
Trr2=(lambdac+2*muc)*Err2+lambdac*(Ephiphi2+Ezz2);
Tphiphi2=(lambdac+2*muc)*Ephiphi2+lambdac*(Ezz2+Err2);
Tzz2=(lambdac+2*muc)*Ezz2+lambdac*(Err2+Ephiphi2);


for kk=1:length(Err1)
u1{kk}=matlabFunction(pi*d1*(Err1(kk)*Trr1(kk)+Ephiphi1(kk)*Tphiphi1(kk)+Ezz1(kk)*Tzz1(kk))*r);
u2{kk}=matlabFunction(pi*d2*(Err2(kk)*Trr2(kk)+Ephiphi2(kk)*Tphiphi2(kk)+Ezz2(kk)*Tzz2(kk))*r);
end

for xx=1:length(u1)
    for n=1:length(a)
        %Energy 
         U1(xx,n)=integral(u1{xx},0,a(n));
         U2(xx,n)=integral(u2{xx},0,a(n));
    end
end
U1=diag(U1);
U2=diag(U2);

%Power spectrum
S_itm=(8*Kb*T./omega*phi1*U1)+(8*Kb*T./omega*phi2*U2);

%Strain
Strain_itm=S_itm/3994.5^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  ETM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma = mirror_sigma;
w = 2*a*sigma;
%w=0.062; % beam radius (m)
d1=182e-9*21; % coating thickness SiO2
d2=131e-9*20; % coating thickness Ta2O5

for ii=1:length(k1)
%constants
p=exp(-k1(ii)^2*w^2/8)/(pi*a^2*J0(ii)^2);
alpha1=p*(lambda+2*mu)/(k1(ii)*mu*(lambda+mu))*(1-Q1(ii)+2*k1(ii)*h.*Q1(ii))/((1-Q1(ii)).^2-4*k1(ii)^2.*h.^2.*Q1(ii));
beta1=p*(lambda+2*mu)*Q1(ii)/(k1(ii)*mu*(lambda+mu)).*(1-Q1(ii)+2*k1(ii)*h)/((1-Q1(ii)).^2-4*k1(ii)^2.*h.^2.*Q1(ii));
gamma1=-p/(2*k1(ii)*mu*(lambda+mu))*((2*k1(ii)^2*h.^2*(lambda+mu)+2*mu*k1(ii)*h).*Q1(ii)+mu*(1-Q1(ii)))/((1-Q1(ii)).^2-4*k1(ii)^2*h.^2.*Q1(ii));
delta1=-p*Q1(ii)/(2*k1(ii)*mu*(lambda+mu)).*(2*k1(ii)^2*h.^2*(lambda+mu)-2*mu*k1(ii)*h-mu*(1-Q1(ii)))/((1-Q1(ii)).^2-4*k1(ii)^2*h.^2.*Q1(ii));
c0=6*a^2./h.^2*J0(ii)*p/z1(ii)^2;
p0=1/(pi*a^2);

%First term of Strain tensors SiO2
A1(ii)=k1(ii)*(gamma1+delta1)/2*(J0r(ii)-J2r(ii));
B1(ii)=k1(ii)*(gamma1+delta1)/2*(J0r(ii)+J2r(ii));
C1(ii)=-1/(lambda+2*mu)*k1(ii)*J0r(ii)*(mu*(alpha1-beta1)+(lambda+2*mu)*(gamma1+delta1));

%First term of Strain tensors Ta
A2(ii)=k1(ii)*(gamma1+delta1)/2*(J0r(ii)-J2r(ii));
B2(ii)=k1(ii)*(gamma1+delta1)/2*(J0r(ii)+J2r(ii));
C2(ii)=-1/(lambdac+2*muc)*k1(ii)*J0r(ii)*(mu*(alpha1-beta1)+(lambdac+2*mu)*(gamma1+delta1));

end

% Strain tensors SiO2
Err1=sum(A1)+((lambda+2*mu)*c0+lambda*p0)/(2*mu*(3*lambda+2*mu));
Ephiphi1=sum(B1)+((lambda+2*mu)*c0+lambda*p0)/(2*mu*(3*lambda+2*mu));
Ezz1=sum(C1)-(lambda*(lambda+2*mu)*c0+(lambda*lambda+3*lambda*mu+2*mu^2)*p0)/(mu*(3*lambda+2*mu)*(lambda+2*mu));

%Stress tensor SiO2
Trr1=(lambda+2*mu)*Err1+lambda*(Ephiphi1+Ezz1);
Tphiphi1=(lambda+2*mu)*Ephiphi1+lambda*(Ezz1+Err1);
Tzz1=(lambda+2*mu)*Ezz1+lambda*(Err1+Ephiphi1);

%Strain tensors Ta
Err2=sum(A2)+((lambda+2*mu)*c0+lambda*p0)/(2*mu*(3*lambda+2*mu));
Ephiphi2=sum(B2)+((lambda+2*mu)*c0+lambda*p0)/(2*mu*(3*lambda+2*mu));
Ezz2=sum(C2)-(lambdac*(lambda+2*mu)*c0+(lambda*lambdac+3*lambda*mu+2*mu^2)*p0)/(mu*(3*lambda+2*mu)*(lambdac+2*muc));

%Stress tensor Ta
Trr2=(lambdac+2*muc)*Err2+lambdac*(Ephiphi2+Ezz2);
Tphiphi2=(lambdac+2*muc)*Ephiphi2+lambdac*(Ezz2+Err2);
Tzz2=(lambdac+2*muc)*Ezz2+lambdac*(Err2+Ephiphi2);



for kk=1:length(Err1)
u1{kk}=matlabFunction(pi*d1*(Err1(kk)*Trr1(kk)+Ephiphi1(kk)*Tphiphi1(kk)+Ezz1(kk)*Tzz1(kk))*r);
u2{kk}=matlabFunction(pi*d2*(Err2(kk)*Trr2(kk)+Ephiphi2(kk)*Tphiphi2(kk)+Ezz2(kk)*Tzz2(kk))*r);
end

for xx=1:length(u1)
    for n=1:length(a)
        %Energy 
         U1(xx,n)=integral(u1{xx},0,a(n));
         U2(xx,n)=integral(u2{xx},0,a(n));
    end
end
U1=diag(U1);
U2=diag(U2);

%Power spectrum
S_etm=(8*Kb*T./omega*phi1*U1)+(8*Kb*T./omega*phi2*U2);


%Strain
Strain_etm=S_etm/3994.5^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Total
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=sqrt(2*Strain_itm+2*Strain_etm);


total_coating_therm_final = S;
end