Select Git revision
RDownScriptalphabeta.m
Forked from
Xisco Jimenez Forteza / RDStackingProject
200 commits behind the upstream repository.
RDownScriptalphabeta.m 6.15 KiB
(* ::Package:: *)
(************************************************************************)
(* This file was generated automatically by the Mathematica front end. *)
(* It contains Initialization cells from a Notebook file, which *)
(* typically will have the same name as this file except ending in *)
(* ".nb" instead of ".m". *)
(* *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs. Doing so is equivalent *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end. *)
(* *)
(* DO NOT EDIT THIS FILE. This entire file is regenerated *)
(* automatically each time the parent Notebook file is saved in the *)
(* Mathematica front end. Any changes you make to this file will be *)
(* overwritten. *)
(************************************************************************)
rootpath="~/git/rdstackingproject/";
SetDirectory[rootpath<>"code"];
<<RDown.m
<<mcmc.m
SetOptions[InputNotebook[], AutoGeneratedPackage -> Automatic]
mysxscase=rootpath<>"/SXS/BBH_SKS_d14.3_q1.22_sA_0_0_0.330_sB_0_0_-0.440";
ntones=2;
tshift=12;
npoints=37000;
(*n=1:19, n=2:12, n=3:8, n=4:5, n=5:3, n=6:1, n=7:0 (not sure; is that really a minimum???) *)
allspectrum = {x0,x1,x2,x3,x4,x5,x6,x7};
allalphas = {\[Alpha]1,\[Alpha]2,\[Alpha]3,\[Alpha]4,\[Alpha]5,\[Alpha]6,\[Alpha]7};
allbetas = {\[Beta]1,\[Beta]2,\[Beta]3,\[Beta]4,\[Beta]5,\[Beta]6,\[Beta]7};
spectrum = allspectrum[[;;ntones+1]]
alphas = allalphas[[;;ntones]]
betas = allbetas[[;;ntones]]
\[Omega]fact=0.1;
\[Tau]fact=0.1;
params = {tshift, ntones, \[Omega]fact, \[Tau]fact}
(* Change them accordingly *)
(*Notebook directory*)
notdir=rootpath<>"code/";
(*Directory of the l's Ringdown data *)
modedir=rootpath;
(* True if you want to export the plots *)
ExportQ=False
ExportDir=notdir
(* PlotStyle *)
markers={"\[FilledCircle]","\[FilledSquare]","\[FilledDiamond]","\[FilledUpTriangle]","\[FilledDownTriangle]","\[EmptyCircle]","\[EmptySquare]","\[EmptyDiamond]","\[EmptyUpTriangle]","\[EmptyDownTriangle]"}
colors=ColorData[97,"ColorList"]
l1=.025;
l2=.001;
s=l1/2;
Options[pdf]={"SNR"->1};
pdf[ansatz_,data_,tend_?NumericQ,x_?(VectorQ[#,NumericQ]&),OptionsPattern[]]:=Block[{ansl,cfit,cfitd,h1red,h2red,norm1,norm2,myTable,snr,tmax},
snr=OptionValue["SNR"];
tmax=data[[1,1]];
ansl=ansatz;
cfitd=Transpose[{data[[All,1]],ansl/.t->data[[All,1]]}];
h1red=Select[data,tmax<=#[[1]]<=tmax+tend&][[All,2]];
h2red=Select[cfitd,tmax<=#[[1]]<=tmax+tend&][[All,2]];
norm1=Total[(Abs@h1red)^2 ];
norm2=Total[(Abs@h2red)^2 ];
myTable=h1red Conjugate@h2red;
-snr^2(1-(Re@Total[myTable]/Sqrt[norm1 norm2]))]
Options[pdfFit]={"SNR"->1};
pdfFit[ansatz_,data_,tend_?NumericQ,x_?(VectorQ[#,NumericQ]&),OptionsPattern[]]:=Block[{ansl,cfit,cfitd,h1red,h2red,norm1,norm2,myTable,snr,tmax},
snr=OptionValue["SNR"];
tmax=data[[1,1]];
ansl=ansatz;
cfitd=Transpose[{data[[All,1]],ansl/.t->data[[All,1]]}];
h1red=Select[data,tmax<=#[[1]]<=tmax+tend&][[All,2]];
h2red=Select[cfitd,tmax<=#[[1]]<=tmax+tend&][[All,2]];
Total[(Re[h1red-h2red]^2+Im[h1red-h2red]^2)]
]
(* Import all the modes you want to use *)
modelist={{2,2,0},{2,2,1},{2,2,2},{2,2,3},{2,2,4},{2,2,5},{2,2,6},{2,2,7}}
(* Files for the modes. Note that in this file notation the fundamental tone corresponds to n=1. However, for the code notation we will always use n=0. *)
Modefiles=Table[modedir<>"l"<>ToString[modelist[[j,1]]]<>"/n"<>ToString[modelist[[j,3]]+1]<>"l"<>ToString[modelist[[j,1]]]<>"m"<>ToString[modelist[[j,2]]]<>".dat",{j,Length@modelist}]
(* Import modes data *)
Modedata=Import/@Modefiles;
modes={{2,2}}
mysxscaserh=Reverse[Select[FileNames["*Lev6/rh*",mysxscase,2],Not@StringMatchQ[#,"*raw*"]&]]
mysxscasemetafile=FileNames["metadata.txt",mysxscase,2][[1]]
metadata=SXSMetaFilesToRules[mysxscasemetafile];
Print["mass1 = ",mass1=("initial-mass1"/.metadata)[[1]]]
Print["mass2 = ",mass2=("initial-mass2"/.metadata)[[1]]]
Print["\[Chi]1 = ",\[Chi]1=Chop[(("initial-dimensionless-spin1"/.metadata)[[-1]])]]
Print["\[Chi]2 = ",\[Chi]2=Chop[(("initial-dimensionless-spin2"/.metadata)[[-1]])]]
Print["m1/m2 = ",q=Max[{mass1/mass2,1}]]
Print["m1\[CenterDot]m2 = ",\[Eta]=q/(1+q)^2*1.]
Print["af = ",af=("remnant-dimensionless-spin"/.metadata)[[-1]]]
Print["mf = ",mf=("remnant-mass"/.metadata)[[1]]]
Print["Tag = ",tag=("alternative-names"/.metadata)[[2]]]
sxsrhs=Flatten[Conjugate/@GetAsymptoticMultiMode[#,2,modes,"ReSample"->True]&/@mysxscaserh,1];
tmax = TimeOfMaximum[sxsrhs[[1]]]
t0 = tmax+tshift
data = Select[sxsrhs[[1]],#[[1]]>= t0&];
(*mysxscase22modeRD=Select[sxsrhs[[1]],#[[1]]\[GreaterEqual] tmax-25&];
data1 = Select[mysxscase22modeRD, #[[1]]\[GreaterEqual]tmax+tshift&]*)
tabvars=Table[
randomvar\[Alpha]s= RandomReal[{-\[Omega]fact, \[Omega]fact}, ntones];
randomvar\[Beta]s = RandomReal[{-\[Omega]fact, \[Omega]fact}, ntones];
ansatz=OvertoneModelV2[ntones,{\[Eta],\[Chi]1,\[Chi]2},t0,"Fit\[Alpha]"->Range[1, ntones],"Fit\[Tau]"->Range[1, ntones],"ModesData"->Modedata]
/.Table[alphas[[i]]->randomvar\[Alpha]s[[i]],{i,ntones}]/.Table[betas[[j]]->randomvar\[Beta]s[[j]],{j,ntones}];
cfit=NonlinearModelFit[data,ansatz,spectrum,t];
cfitd=Transpose[{data[[All,1]],Normal[cfit]/.t->data[[All,1]]}];
alphasandmis = Join[alphas/.Table[alphas[[i]]->randomvar\[Alpha]s[[i]],{i,ntones}], {1-EasyMatchT[data,cfitd,t0,t0+90]}];
(* This extra step is to make sure everything's correct, and also to make the output a uniform shape *)
betasandtfit = Join[betas/.Table[betas[[j]]->randomvar\[Beta]s[[j]],{j,ntones}], {tshift}];
complexamps = spectrum/.cfit["BestFitParameters"];
amplitudes = Re[Sqrt[complexamps * Conjugate[complexamps]]];
phases = Mod[Arg[complexamps], 2 Pi];
{alphasandmis, betasandtfit, amplitudes, phases}, {j,npoints}];
Export[rootpath<>"plots/n="<>ToString[ntones]<>"_params.fits", params]
Export[rootpath<>"plots/n="<>ToString[ntones]<>"_t0="<>ToString[tshift]<>"M_data.fits", tabvars]