% Sensitivity and inspiral range of KAGRA (and KAGRA+)
% Written by K. Komori, based on Mathematica code written by K. Somiya
% Reviewed by K. Yamamoto, T. Yamamoto, Y. Enomoto and approved

clear
close all

%% DEFINE PARAMETERS %%%%%%%%%%
%RSE configuration

RSE='DRSE';     % 'BRSE', 'BRSEhomo', or 'DRSE' (BRSEhomo for optimized homodyne readout phase)

%%%%%physical constants and other basic parameters

%from Rev. Mod. Phys. 88, 035009 (2016)
c=2.99792458e8;        %speed of light (m/s)
hbar=1.054571800e-34;  %Dirac constant (J*s)
kb=1.38064852e-23;     %Boltzmann constant (J/K)
e=1.6021766208e-19;    %electric charge (C)
G=6.67384e-11;         %gravitational constant (m^3/kg/s^2)
g=9.80665;             %gravitational acceleration (m/s^2)

lambda=1.064e-6;        %laser wave length (m)
omegal=2*pi*c/lambda;   %laser frequency (Hz)

L=3e3;     %arm length (m)

fres=1000;               %frequency resolution
f=logspace(0,4,fres)';
omega=2*pi.*f;
% Frequency resolution could be insufficient for violin modes. Change the resolution if necessary.


%%%%%interferometer parameters

%Detune phase and homodyne phase are optimized based on BNS range
%Re-optimization is necessary if other parameters are changed
if strcmp(RSE,'DRSE')     %detuned RSE
    phidet=pi/180*86.5;   %detune phase
    xi=pi/180*135.1;      %homodyne phase
elseif strcmp(RSE,'BRSEhomo')    %broadband RSE with homodyne
    phidet=pi/180*90;
    xi=pi/180*119.1;
elseif strcmp(RSE,'BRSE')    %broadband RSE without homodyne
    phidet=pi/180*90;
    xi=pi/180*90;
end

trans=0.004;          %power transmittivity of ITM
rSRM=0.92;            %amplitude reflectivity of SRM
roundtriploss=100e-6; %round trip loss of arm
lambdaSR=0.002;       %power loss at SRM
lambdaPD=0.1;         %power loss at PD
Sadd=0*(pi^2/8-1);

%%%%%mirror parameters

rhom=4.0e3;      %density of mirror (kg/m^3) %sapphire
radius=11e-2;    %radius of mirror (m)
height=15e-2;    %thickness of mirror (m)
mass=22.8;       %mirror mass (kg)
w1=3.5e-2;       %beam radius at ITM (m)
w2=3.5e-2;       %beam radius at ETM (m)
Tm=22;           %temperature at TM (K)
E0=400e9;        %Youngs modulus of mirror (Pa), see Phys. Lett. A 264, 1 (1999)
sigma=0.29;      %poisson ratio of mirror, see Phys. Lett. A 264, 1 (1999)
Qm=1e8;          %Q-value of mirror


%From Uchiyama-san. Origninally from Touloukian and Ho, Thermophysical Properties of Matter: The TRPC Data Series, Plenum Press (1970).
%see PRL 108, 141101 (2012).
kappa20=(15700/15880)*((5*Tm^2.75)^(-4/5)+(10^10.25*Tm^(-3.8))^(-4/5))^(-5/4); %thermal conductivity
Cth20=(0.69/0.80)*Tm^3.14/3.8/rhom;  %specific heat per unit mass
alpha20=(5.6/5.496)*Tm^2.99/10^12.15;  %thermal linear expansion


%%%%%coating parameters

%number of coating layers
Ncoa1=11;       %ITM
Ncoa2=20;       %ETM

%Youngs modulus (Pa)
Ys=73.2e9;       %silica, see PRD 93, 012007 (2016)
Yt=140e9;        %tantala, see PRD 70, 082003 (2004)

%poisson ratio
sigmas=0.164;    %silica, see PRD 93, 012007 (2016)
sigmat=0.23;     %tantala, see PRD 70, 082003 (2004)

%refractive index
n1=1.442;        %silica, see PRD 93, 012007 (2016)
n2=2.035;        %tantala, see CQG 27, 084030 (2010)

