Select Git revision
Forked from
finesse / pykat
Source project has a limited visibility.
RDown.m 40.93 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. *)
(************************************************************************)
(* Probably not all the extra packages are really needed *)
BeginPackage["RDown`",{"HypothesisTesting`"}];
EasyMatchT::usage="EasyMatchT[h1_,h2_,tMin_,tMax_]. Time domain version of the match as in eq. (2) of 1903.08284. tMin and tMax are the time interval where to compute the match.";
PlanckTaper::usage="PlanckTaper[t,t1,t2] is the one-sided version of the tapering function described in CQG27:084020, arXiv:1003.2939 [gr-qc].";
TimeOfMaximum::usage="TimeOfMaximum[list_?ListQ] determines the time at which a time series reaches the absolute maximum of the modulus.";
GetAsymptoticMultiMode::usage="GetAsymptoticMode[h5file_,order_,modes_list,OptionsPattern[{\"Verbose\"\[Rule]False}]].";
EradUIB2017::usage="EradUIB2017[\[Eta]_,\[Chi]1_,\[Chi]2_] is fit for the energy radiated of the coalescence of 2 non-precessing quasicurcular BHs as in PRD:DY11536";
FinalSpinUIB2017::usage="FinalSpinUIB2017[\[Eta]_,\[Chi]1_,\[Chi]2_] is fit for the final spin of the coalescence of 2 non-precessing quasicurcular BHs as in PRD:DY11536";
IsFPNumberQ::usage="IsFPNumberQ[str_String] checks whether a string represents a real Number.";
ReadFloatingPointNumbers::usage="ReadFloatingPointNumbers[line_String, matchAnyFP_:True] converts string that contains a Fortran or C-style floating point number to a Mathematica Number.";
StringToNumber::usage="StringToNumber[x_] converts a String x to a number, allowing more general number formats than ToExpression.";
SXSMetaFilesToRules::usage="SXSMetaFilesToRules[filename_String] converts a SXS metadata.txt file to a list of rules.";
SXSParClassification::usage="SXSParClassification[sxsdir_,ClassStr_]. Given a list of SXS NR. data folders 'sxsdir', it returns all the cases that match a certain criterion 'ClassStr' (MassRatio range, Precessing or not, Initial Distance, Orbits Number)taking as reference the SXS metadata.txt files. If it is used iteratively, one could do different classifications ";
TakeColumn::usage="TakeColumn[list1_?ListQ,columns]. Given a list, take the columns specified by 'columns'. ";
ReSampleTD::usage="ReSampleTD[data_,sample_]. Resample your data with step-size 'sample'.";
AtomsList::usage="AtomsList[expr_]. Split your expression in elements";
InterpolationDomain::usage="InterpolationDomain[IntFunction]. Outputs the domain of your interpolated function IntFunction";
positionDuplicates::usage="positionDuplicates[list_]. Finds the duplicate elements in your list";
\[Omega]lmn::usage="\[Omega]lmn[l_,m_,n_,\[Eta]_,\[Chi]1_,\[Chi]2_,ModesData\[Rule]{}]. It computes the lmn frequency and damping time for a BBH with parameters {\[Eta],\[Chi]1,\[Chi]2}. ModesData can be feed with a datafile containing the modes data, which it speeds up the computation. ";
DetectConvergence::usage="DetectConvergence[dampsignal_,\[Tau]_,OptionsPattern[{Test\[Rule]FindMaximum}]]. Given a damped signal and a damping time, it estimaes a time t0 where the decay is order 1%";
OvertoneModel::usage="OvertoneModelV[overtones_,pars_,OptionsPattern[]]. It provides a RDown ansatz for the desired number 'overtones' with parameters {\[Eta],\[Chi]1,\[Chi]2}.";
OvertoneModelV2::usage="OvertoneModelV2[overtones_,pars_,ti_,OptionsPattern[]]. Upgraged version of OvertoneModelV2";
AICcRes::usage="AICRes[residuals_,coeff_]. It computes the AICc value from the residuals given the number of coefficients 'coeff' ";
BICRes::usage="BICRes[residuals_,coeff_]. It computes the BIC value from the residuals given the number of coefficients 'coeff' ";
PValueTest::usage="PValueTest[residuals_,coeff_]. It computes the P value from the residuals given the number of coefficients 'coeff' ";
CNMinimize::usage="CNMinimize[data_,ansatz_,coeffs_,{x_,xlist_},OptionsPattern[]]. Minimizator for complex data where coeffs are the free coefficients to adjust, x is coordinate variable and x_list its values";
Arg0DStructs::usage="Arg0DStructs[re_,im_] and Arg0DStructs[reim_] computes the phase angle of complex data. the two-argument version is significantly slower, the one-argument version is just a backwards compatibility interface to UnwrappedPhase.";
Abs0DStructs::usage="Abs0DStructs[re_,im_] computes the absolute value of complex data.";
UnwrapPhaseVector::usage="UnwrapPhaseVector[x_?VectorQ] unwraps the phase of a real vector representing a phase angle.";
UnwrapPhaseMatrix::usage="UnwrapPhaseMatrix[x_?MatrixQ] unwraps the phase of a real time series where the values represent a phase angle.";
UnwrappedPhase::usage="UnwrapPhase[x_?VectorQ] computes the unwrapped phase of a complex vector representing time series values; UnwrapPhase[x_?MatrixQ] computes the unwrapped phase of a complex time series.";
CombineColumns::usage="CombineColumns[list1_,list2_] combines columns from two lists in into a single list, e.g. CombineColumns[{1,2,3},{4,5,6}] will yield {{1,4},{2,5},{3,6}}.";
fH::usage="fH[f,M]";
fNU::usage="fNU[f,M]";
HessianH::usage="HessianH[f_,t_List?VectorQ]. f is the function you want to compute the Hessian and t the list of parameters."
CovarianceMatrix::usage="CovarianceMatrix[res_,vars_,pars_]. It computes the convarianc:e matrix from the residuals and the list of parameters."
CorrelationMatrix::usage="CovarianceMatrix[res_,vars_,pars]. It computes the correlation matrix from the residuals and the list of parameters."
LogLikelihoodDist::usage="Likelihood[data_,ansatz_,vars_,x_]. It computes the likelihood."
ComputeFitBands::usage="ComputeFitBands[data_,ansatz_,vars_,pars_,ConfidenceLevel\[Rule]val]. It computes the fit ConfidentBands at the val confidence level."
tPhys::usage="tPhys[t,M]";
tNR::usage="tNR[t,M]";
Begin["`Private`"]
Options[EasyMatchT]={PadRight->7, Listable->True, FreqFormat->{"Mf","Index","Index"},DisableChecks->False};
Options[EasyMatchT]=Options[EasyMatchT];
EasyMatchT[h1_?ListQ,h2_?ListQ,tMin_,tMax_,OptionsPattern[]]:=Module[{t,fHz,scale,h1red,h2red,tableS,norm1,norm2,myTable,myTable2,prod,psddat,fMinInd,fMaxInd,dt},
h1red=Select[h1,tMin<=#[[1]]<=tMax&][[All,2]];
h2red=Select[h2,tMin<=#[[1]]<=tMax&][[All,2]];
norm1=Total[(Abs@h1red)^2];
norm2=Total[(Abs@h2red)^2];
myTable=h1red Conjugate@h2red;
Re@Total[myTable]/Sqrt[norm1 norm2]
];
IsFPNumberQ[str_String]:=StringMatchQ[StringTrim@str,RegularExpression["[+-]?\\d+\\.?\\d*$"]] ||
StringMatchQ[StringTrim@str,RegularExpression["[-+]?\\d*\\.\\d+([eE][-+]?\\d+)?"]]
TimeOfMaximum[list_?ListQ,OptionsPattern[{"MinTime"->""}]]:=Module[{reducedList,vals,maxval,maxpos,mintime},
mintime=OptionValue["MinTime"];
If[NumberQ@mintime,
reducedList=Select[list,#[[1]]>= mintime&];,
reducedList=list;
];
vals=Abs@TakeColumn[reducedList,2];
maxval=Max@vals;
If[Length@Union@vals>1,
maxpos=First@First@Position[vals,maxval];
First[reducedList[[maxpos]]],
"NOMAX"]
];
PlanckTaper[t_?VectorQ,t1_,t2_]:=Module[{tmid=Select[t,t1<#<t2&]},Join[
ConstantArray[0,Length@Select[t,#<=t1&]],
1/(Exp[(t2-t1)/(tmid-t1)+(t2-t1)/(tmid-t2)]+1),
ConstantArray[1,Length@Select[t,#>=t2&]]
]];
PlanckTaper[t_?NumberQ,t1_,t2_]:=
Piecewise[{{0,t<=t1},
{1/(Exp[(t2-t1)/(t-t1)+(t2-t1)/(t-t2)]+1),t1<t<t2},
{1,t>= t2}
}];
StringToNumber[y_]:=Module[{temp,x},
x = ToString@y;
temp = StringReplace[x," "->""];
If[StringMatchQ[temp,NumberString],
ToExpression@temp,
temp = ReadFloatingPointNumbers[temp];
If[ListQ@temp && Length@temp > 0,
First@temp,
Indeterminate
]]
]
ReadFloatingPointNumbers[lineIn_String,matchAnyFP_:True]:=Module[{line,fpregex,simpleHackRule,temp},
temp=Map[ToString,Range[0,9]];
simpleHackRule=temp/.xx_/;StringQ@xx:>xx->xx<>".0";line=lineIn/.simpleHackRule;
simpleHackRule=temp/.xx_/;StringQ@xx:>xx<>"."->xx<>".0";
line=line/.simpleHackRule;
If[matchAnyFP,fpregex="[-+]?\\d*\\.\\d+([eE][-+]?\\d+)?";(*matches any floating point number*),fpregex="[-+]?\\d*\\.?\\d+[eE][-+]?\\d+";(*matches only floating point numbers in exponential format*)];
Map[Read[StringToStream[#],Number]&,StringCases[line,RegularExpression[fpregex]]]
];
ReSampleTD[data_,sample_]:=Module[{dataInt,domain,times},
dataInt=Interpolation@data;
domain=InterpolationDomain@dataInt;
times=Range[domain[[1,1]],domain[[1,2]],sample];
Transpose[{times,dataInt@times}]
]
(* ::Code::Initialization:: *)
AtomsList[expr_]:=Union@Select[Level[expr,{0,Infinity}],AtomQ];
InterpolationDomain[fun_]:=Module[{min,max},fun[[1]]];
TakeColumn[list1_?ListQ,list2_?ListQ]:=Map[Part[#,list2]&,list1];
TakeColumn[list1_?ListQ,n_?IntegerQ]:=(list1//Transpose)[[n]];
positionDuplicates[listaux_]:=GatherBy[Range@Length[listaux],listaux[[#]]&];
Arg0DStructs[reim_?MatrixQ]:=UnwrappedPhase[reim]
Arg0DStructs[reim_?MatrixQ,thresh_?NumberQ]:=Arg0DStructs[reim]
Arg0DStructs[re_?MatrixQ,im_?MatrixQ]:=Module[{revals,imvals,reim,reimvals,times,i,j,diff,
tmp,twopi,zero,test,thresh=N@\[Pi]},
times=TakeColumn[re,1];
revals=TakeColumn[re,2];
imvals=TakeColumn[im,2];
reim = revals + I imvals;
reimvals=Arg@TakeColumn[reim,2];
tmp=N[2 \[Pi]];
twopi=ConstantArray[tmp,Length@times];
zero=ConstantArray[0.,Length@times];
Do[
test=(reimvals[[i+1]]-reimvals[[i]]);
If[(test>thresh),
reimvals=reimvals-Join[zero[[;;i]],twopi[[i+1;;]]];
];
If[(test<-thresh),
reimvals=reimvals+Join[zero[[;;i]],twopi[[i+1;;]]];
];,
{i,1,Length@reimvals-1}
];
CombineColumns[times,reimvals]
];
Abs0DStructs[re_,im_]:=Module[{revals,imvals,reimvals,times,dt,result},
times=Map[First,re];
revals=Map[Last,re];
imvals=Map[Last,im];
reimvals=revals + I imvals;
reimvals=Abs@reimvals;
result=CombineColumns[times,reimvals]
];
Abs0DStructs[reim_]:=Module[{reimvals,times,dt,result},
times=Map[First,reim];
reimvals=Map[Last,reim];
reimvals=Abs@reimvals;
result=CombineColumns[times,reimvals]
];
UnwrapPhaseVector=Compile[{{data,_Real,1}},
(* Juergen Tischer,
http://forums.wolfram.com/mathgroup/archive/1998/May/msg00105.html*)
FoldList[Round[(#1-#2)/(2Pi)] 2Pi+#2&,First[data],Rest[data]]
];
UnwrappedPhase[data_?VectorQ]:=UnwrapPhaseVector[Arg@data]
UnwrappedPhase[data_?MatrixQ]:=CombineColumns[TakeColumn[data,1],UnwrapPhaseVector[Arg@TakeColumn[data,2]]]
CombineColumns[list1_,list2_]:=If[Length@list1==Length@list2, {list1,list2}//Transpose,0{list1,list1}//Transpose]
CombineColumns[list_]:=If[Length@Union@(Length /@ list) == 1, list//Transpose,{}]
fH[f_,M_]:=Module[{c,G,MS},
c=2.99792458 10^8;G=6.67259*10^-11;MS=1.9885*10^30;
(c)^3/((M)MS(G)) f];(* in Solar mass units *)
fNU[f_,M_]:=Module[{c,G,MS},c=2.99792458 10^8;G=6.67259*10^-11;MS=1.9885*10^30;((M)MS(G))/(c)^3 f]
tPhys[t_,M_]:=Module[{c,G,MS},c=2.99792458 10^8;G=6.67259*10^-11;MS=1.9885*10^30;
((M)MS(G))/(c)^3 t]
tNR[t_,M_]:=Module[{c,G,MS},c=2.99792458 10^8;G=6.67259*10^-11;MS=1.9885*10^30;
1/(((M)MS(G))/(c)^3) t]
Options[CNMinimize]=Join[Options@NMinimize,{"Method"->Automatic}];
CNMinimize[data_,ansatz_,vars_,{x_,xlist_},OptionsPattern[]]:=Module[{likelihood,ansatzlist,method},
method=OptionValue["Method"];
ansatzlist=ansatz/.x->xlist;
likelihood=Total[(Re[ansatzlist-data]^2+Im[ansatzlist-data]^2)];
NMinimize[likelihood,vars,Method->method]
];
GetAsymptoticMultiMode[h5file_,order_,list_,OptionsPattern[{"Verbose"->False,"TiReIm"->True,"ZeroAlign"->False,"ReSample"->True,"Sampling"->0.5}]]:=Module[{datasets,targetstring,pos,VPrint,
l,m,positions,timreim,tt,re,im,import,zeroalign,tmax,resample,wave,sampling},
VPrint[expr___]:=If[OptionValue["Verbose"],Print[expr]];
timreim=OptionValue["TiReIm"];
zeroalign=OptionValue["ZeroAlign"];
resample=OptionValue["ReSample"];
sampling=OptionValue["Sampling"];
datasets=Import[h5file,"Datasets"];
positions=Table[
If[NumberQ[order],l=list[[i,1]]; m=list[[i,2]];
targetstring="/Extrapolated_N"<>ToString[order]<>".dir/Y_l"<>ToString[l]<>"_m"<>ToString[m]<>".dat";
,
If[order==="outermost",targetstring="/OutermostExtraction.dir/Y_l"<>ToString[l]<>"_m"<>ToString[m]<>".dat";]
];
VPrint["Targeting entry: ",targetstring];
pos=Flatten[Position[datasets,targetstring]];Switch[Length[pos],
0,Print["Data not found"];Return[Null];,
1,pos=pos[[1]];,
_,Print["Error, two entries with this name"];
Return[Null];];
VPrint["Annotations: ",Import[h5file,{"Annotations",pos}]];
VPrint["Data Format: ",Import[h5file,{"DataFormat",pos}]];
VPrint["Data dimensions: ",Import[h5file,{"Dimensions",pos}]];
VPrint["Data encoding: ",Import[h5file,{"DataEncoding",pos}]];
pos,{i,1,Length[list]}];
If[Length[list]==1,
import={Import[h5file,{"Data",positions}]};,
import=Import[h5file,{"Data",positions}];
];
If[timreim,wave=(#1/. {tt_,re_,im_}->{tt,re+I im}&)/@import;,wave=import];
If[resample,wave=ReSampleTD[#,sampling]&/@wave;];
If[zeroalign,tmax=TimeOfMaximum[wave[[1,-Round[(Length@wave[[1]])/2];;-1]]];Table[wave[[i,All,1]]=wave[[i,All,1]]-tmax,{i,Length@import}]];
wave
]
SXSMetaFilesToRules[filePath_]:=Module[{filepath,fileList,meta1,pos1,meta2,meta3,meta4,value,var,list},
If[ListQ@filePath,filepath=filePath[[1]],filepath=filePath];
If[Not@FileExistsQ@filepath,Print["File not found"];Return[]];
(*Reading the file*)
fileList=ReadList[filepath,String];
(*Delete comments*)
meta1=Delete[fileList,Position[StringMatchQ[fileList,"#*"],True]];
(*Fix eccentricity*)
meta1[[Flatten@Position[StringMatchQ[meta1,"*<*e-*"],True]]]=StringReplace[meta1[[Flatten@Position[StringMatchQ[meta1,"*<*e-*"],True]]],"<"->""];
(*Find = *)
meta2=meta1[[Flatten@Position[StringMatchQ[meta1,"*=*"],True]]];
meta3=StringSplit[meta2,"="];
(* Select non-empty fields*)
meta3=Select[meta3,Length@#>1&];
(*Delete spaces*)
meta4=Transpose[{StringReplace[meta3[[All,1]]," "->""],StringReplace[meta3[[All,2]]," "->""]}];
var=meta4[[All,1]];
value=meta4[[All,2]];
value=StringSplit[#,","]&/@value;
(*Delete empty elements*)
(*value=Select[value, UnsameQ[#, {}] &];*)
list={};
Do[
If[Length@value[[i]]==0,
{}
,
If[IsFPNumberQ@value[[i,1]],
value[[i]]=(StringToNumber@#)&/@value[[i]];
]; list=AppendTo[list,{var[[i]]->value[[i]]}];
];
,
{i,1,Length@value}];
Flatten@list
]
SXSParClassification[sxsdir_?ListQ,ClassStr_?ListQ,OptionsPattern[{"\[Epsilon]"->0.001,"HighSpin"->0.8,"UnRepeated"->False,"Verbose"->False,"Mass1-Str"->"initial-mass1","Mass2-Str"->"initial-mass2"}]]:=Module[{metafiles,metadata,
orbitStr="number-of-orbits",dStr="initial-separation",mass1Str,mass2Str,spin1Str="initial-dimensionless-spin1",
spin2Str="initial-dimensionless-spin2",eccStr="relaxed-eccentricity",spin1Dim,spin2Dim,mass1,mass2,massratio,eccentricity,dist,orbit,select,pos,condition,A1,A2,precvalue,precvalueNorm,\[Epsilon],
spin1Norm,spin2Norm,highspin,spintest,\[Chi]eff,sxsdirout,spinz,spinzDiff,auxDist,posdup,posdupDist,posdistecc,unrepeated,verbose,sxsdiroutaux,precvalue1,precvalue2,precvalueNorm1,precvalueNorm2},
Print["Classification Input Variables. Examples: {{'MassRatio', '0.99<#<1.1'}},{{'Distance', '#>16'}},{{'Orbits', '#>25'}},{{'Precessing'}},
{{'Non-Precessing'}},{{'High-Spin'}},{{'\[Chi]eff','#>0.6'}},{{'\[Chi]1','#>0.6'}},{{'\[Chi]2','#>0.6'}},{{'Unequal'}}"];
Print["Take care! Some of the sxs file names are not consistent with the metadata files"];
Print["The spin definition is consitent with 'initial-spin' values and not relaxed ones"];
Print["The mass definition is consitent with 'initial-mass' values and not relaxed ones"];
mass1Str=OptionValue["Mass1-Str"];
mass2Str=OptionValue["Mass2-Str"];
(* Kill the loop if the root directory is wrong *)
If[And@@(Not/@DirectoryQ/@sxsdir),Print[Style["Directory not found",Red]];Return[{}]];
sxsdirout=sxsdir;
\[Epsilon]=OptionValue["\[Epsilon]"];
highspin=OptionValue["HighSpin"];
unrepeated=OptionValue["UnRepeated"];
verbose=OptionValue["Verbose"];
metafiles=Flatten[FileNames["metadata.txt",#,4]&/@sxsdirout,1];
metadata=SXSMetaFilesToRules[#]&/@metafiles;
mass1=Flatten@((mass1Str/.#)&/@metadata);
mass2=Flatten@((mass2Str/.#)&/@metadata);
massratio=mass1/mass2;
dist=Flatten@((dStr/.#)&/@metadata);
orbit=Flatten@((orbitStr/.#)&/@metadata);
spin1Dim=((spin1Str/.#)&/@metadata);
spin2Dim=((spin2Str/.#)&/@metadata);
A1=(1+3 massratio/(4.));
A2=(1+3 /(4.*massratio) );
\[Chi]eff=massratio/(1.+massratio)*spin1Dim +1./(1+massratio)*spin2Dim;
eccentricity=Flatten[(eccStr/.#)&/@metadata];
spinz=Chop/@Transpose[{TakeColumn[spin1Dim,3],TakeColumn[spin2Dim,3]}];
spinzDiff=Abs[(#[[2]]-#[[1]])&/@spinz];
condition=ClassStr[[All,1]];
select=Table[If[Length@ClassStr[[i]]==2,ToExpression@(ClassStr[[i,2]]),"Null"],{i,1,Length@ClassStr}];
pos=Table[i,{i,Length@sxsdirout}];
Do[
Which[condition[[i]]=="MassRatio",
If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
pos=Flatten@Position[massratio,_?(Evaluate[select[[i]]]&)];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];
A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];sxsdirout=sxsdirout[[pos]];,
condition[[i]]== "Distance",
If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
pos=Flatten@Position[sxsdirout=sxsdirout[[pos]];,_?(Evaluate[select[[i]]]&)];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];
A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];
sxsdirout=sxsdirout[[pos]];,
condition[[i]]== "Orbits",
If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
pos=Flatten@Position[orbit,_?(Evaluate[select[[i]]]&)];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];
A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];
sxsdirout=sxsdirout[[pos]];,
condition[[i]]== "Non-Precessing",
(*precvalue=(#)\[Cross]{0,0,1}&/@((spin1Dim)A1+(spin2Dim)A2);*)
precvalue1=(#)\[Cross]{0,0,1}&/@((spin1Dim)A1);
precvalue2=(#)\[Cross]{0,0,1}&/@((spin2Dim)A2);
precvalueNorm1=Norm[#]&/@precvalue1;
precvalueNorm2=Norm[#]&/@precvalue2;
precvalueNorm=precvalueNorm1^2+precvalueNorm2^2;
pos=Flatten@Position[precvalueNorm,_?(#<\[Epsilon] &)];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];
A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];
sxsdirout=sxsdirout[[pos]];,
condition[[i]]== "Unequal",
pos=Flatten@Position[spinzDiff,_?(#>\[Epsilon] &)];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];
sxsdirout=sxsdirout[[pos]];,
condition[[i]]=="Precessing",
precvalue1=(#)\[Cross]{0,0,1}&/@((spin1Dim)A1);
precvalue2=(#)\[Cross]{0,0,1}&/@((spin2Dim)A2);
precvalueNorm1=Norm[#]&/@precvalue1;
precvalueNorm2=Norm[#]&/@precvalue2;
precvalueNorm=precvalueNorm1^2+precvalueNorm2^2;
(*precvalue=(#)\[Cross]{0,0,1}&/@((spin1Dim)A1+(spin2Dim)A2);*)
pos=Flatten@Position[precvalueNorm,_?(#>\[Epsilon] &)];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];
A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];
sxsdirout=sxsdirout[[pos]];,
condition[[i]]== "High-Spin",
spin1Norm=Norm[#]&/@spin1Dim;
spin2Norm=Norm[#]&/@spin2Dim;
spintest=Table[If[Abs@spin1Norm[[i]]>=highspin || Abs@spin2Norm[[i]]>=highspin,True,False],{i,1,Length@spin1Norm}];
pos=Flatten@Position[spintest,True];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];
A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];
sxsdirout=sxsdirout[[pos]];,
condition[[i]]== "\[Chi]eff",
If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
pos=Flatten@Position[Norm[#]&/@\[Chi]eff,_?(Evaluate[select[[i]]] &)];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];
A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];
sxsdirout=sxsdirout[[pos]];,
condition[[i]]== "\[Chi]1",
If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
pos=Flatten@Position[Norm[#]&/@spin1Dim,_?(Evaluate[select[[i]]] &)];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];
A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];
sxsdirout=sxsdirout[[pos]];,
condition[[i]]== "\[Chi]2",
If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
pos=Flatten@Position[Norm[#]&/@spin2Dim,_?(Evaluate[select[[i]]] &)];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];
A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];
sxsdirout=sxsdirout[[pos]];,
condition[[i]]== "relaxed-eccentricity",
If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
pos=Flatten@Position[eccentricity,_?(Evaluate[select[[i]]] &)];
massratio=massratio[[pos]];
dist=dist[[pos]];
orbit=orbit[[pos]];
spin1Dim=spin1Dim[[pos]];
spin2Dim=spin2Dim[[pos]];
A1=A1[[pos]];
A2=A2[[pos]];
spinzDiff=spinzDiff[[pos]];
eccentricity=eccentricity[[pos]];
\[Chi]eff=\[Chi]eff[[pos]];
sxsdirout=sxsdirout[[pos]];,
True,
Print[Style["Wrong input",Bold,Red,16]];
Break[];
];
,{i,1,Length@ClassStr}];
If[unrepeated,
Print["Taking among the repeated cases only those with lower eccentricity and larger D (just in case ei=ej)"];
(* Selecting Case with lower e *)
auxDist=Transpose[{Round[#&/@massratio,0.1],Round[Chop[#,10^(-2)],0.01]&/@spin1Dim,Round[Chop[#,10^(-2)],0.01]&/@spin2Dim,dist,eccentricity}];
posdup=positionDuplicates@auxDist[[All,1;;3]];
posdistecc=Flatten[#,1]&/@Table[Position[auxDist[[posdup[[i]],5]],Min@auxDist[[posdup[[i]],5]]],{i,1,Length@posdup}];
posdistecc=Flatten@Table[posdup[[i,posdistecc[[i]]]],{i,1,Length@posdup}];sxsdirout=sxsdirout[[posdistecc]];
(* Selecting Case with larger D *)
auxDist=auxDist[[posdistecc]];
posdup=positionDuplicates@auxDist[[All,1;;3]];
posdistecc=Flatten[#,1]&/@Table[Position[auxDist[[posdup[[i]],4]],Max@auxDist[[posdup[[i]],4]]],{i,1,Length@posdup}];
posdistecc=Flatten@Table[posdup[[i,posdistecc[[i]]]],{i,1,Length@posdup}];
auxDist=auxDist[[posdistecc]];
sxsdirout=sxsdirout[[posdistecc]];
sxsdiroutaux=SortBy[Table[Join[{sxsdirout[[i]]},auxDist[[i]]],{i,1,Length@sxsdirout}],#[[2]]&];
If[verbose,
Print[Prepend[Table[ToString@#&/@sxsdiroutaux[[i]],{i,1,Length@sxsdiroutaux}],{"Case","q","\[Chi]1","\[Chi]2","D","e"}]//TableForm]];
sxsdiroutaux[[All,1]],
auxDist=Transpose[{Round[#&/@massratio,0.1],Round[Chop[#,10^(-2)],0.01]&/@spin1Dim,Round[Chop[#,10^(-2)],0.01]&/@spin2Dim,dist,eccentricity}];
(*sxsdiroutaux=Table[Join[{sxsdirout[[i]]},auxDist[[i]]],{i,1,Length@sxsdirout}];*)
sxsdiroutaux=SortBy[Table[Join[{sxsdirout[[i]]},auxDist[[i]]],{i,1,Length@sxsdirout}],#[[2]]&];
If[verbose,
Print[Prepend[Table[ToString@#&/@sxsdiroutaux[[i]],{i,1,Length@sxsdiroutaux}],{"Case","q","\[Chi]1","\[Chi]2","D","e"}]//TableForm]];
sxsdiroutaux[[All,1]]]
]
Options[\[Omega]lmn]={"ModesData"->{},"ModesFile"->""};
\[Omega]lmn[l_,m_,n_,\[Eta]_,\[Chi]1_,\[Chi]2_,OptionsPattern[]]:=Module[{file,data,modesdata,modesfile,mstr,sp,fmass,f\[Chi],\[Chi]\[Omega]\[Tau],\[Omega],\[Tau],llim,ulim},
modesdata=OptionValue["ModesData"];
modesfile=OptionValue["ModesFile"];
If[Length[modesdata]==0,
If[n<=7,
If[m>=0,mstr="m",mstr="mm"];
file=modesfile<>"l"<>ToString[l]<>"/n"<>ToString[n+1]<>"l"<>ToString[l]<>mstr<>ToString[Abs[m]]<>".dat";,
file=modesfile<>"l"<>ToString[l]<>"/n7l"<>ToString[l]<>"m"<>ToString[Abs[m]]<>".dat"];
If[Not@FileExistsQ[file],Print["File not found"];Return[];];
data=Import[file];,
data=modesdata;
];
\[Chi]\[Omega]\[Tau]=TakeColumn[data,{1,2,3}];
fmass=1-EradUIB2017[\[Eta],\[Chi]1,\[Chi]2];
f\[Chi]=FinalSpinUIB2017[\[Eta],\[Chi]1,\[Chi]2];
ulim=Flatten[Select[\[Chi]\[Omega]\[Tau],#[[1]]>=f\[Chi]&,1]];
llim=Flatten[Select[\[Chi]\[Omega]\[Tau],#[[1]]-ulim[[1]]<0&][[-1]]];
\[Omega]=Mean[Join[{ulim},{llim}]][[2]]/fmass;
\[Tau]=Mean[(Abs[1/Join[{ulim},{llim}]]*fmass)][[-1]];
{\[Omega],\[Tau]}
]
EradUIB2017[\[Eta]_,\[Chi]1_,\[Chi]2_]:=Module[{m1,m2,S},
m1=1/2 (1+Sqrt[1-4 \[Eta]]);
m2=1/2 (1-Sqrt[1-4 \[Eta]]);
S= (m1^2 \[Chi]1 + m2^2 \[Chi]2)/(m1^2 + m2^2);
(((1-(2 Sqrt[2])/3) \[Eta]+0.5609904135313374` \[Eta]^2-0.84667563764404` \[Eta]^3+3.145145224278187` \[Eta]^4) (1+S^3 (-0.6320191645391563`+ 4.952698546796005` \[Eta]-10.023747993978121` \[Eta]^2)+S^2 (-0.17762802148331427`+ 2.176667900182948` \[Eta]^2)+S (-0.13084389181783257`- 1.1387311580238488` \[Eta]+5.49074464410971` \[Eta]^2)))/(1+S (-0.9919475346968611`+ 0.367620218664352` \[Eta]+4.274567337924067` \[Eta]^2))
-0.01978238971523653` S (1-4.91667749015812` \[Eta]) (1-4 \[Eta])^0.5` \[Eta] (\[Chi]1-\[Chi]2)-0.09803730445895877` (1-4 \[Eta])^0.5` (1-3.2283713377939134` \[Eta]) \[Eta]^2 (\[Chi]1-\[Chi]2)+0.01118530335431078` \[Eta]^3 (\[Chi]1-\[Chi]2)^2]
FinalSpinUIB2017[\[Eta]_,\[Chi]1_,\[Chi]2_]:=Module[{m1,m2,S},
m1=1/2 (1+Sqrt[1-4 \[Eta]]);
m2=1/2 (1-Sqrt[1-4 \[Eta]]);
S= (m1^2 \[Chi]1 + m2^2 \[Chi]2)/(m1^2 + m2^2);
(m1^2 \[Chi]1 + m2^2 \[Chi]2) +(2 Sqrt[3] \[Eta]+20.0830030082033` \[Eta]^2-12.333573402277912` \[Eta]^3)/(1+7.2388440419467335` \[Eta])+
0.3223660562764661` (1-4 \[Eta])^0.5` \[Eta]^2 (1+9.332575956437443` \[Eta]) (\[Chi]1-\[Chi]2)-0.059808322561702126` \[Eta]^3 (\[Chi]1-\[Chi]2)^2+
(2.3170397514509933` (1-4 \[Eta])^0.5` (1-3.2624649875884852` \[Eta]) \[Eta]^3 (\[Chi]1-\[Chi]2) (1/4 (1+Sqrt[1-4 \[Eta]])^2 \[Chi]1+1/4 (1-Sqrt[1-4 \[Eta]])^2 \[Chi]2))/(1/4 (1-Sqrt[1-4 \[Eta]])^2+1/4 (1+Sqrt[1-4 \[Eta]])^2)+
(((0.` -0.8561951310209386` \[Eta]-0.09939065676370885` \[Eta]^2+1.668810429851045` \[Eta]^3) (1/4 (1+Sqrt[1-4 \[Eta]])^2 \[Chi]1+1/4 (1-Sqrt[1-4 \[Eta]])^2 \[Chi]2))/(1/4 (1-Sqrt[1-4 \[Eta]])^2+1/4 (1+Sqrt[1-4 \[Eta]])^2)+((0.` +0.5881660363307388` \[Eta]-2.149269067519131` \[Eta]^2+3.4768263932898678` \[Eta]^3) (1/4 (1+Sqrt[1-4 \[Eta]])^2 \[Chi]1+1/4 (1-Sqrt[1-4 \[Eta]])^2 \[Chi]2)^2)/(1/4 (1-Sqrt[1-4 \[Eta]])^2+1/4 (1+Sqrt[1-4 \[Eta]])^2)^2+((0.` +0.142443244743048` \[Eta]-0.9598353840147513` \[Eta]^2+1.9595643107593743` \[Eta]^3) (1/4 (1+Sqrt[1-4 \[Eta]])^2 \[Chi]1+1/4 (1-Sqrt[1-4 \[Eta]])^2 \[Chi]2)^3)/(1/4 (1-Sqrt[1-4 \[Eta]])^2+1/4 (1+Sqrt[1-4 \[Eta]])^2)^3)/(1+((-0.9142232693081653`+ 2.3191363426522633` \[Eta]-9.710576749140989` \[Eta]^3) (1/4 (1+Sqrt[1-4 \[Eta]])^2 \[Chi]1+1/4 (1-Sqrt[1-4 \[Eta]])^2 \[Chi]2))/(1/4 (1-Sqrt[1-4 \[Eta]])^2+1/4 (1+Sqrt[1-4 \[Eta]])^2))]
Eth[n_,f_] := - Sin[t]^n (D[f/Sin[t]^n,t]+I D[f/Sin[t]^n,p]/Sin[t])
Yp1[l_,m_]:= FullSimplify[Sqrt[ (l-1)!/(l+1)!]Eth[0,SphericalHarmonicY[l,m,t,p]]]
Yp2[l_,m_]:= FullSimplify[Sqrt[ (l-2)!/(l+2)!]Eth[1,Eth[0,SphericalHarmonicY[l,m,t,p]]]]
Ethbar[n_,f_] := - Sin[t]^(-n) (D[f Sin[t]^n,t]-I D[f Sin[t]^n,p]/Sin[t])
Ym1[l_,m_]:= FullSimplify[Sqrt[ (l-1)!/(l+1)!](-1)^(-1) Ethbar[0,SphericalHarmonicY[l,m,t,p]]]
Ym2[l_,m_]:= FullSimplify[Sqrt[ (l-2)!/(l+2)!](-1)^(-2)Ethbar[-1,Ethbar[0,SphericalHarmonicY[l,m,t,p]]]]
(* same source, direct formula using Wigner d-functions *)
wd[n_,l_,m_]:=Sum[
(-1)^i Sqrt[(l+m)! (l-m)! (l+n)! (l-n)!] / ((l+m-i)!(l-n-i)!i!(i+n-m)!)
Cos[t/2]^(2l+m-n-2i)Sin[t/2]^(2i+n-m),
{i,Max[0,m-n],Min[l+m,l-n]}
]
Ydirect[n_,l_,m_]:=FullSimplify[(-1)^n Sqrt[(2l+1)/(4Pi)] wd[-n,l,m]E^(I m p)]
Options[DetectConvergence]={"Test"->FindMaximum};
DetectConvergence[res_,\[Tau]_,OptionsPattern[]]:=Module[{intres,mins,tmax,t,fits,a,tol,tt,yy,test},
test=OptionValue["Test"];
tmax= 4 * \[Tau];
intres=Interpolation@res;
mins=Quiet[Table[test[intres@t,{t,i}],{i,res[[1,1]]+5,res[[1,1]]+5 +tmax,10}]/.{yy_,tt_}->{t/.tt,yy}];
mins=Select[mins,res[[1,1]]<#[[1]]<res[[-1,1]]&];
SortBy[mins,First]
]
Options[OvertoneModel]=Join[Options[\[Omega]lmn],{"FixExtra"->True,"Fit\[Alpha]"->{},"Fit\[Tau]"->{},"Mode"->{2,2},"Vary\[Omega]"->False,"\[Omega]val"->{-0.05,0.05}}];
OvertoneModel[overtones_,pars_,OptionsPattern[]]:=Block[{ansatz,fixetra,fit\[Alpha],fit\[Tau],l,m,mode,modesdata,modesfile,modto0,var,vary\[Omega],Global`t,ti,\[Tau]s,\[Alpha],\[Beta],\[Omega]s,\[Eta],\[Chi]1,\[Chi]2,\[Omega]val},
fixetra=OptionValue["FixExtra"];
fit\[Alpha]=OptionValue["Fit\[Alpha]"];
fit\[Tau]=OptionValue["Fit\[Tau]"];
mode=OptionValue["Mode"];
vary\[Omega]=OptionValue["Vary\[Omega]"];
\[Omega]val=OptionValue["\[Omega]val"];
modesdata=OptionValue["ModesData"];
modesfile=OptionValue["ModesFile"];
l=mode[[1]];
m=mode[[2]];
{\[Eta],\[Chi]1,\[Chi]2}=pars[[1;;3]];
(* Freqs. & damping times from Vitor *)
If[Length@modesdata==0,
\[Omega]s=Table[\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesFile"->modesfile][[1]],{n,0,overtones}];
var=Range[0,7];
If[vary\[Omega],\[Omega]s=Table[If[n>0,var[[n+1]]=RandomReal[{1+\[Omega]val[[1]],1+\[Omega]val[[2]]}],var[[1]]=1];var[[n+1]]*\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesFile"->modesfile][[1]],{n,0,overtones}]];
\[Tau]s=Table[\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesFile"->modesfile][[2]],{n,0,overtones}];
,
\[Omega]s=Table[\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesData"->modesdata[[n+1]]][[1]],{n,0,overtones}];
var=Range[0,7];
If[vary\[Omega],\[Omega]s=Table[If[n>0,var[[n+1]]=RandomReal[{1+\[Omega]val[[1]],1+\[Omega]val[[2]]}],var[[1]]=1];var[[n+1]]*\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesData"->modesdata[[n+1]]][[1]],{n,0,overtones}]];
\[Tau]s=Table[\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesData"->modesdata[[n+1]]][[2]],{n,0,overtones}];
];
ansatz=Sum[ToExpression["x"<>ToString[n]] Exp[-(Global`t-ti)/(\[Tau]s[[n+1]](1+ToExpression["\[Beta]"<>ToString[n]]))] (Cos[\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]]) (Global`t)]+I Sin[ \[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]]) (t)] ),{n,0,overtones}];modto0=Complement[Table[i,{i,0,overtones}],fit\[Alpha]];
modto0=Table[ToExpression["\[Alpha]"<>ToString[modto0[[i]]]],{i,Length@modto0}];
ansatz=ansatz/.(Table[modto0[[i]]->0,{i,Length@modto0}]);
modto0=Complement[Table[i,{i,0,overtones}],fit\[Tau]];
modto0=Table[ToExpression["\[Beta]"<>ToString[modto0[[i]]]],{i,Length@modto0}];
ansatz=ansatz/.(Table[modto0[[i]]->0,{i,Length@modto0}])
]
Options[OvertoneModelV2]=Join[Options[\[Omega]lmn],{"FixExtra"->True,"Fit\[Alpha]"->{},"Fit\[Tau]"->{},"Mode"->{2,2},"Vary\[Omega]"->False,"Mixing"->{False,{2,2}},"ReIm"->False,"AmpPhase"->False,"AmpPhaseShift"->False,"QualityFactor"->False,"ParameterValues"->{},"\[Omega]val"->{-0.05,0.05},"Export_\[Omega]val"->False}];
OvertoneModelV2[overtones_,pars_,ti_,OptionsPattern[]]:=Block[{ansatz,ampansatz,amphase,amphaseshift,ex\[Omega]val,fixetra,fit\[Alpha],fit\[Tau],im,imm,l,m,lm,modesdata,modesfile,mm,mixing,mode,modto0,parvals,phaseansatz,qfact,qfactm,qualfactorQ,re,reim,rem,Global`t,var,vary\[Omega],x,y,\[Tau]s,\[Tau]sm,\[Alpha],\[Beta],\[Omega]m,\[Omega]s,\[Omega]m\[Tau]sm,\[Eta],\[Chi]1,\[Chi]2,mixmode,\[Omega]val},
fixetra=OptionValue["FixExtra"];
fit\[Alpha]=OptionValue["Fit\[Alpha]"];
fit\[Tau]=OptionValue["Fit\[Tau]"];
mode=OptionValue["Mode"];
vary\[Omega]=OptionValue["Vary\[Omega]"];
mixing=OptionValue["Mixing"];
reim=OptionValue["ReIm"];
amphase=OptionValue["AmpPhase"];
amphaseshift=OptionValue["AmpPhaseShift"];
parvals=OptionValue["ParameterValues"];
\[Omega]val=OptionValue["\[Omega]val"];
ex\[Omega]val=OptionValue["Export_\[Omega]val"];
modesdata=OptionValue["ModesData"];
modesfile=OptionValue["ModesFile"];
qualfactorQ=OptionValue["QualityFactor"];
{lm,mm}=mixing[[2]];
l=mode[[1]];
m=mode[[2]];
{\[Eta],\[Chi]1,\[Chi]2}=pars[[1;;3]];
(* Freqs. & damping times from Vitor *)
If[Length@modesdata==0,
\[Omega]s=Table[\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesFile"->modesfile][[1]],{n,0,overtones}];
var=Range[0,7];
If[vary\[Omega],\[Omega]s=Table[If[n>0,var[[n+1]]=RandomReal[{1+\[Omega]val[[1]],1+\[Omega]val[[2]]}],var[[1]]=1];var[[n+1]]*\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesFile"->modesfile][[1]],{n,0,overtones}]];
\[Tau]s=Table[\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesFile"->modesfile][[2]],{n,0,overtones}];
qfact=\[Omega]s*\[Tau]s;
,
\[Omega]s=Table[\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesData"->modesdata[[n+1]]][[1]],{n,0,overtones}];
var=Range[0,7];
If[vary\[Omega],\[Omega]s=Table[If[n>0,var[[n+1]]=RandomReal[{1+\[Omega]val[[1]],1+\[Omega]val[[2]]}],var[[1]]=1];var[[n+1]]*\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesData"->modesdata[[n+1]]][[1]],{n,0,overtones}]];
\[Tau]s=Table[\[Omega]lmn[l,m,n,\[Eta],\[Chi]1,\[Chi]2,"ModesData"->modesdata[[n+1]]][[2]],{n,0,overtones}];
qfact=\[Omega]s*\[Tau]s;
];
If[Not@amphase,
If[qualfactorQ,
ansatz=Sum[If[Not@reim,x=ToExpression["x"<>ToString[n]],x=ToExpression[("x"<>ToString[n])]+I ToExpression[("y"<>ToString[n])]]; x Exp[-(Global`t-ti)(\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]])/(qfact[[n+1]](1+ToExpression["\[Beta]"<>ToString[n]])))] (Cos[\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]]) (Global`t)]+I Sin[ \[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]]) (Global`t)] ),{n,0,overtones}];
If[mixing[[1]],
\[Omega]m\[Tau]sm=Table[\[Omega]lmn[lm,mm,n,\[Eta],\[Chi]1,\[Chi]2,"ModesFile"->modesfile],{n,0,overtones}];
\[Omega]m=\[Omega]m\[Tau]sm[[All,1]];
\[Tau]sm=\[Omega]m\[Tau]sm[[All,2]];
qfactm=\[Omega]m*\[Tau]sm;
ansatz=ansatz+Sum[If[Not@reim,x=ToExpression["xm"<>ToString[n]],x=ToExpression["xm"<>ToString[n]]+I ToExpression["ym"<>ToString[n]]];x Exp[-(Global`t-ti)(\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]])/(qfactm[[n+1]](1+ToExpression["\[Beta]"<>ToString[n]])))] (Cos[\[Omega]m[[n+1]] (Global`t)]+I Sin[ \[Omega]m[[n+1]] (Global`t)] ),{n,0,overtones}]];
,
ansatz=Sum[If[Not@reim,If[amphaseshift,x=ToExpression["x"<>ToString[n]]*Exp[I ToExpression["y"<>ToString[n]]],x=ToExpression["x"<>ToString[n]]],x=ToExpression[("x"<>ToString[n])]+I ToExpression[("y"<>ToString[n])]]; x Exp[-(Global`t-ti)(1/(\[Tau]s[[n+1]](1+ToExpression["\[Beta]"<>ToString[n]])))] (Cos[\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]]) (Global`t)]+I Sin[ \[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]]) (Global`t)] ),{n,0,overtones}];
If[mixing[[1]],
\[Omega]m\[Tau]sm=Table[\[Omega]lmn[lm,mm,n,\[Eta],\[Chi]1,\[Chi]2,"ModesFile"->modesfile],{n,0,overtones}];
\[Omega]m=\[Omega]m\[Tau]sm[[All,1]];
\[Tau]sm=\[Omega]m\[Tau]sm[[All,2]];
ansatz=ansatz+Sum[If[Not@reim,x=ToExpression["xm"<>ToString[n]],x=ToExpression["xm"<>ToString[n]]+I ToExpression["ym"<>ToString[n]]];x Exp[-(Global`t-ti)/(\[Tau]sm[[n+1]])] (Cos[\[Omega]m[[n+1]] (Global`t)]+I Sin[ \[Omega]m[[n+1]] (Global`t)] ),{n,0,overtones}]];
];
,
If[qualfactorQ,
re=Sum[ToExpression["x"<>ToString[n]] Exp[-(Global`t-ti)(\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]])/(qfact[[n+1]](1+ToExpression["\[Beta]"<>ToString[n]])))]Cos[\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]]) (Global`t)+ToExpression["\[Phi]"<>ToString[n]]],{n,0,overtones}]; im=Sum[ToExpression["x"<>ToString[n]] Exp[-(Global`t-ti)(\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]])/(qfact[[n+1]](1+ToExpression["\[Beta]"<>ToString[n]])))]Sin[\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]]) (Global`t)+ToExpression["\[Phi]"<>ToString[n]]],{n,0,overtones}];
If[mixing[[1]],
\[Omega]m\[Tau]sm=Table[\[Omega]lmn[lm,mm,n,\[Eta],\[Chi]1,\[Chi]2,"ModesFile"->modesfile],{n,0,overtones}];
\[Omega]m=\[Omega]m\[Tau]sm[[All,1]];
\[Tau]sm=\[Omega]m\[Tau]sm[[All,2]];
qfactm=\[Omega]m*\[Tau]sm;
rem=Sum[ToExpression["x"<>ToString[n]] Exp[-(Global`t-ti)(\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]])/(qfactm[[n+1]](1+ToExpression["\[Beta]"<>ToString[n]])))] (Cos[\[Omega]m[[n+1]] (Global`t)]+I Sin[ \[Omega]m[[n+1]] (Global`t)] ),{n,0,overtones}];
imm=Sum[ToExpression["x"<>ToString[n]] Exp[-(Global`t-ti)(\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]])/(qfactm[[n+1]](1+ToExpression["\[Beta]"<>ToString[n]])))] (Sin[\[Omega]m[[n+1]] (Global`t)]+I Sin[ \[Omega]m[[n+1]] (Global`t)] ),{n,0,overtones}];
,
rem=0;
imm=0;
];
,
re=Sum[ToExpression["x"<>ToString[n]] Exp[-(Global`t-ti)/(\[Tau]s[[n+1]](1+ToExpression["\[Beta]"<>ToString[n]]))]Cos[\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]]) (Global`t)+ToExpression["\[Phi]"<>ToString[n]]],{n,0,overtones}];
im=Sum[ToExpression["x"<>ToString[n]] Exp[-(Global`t-ti)/(\[Tau]s[[n+1]](1+ToExpression["\[Beta]"<>ToString[n]]))]Sin[\[Omega]s[[n+1]](1+ToExpression["\[Alpha]"<>ToString[n]]) (Global`t)+ToExpression["\[Phi]"<>ToString[n]]],{n,0,overtones}];
If[mixing[[1]],
\[Omega]m\[Tau]sm=Table[\[Omega]lmn[lm,mm,n,\[Eta],\[Chi]1,\[Chi]2,"ModesFile"->modesfile],{n,0,overtones}];
\[Omega]m=\[Omega]m\[Tau]sm[[All,1]];
\[Tau]sm=\[Omega]m\[Tau]sm[[All,2]];
rem=Sum[ToExpression["x"<>ToString[n]] Exp[-(Global`t-ti)/(\[Tau]sm[[n+1]])] (Cos[\[Omega]m[[n+1]] (Global`t)]+I Sin[ \[Omega]m[[n+1]] (Global`t)] ),{n,0,overtones}];
imm=Sum[ToExpression["x"<>ToString[n]] Exp[-(Global`t-ti)/(\[Tau]sm[[n+1]])] (Sin[\[Omega]m[[n+1]] (Global`t)]+I Sin[ \[Omega]m[[n+1]] (Global`t)] ),{n,0,overtones}];
,
rem=0;
imm=0;
];
];
ampansatz=Sqrt[(re+rem)^2+(im+imm)^2];
phaseansatz=ArcTan[im/re];
ansatz={ampansatz,phaseansatz};
];
modto0=Complement[Table[i,{i,0,overtones}],fit\[Alpha]];
modto0=Table[ToExpression["\[Alpha]"<>ToString[modto0[[i]]]],{i,Length@modto0}];
ansatz=ansatz/.(Table[modto0[[i]]->0,{i,Length@modto0}]);
modto0=Complement[Table[i,{i,0,overtones}],fit\[Tau]];
modto0=Table[ToExpression["\[Beta]"<>ToString[modto0[[i]]]],{i,Length@modto0}];
ansatz=ansatz/.(Table[modto0[[i]]->0,{i,Length@modto0}]);
If[ex\[Omega]val,{ansatz,var},ansatz]
]
Options[AICcRes]=Options[AICcRes];
AICcRes[res_,coeff_]:=Module[{err,n,newfit,bracketedvars,ress,weigths},
err=1;
n=Length@res;
ress=Total[err (Abs[res])^2];
n Log[ress/n]+2(coeff+1)+(2*coeff(coeff+1) )/(n - coeff -1)
]
Options[BICRes]=Options[AICcRes];
BICRes[res_,coeff_]:=Module[{err,fact1,fact2,n,newfit,bracketedvars,ress,weigths},
err=1;
n=Length@res;
ress=Total[err (Abs[res])^2];
fact1=n Log[ress/n];
fact2=Log[n](coeff+1.);
fact1+fact2
]
Options[PValueTest]={"IRoot"->10};
PValueTest[dfnum_,dfden_,signl_,OptionsPattern[]]:=Module[{xmax,iroot},
iroot=OptionValue["IRoot"];xmax/.FindRoot[NIntegrate[PDF[FRatioDistribution[dfnum,dfden],x],{x,0,xmax}]==1-signl,{xmax,iroot}]]
HessianH[f_,t_List?VectorQ]:=D[f,{t,2}]
Options[CovarianceMatrix]={"Tolerance"->10^(-32)};
CovarianceMatrix[Xx_,length_,pars_,truepars_,OptionsPattern[]]:=Module[{hessian,tol,resreim},
tol=OptionValue["Tolerance"];
hessian=Rationalize[HessianH[Log[Xx],pars]/.truepars,tol];
(2Inverse[hessian])/(length-2)*1.]
Options[CorrelationMatrix]=Options[CovarianceMatrix];
CorrelationMatrix[Xx_,length_,pars_,truepars_,OptionsPattern[]]:=Module[{Ds,Dsi,covmat,hessian,tol},
tol=OptionValue["Tolerance"];
covmat=CovarianceMatrix[Xx,length,pars,truepars];
Ds=Rationalize[DiagonalMatrix[Diagonal[Sqrt[Re@covmat]]],tol];
Dsi=Inverse[Ds];
Dsi.covmat.Dsi
]
LogLikelihoodDist[data_,ansatz_,vars_,x_]:=Module[{ansatzlist},
ansatzlist=ansatz/.vars/.x->data[[All,1]];
-Total[(Re[ansatzlist-data[[All,2]]]^2+Im[ansatzlist-data[[All,2]]]^2)]]
Options[ComputeFitBands]={"ConfidenceLevel"->0.9};
ComputeFitBands[data_,ansatz_,vars_,pars_,OptionsPattern[];]:=Module[{ansatzv,colevel,covmat,dvec,lower,res,resce,resre,resim,times,tval,Xx,upper,a1},
colevel=OptionValue["ConfidenceLevel"];
tval=StudentTCI[0,1,Length@vars,ConfidenceLevel->colevel][[2]];
res=(ansatz/.t->data[[All,1]])-data[[All,2]];
resce=ComplexExpand@Total[res];
(resre=Simplify[Re[ComplexExpand[resce]],Element[vars,Reals]]);
(resim=Simplify[Im[ComplexExpand[resce]],Element[vars,Reals]]);
Xx=resre^2+resim^2;
covmat=CovarianceMatrix[Xx,Length@res,vars,pars];
times=data[[All,1]];
lower=Table[{times[[i]],ansatzv=ansatz/.t->times[[i]];dvec=Table[D[ansatzv,vars[[j]]],{j,Length@vars}]/.pars;(ansatzv/.pars)-tval Sqrt[dvec.covmat.dvec]},{i,Length@times}];
upper=Table[{times[[i]],ansatzv=ansatz/.t->times[[i]];dvec=Table[D[ansatzv,vars[[j]]],{j,Length@vars}]/.pars;(ansatzv/.pars)+tval Sqrt[dvec.covmat.dvec]},{i,Length@times}];
{lower,upper}
]
(* ::Code::Initialization:: *)
End[];
EndPackage[];