%%bLCGT

clear all;
clf 

dataSQL = load('hSQL.dat');
f = dataSQL(:,1);
univec = ones(size(f));

hb = 1.05457266912510183 * 10^-34; %(* reduced Planck constant (J*s) *)
Omega0 = 1.8 * 10^15;
c = 2.99792458000000117 * 10^8; %(* speed of light (m/s) *)
CapitalOmega = 2 * pi * f;
kB = 1.38 * 10^(-23); %(* Boltzmann's constant (J/K) *)
g = 9.8; %(* acceleration of the gravity (m/s^2) *)
L = 3000;

lambda = 1064 * 10^(-9); %(* wavelength (m) *)

l = 3*10^3;  %(* length of cavity (m) *)
w1 = 3.43*10^(-2);  %(* ITM beam radius at mirrors (m) *)
w2 = 4.53*10^(-2);  %(* ETM beam radius at mirrors (m) *)

rho = 4.0*10^3; %(* density (kg/m^3) *)
radius = 12.5*10^(-2); %(* radius (m); used to be 12.5cm *)
height = 15*10^(-2); %(* height (m) *)
mass = rho*pi*radius^2*height; %(* mass (kg) *)
Tm = 20; %(* temperature of mirror (K) *)
E0 = 4.0*10^(11); %(* Young's modulus (Pa) *)
sigma = 0.29; %(* Poisson ratio *)
Qmirror = 10^(8); %(* loss angle of material *)
alpha20 = 5.6*10^(-9); %(* thermal expansion coefficient at 20K (/K) *)
Cth20 = 0.69; %(* specific heat at 20 K (J/kg*K) *)
kappa20 = 1.57*10^4; %(* thermal conductivity at 20 K (W/m*K) *)

Ncoa1 = 9;
Ncoa2 = 18;

Ys = 7.2*10^10; %(* Young's modulus for silica (Pa=N/m^2) *)
Yt = 1.4*10^11; %(* Young's modulus for Ta2O5 (Pa=N/m^2) *)
sigmas = 0.17; %(* Poisson ratio for silica *)
sigmat = 0.23; %(* Poisson ratio for Ta2O5 or sapphire *)
n1 = 1.45;  %(* Refraction index of silica *)
n2 = 2.07;  %(* Refraction index of tantala *)
dcoa1s = (Ncoa1 + 1)/4*10^(-6)/n1; %(* thickness of ITM silica coating (m) *)
dcoa1t = Ncoa1/4*10^(-6)/n2; %(* thickness of ITM tantala coating (m) *)
dcoa2s = (Ncoa2 + 1)/4*10^(-6)/n1; %(* thickness of ETM silica coating (m) *)
dcoa2t = Ncoa2/4*10^(-6)/n2; %(* thickness of ETM tantala coating (m) *)
Phic1 = 3.0*10^(-4); %(* loss angle of silica coating *)
Phic2 = 5.0*10^(-4); %(* loss angle of tantala coating *)

Cs = 1.64*10^6;
Kappas = 1.38/Cs;
Alphaexps = 5.1*10^-7;
Ct = 2.1*10^6;
Kappat = 33/Ct;
Alphaexpt = 3.6* 10^-6;

%%1.6 Suspension

Omega = CapitalOmega;

m1 = 60;
m2 = mass; 
m3 = 30;

mjoint = 0.1;  %(* mini GAS joint *)
fmg = 3;  %(* mini GAS freq *)

Mu21 = m2/m1;
Mu31 = m3/m1;

r1 = 0.0006/2;  %(* This doesn't matter much *)
r2 = 0.0016/2;  %(* 1W heat extraction *)
r3 = 0.0004/2; 

I1 = pi/4*r1^4;
I2 = pi/4*r2^4;
I3 = pi/4*r3^4;

Y1 = 161*10^9;  %(* Tungsten *)
Y2 = 400*10^9;  %(* Sapphire *)
Y3 = 130*10^9;  %(* BeCu *)

N1 = 4;
N2 = 4;
N3 = 4;

lsus1 = 0.5; 
lsus2 = 0.3; 
lsus3 = 0.3;

b1 = 2;
b2 = 2;
b3 = 2;

D1 = (b1*N1*sqrt((m1*g/N1)*Y1*I1))/(2*m1*g*lsus1);  %(* Saulson(1990); b is added. *)
D2 = (b2*N2*sqrt((m2*g/N2)*Y2*I2))/(2*m2*g*lsus2);
D3 = (b3*N3*sqrt((m3*g/N3)*Y3*I3))/(2*m3*g*lsus3);

Omega1 = sqrt(g/lsus1)*(1 + D1);
Omega2 = sqrt(g/lsus2)*(1 + D2);
Omega3 = sqrt(g/lsus3)*(1 + D3);

T1 = 10;
T2 = 16;  %(* temp *)
T3 = 16;  %(* temp *)

Phip = 2*10^-7;  %(* Sapphire *)
Phip2 = 10^-4;  %(* Tungsten *)
Phip3 = 5*10^-6;  %(* BeCu *)
Phicl = 10^-3;
Phicl2 = 10^-3;

Q1 = (Omega1./Omega* Phip2*D1 + Phicl*univec).^-1;
Q2 = (Omega2./Omega* Phip*D2 + 0*univec).^-1;
Q3 = (Omega3./Omega* Phip3*D3 + 0*10^-7*univec).^-1;


VHC = 1/200; 

lsus = lsus2; 
dwire = 2*r2;
Ewire = Y2;

n = N2;

Icwire = I2; %(* moment of cross section (m^4) *)

ksus = sqrt((-mass*g/n +  sqrt((mass*g/n)^2 + 4*Ewire*(1 + 1i*Phip)*Icwire*rho*pi*(dwire/2)^2.*Omega.^2))/(2*Ewire*(1 + 1i*Phip)*Icwire));
kele = sqrt((mass*g/n + sqrt((mass*g/n)^2 + 4*Ewire*(1 + 1i*Phip)*Icwire*rho*pi*(dwire/2)^2.*Omega.^2))/(2*Ewire*(1 + 1i*Phip)*Icwire));
Ksusp = n*Ewire*(1 + 1i*Phip)*Icwire.*(kele.*ksus.*(kele.^2 + ksus.^2).*(sinh(kele*lsus).*kele.*cos(ksus*lsus) + cosh(kele*lsus).*ksus.*sin(ksus*lsus)))./(2*kele.*ksus.*(1 - cos(ksus*lsus).*cosh(kele*lsus)) + (kele.^2 - ksus.^2).*sin(ksus*lsus).*sinh(kele*lsus));

%%Impedance matrix

D11 = -Omega.^2 + Omega1^2 + Mu21 * Ksusp/m2 + Mu31 * Omega2^2 + 1i*Omega.*(Omega1./Q1 + Mu31*Omega3./Q3);  
D12 = -Mu21*Ksusp/m2;
D13 = -Mu31*Omega3^2 - 1i*Omega.*(Mu31*Omega3./Q3);
D21 = -Ksusp/m2;
D22 = -Omega.^2 + Ksusp/m2;
D23 = 0;
D31 = -Omega3^2 - 1i*Omega.*(Omega3./Q3);
D32 = 0;
D33 = -Omega.^2 + Omega3^2 + 1i*Omega.*(Omega3./Q3);

Z=zeros(3,3,length(f));

Z(1,1,:)=m1*D11./(1i.*Omega);
Z(1,2,:)=m1*D12./(1i.*Omega);
Z(1,3,:)=m1*D13./(1i.*Omega);
Z(2,1,:)=m2*D21./(1i.*Omega);
Z(2,2,:)=m2*D22./(1i.*Omega);
Z(2,3,:)=m2*D23./(1i.*Omega);
Z(3,1,:)=m3*D31./(1i.*Omega);
Z(3,2,:)=m3*D32./(1i.*Omega);
Z(3,3,:)=m3*D33./(1i.*Omega);

%%Thermal noise spectrum (PPP Eq.(69))

Tmat=zeros(3,3,length(f));

Tmat(1,1,:)=T1;
Tmat(2,2,:)=T2;
Tmat(3,3,:)=T3;

Fsq = 4*kB*real(Tmat.*Z);

Dmat=zeros(3,3,length(f));
Dmat(1,1,:)=D11;
Dmat(1,2,:)=D12;
Dmat(1,3,:)=D13;
Dmat(2,1,:)=D21;
Dmat(2,2,:)=D22;
Dmat(2,3,:)=D23;
Dmat(3,1,:)=D31;
Dmat(3,2,:)=D32;
Dmat(3,3,:)=D33;

invD=zeros(3,3,length(f));
for ff=1:length(f)
    fDmat=Dmat(:,:,ff);
    finvD=inv(fDmat);
    invD(:,:,ff)=finvD(:,:);
end

Sth=zeros(1,length(f));
Sth(:)=Fsq(1,1,:).*abs(invD(2,1,:)).^2/m1^2 + Fsq(2,2,:).*abs(invD(2,2,:)).^2/m2^2 + Fsq(3,3,:).*abs(invD(2,3,:)).^2/m3^2 + Fsq(1,2,:)/m1/m2*2.*real(invD(2,1,:).*conj(invD(2,2,:))) + Fsq(2,3,:)/m2/m3*2.*real(invD(2,2,:).*conj(invD(2,3,:))) + Fsq(1,3,:)/m1/m3*2.*real(invD(2,1,:).*conj(invD(2,3,:)));

%%Vertical thermal noise

Mut = m1/(m1 + m2 + m3);

Omegav1 = sqrt((2*pi*fmg)^2/Mut);
Omegav2 = sqrt((pi*r2^2*Y2/lsus2)/(m2/4));
Omegav3 = sqrt((pi*r3^2*Y3/lsus3)/(m3/4));

Qv1 = (sqrt(m1/mjoint)./(2*Omegav2*Omegav1^2*Omega).*(Omegav1./Omega*Phip2 + univec*Phicl2) + univec*Phicl./(Omega*sqrt(m1/mjoint))).^-1;  %(* see minigas2.nb *)
Qv2 = (Omegav2./Omega*Phip).^-1;  %(* no visvous damping *)
Qv3 = (Omegav3./Omega*Phip3 + 0*10^-7).^-1; 

D11v = -Omega.^2 + Omegav1^2 + Mu21*Omegav2^2 + Mu31*Omegav3^2 + 1i*Omega.*(Omegav1./Qv1 + Mu21*Omegav2./Qv2 + Mu31*Omegav3./Qv3);
D12v = -Mu21*Omegav2^2 - 1i*Omega.*(Mu21*Omegav2./Qv2);
D13v = -Mu31*Omegav3^2 - 1i*Omega.*(Mu31*Omegav3./Qv3);
D21v = -Omegav2^2 - 1i*Omega.*(Omegav2./Qv2);
D22v = -Omega.^2 + Omegav2^2 + 1i*Omega.*(Omegav2./Qv2);
D23v = 0;
D31v = -Omegav3^2 - 1i*Omega.*(Omegav3./Qv3);
D32v = 0;
D33v = -Omega.^2 + Omegav3^2 + 1i*Omega.*(Omegav3./Qv3);

Zv=zeros(3,3,length(f));

Zv(1,1,:)=m1*D11v./(1i*Omega);
Zv(1,2,:)=m1*D12v./(1i*Omega);
Zv(1,3,:)=m1*D13v./(1i*Omega);
Zv(2,1,:)=m2*D21v./(1i*Omega);
Zv(2,2,:)=m2*D22v./(1i*Omega);
Zv(2,3,:)=m2*D23v./(1i*Omega);
Zv(3,1,:)=m3*D31v./(1i*Omega);
Zv(3,2,:)=m3*D32v./(1i*Omega);
Zv(3,3,:)=m3*D33v./(1i*Omega);

%%Thermal noise spectrum 

Fsqv = 4*kB*real(Tmat.*Zv);

Dvmat=zeros(3,3,length(f));
Dvmat(1,1,:)=D11v;
Dvmat(1,2,:)=D12v;
Dvmat(1,3,:)=D13v;
Dvmat(2,1,:)=D21v;
Dvmat(2,2,:)=D22v;
Dvmat(2,3,:)=D23v;
Dvmat(3,1,:)=D31v;
Dvmat(3,2,:)=D32v;
Dvmat(3,3,:)=D33v;

invDv=zeros(3,3,length(f));
for ff=1:length(f)
    fDvmat=Dvmat(:,:,ff);
    finvDv=inv(fDvmat);
    invDv(:,:,ff)=finvDv(:,:);
end

Sthv=zeros(1,length(f));
Sthv(:)=Fsqv(1,1,:).*abs(invDv(2,1,:)).^2/m1^2 + Fsqv(2,2,:).*abs(invDv(2,2,:)).^2/m2^2 + Fsqv(3,3,:).*abs(invDv(2,3,:)).^2/m3^2 + Fsqv(1,2,:)/m1/m2*2.*real(invDv(2,1,:).*conj(invDv(2,2,:))) + Fsqv(2,3,:)/m2/m3*2.*real(invDv(2,2,:).*conj(invDv(2,3,:))) + Fsqv(1,3,:)/m1/m3*2.*real(invDv(2,1,:).*conj(invDv(2,3,:)));

%%Seismic noise

hseis=zeros(1,length(f));
hseis(:) = 2*sqrt((2*10^-17*f.^-6.5).^2 + (2*10^-16* f.^-8).^2);

%%Suspension Thermal noise

hsusp=zeros(1,length(f));
hsusp(:) = sqrt(abs(Sth) + abs(Sthv)*VHC^2)*2/L;

%%Mirror Thermal noise

hmirrorstr=zeros(1,length(f));
hmirrorstr(:) = sqrt(2)/L*sqrt(4*kB*Tm*(1-sigma^2)./(sqrt(pi)*E0*w1*Qmirror*Omega)+4*kB*Tm*(1-sigma^2)./(sqrt(pi)*E0*w2*Qmirror*Omega));

%%Mirror Thermoelastic noise

hmirrorthermo=zeros(1,length(f));

Omegac1 = pi*Cth20*rho*w1^2*f/kappa20;
Omegac2 = pi*Cth20*rho*w2^2*f/kappa20;
Jhyper1 = 1./(Omegac1.^2).*(1-real(exp((Omegac1*1i)/2).*((1-Omegac1*1i).*(1-erfz((sqrt(Omegac1)*(1+1i))/2)))))-sqrt(1./(pi*Omegac1.^3));
Jhyper2 = 1./(Omegac2.^2).*(1-real(exp((Omegac2*1i)/2).*((1-Omegac2*1i).*(1-erfz((sqrt(Omegac2)*(1+1i))/2)))))-sqrt(1./(pi*Omegac2.^3));

hmirrorthermo(:)=sqrt(2)/L*sqrt((4*kB*Tm^2*alpha20^2*(1+sigma)^2*w1)/(sqrt(pi)*kappa20)*Jhyper1+(4*kB*Tm^2*alpha20^2*(1+sigma)^2*w2)/(sqrt(pi)*kappa20)*Jhyper2);

%%Coating thermal noise
hmirrorcoa=zeros(1,length(f));
hmirrorcoa(:) = sqrt(2)/L*sqrt((4*kB*Tm*dcoa1s*Phic1)./(pi*w1^2*Omega)*(Ys^2*(1+sigma)^2*(1-2*sigma)^2+E0^2*(1+sigmas)^2*(1-2*sigmas))/(E0^2*Ys*(1-sigmas^2))+(4*kB*Tm*dcoa1t*Phic2)./(pi*w1^2*Omega)*(Yt^2*(1+sigma)^2*(1-2*sigma)^2+E0^2*(1+sigmat)^2*(1-2*sigmat))/(E0^2*Yt*(1-sigmat^2))+(4*kB*Tm*dcoa2s*Phic1)./(pi*w2^2*Omega)*(Ys^2*(1+sigma)^2*(1-2*sigma)^2+E0^2*(1+sigmas)^2*(1-2*sigmas))/(E0^2*Ys*(1-sigmas^2))+(4*kB*Tm*dcoa2t*Phic2)./(pi*w2^2*Omega)*(Yt^2*(1+sigma)^2*(1-2*sigma)^2+E0^2*(1+sigmat)^2*(1-2*sigmat))/(E0^2*Yt*(1-sigmat^2)));

%%Thermal noise total
hmirror=sqrt(hmirrorstr.^2+hmirrorthermo.^2+hmirrorcoa.^2);

%%Quantum noise

T = 0.004;
rho = 0.92; % SRM amplitude reflectivity
%phi = pi/2;
%zeta = pi/2;
phi = 1.50946;
zeta = -0.792759;
m = mass; % test mass (kg) 
P0 = 75; % power into interderometer(W)
Gpower = 11; % power recycling gain
I0 = P0*Gpower; % incident power to the beamsplitter

hSQL = sqrt((8*hb)./(m.*Omega.^2*L^2));
gamma = T*c/(4*L); % cavity pole 
tau = sqrt(1-rho^2); % SRM transmittance 
beta = atan(Omega/gamma); % phase shift in the arm cavity 

ISQL = (m*L^2*gamma^4)/(4*Omega0); % the power with which the quantum noise reaches the SQL at the cavity pole frequency
kappa = (2*(I0/ISQL)*gamma^4)./(Omega.^2.*(gamma^2+Omega.^2)); % coefficient that shows how much the mirror moves

roundtriploss = 100*10^(-6); 
epsilon = 2/T*roundtriploss; % loss parameter of arms defined in KLMTV paper
lambdaSR = 0.02; % loss in power at the signal recycling mirror 
lambdaPD = 0.1;  % 1 - PD's quantum efficiency in power 
Sadd = 0*(pi^2/8 - 1);

M = 1+rho^2.*exp(4*1i.*beta)-2*rho*(cos(2*phi)+kappa/2*sin(2*phi)).*exp(2*1i*beta)+lambdaSR*rho*(-rho*exp(2*1i*beta)+cos(2*phi)+kappa/2*sin(2*phi)).*exp(2*1i*beta)+epsilon*rho*(2*cos(beta).^2.*(-rho*exp(2*1i*beta)+cos(2*phi))+kappa/2.*(3+exp(2*1i*beta))*sin(2*phi)).*exp(2*1i*beta);
C11 = sqrt(1-lambdaPD)*((1+rho^2)*(cos(2*phi)+kappa/2*sin(2*phi))-2*rho*cos(2*beta)-1/4*epsilon*(-2*(1+exp(2*1i*beta)).^2*rho+4*(1+rho^2)*cos(beta).^2*cos(2*phi)+(3+exp(2*1i*beta)).*kappa*(1+rho^2)*sin(2*phi))+lambdaSR*(exp(2*1i*beta)*rho-1/2*(1+rho^2)*(cos(2*phi)+kappa/2*sin(2*phi))));
C12 = sqrt(1-lambdaPD)*tau^2*(-sin(2*phi)-kappa*sin(phi)^2+1/2*epsilon*sin(phi)*((3+exp(2*1i*beta)).*kappa*sin(phi)+4*cos(beta).^2*cos(phi))+1/2*lambdaSR*(sin(2*phi)+kappa*sin(phi)^2));
C21 = sqrt(1-lambdaPD)*tau^2*(sin(2*phi)-kappa*cos(phi)^2+1/2*epsilon*cos(phi)*((3+exp(2*1i*beta)).*kappa*cos(phi)+4*cos(beta).^2*sin(phi))+1/2*lambdaSR*(-sin(2*phi)+kappa*cos(phi)^2));
C22 = C11;
D1sig = sqrt(1-lambdaPD)*(-(1+rho*exp(2*1i*beta))*sin(phi)+1/4*epsilon*(3+rho+2*rho*exp(4*1i*beta)+(1+5*rho)*exp(2*1i*beta))*sin(phi)+1/2*lambdaSR*exp(2*1i*beta)*rho*sin(phi));
D2sig = sqrt(1-lambdaPD)*(-(-1+rho*exp(2*1i*beta))*cos(phi)+1/4*epsilon*(-3+rho+2*rho*exp(4*1i*beta)+(-1+5*rho)*exp(2*1i*beta))*cos(phi)+1/2*lambdaSR*exp(2*1i*beta)*rho*cos(phi));
P11 = 1/2*sqrt(1-lambdaPD)*sqrt(lambdaSR)*tau*(-2*rho*exp(2*1i*beta)+2*cos(2*phi)+kappa*sin(2*phi));
P12 = -sqrt(1-lambdaPD)*sqrt(lambdaSR)*tau*sin(phi)*(2*cos(phi)+kappa*sin(phi));
P21 = sqrt(1-lambdaPD)*sqrt(lambdaSR)*tau*cos(phi)*(2*sin(phi)-kappa*cos(phi));
P22 = P11;
Q11 = sqrt(lambdaPD)*(exp(-2*1i*beta)+rho^2*exp(2*1i*beta)-rho*(2*cos(2*phi)+kappa*sin(2*phi))+1/2*epsilon*rho*(exp(-2*1i*beta)*cos(2*phi)+exp(2*1i*beta).*(-2*rho-2*rho*cos(2*beta)+cos(2*phi)+kappa*sin(2*phi))+2*cos(2*phi)+3*kappa*sin(2*phi))-1/2*lambdaSR*rho*(2*rho*exp(2*1i*beta)-2*cos(2*phi)-kappa*sin(2*phi)));
Q12 = 0;
Q21 = 0;
Q22 = Q11;
N11 = sqrt(1-lambdaPD)*sqrt(epsilon/2)*tau*(kappa.*(1+rho*exp(2*1i*beta))*sin(phi)+2*cos(beta).*(exp(-1i*beta)*cos(phi)-rho*exp(1i*beta).*(cos(phi)+kappa*sin(phi))));
N12 = -sqrt(1-lambdaPD)*sqrt(2*epsilon)*tau*(exp(-1i*beta)+rho*exp(1i*beta)).*cos(beta)*sin(phi);
N22 = -sqrt(1-lambdaPD)*sqrt(2*epsilon)*tau*(-exp(-1i*beta)+rho*exp(1i*beta)).*cos(beta)*cos(phi);
N21 = sqrt(1-lambdaPD)*sqrt(epsilon/2)*tau*(-kappa.*(1+rho)*cos(phi)+2*cos(beta).*(exp(-1i*beta)+rho*exp(1i*beta)).*cos(beta)*sin(phi));
      
Sh=zeros(1,length(f));
Sh(:) = hSQL.^2./(2*kappa).*(abs(C11*sin(zeta)+C21*cos(zeta)).^2+abs(C12*sin(zeta)+C22*cos(zeta)).^2+abs(P12*sin(zeta)+P22*cos(zeta)).^2+abs(P11*sin(zeta)+P21*cos(zeta)).^2+abs(Q12*sin(zeta)+Q22*cos(zeta)).^2+abs(Q11*sin(zeta)+Q21*cos(zeta)).^2+abs(N11*sin(zeta)+N21*cos(zeta)).^2+abs(N12*sin(zeta)+N22*cos(zeta)).^2+(abs(M)*sqrt(Sadd)).^2)./(tau^2*abs(D1sig*sin(zeta)+D2sig*cos(zeta)).^2);
%Sh(:) = hSQL.^2./(2*kappa).*(abs(C11*sin(zeta)+C21*cos(zeta)).^2+abs(C12*sin(zeta)+C22*cos(zeta)).^2)./(tau^2*abs(D1sig*sin(zeta)+D2sig*cos(zeta)).^2);

loglog(f,sqrt(Sh),'-k');
hold on
loglog(f,hmirror,'-r');
loglog(f,hsusp,'-b');
loglog(f,hseis,'-c');
loglog(f,hSQL,'-g');
xlabel('Frequency [Hz]','FontSize',12)
ylabel('Sensitivity [1/rtHz]','FontSize',12)
title('bLCGT','FontSize',14)
legend('QN','mirror','suspension','seismic','SQL',1);
axis([5,5e3,1e-25,1e-20]);
saveas(gcf,'bLCGT.png');
hold off

dlmwrite('QN.dat', [transpose(f);sqrt(Sh)]', '\t')
dlmwrite('hmirror.dat', [transpose(f);hmirror]', '\t')
dlmwrite('hsusp.dat', [transpose(f);hsusp]', '\t')
dlmwrite('hseis.dat', [transpose(f);hseis]', '\t')