%thickness of layers (m)
dcoa1s=(Ncoa1+1)*lambda/(4*n1);  %silica, ITM
dcoa1t=Ncoa1*lambda/(4*n2);      %tantala, ITM
dcoa2s=(Ncoa2+1)*lambda/(4*n1);  %silica, ETM
dcoa2t=Ncoa2*lambda/(4*n2);      %tantala, ITM

%loss angle
phic1=3e-4;       %silica coating
phic2=5e-4;       %tantala coating


%%%%%suspension parameters

rhos=4.0e3;       %density of suspension wire (kg/m^3) %sapphire

%mass (kg)
m1=20.5;          %IM
m2=0.055*4;       %blade spring
m3=mass;          %TM

%resonant frequency of blade spring (Hz)
wh2=2*pi*2e3;     %horizontal
fmg=14.5;         %vertical

%radius of wire (m)
r1=0.6e-3/2;      %wire1 between IM and MN
r2=1.6e-3/2;      %wire2 between blade spring and TM

%Youngs modulus (Pa)
Y1=1.3e11;        %wire1 %CuBe
Y2=4.0e11;        %blade spring %sapphire
Y3=4.0e11;        %wire2 %sapphire

%length of wire (m)
l1=26.1e-2;       %wire1
l2=0.5e-2;        %blade spring
l3=35e-2;         %wire2

%loss angle of wire
phi1=5e-6;        %wire1
phi2=7e-7;        %blade spring
phi3=2e-7;        %wire2

%temperature (K)
T1=16;            %wire1
T2=16;            %blade spring
T3=(T2+Tm)/2;     %wire2

%loss angle normalized by T3
phih1=phi1*T1/T3; %horizontal wire1
phih2=phi2*T2/T3; %horizontal blade spring
phih3=phi3;       %horizontal wire2
phiv1=phi1*T1/T3; %vertical wire1
phiv2=phi2*T2/T3; %vertical blade spring
phiv3=phi3;       %vertical wire2

%coupling from vertical mode to horizontal mode
VHC=1/200;

%number of wires
N1=4;
N2=4;
N3=4;

%number of bending points
b1=2;
b2=2;
b3=2;

%% DETERMINATION OF LASER POWER AT BS %%%%%%%%%%
%Power at BS is determined from sapphre mirror temperature and absorption.

subloss=50e-6*1e2;       %energy loss rate of substrate (/m)
%see JGW-G1706429-v1 (2017)
coaloss=0.5e-6;          %energy loss rate of coating
outside=0.05;            %radiation energy from outside (W)
%see CQG 31, 224003 (2014)

alphak=7.98;
betak=2.2;      %kappa=alphak*T^betak ,thermal conductivity (W/K/m)
%see CQG 31, 105004 (2014)

Kabs=N3*alphak/(betak+1)*pi*r2^2/l3*(Tm^(betak+1)-T2^(betak+1)); %total heatflow
Pmich=(Kabs-outside)/(2*subloss*height+coaloss*4/trans); %power between BS and ITM
I0=2*Pmich;          %power at BS
fprintf('Power at BS: %g W\n', I0);

%% NOISE CALCULATIONS %%%%%%%%%%
%%%%%seismic noise%%%%%
% fitting function of the seismic noise calculated by R. Takahashi.
hseisorig=2*2.*sqrt((4.1e-17.*f.^(-6.8)).^2+(1.9e-16.*f.^(-8)).^2+(1.8e-18.*f.^(-5.7)).^2);

% gravity gradient noise
% fitting function of the GGN noise calculated by J. Harms (see CQG 29, 124007 (2012)).
hGGN=1.2e-16./f.^4./L;

hseis=sqrt(hseisorig.^2+hGGN.^2);


%%%%%suspension thermal noise%%%%%
%see VIR-015E-09 (2009)
%violin mode parameters
lsus=l3;
Ewire=Y3*(1+i*phi3);   %complex Youngs modulus (Pa)
Icwire=pi*r2^4/4;      %area moment of inertia (m^2)
tension=mass*g/N3;     %tension per wire (N)
rholine=rhos*pi*r2^2;  %linear density (kg/m)

% wave number of elastic fibers (1/m)
ks=sqrt((-tension+sqrt(tension^2+4.*Ewire.*Icwire.*rholine.*omega.^2))./(2.*Ewire.*Icwire));

% wave number associated with the flexural stiffness
ke=sqrt((tension+sqrt(tension^2+4.*Ewire.*Icwire.*rholine.*omega.^2))./(2.*Ewire.*Icwire));

