Exit[]

SetDirectory[NotebookDirectory[]];
path=SetDirectory[NotebookDirectory[]];

Calculations of the chirp mass from the frequency-time map of GW170817

Import["GW170817.dat"];
{time,freq}={-%[[All,1]],%[[All,2]]};

(* --- This is the Time - Frequency map --- *)

xfreq=Join[Table[{10i,,{0.015,0}},{i,1,79}],Table[{50i,ToString[50i],{0.03,0}},{i,1,20}]];
xfreqn=Join[Table[{10i,,{0.015,0}},{i,1,79}],Table[{50i,,{0.03,0}},{i,1,20}]];

xt=Join[Table[{1i,,{0.015,0}},{i,1,100}],Table[{5i,ToString[5i],{0.03,0}},{i,1,100}]];
xtn=Join[Table[{1i,,{0.015,0}},{i,1,100}],Table[{5i,,{0.03,0}},{i,1,100}]];

ListPlot[Thread[{time,freq}],Frame->True,FrameStyle->Directive[16,GrayLevel[0],FontFamily->"American Typewriter",Italic,AbsoluteThickness[1.2]],FrameLabel->{{Style["𝒻 [Hz]",22],None},{Style["t-(Subscript[t, 0]) [s]",22],None}},PlotMarkers->{"●",14},PlotStyle->Black,PlotRange->All,FrameTicks->{{xfreq,xfreqn},{xt,xtn}}]

(* --- Fit of (f^(-8/3)) vs time to extact the chirp mass --- *)

fManalytic[ℳ_,t_]:=With[{Msun=1.98847 (10^30),c=299792,G=6.67408 (10^(-11)) (10^(-9))},((((8π))^(8/3))/5)(((((G ℳ Msun)/(c^3))))^(5/3))t]

Show[ListPlot[Thread[{time,(freq^(-8/3))}],Frame->True,FrameStyle->Directive[16,GrayLevel[0],FontFamily->"American Typewriter",Italic,AbsoluteThickness[1.2]],FrameLabel->{{Style["(𝒻^(-8/3))",22],None},{Style["t-(Subscript[t, 0]) [s]",22],None}},PlotStyle->Black,PlotRange->All,ImageSize->400,PlotMarkers->{"◦",22},PlotLegends->Placed[{Style["data"]},Scaled[{0.85,0.12}]]],Plot[{fManalytic[1.4,t],fManalytic[1.2,t],fManalytic[1,t]},{t,0.5,30},PlotStyle->{Directive[AbsoluteThickness[2]],Directive[Dashing[{0.03,0.03}],AbsoluteThickness[2]],Directive[Dashing[{0.05,0.05,0.02}],AbsoluteThickness[2]]},PlotLegends->Placed[{Style["ℳ = 1.4 (Subscript[M, ⊙])"],Style["ℳ = 1.2 (Subscript[M, ⊙])"],Style["ℳ = 1 (Subscript[M, ⊙])"]},Scaled[{0.23,0.78}]]]]

(* --- Compute the chirp mass point-by-point in the frequency-time map --- *)

νGW[Δt_]:=With[{c=299792,msun=1.47678,G=6.67408 (10^(-11)) (10^(-9))},((5^(3/8))/(8π))(((((c^3)/(G ℳ ))))^(5/8))(1/((Δt)^(3/8)))]

chripfreqs=Table[With[{Msun=1.98847 (10^30)},{freq[[i]],Solve[νGW[time[[i]]]==freq[[i]],ℳ][[1,1,2]]/Msun}],{i,1,Length[freq]}];

xM=Join[Table[{1.1+0.001i,,{0.015,0}},{i,1,300}],Table[{1+0.01i,ToString[1+0.01i],{0.03,0}},{i,1,100}]];
xMn=Join[Table[{1.1+0.001i,,{0.015,0}},{i,1,300}],Table[{1+0.01i,,{0.03,0}},{i,1,100}]];

xfreq=Join[Table[{40+0.5i,,{0.015,0}},{i,1,200}],Table[{2i,ToString[2i],{0.03,0}},{i,1,30}]];
xfreqn=Join[Table[{40+0.5i,,{0.015,0}},{i,1,200}],Table[{2i,,{0.03,0}},{i,1,30}]];

ListPlot[chripfreqs,Filling->Bottom,Frame->True,FrameStyle->Directive[16,GrayLevel[0],FontFamily->"American Typewriter",Italic,AbsoluteThickness[1.2]],FrameLabel->{{Style["ℳ [(Subscript[M, ⊙])]",20],None},{Style["Hz",22],None}},PlotMarkers->{"●",14},PlotStyle->Black,PlotRange->{{40,52.5},{1.19,1.23}},ImageSize->400,FrameTicks->{{xM,xMn},{xfreq,xfreqn}}]

