Skip to content
Snippets Groups Projects
Select Git revision
  • 64ce284db98c7ce9d66b278cba720b0ccbf03fac
  • master default protected
2 results

RDownScript.m

Blame
  • Forked from Xisco Jimenez Forteza / RDStackingProject
    Source project has a limited visibility.
    RDownScript.m 4.44 KiB
    (* ::Package:: *)
    
    (* ::Title:: *)
    (*RDown excitation time*)
    
    
    (* ::Input:: *)
    (*Quit[]*)
    
    
    (* ::Section:: *)
    (*Setup*)
    
    
    (* ::Subsection::Closed:: *)
    (*Load package and set root paths*)
    
    
    rootpath="/work/francisco.jimenez/sio/git/rdstackingproject/";
    
    
    SetDirectory[rootpath<>"code"];
    <<RDown.m
    <<mcmc.m
    
    
    mysxscase=rootpath<>"SXS/BBH_SKS_d14.3_q1.22_sA_0_0_0.330_sB_0_0_-0.440";
    tshift=0;
    npoints=10;
    ntones=1;
    \[Omega]fact=0.05;
    
    
    (* ::Subsection::Closed:: *)
    (*Setup paths and other notebook folders*)
    
    
    (* 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
    
    
    (* ::Subsection::Closed:: *)
    (*Plot Formatting*)
    
    
    (* PlotStyle *)
    markers={"\[FilledCircle]","\[FilledSquare]","\[FilledDiamond]","\[FilledUpTriangle]","\[FilledDownTriangle]","\[EmptyCircle]","\[EmptySquare]","\[EmptyDiamond]","\[EmptyUpTriangle]","\[EmptyDownTriangle]"}
    colors=ColorData[97,"ColorList"]
    l1=.025;
    l2=.001;
    s=l1/2;
    
    
    (* ::Subsection::Closed:: *)
    (*Other functions*)
    
    
    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)]
    ]
    
    
    (* ::Subsection:: *)
    (*Imports*)
    
    
    (* ::Subsubsection::Closed:: *)
    (*RDown data*)
    
    
    (* ::Input:: *)
    (*(* The QNM spectrum is obtained from https://pages.jh.edu/~eberti2/ringdown/. Up to now, we can load any lmn combination up to l=4 and n=7. To extend the l number you need to download the l>4 files. *)*)
    
    
    (* 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;
    
    
    (* ::Subsubsection::Closed:: *)
    (*NR data + EOBFits*)
    
    
    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];
    
    
    (* ::Section:: *)
    (*Results for the paper*)
    
    
    (* ::Subsubsection::Closed:: *)
    (*Vary the fundamental frequency*)
    
    
    t0=TimeOfMaximum[sxsrhs[[1]]]+tshift+20
    data=Select[sxsrhs[[1]],#[[1]]>= t0&];
    
    
    tab=Table[
    ansatz=OvertoneModelV2[ntones,{\[Eta],\[Chi]1,\[Chi]2},t0,"Fit\[Alpha]"->{},"Vary\[Omega]"->True,"\[Omega]val"->{-\[Omega]fact,\[Omega]fact},"Export_\[Omega]val"->True,"ModesData"->Modedata];
    cfit=NonlinearModelFit[data,ansatz[[1]],{x0,x1,x2,x3,x4,x5,x6,x7},t];
    cfitd=Transpose[{data[[All,1]],Normal[cfit]/.t->data[[All,1]]}];
    {ansatz[[2]],{x0,x1,x2,x3,x4,x5,x6,x7}/.cfit["BestFitParameters"],1-EasyMatchT[data,cfitd,t0,t0+90]},{j,npoints}];
    
    
    Export[rootpath<>"results_data/freq_bfitpars_mmatch.dat",tab]