% effective spring constant (N/m)
% Calculated by K. Somiya and confirmed by K. Komori. The spring constant Eq.(48) in PPPnote (VIR-015E-09) is the limitting value of ke*lsus -> infinity in our case.
Ksusp=N3.*Ewire.*Icwire.*ke.*ks.*(ke.^2+ks.^2).*(ke.*cos(ks.*lsus).*sinh(ke.*lsus)+ks.*sin(ks.*lsus).*cosh(ke.*lsus))./(2.*ke.*ks.*(1-cos(ks.*lsus).*cosh(ke.*lsus))+(ke.^2-ks.^2).*sin(ks.*lsus).*sinh(ke.*lsus));

%spring constant and resonant frequency
Ecube=Y1*(1+i*phih1);          %complex Youngs modulus (Pa)
I1=pi*r1^4/4;                  %area moment of inertia (m^2)
tensioncube=(m1+m2+m3)*g/N1;   %tension per CuBe wire (N)

kh1=N1*tensioncube./(l1-2.*sqrt(Ecube.*I1./tensioncube));      %hor. spring constatnt of wire1 (N/m)
kh2=m2*wh2^2*(1+i*phih2);      %hor. spring constatnt of blade spring (N/m)
kh3=Ksusp;

wv1=sqrt(pi*r1^2*Y1*N1/(l1*m1));  %ver. wire1 (Hz)
kv1=m1*wv1^2*(1+i*phiv1);         %ver. spring constatnt of wire1 (N/m)

wv2=2*pi*fmg*sqrt((m2+m3)/m2);         %ver. blade spring (Hz)
kv2=(m2+m3)*(2*pi*fmg)^2*(1+i*phiv2);  %ver. spring constatnt of blade spring (N/m)

wv3=sqrt(pi*r2^2*Y3*N3/(l3*m3));  %ver. wire2 (Hz)
kv3=m3*wv3^2*(1+i*phiv3);         %ver. spring constatnt of wire2 (N/m)


%matrix of horizontal spring constatnt (N/m)
Dh11=-m1.*omega.^2+kh1+kh2;
Dh12=-kh2*ones(fres,1);
Dh13=0*ones(fres,1);
Dh21=Dh12;
Dh22=-m2.*omega.^2+kh2+kh3;
Dh23=-kh3;
Dh31=Dh13;
Dh32=Dh23;
Dh33=-m3.*omega.^2+kh3;

%matrix of horizontal impedance (N/m/Hz)
Zh11=Dh11./(i.*omega);
Zh12=Dh12./(i.*omega);
Zh13=Dh13./(i.*omega);
Zh21=Dh21./(i.*omega);
Zh22=Dh22./(i.*omega);
Zh23=Dh23./(i.*omega);
Zh31=Dh31./(i.*omega);
Zh32=Dh32./(i.*omega);
Zh33=Dh33./(i.*omega);

% matrix of vertical spring constatnt (N/m)
Dv11=-m1.*omega.^2+kv1+kv2;
Dv12=-kv2*ones(fres,1);
Dv13=0*ones(fres,1);
Dv21=Dv12;
Dv22=-m2.*omega.^2+kv2+kv3;
Dv23=-kv3*ones(fres,1);
Dv31=Dv13;
Dv32=Dv23;
Dv33=-m3.*omega.^2+kv3;

%matrix of vertical impedance (N/m/Hz)
Zv11=Dv11./(i.*omega);
Zv12=Dv12./(i.*omega);
Zv13=Dv13./(i.*omega);
Zv21=Dv21./(i.*omega);
Zv22=Dv22./(i.*omega);
Zv23=Dv23./(i.*omega);
Zv31=Dv31./(i.*omega);
Zv32=Dv32./(i.*omega);
Zv33=Dv33./(i.*omega);

