Back to index.
%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