Print[Style["Phase noise Advanced Virgo","Section"]]
Directory[]
SetOptions[Plot,BaseStyle->{FontFamily->"Times",FontSize->14}];
kB=1.3806503 (10^(-23)); (* m2 kg s-2 K-1 *)
AMU=1.660538782 (10^(-27));  (* kg *)
NA = 6.0221415 (10^23);  (* mol-1 *)
R = 8.314472 (10^(-2)); (* m3 mbar-1 mol-1 K-1 *)
mN2 = 28 AMU;
mH2O = 18 AMU;
mH2 = 2 AMU;
mC = 44 AMU;
PH2O = 1.5 (10^(-10)); (* mbar *)
PN2 = 5.0 (10^(-10)); (* mbar *)
PH2 = 1.0 (10^(-9)); (* mbar *)
PC = 1.0 (10^(-13)); (* mbar *)
P=PH2O+PN2+PH2+PC;
T = 300;(* K *)
L = 3000; (* meters *)
AlphaN2 =  1.6 (10^(-30)); (* m3 *)
AlphaH2O =  1.45 (10^(-30)); (* m3 *)
AlphaH2 =  0.37 (10^(-30)); (* m3 *)
AlphaC =  5.3 (10^(-30)); (* m3 *)
Diam = 1.2 ;(* meters *)
v0N2=Sqrt[ 2kB T / mN2 ];
v0H2O=Sqrt[ 2kB T / mH2O ];
v0H2=Sqrt[ 2kB T / mH2 ];
v0C=Sqrt[ 2kB T / mC ];
rhoH2O = NA PH2O/(R T);
rhoN2 = NA PN2/(R T);
rhoH2 = NA PH2/(R T);
rhoC = NA PC/(R T);
wZ = 50 (10^(-3));(* m *)
w0 = 8.5 (10^(-3));(* m *)
z0=1385;(* m *)
zR=213;(* m *)
Print["Partial pressure H2O [ mbar ]: ", PH2O]
Print["Partial pressure N2 [ mbar ]: ", PN2]
Print["Partial pressure H2 [ mbar ]: ", PH2]
Print["Partial pressure hydrocarbons [ mbar ]: ", PC]
Print["Total pressure [ mbar ]: ", P]
Print["Temperature [ K ]: ", T]
Print["Molecular mass H2O [ amu ]: ", mH2O/AMU]
Print["Molecular mass N2 [ amu ]: ", mN2/AMU]
Print["Molecular mass H2 [ amu ]: ", mH2/AMU]
Print["Molecular mass hydrocarbons [ amu ]: ", mC/AMU]
Print["Molecular polarizability H2O [ m3 ]: ", AlphaH2O]
Print["Molecular polarizability N2 [ m3 ]: ", AlphaN2]
Print["Molecular polarizability H2 [ m3 ]: ", AlphaH2]
Print["Molecular polarizability hydrocarbons [ m3 ]: ", AlphaC]
Print["Cavity length [m ]: ", L]
Print["Average beam waist [ mm ]: ", 1000 wZ]
Print["Rayleigh length [ m ]: ", zR]
Print["Beam waist [ mm ]: ", 1000 w0]
Print["Beam waist position [ m from IM ]: ", z0]
Print["Beam radius on IM [ mm ]: ",1000 w0 Sqrt[1+ z0^2/zR^2]]
Print["Beam radius on EM [ mm ]: ",1000w0 Sqrt[1+ (L-z0)^2/zR^2]]
Print["Molecular velocity H2O [ m/s ]: ", v0H2O]
Print["Molecular velocity N2 [ m/s ]: ", v0N2]
Print["Molecular velocity H2 [ m/s ]: ", v0H2]
Print["Molecular velocity hydrocarbons [ m/s ]: ", v0C]
Print["Particle density H2O [ #/cm3 ]: ", rhoH2O/1000000]
Print["Particle density N2 [ #/cm3 ]: ", rhoN2/1000000]
Print["Particle density H2 [ #/cm3 ]: ", rhoH2/1000000]
Print["Particle density hydrocarbons [ #/cm3 ]: ", rhoC/1000000]

fstart = 10;
fend = 10000;
LowBound = 2.0*(10^(-26));
HighBound = 1*(10^(-21));
(*
Ntotal=10;
*)
AdV=ReadList["D:\\nikhef\\Experiments\\APP\\Virgo\\Technical projects\\Vacuum\\modeling\\phase noise\\AdV_refsens.txt",{Number,Number}];
etB=ReadList["D:\\nikhef\\Experiments\\APP\\Virgo\\Technical projects\\Vacuum\\modeling\\phase noise\\ET_B_data.txt",{Number,Number}];
etC=ReadList["D:\\nikhef\\Experiments\\APP\\Virgo\\Technical projects\\Vacuum\\modeling\\phase noise\\ET_C_data.txt",{Number,Number}];

freq=300;
Sqrt[8 rhoH2O (2 Pi AlphaH2O )^2 NIntegrate[Exp[-2 Pi freq w0 Sqrt[1+(z-z0)^2/(zR)^2] / v0H2O ]/ (v0H2O w0 Sqrt[1+(z-z0)^2/(zR)^2]),{z,0,3000}]+ 8 rhoN2 (2 Pi AlphaN2 )^2 NIntegrate[Exp[-2 Pi freq w0 Sqrt[1+(z-z0)^2/(zR)^2] / v0N2 ]/ (v0N2 w0 Sqrt[1+(z-z0)^2/(zR)^2]),{z,0,3000}]+ 8 rhoH2 (2 Pi AlphaH2 )^2 NIntegrate[Exp[-2 Pi freq w0 Sqrt[1+(z-z0)^2/(zR)^2] / v0H2 ]/ (v0H2 w0 Sqrt[1+(z-z0)^2/(zR)^2]),{z,0,3000}]+ 8 rhoC (2 Pi AlphaC )^2 NIntegrate[Exp[-2 Pi freq w0 Sqrt[1+(z-z0)^2/(zR)^2] / v0C ]/ (v0C w0 Sqrt[1+(z-z0)^2/(zR)^2]),{z,0,3000}]]/L
Sqrt[8 rhoH2O (2 Pi AlphaH2O )^2 NIntegrate[Exp[-2 Pi freq w0 Sqrt[1+(z-z0)^2/(zR)^2] / v0H2O ]/ (v0H2O w0 Sqrt[1+(z-z0)^2/(zR)^2]),{z,0,3000}]]/L
Sqrt[8 rhoN2 (2 Pi AlphaN2 )^2 NIntegrate[Exp[-2 Pi freq w0 Sqrt[1+(z-z0)^2/(zR)^2] / v0N2 ]/ (v0N2 w0 Sqrt[1+(z-z0)^2/(zR)^2]),{z,0,3000}]]/L
Sqrt[ 8 rhoH2 (2 Pi AlphaH2 )^2 NIntegrate[Exp[-2 Pi freq w0 Sqrt[1+(z-z0)^2/(zR)^2] / v0H2 ]/ (v0H2 w0 Sqrt[1+(z-z0)^2/(zR)^2]),{z,0,3000}]]/L
Sqrt[8 rhoC (2 Pi AlphaC )^2 NIntegrate[Exp[-2 Pi freq w0 Sqrt[1+(z-z0)^2/(zR)^2] / v0C ]/ (v0C w0 Sqrt[1+(z-z0)^2/(zR)^2]),{z,0,3000}]]/L