for k=1:fres;
    %%calculation of inverse
    Zh(:,:,k)=[Zh11(k) Zh12(k) Zh13(k); Zh21(k) Zh22(k) Zh23(k); Zh31(k) Zh32(k) Zh33(k)];
    invZh(:,:,k)=Zh(:,:,k)^(-1);
    Zv(:,:,k)=[Zv11(k) Zv12(k) Zv13(k); Zv21(k) Zv22(k) Zv23(k); Zv31(k) Zv32(k) Zv33(k)];
    invZv(:,:,k)=Zv(:,:,k)^(-1);
    %%%calculation of eq.(63) in PPPnote (VIR-015E-09)%%%
    % The fluctuation of TM is derived from 3*3 component of the inverse matrix of impedance.
    Fsqh(:,:,k)=4*kb*T3.*real(invZh(:,:,k))./omega(k).^2;
    Sthh(k,1)=Fsqh(3,3,k);
    Fsqv(:,:,k)=4*kb*T3.*real(invZv(:,:,k))./omega(k).^2;
    Sthv(k,1)=Fsqv(3,3,k);
end

%strain of thermal noise with 4 TMs
hsusph=sqrt(abs(Sthh)).*2/L;
hsuspv=sqrt(abs(Sthv)).*VHC*2/L;
hsusp=sqrt(hsusph.^2+hsuspv.^2);


%%%%%%mirror thermal noise%%%%%
%strain of substrate thermal noise for 2 ITMs and ETMs
%see CQG, 29, 124007 (2012)
hmirsub=sqrt(2)/L.*sqrt(4*kb*Tm*(1-sigma^2)./(sqrt(pi)*E0*Qm.*omega).*(1/w1+1/w2));

%thermoelastic noise calculations
%see PRD, 82, 127101 (2010)
%normalized frequency
omegac1=omega.*Cth20*rhom*w1^2./(2*kappa20);
omegac2=omega.*Cth20*rhom*w2^2./(2*kappa20);

%calculation of Jhyper
%It is confirmed by K. Komori that this integration is the same as the analytical calculation using error function.
resolution=0.01;
up=10;    %upper limit of integration instead of infinity
bin=floor(up/resolution);
u=linspace(0,up,bin+1);

for k=1:fres;
    integral11(k,:)=u.^2.*exp(-u.^2./2).*(1-u./sqrt(2).*sqrt((u.^2+sqrt(u.^4+omegac1(k).^2))./(u.^4+omegac1(k).^2)));
    integral12(k,:)=integral11(k,:);
    integral21(k,:)=u.^2.*exp(-u.^2./2).*(1-u./sqrt(2).*sqrt((u.^2+sqrt(u.^4+omegac2(k).^2))./(u.^4+omegac2(k).^2)));
    integral22(k,:)=integral21(k,:);
end

integral11(:,1)=[];
integral12(:,bin)=[];
integral1=sum(integral11+integral12,2)./2.*resolution;

integral21(:,1)=[];
integral22(:,bin)=[];
integral2=sum(integral21+integral22,2)./2.*resolution;

Jhyper1int=sqrt(2/pi).*integral1./omegac1.^2;
Jhyper2int=sqrt(2/pi).*integral2./omegac2.^2;

%strain of thermoelastic noise for 2 ITMs and ETMs
hmirelaint=sqrt(2)/L.*sqrt(4*kb*Tm^2*alpha20^2*(1+sigma)^2/(sqrt(pi)*kappa20).*(w1.*Jhyper1int+w2.*Jhyper2int));

%coating thermal noise calculations
%see CQG, 19, 897 (2002)
%thermal noise of silica coating
hmircoas=sqrt(4*kb*Tm./(pi*E0^2.*omega)).*sqrt(phic1*(Ys^2*(1+sigma)^2*(1-2*sigma)^2+E0^2*(1+sigmas)^2*(1-2*sigmas))/(Ys*(1-sigmas^2))*(dcoa1s/w1^2+dcoa2s/w2^2));

%thermal noise of tantala coating
hmircoat=sqrt(4*kb*Tm./(pi*E0^2.*omega)).*sqrt(phic2*(Yt^2*(1+sigma)^2*(1-2*sigma)^2+E0^2*(1+sigmat)^2*(1-2*sigmat))/(Yt*(1-sigmat^2))*(dcoa1t/w1^2+dcoa2t/w2^2));

%strain of coating thermal noise for 2 ITMs and ETMs
hmircoa=sqrt(2)/L.*sqrt(hmircoas.^2+hmircoat.^2);

%total mirror thermal noise
hmir=sqrt(hmirsub.^2+hmirelaint.^2+hmircoa.^2);


