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

components.py

Blame
  • 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[];