Modified for alpha and beta

code/RDownScriptalphabeta.m

+(* ::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="/Users/RayneLiu/Desktop/Research/AEI/Simulations/SXS/BBH_SKS_d14.3_q1.22_sA_0_0_0.330_sB_0_0_-0.440";
+tshift=0;
+npoints=100;
+ntones=7;
+allspectrum = {x0,x1,x2,x3,x4,x5,x6,x7};
+spectrum = allspectrum[[;;ntones+1]];
+\[Omega]fact=0.1;
+\[Tau]fact=0.1;
+params = {tshift, ntones, \[Omega]fact, \[Tau]fact}
+
+
+combinations = DeleteCases[Tuples[Range[0, ntones],2],Alternatives@@Subsets[Reverse[Range[0, ntones]],{2}]]
+
+
+(* 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];
+
+
+t0=TimeOfMaximum[sxsrhs[[1]]]+tshift
+data=Select[sxsrhs[[1]],#[[1]]>= t0&];
+
+
+tab\[Alpha]s=Table[Table[randomvar1= RandomReal[{-\[Omega]fact, \[Omega]fact}];randomvar2 = RandomReal[{-\[Omega]fact, \[Omega]fact}];
+ansatz=OvertoneModelV2[ntones,{\[Eta],\[Chi]1,\[Chi]2},t0,"Fit\[Alpha]"->combinations[[k]],"Fit\[Tau]"->{},"ModesData"->Modedata]/.ToExpression["\[Alpha]"<>ToString[combinations[[k]][[1]]]]->randomvar1/.ToExpression["\[Alpha]"<>ToString[combinations[[k]][[2]]]] -> randomvar2;
+cfit=NonlinearModelFit[data,ansatz,spectrum,t];
+cfitd=Transpose[{data[[All,1]],Normal[cfit]/.t->data[[All,1]]}];
+comp\[Alpha]s = spectrum/.cfit["BestFitParameters"];
+{{randomvar1, randomvar2}, combinations[[k]], (*ansatz, cfit["BestFitParameters"], comp\[Alpha]s,*) Re[Sqrt[comp\[Alpha]s * Conjugate[comp\[Alpha]s]]], Mod[Arg[comp\[Alpha]s], 2 Pi], {tshift, 1-EasyMatchT[data,cfitd,t0,t0+90]}},{j,npoints}], {k, Length[combinations]}];
+
+
+tab\[Alpha]\[Beta]s=Flatten[Table[Table[randomvar1= RandomReal[{-\[Omega]fact, \[Omega]fact}];randomvar2 = RandomReal[{-\[Tau]fact, \[Tau]fact}];
+ansatz=OvertoneModelV2[ntones,{\[Eta],\[Chi]1,\[Chi]2},t0,"Fit\[Alpha]"->{k},"Fit\[Tau]"->{l},"ModesData"->Modedata]/.ToExpression["\[Alpha]"<>ToString[k]]->randomvar1/.ToExpression["\[Beta]"<>ToString[l]] -> randomvar2;
+cfit=NonlinearModelFit[data,ansatz,spectrum,t];
+cfitd=Transpose[{data[[All,1]],Normal[cfit]/.t->data[[All,1]]}];
+comp\[Alpha]\[Beta]s = spectrum/.cfit["BestFitParameters"];
+{{randomvar1, randomvar2}, {k, l}, (*ansatz, cfit["BestFitParameters"], comp\[Alpha]\[Beta]s,*) Re[Sqrt[comp\[Alpha]\[Beta]s * Conjugate[comp\[Alpha]\[Beta]s]]], Mod[Arg[comp\[Alpha]\[Beta]s], 2 Pi], {tshift, 1-EasyMatchT[data,cfitd,t0,t0+90]}}, {j,npoints}], {k, Range[0, ntones]}, {l, Range[0, ntones]}], 1];
+
+
+tab\[Beta]s=Table[Table[randomvar1= RandomReal[{-\[Tau]fact, \[Tau]fact}];randomvar2 = RandomReal[{-\[Tau]fact, \[Tau]fact}];
+ansatz=OvertoneModelV2[ntones,{\[Eta],\[Chi]1,\[Chi]2},t0,"Fit\[Alpha]"->{},"Fit\[Tau]"->combinations[[k]],"ModesData"->Modedata]/.ToExpression["\[Beta]"<>ToString[combinations[[k]][[1]]]]->randomvar1/.ToExpression["\[Beta]"<>ToString[combinations[[k]][[2]]]] -> randomvar2;
+cfit=NonlinearModelFit[data,ansatz,spectrum,t];
+cfitd=Transpose[{data[[All,1]],Normal[cfit]/.t->data[[All,1]]}];
+comp\[Beta]s = spectrum/.cfit["BestFitParameters"];
+{{randomvar1, randomvar2}, combinations[[k]], (*ansatz, cfit["BestFitParameters"], comp\[Beta]s,*) Re[Sqrt[comp\[Beta]s * Conjugate[comp\[Beta]s]]], Mod[Arg[comp\[Beta]s], 2 Pi], {tshift, 1-EasyMatchT[data,cfitd,t0,t0+90]}},{j,npoints}], {k, Length[combinations]}];
+
+
+Export[rootpath<>"plots/params.fits", params]
+Export[rootpath<>"plots/freq_bfitpars_mmatch_as.fits", tab\[Alpha]s]
+Export[rootpath<>"plots/freq_bfitpars_mmatch_abs.fits", tab\[Alpha]\[Beta]s]
+Export[rootpath<>"plots/freq_bfitpars_mmatch_bs.fits", tab\[Beta]s]
+
+
+