%%%%%%quantum noise%%%%%
%see PRD, 64, 042006 (2001)
%parameters
gamma=trans*c/(4*L);       %cavity pole (Hz)
tSRM=sqrt(1-rSRM^2);       %SRM transmittance
beta=atan(omega./gamma);   %phase shift in the arm cavity

%laser power quantum noise reaches the SQL at the cavity pole
ISQL=mass*L^2*gamma^4/(4*omegal);

%coefficient showing how much the mirror moves
kappa=I0/ISQL*2*gamma^4./(omega.^2.*(gamma^2+omega.^2));

%loss parameter of arms defined in KLMTV paper
epsilon=2*roundtriploss/trans;


%quantum noise calculations
% strain of SQL
hsql=sqrt(8*hbar./(mass*omega.^2.*L^2));

%components of eq.(5.13) in PRD, 64, 042006 (2001)
C11=sqrt(1-lambdaPD).*((1+rSRM^2).*(cos(2*phidet)+kappa.*sin(2*phidet)/2)-2*rSRM.*cos(2.*beta)-epsilon/4.*(-2.*(1+exp(2*i.*beta)).^2.*rSRM+4*(1+rSRM^2).*cos(beta).^2.*cos(2*phidet)+(3+exp(2*i.*beta)).*kappa.*(1+rSRM^2)*sin(2*phidet))+lambdaSR.*(exp(2*i.*beta).*rSRM-(1+rSRM^2)/2.*(cos(2*phidet)+kappa.*sin(2*phidet)/2)));
C22=C11;
C12=sqrt(1-lambdaPD)*tSRM^2.*(-(sin(2*phidet)+kappa.*sin(phidet)^2)+epsilon*sin(phidet)/2.*((3+exp(2*i.*beta)).*kappa.*sin(phidet)+4.*cos(beta).^2.*cos(phidet))+lambdaSR/2.*(sin(2*phidet)+kappa.*sin(phidet)^2));
C21=sqrt(1-lambdaPD)*tSRM^2.*((sin(2*phidet)-kappa.*cos(phidet)^2)+epsilon*cos(phidet)/2.*((3+exp(2*i.*beta)).*kappa.*cos(phidet)-4.*cos(beta).^2.*sin(phidet))+lambdaSR/2.*(-sin(2*phidet)+kappa.*cos(phidet)^2));
D1sig=sqrt(1-lambdaPD).*(-(1+rSRM.*exp(2*i.*beta)).*sin(phidet)+epsilon/4.*(3+rSRM+2*rSRM.*exp(4*i.*beta)+exp(2*i.*beta).*(1+5*rSRM)).*sin(phidet)+lambdaSR/2.*exp(2*i.*beta).*rSRM*sin(phidet));
D2sig=sqrt(1-lambdaPD).*(-(-1+rSRM.*exp(2*i.*beta)).*cos(phidet)+epsilon/4.*(-3+rSRM+2*rSRM.*exp(4*i.*beta)+exp(2*i.*beta).*(-1+5*rSRM)).*cos(phidet)+lambdaSR/2.*exp(2*i.*beta).*rSRM*cos(phidet));
P11=sqrt(1-lambdaPD)/2*sqrt(lambdaSR)*tSRM.*(-2*rSRM.*exp(2*i.*beta)+2*cos(2*phidet)+kappa.*sin(2*phidet));
P22=P11;
P12=-sqrt(1-lambdaPD)*sqrt(lambdaSR)*tSRM*sin(phidet).*(2*cos(phidet)+kappa.*sin(phidet));
P21=sqrt(1-lambdaPD)*sqrt(lambdaSR)*tSRM*cos(phidet).*(2*sin(phidet)-kappa.*cos(phidet));
Q11=sqrt(lambdaPD).*(exp(-2*i.*beta)+rSRM^2.*exp(2*i.*beta)-rSRM.*(2*cos(2*phidet)+kappa.*sin(2*phidet))+epsilon/2*rSRM.*(exp(-2*i.*beta).*cos(2*phidet)+exp(2*i.*beta).*(-2*rSRM-2*rSRM.*cos(2.*beta)+cos(2*phidet)+kappa.*sin(2*phidet))+2*cos(2*phidet)+3.*kappa.*sin(2*phidet))-lambdaSR/2*rSRM.*(2*rSRM.*exp(2*i.*beta)-2*cos(2*phidet)-kappa.*sin(2*phidet)));
Q22=Q11;
Q12=0;
Q21=0;
N11=sqrt(1-lambdaPD)*sqrt(epsilon/2)*tSRM.*(kappa.*(1+rSRM.*exp(2*i.*beta)).*sin(phidet)+2.*cos(beta).*(exp(-i.*beta).*cos(phidet)-rSRM.*exp(i.*beta).*(cos(phidet)+kappa.*sin(phidet))));
N22=-sqrt(1-lambdaPD)*sqrt(2*epsilon)*tSRM.*(-exp(-i.*beta)+rSRM.*exp(i.*beta)).*cos(beta).*cos(phidet);
N12=-sqrt(1-lambdaPD)*sqrt(2*epsilon)*tSRM.*(exp(-i.*beta)+rSRM.*exp(i.*beta)).*cos(beta).*sin(phidet);
N21=sqrt(1-lambdaPD)*sqrt(epsilon/2)*tSRM.*(-kappa.*(1+rSRM).*cos(phidet)+2.*cos(beta).*(exp(-i.*beta)+rSRM.*exp(i.*beta)).*cos(beta).*sin(phidet));
M=0;