(*--- This is to show a larger frequency band with ℳ varying --- *)
xM=Join[Table[{0.1i,,{0.015,0}},{i,1,300}],Table[{0.5i,ToString[0.5i],{0.03,0}},{i,1,100}]];
xMn=Join[Table[{0.1i,,{0.015,0}},{i,1,300}],Table[{0.5i,,{0.03,0}},{i,1,100}]];

xfreq=Join[Table[{10i,,{0.015,0}},{i,1,79}],Table[{50i,ToString[50i],{0.03,0}},{i,1,20}]];
xfreqn=Join[Table[{10i,,{0.015,0}},{i,1,79}],Table[{50i,,{0.03,0}},{i,1,20}]];

ListPlot[chripfreqs,Filling->Bottom,Frame->True,FrameStyle->Directive[16,GrayLevel[0],FontFamily->"American Typewriter",Italic,AbsoluteThickness[1.2]],FrameLabel->{{Style["ℳ [(Subscript[M, ⊙])]",20],None},{Style["Hz",22],None}},PlotMarkers->{"●",14},PlotStyle->Black,PlotRange->All,ImageSize->400,FrameTicks->{{xM,xMn},{xfreq,xfreqn}}]

(* --- Find the lower bound on the total mass of the system --- *)

Mm2[m2_]:=Solve[((((m1 m2))^(3/5))(M^(-1/5))//.m1->M-m2)==chripfreqs[[1,2]],M][[1,1,2]]

Mmin=Table[Mm2[m2],{m2,0.5,2,0.01}]//Min

xm=Join[Table[{0.05i,,{0.015,0}},{i,1,200}],Table[{0.2i,ToString[0.2i],{0.03,0}},{i,1,30}]];
xmn=Join[Table[{0.05i,,{0.015,0}},{i,1,200}],Table[{0.2i,,{0.03,0}},{i,1,30}]];

xl=Join[Table[{0.02i,,{0.015,0}},{i,1,200}],Table[{0.1i,ToString[0.1i],{0.03,0}},{i,1,30}]];
xln=Join[Table[{0.02i,,{0.015,0}},{i,1,200}],Table[{0.1i,,{0.03,0}},{i,1,30}]];

Show[Plot[Mm2[x],{x,0.5,2.2},PlotRange->{{0.9,2.2},{2.77,3}},Frame->True,FrameStyle->Directive[16,GrayLevel[0],FontFamily->"American Typewriter",Italic,AbsoluteThickness[1.2]],PlotStyle->Directive[Red,AbsoluteThickness[2]],FrameTicks->{{xl,xln},{xm,xmn}},Epilog->Inset[Style["(Subscript[M, min]) = "<>ToString[NumberForm[Mmin,3]],18,Black,FontFamily->"American Typewriter",Italic],Scaled[{0.78,.23}]],FrameLabel->{{Style["M [(Subscript[M, ⊙])]",22],None},{Style["(Subscript[m, 2]) [(Subscript[M, ⊙])]",22],None}}],ListPlot[{{0.5,Mmin},{2.4,Mmin}},Joined->True,PlotStyle->Directive[Black,Dashing[{0.03,0.03}]]]]

(* --- Find the orbital distance between the bodies and the velocity, assuming M=2.8(Subscript[M, ⊙]) --- *)

compct=With[{msun=1.4768,cc=299792,M=Mmin},Thread[{freq,(((((M msun)/(((π freq/cc))^2))))^(1/3))}]];

vel=With[{cc=299792},Thread[{freq,π freq compct[[All,2]]/2/cc}]];

xfreq=Join[Table[{10i,,{0.015,0}},{i,0,100}],Table[{50i,ToString[50i],{0.03,0}},{i,1,20}]];
xfreqn=Join[Table[{10i,,{0.015,0}},{i,0,100}],Table[{50i,,{0.03,0}},{i,1,20}]];

xv=Join[Table[{0.0025i,,{0.015,0}},{i,0,100}],Table[{0.01i,ToString[0.01i],{0.03,0}},{i,1,20}]];
xvn=Join[Table[{0.0025i,,{0.015,0}},{i,0,100}],Table[{0.01i,,{0.03,0}},{i,1,20}]];

GraphicsColumn[{ListPlot[compct,Frame->True,PlotRange->{{40,160},{Automatic,Automatic}},FrameStyle->Directive[16,GrayLevel[0],FontFamily->"American Typewriter",Italic,AbsoluteThickness[1.2]],FrameTicks->{{xfreq,xfreqn},{xfreqn,xfreqn}},PlotMarkers->{"●",8},PlotStyle->Directive[Black,AbsoluteThickness[2]],FrameLabel->{{Style["(Subscript[ℓ, 0]) [km]",22],None},{None,None}},ImagePadding->{{70,10},{20,10}},AspectRatio->0.5,ImageSize->400],

ListPlot[vel,Frame->True,PlotRange->{{40,160},{Automatic,Automatic}},FrameStyle->Directive[16,GrayLevel[0],FontFamily->"American Typewriter",Italic,AbsoluteThickness[1.2]],FrameTicks->{{xv,xvn},{xfreq,xfreqn}},PlotMarkers->{"●",8},PlotStyle->Directive[Black,AbsoluteThickness[2]],FrameLabel->{{Style["v/c",22],None},{Style["𝒻 [Hz]",22],None}},ImagePadding->{{70,10},{20,10}},AspectRatio->0.5,ImageSize->400]},Spacings->-30]

Binary parameters for GW150914

(* Adaptded from International Physics Olimpiad 2018 *)

(* --- The GW freqeuncy is (Subscript[ν, GW]) = 1/T = 1/(2Δt) where T is the period. We have 2 couples of data ((Subscript[Δt, i]),(Subscript[f, i])) for 2 unknown variables given by ((Subscript[t, 0]),ℳ). We can then solve the system for the chirp mass. --- *)

t={0,0.009,0.02,0.0275,0.034,0.04};

{tAB,tCD,tDE}={t[[2]]+t[[1]],t[[3]]+t[[4]],t[[5]]+t[[6]]}0.5;
{ΔtAB,ΔtCD,ΔtDE}={t[[2]]-t[[1]],t[[4]]-t[[3]],t[[6]]-t[[5]]};
{fAB,fCD,fDE}={0.5/ΔtAB,0.5/ΔtCD,0.5/ΔtDE};

eq[t_,f_]:=With[{Msun=1.98847 (10^30),c=299792,G=6.67408 (10^(-11)) (10^(-9))},((((8π))^(8/3))/5)(((((G ℳ Msun)/(c^3))))^(5/3))(t0-t)-(f^(-8/3))]

(* --- See the three cases --- *)

{ℳest,t0est}=NSolve[{eq[tAB,fAB],eq[tDE,fDE]},{ℳ,t0},Reals][[1,All,2]]

{ℳest,t0est}=NSolve[{eq[tCD,fCD],eq[tDE,fDE]},{ℳ,t0},Reals][[1,All,2]]

{ℳest,t0est}=NSolve[{eq[tAB,fAB],eq[tCD,fCD]},{ℳ,t0},Reals][[1,All,2]]

(* --- Find the lower bound on the total mass of the system --- *)

Mm2[m2_]:=Solve[((((m1 m2))^(3/5))(M^(-1/5))//.m1->M-m2)==ℳest,M][[1,1,2]]

Mmin=Table[Mm2[m2],{m2,10,80,0.1}]//Min

xm=Join[Table[{2i,,{0.015,0}},{i,1,200}],Table[{10i,ToString[10i],{0.03,0}},{i,1,30}]];
xmn=Join[Table[{2i,,{0.015,0}},{i,1,200}],Table[{10i,,{0.03,0}},{i,1,30}]];

xl=Join[Table[{2i,,{0.015,0}},{i,1,200}],Table[{10i,ToString[10i],{0.03,0}},{i,1,30}]];
xln=Join[Table[{2i,,{0.015,0}},{i,1,200}],Table[{10i,,{0.03,0}},{i,1,30}]];

Show[Plot[Mm2[x],{x,15,80},PlotRange->{{15,80},{64,98}},Frame->True,FrameStyle->Directive[16,GrayLevel[0],FontFamily->"American Typewriter",Italic,AbsoluteThickness[1.2]],PlotStyle->Directive[Red,AbsoluteThickness[2]],Epilog->Inset[Style["(Subscript[M, min]) = "<>ToString[NumberForm[Mmin,3]],18,Black,FontFamily->"American Typewriter",Italic],Scaled[{0.78,.22}]],FrameLabel->{{Style["M [(Subscript[M, ⊙])]",22],None},{Style["(Subscript[m, 2]) [(Subscript[M, ⊙])]",22],None}},FrameTicks->{{xm,xmn},{xl,xln}}],ListPlot[{{10,Mmin},{90,Mmin}},Joined->True,PlotStyle->Directive[Black,Dashing[{0.03,0.03}]]]]

(* --- Separation in km at the two freqeuncies we have computed --- *)
D1=With[{msun=1.4768,cc=299792,M=Mmin},(((((M msun)/(((π fAB/cc))^2))))^(1/3))]

D2=With[{msun=1.4768,cc=299792,M=Mmin},(((((M msun)/(((π fCD/cc))^2))))^(1/3))]

D3=With[{msun=1.4768,cc=299792,M=Mmin},(((((M msun)/(((π fDE/cc))^2))))^(1/3))]

(* --- Compute the orbital velocity at the specific frequecy [km/s] --- *)

v1=(π fAB)D1/2
%/299792 (* ratio with the speed of light *)

v2=(π fCD)D2/2
%/299792 (* ratio with the speed of light *)

v3=(π fDE)D3/2
%/299792 (* ratio with the speed of light *)