%power of shot noise
Ssh=hsql.^2./(2*tSRM^2.*kappa.*abs(D1sig.*sin(xi)+D2sig.*cos(xi)).^2).*(abs(C11.*sin(xi)+C21.*cos(xi)).^2+abs(C12.*sin(xi)+C22.*cos(xi)).^2+abs(P11.*sin(xi)+P21.*cos(xi)).^2+abs(P12.*sin(xi)+P22.*cos(xi)).^2+abs(Q11.*sin(xi)+Q21.*cos(xi)).^2+abs(Q12.*sin(xi)+Q22.*cos(xi)).^2+abs(N11.*sin(xi)+N21.*cos(xi)).^2+abs(N12.*sin(xi)+N22.*cos(xi)).^2+(abs(M).*sqrt(Sadd)).^2);

%These calculations are consistence with the result from Haino.

%strain of shot noise without squeezing
hsh=sqrt(Ssh);

%squeezing related parameters
Lfilter=100;               %filter cavity length (m)
kappafilter=2*pi*45;       %line width of filter cavity (Hz)
detune=kappafilter;       %detuning of filter cavity (Hz)
Ttot=kappafilter*4*Lfilter/c; %transmittance of filter cavity
squeezelevel=6;            %dB
rsq=10^(squeezelevel/20);
sqangle=90*pi/180;         %squeezing angle

%covariance matrix of vacuum fluctuations just after the squeezer
sq11=cosh(2*rsq)+sinh(2*rsq)*cos(2*sqangle);
sq12=sinh(2*rsq)*sin(2*sqangle);
sq22=cosh(2*rsq)-sinh(2*rsq)*cos(2*sqangle);

%loss
lambdasq1=0.05;   %between squeezer and filter cavity
lambdasq2=30e-6;  %inside the filter cavity
lambdasq3=0.05;   %between filter cavity and SRM

%quantum noise calculation with squeezing (introduced by Yutaro Enomoto)
%amplitude reflectivity for sidebands of the filter cavity
reflp=-sqrt(1-Ttot+lambdasq2)+((Ttot-lambdasq2)*sqrt(1-lambdasq2).*exp(i*2*Lfilter.*(detune+omega)./c))./(1-sqrt(1-Ttot).*exp(i*2*Lfilter.*(detune+omega)./c)); %upper
reflm=-sqrt(1-Ttot+lambdasq2)+((Ttot-lambdasq2)*sqrt(1-lambdasq2).*exp(i*2*Lfilter.*(detune-omega)./c))./(1-sqrt(1-Ttot).*exp(i*2*Lfilter.*(detune-omega)./c)); %lower

cov11=(1-lambdasq1)*(1-lambdasq3).*(abs(reflp+conj(reflm)).^2.*sq11/4-imag(reflp.*reflm).*sq12+abs(-reflp+conj(reflm)).^2.*sq22/4)+lambdasq1*(1-lambdasq3).*(abs(reflp).^2+abs(reflm).^2)./2+(1-lambdasq3).*(2-abs(reflp).^2-abs(reflm).^2)./2+lambdasq3;
cov12=(1-lambdasq1)*(1-lambdasq3).*(-(conj(reflp)+reflm).*(-reflp+conj(reflm)).*sq11/(4*i)+(abs(reflp+conj(reflm)).^2./4-abs(-reflp+conj(reflm)).^2./4).*sq12-(-conj(reflp)+reflm).*(reflp+conj(reflm)).*sq22/(4*i))+lambdasq1*(1-lambdasq3).*(abs(reflp).^2-abs(reflm).^2)./(2*i)+(1-lambdasq3).*(-abs(reflp).^2+abs(reflm).^2)./(2*i);
cov22=(1-lambdasq1)*(1-lambdasq3).*(abs(-reflp+conj(reflm)).^2.*sq11/4+imag(reflp.*reflm).*sq12+abs(reflp+conj(reflm)).^2.*sq22/4)+lambdasq1*(1-lambdasq3).*(abs(reflp).^2+abs(reflm).^2)./2+(1-lambdasq3).*(2-abs(reflp).^2-abs(reflm).^2)./2+lambdasq3;

%power of shot noise
Sshsq=hsql.^2./(2*tSRM^2.*kappa.*abs(D1sig.*sin(xi)+D2sig.*cos(xi)).^2).*(abs(C11.*sin(xi)+C21.*cos(xi)).^2.*cov11+abs(C12.*sin(xi)+C22.*cos(xi)).^2.*cov22+2.*real(conj(C11.*sin(xi)+C21.*cos(xi)).*(C12.*sin(xi)+C22.*cos(xi)).*cov12)+abs(P11.*sin(xi)+P21.*cos(xi)).^2+abs(P12.*sin(xi)+P22.*cos(xi)).^2+abs(Q11.*sin(xi)+Q21.*cos(xi)).^2+abs(Q12.*sin(xi)+Q22.*cos(xi)).^2+abs(N11.*sin(xi)+N21.*cos(xi)).^2+abs(N12.*sin(xi)+N22.*cos(xi)).^2+(abs(M).*sqrt(Sadd)).^2);

%strain of shot noise with squeezing
hshsq=sqrt(Sshsq);


%%%%%total noise%%%%%
htotal=sqrt(hseis.^2+hsusp.^2+hmir.^2+hsh.^2);


%% INSPIRAL RANGE %%%%%%%%%%
%parameters
%see http://asa.usno.navy.mil/SecK/Constants.html
msolar=1.9884e30;                               %mass of the sun (kg)
Gpc=3.085677581491367278913937957796471611e25;  %1Gpc (m)

%lower limit of frequency integration
fminset=10;
[resfmin,fminind]=min(abs(f-fminset));
fmin=f(fminind);

%threshold of signal to noise ratio
SNR=8;

%sky average constant given by Kanda-san and see PRD, 47, 2198 (1993).
averageconst=0.442;

%resoliution of mass for inspiral range
massreso=200;
mrange=cat(1,logspace(0,3,massreso)',1.4);  %last one for 1.4Msun

%integrated function for inspiral range
integfunc=f.^(-7/3)./htotal.^2;

%upper limit of frequency integration normalized by ISCO frequency
iscofactor=1;

%%integration
for m=1:massreso+1
    massstar1(m,1)=mrange(m).*msolar; %mass of star1 (kg)
    massstar2(m,1)=massstar1(m,1); %mass of star2 (kg)  
    mtotal(m,1)=massstar1(m,1)+massstar2(m,1); %total mass (kg)
    mchirp(m,1)=(massstar1(m,1).*massstar2(m,1)).^(3/5)./(massstar1(m,1)+massstar2(m,1)).^(1/5); %chirp mass (kg)

    fisco(m,1)=c^3./(6^(3/2)*pi*G.*mtotal(m,1)); %ISCO frequency
    [resfmax(m,1),fmaxind(m,1)]=min(abs(f-iscofactor*fisco(m,1)));
    fmax(m,1)=f(fmaxind(m,1));
    
    sensitivityinteg1=integfunc(fminind:fmaxind(m,1)-1,1);
    sensitivityinteg2=integfunc(fminind+1:fmaxind(m,1),1);
    fbinwidth1=f(fminind:fmaxind(m,1)-1,1);
    fbinwidth2=f(fminind+1:fmaxind(m,1),1);
    sensitivityinteg(m,1)=sum((sensitivityinteg1+sensitivityinteg2)./2.*(fbinwidth2-fbinwidth1));

    clear sensitivityinteg1 sensitivityinteg2 fbinwidth1 fbinwidth2

    % inspiral range (Gpc)
    %see eq 7.182 of M. Maggiore 'Gravitational Waves Volume 1', Oxford University, Press Inc. (2007)
    ir(m,1)=4*averageconst/SNR.*sqrt(sensitivityinteg(m,1)).*sqrt(5/96)*pi^(-2/3)*c.*(G.*mchirp(m,1)./c^3).^(5/6)./Gpc;
end

BNSIR=ir(massreso+1,1)*1e3;  % inspiral range for 1.4-1.4Msun BNS (Mpc)
fprintf('Inspiral range for BNS: %g Mpc\n', BNSIR);

%% PLOT GRAPHS AND SAVE DATA %%%%%%%%%%
set(0,'defaultAxesFontName','helvetica');
set(0,'defaultTextFontName','helvetica');
set(0,'DefaultLineMarkerSize',15);
set(0,'DefaultLineLineWidth',2);
set(0,'defaultAxesFontSize',15);
set(0,'defaultTextFontSize',15);

loglog(f,hseis,'m-')
hold on
%grid on
loglog(f,hsusp,'g-')
loglog(f,hmir,'b-')
loglog(f,hsh,'r-')
%loglog(f,hshsq,'r-')
loglog(f,htotal,'k-')
legend('seismic+gravity','suspension','mirror','quantum','total')
xlabel('Frequency [Hz]')
ylabel('Strain [/rHz]')
xlim([5,5000])
ylim([5e-25,5e-21])
grid on
saveas(gcf,['sensitivity_',RSE,'.png'])

%{
figure(3)
loglog(f,hmirsub,'m-')
hold on
loglog(f,hmirelaint,'g-')
loglog(f,hmircoa,'b')
set(gca,'FontSize',18)
legend('substrate','elastic','coating')
xlabel('Frequency [Hz]')
ylabel('Strain [/rHz]')
xlim([5,5000])
ylim([5e-25,5e-21])
%}

figure(2)
loglog(mchirp(1:massreso)./msolar,ir(1:massreso),'b-')
xlabel('Chirp mass [Msolar]')
ylabel('Inspiral range [Gpc]')
xlim([1,1000])
ylim([0.1,10])
title(['\fontsize{10} SNR=',num2str(SNR),', fmin=',num2str(fmin,'%.1f'),' Hz, skyaverage=',num2str(averageconst,'%.3f')])
grid on

%data output
data(:,1)=f;
data(:,2)=hseis;
data(:,3)=hsusp;
data(:,4)=hmir;
data(:,5)=hsh;
data(:,6)=hshsq;
data(:,7)=htotal;

fileID = fopen(['spectrum_',RSE,'.txt'],'w');
fprintf(fileID,['#',RSE,' ',num2str(BNSIR,'%.0f'), ' Mpc\n']);
fprintf(fileID,'#frequency seismic suspension_thermal mirror_thermal quantum quantum_squeezing total\n');
fprintf(fileID,'%1.3e %1.3e %1.3e %1.3e %1.3e %1.3e %1.3e\n',data');
fclose(fileID);

%% SENSITIVITY COMPARISON PLOT %%%%%%%%%%
%plot sensitivities of BRSE/BRSEhomo/DRSE if sensitivity files exist in current folder
plotflag=1;  % 1 for plot
config={'BRSE','BRSEhomo','DRSE'};
if plotflag && exist('spectrum_BRSE.txt', 'file') == 2 && exist('spectrum_BRSEhomo.txt', 'file') == 2 && exist('spectrum_DRSE.txt', 'file') == 2
    figure(100)
    leg={}; % legends
    for kk=1:3
        sens=importdata(['spectrum_',config{kk},'.txt']);
        d=textscan(sens.textdata{1},'%s');    % extract inspiral range info from comment line
        loglog(sens.data(:,1),sens.data(:,7))
        leg=cat(1,leg,[config{kk},' ',d{1}{2},' Mpc']);
        hold on
    xlabel('Frequency [Hz]')
    ylabel('Strain [/rHz]')
    xlim([5,5000])
    ylim([5e-25,5e-21])
    legend(leg)
    grid on
    end
    saveas(gcf,['sensitivity_all.png'])
    saveas(gcf,['sensitivity_all.eps'],'epsc')
end