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

test_hom.py

Blame
  • BBHReduce.m 131.20 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.                                                         *)
    (************************************************************************)
    
    
    
    BeginPackage["BBHReduce`",{"NRLists`","NRFiles`","NRStrings`","NRWaves`","NinjaBase`","SXSFormat`","NRTimeSeries`"}];
    
    
    ValPrint::usage="ValPrint[x_?StringQ]";
    
    
    RunsFromRunsFile::usage="RunsFromRunsFile[file_?StringQ]";
    ContainsRun::usage="ContainsRun[directory_, runIdentifierStrings_]";
    IsBAMObsoleteDirectory::usage="IsBAMObsoleteDirectory[dirname_?StringQ]";
    IsBAMEvolutionParfile::usage="IsBAMEvolutionParfile[ filecontent_]";
    IsBAMInitialDataParfile::usage="IsBAMInitialDataParfile[ filecontent_]";
    IsBAMEvolutionDirectory::usage="IsBAMEvolutionDirectory[dirname_?StringQ]";
    ParfileInDirectory::usage="ParfileInDirectory[dirname_?StringQ]";
    levelFun::usage="levelFun[str_]";
    HasModesDirectory::usage="HasModesDirectory[dirname_?StringQ]";
    HasModesFiles::usage="HasModesFiles[dirname_?StringQ,pattern_?StringQ], pattern can e.g. be  \"hmod*\"";
    ModesDirectory::usage="ModesDirectory[dirname_?StringQ]";
    
    BAMDataDirectories::usage="BAMEvolutionDataDirectories[runName_?StringQ,
       OptionsPattern[{
    \"RunsRoot\"       \[Rule] HomeDirectory,
    \"ReducedRoot\"    \[Rule] HomeDirectory ,
    \"TraverseLevels\" \[Rule] 4}]] searches for the evolution and reduced-data directories for a specific run";
    
    
    LocateInitialDataDirectory::usage="LocateInitialDataDirectory[rootDir_, psidfile_]";
    InitialDataParameters::usage="InitialDataParameters[idDir_, psidfile_]";
    
    ParfileToRules::usage="ParfileToRules[filename_String] converts a parameter file to a list of rules.";
    BAMParfileToRules::usage="BAMParfileToRules[filename_String] coverts a initial data parameter file to a list of rules.";
    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 ";
    BAMMetaFilesToRules::usage="BAMMetaFilesToRules[filename_String] converts a BAM .bbh file to a list of rules."
    RITMetaFilesToRules::usage="RITMetaFilesToRules[filename_String] converts a RIT Metadata file to a list of rules.";
    RITParClassification::usage="RITParClassification[ritdir_,ClassStr_]. Given a list of RIT NR. data folders 'ritdir', it returns all the cases that match a certain criterion 'ClassStr' (MassRatio range, Precessing or not, Initial Distance, Orbits Number)taking as reference the RIT metadata.txt files. If it is used iteratively, one could do different classifications ";
    
    
    BAMStringParameter::usage="BAMStringParameter[directory_, parametername_]";
    BAMNumberParameter::usage="BAMNumberParameter[directory_, parametername_]";
    BAMNumberParameterInFile::usage="BAMNumberParameterInFile[file_, parametername_]";
    BAMNumberParametersInFile::usage="BAMNumberParametersInFile[file_, parametername_]";
    BAMNumberParameters::usage="BAMNumberParameters[directory_, parametername_]";
    PSIDHashedNumberParameter::usage="PSIDHashedNumberParameter[file_, parametername_]";
    PSIDNumberParameter::usage="PSIDNumberParameter[file_, parametername_]";
    PSIDReadHeader::usage="PSIDReadHeader[file_]";
    PSIDReadData::usage="PSIDReadData[file_] reads initial data for function \!\(\*FormBox[\(\(\*SubscriptBox[\(U\), \(ijk\)] = \\\ \(TraditionalForm\`\(\(U\)\((\)\*SubscriptBox[\(A\), \(i\)]\)\), \(TraditionalForm\`\*SubscriptBox[\(B\), \(j\)]\), \(TraditionalForm\`\*SubscriptBox[\(\[CurlyPhi]\), \(k\)]\)\)
    StyleBox[\")\",\nFontSize->10]\),
    TraditionalForm]\). The function which appears in the conformal factor is u = (A-1)U";
    
    PSID2Rules::usage="PSID2Rules[filename_?StringQ] converts a BAM PSID file to a list of rules.";
    
    NMovingLevels::usage="NMovingLevels[parRules_] computes the number of moving levels in a BAM parameter file from a list of rules parRules, which corresponds to the content of the parameter file.";
    
    
    CreateDataReduceDirectory::usage="CreateDataReduceDirectory[dirname_],
     CreateDataReduceDirectory[reduceroot_, dirname_]";
    LocateMode::usage="LocateMode[modeDir_, lmode_, mmode_]";
    LocateModes::usage="LocateModes[modeDir_, lmode_, mmode_]";
    CopyL2mode::usage="CopyL2mode[configStr_, modesdir_, reducedir_]";
    CopyL2modes::usage="CopyL2modes[configStr_, modesdir_, reducedir_]";
    CopyModes::usage="CopyModes[configStr_, modesdir_, reducedir_, Lmode_, Mmode_] ";
    CopyFiles::usage="CopyFiles[configStr_, modesdir_, reducedir_, patterns_]";
    
    
    FormatPunctureData::usage="FormatPunctureData[mp_,string_,M_].";
    SafeFormatPunctureData2::usage="SafeFormatPunctureData2[mp1_,mp2_,string_,m1_,m2_].";
    SafeFormatPunctureDataCactus::usage="SafeFormatPunctureDataCactus[mp_,string_,m1_,m2_].";
    
    BAMModesFilesTo3Col::usage="BAMModesFilesTo3Col[modesFile_] reads BAM style {r,i}psi4modes_* files and saves the content in 3-colums
    format as {time,Re@values,Im@values};"
    
    BAMHorizonFilesToNRARFiles::usage="BAMHorizonFilesToNRARFiles[File_,OptionsPattern[{\"DeleteSourceFiles\" \[Rule] False}]]
     converts a BAM-style horizon file to NRAR-style horizon mass, spin, and horizon trajectory files.";
    
    BAMTrajectoryFileTo4Col::usage="BAMTrajectoryFileTo4Col[File_,OptionsPattern[{\"DeleteSourceFiles\" \[Rule] False}]]
    converts a BAM-style puncture trajectory file to NRAR-style trajectory file.";
    
    
    BAMExtractionRadii::usage="BAMExtractionRadii[modesdir_] lists all BAM extraction radii for psi4.";
    
    SymmetriesFromParfile::usage="SymmetriesFromParfile[parfile_].";
    
    WriteModeDecompConfigFile::usage="WriteModeDecompConfigFile[DirectoryRules_?ListQ]";
    
    WaveExtractionRadii::usage="WaveExtractionRadii[rootDir_, modesDir_]";
    ADMReduce::usage="ADMReduce[rootDir_, reduceDir_]";
    BBHDataReduce::usage="BBHDataReduce[modesdir_, IDroot_, ReduceRoot_]";
    
    CurateBBHData::usage="CurateBBHData[modesdir_, IDroot_, ReduceRoot_] is a new development version of BBHDataReduce.";
    
    NRARPsi4ToStrain::usage="NRARPsi4ToStrain[metaFile_,LMAX_,\[Omega]GWStart_] reads a NRAR-style metadata file and converts 
    psi4 modes to strain using the FFI algorithm, exporting the results to a directory FFIStrainModes.";
    
    
    CactusThornsAvailable::usage="CactusThornsAvailable[cactusdir_] list available Cactus thorns.";
    NormalizeCactusParfile::usage="NormalizeCactusParfile[parfile_,cactusDir_,OptionsPattern[
    {\"ThornListOutputFile\"\[Rule] \"\",\"ParameterOutputFile\"\[Rule] \"\"}]] normalizes a Cactus parfile for easier comparison or parsing.";
    
    CompleteCactusParfile::usage="CompleteCactusParfile[parfile_,cactusDir_,OptionsPattern[
    {\"ThornListOutputFile\"\[Rule] \"\",\"ParameterOutputFile\"\[Rule] \"\"}]]";
    
    SetParfileEntryValue::usage="SetParfileEntryValue[text_,key_,value_]";
    SetParfileValue::usage="SetParfileValue[text_,key_,value_]";
    SetParfileVectorValue::usage="SetParfileVectorValue[text_,key_,component_,value_]";
    
    
    SXSLuminosityFromMetaFiles::usage="SXSLuminosityFromMetaFiles:[metaFile_,modes_]; Computes Luminosity from the SXS metadata given a list of modes";
    
    
    AHBAMsmall::usage="AHBAMsmall[q,s]. Computes the aparent horizon (horizon radius) for the smallest black hole in the BAM coordinates. See notebook /BBHReduce/AparentHorizonFit.nb"
    AHBAMbig::usage="AHBAMbig[q,s]. Computes the aparent horizon (horizon radius) for the biggest black hole in the BAM coordinates. See notebook /BBHReduce/AparentHorizonFit.nb"
    AHBAMsmall2017::usage="AHBAMsmall2017[q,s]. Computes the aparent horizon (horizon radius) for the smallest black hole in the BAM coordinates. See notebook /BBHReduce/AparentHorizonFit.nb"
    AHBAMbig2017::usage="AHBAMbig2017[q,s]. Computes the aparent horizon (horizon radius) for the biggest black hole in the BAM coordinates. See notebook /BBHReduce/AparentHorizonFit.nb"
    
    
    Begin["`Private`"];
    
    
    ValPrint[x_?StringQ]:=Print[First@StringSplit[x,"$"]<>" = ",ToExpression@x]
    
    
    RunsFromRunsFile[file_?StringQ]:=Module[{content,runs},
    content = StringSplit[Import[file,"String"],EndOfLine];
    
    runs=Flatten@Map[StringCases[#,"*"~~x:Except[WhitespaceCharacter]...~~  WhitespaceCharacter... -> x]&, content];
    StringTrim[#,RegularExpression["/$"]]&/@runs
    ];
    
    
    ContainsRun[directory_,runIdentifierStrings_]:=Length@Intersection[{LastInPath@directory},runIdentifierStrings]>0
    
    
    IsBAMObsoleteDirectory[dirname_?StringQ]:=StringMatchQ[dirname, "*_old"]||StringMatchQ[dirname, "*_previous"]
    
    
    IsBAMEvolutionParfile[fileORfilecontent_]:=Module[{joined},
    
    If[ListQ@fileORfilecontent,
    joined=StringJoin@fileORfilecontent;,
    
    If[FileType@fileORfilecontent == File,
      joined = StringJoin@Import[fileORfilecontent,"String"];
    ];
    ];
    
    StringMatchQ[joined,"*bampi_*"]&&
    StringMatchQ[joined,"*amr*"]&&StringMatchQ[joined,"*physics*"]
    ];
    
    
    IsBAMInitialDataParfile[ fileORfilecontent_]:=Module[{joined},
    
    If[FileType@fileORfilecontent == File,
    joined = StringJoin@Import[fileORfilecontent,"String"];
    ];
    
    If[ListQ@fileORfilecontent,
    joined=StringJoin@fileORfilecontent;
    ];
    
     StringMatchQ[joined,"*nx*"]&&
    StringMatchQ[joined,"*iterate*"]&&StringMatchQ[joined,"*physics*"]
    ];
    
    
    IsBAMEvolutionDirectory[dirname_?StringQ]:=Module[{parfiles,content,i},
    
    parfiles=Join[
    FileNames["*.par",dirname],
    FileNames["*.par.gz",dirname],
    FileNames["*.par.bz2",dirname]];
    
    content=Table[Import[parfiles[[i]],"String"],{i,1,Length@parfiles}];
    
    IsBAMEvolutionParfile@content
    ]
    
    
    ParfileInDirectory[dirname_?StringQ,OptionsPattern[{"Style"-> "BAM"}]]:=Module[{guess,parfiles,sel,style},
    
    style = OptionValue["Style"];
    
    guess=dirname<>"/"<>LastInPath@dirname<>".par";
    If[FileType@guess == File, 
    
     sel = guess,
    
     parfiles=Join[
      FileNames["*.par",dirname],
      FileNames["*.par.gz",dirname],
      FileNames["*.par.bz2",dirname]
      ];
    
    Print[parfiles];
    
    sel=Switch[style,
    "BAM",First@Select[parfiles, IsBAMEvolutionParfile@# || IsBAMInitialDataParfile@# &],
    "Cactus",First@parfiles,
    _,First@parfiles
    ];
    ];
    
    Print["ParfileInDirectory identifies parameter file ", sel];
    
    sel
    ];
    
    
    levelFun[str_]:=ToExpression@First@StringCases[str,"hmod.r"~~r:NumberString..~~".l"~~l:NumberString-> {r,l}];
    
    
    HasModesDirectory[dirname_?StringQ]:=
    safeDirectoryQ[StringReplace[dirname<>"/Modes","//"->"/"]]||
    safeDirectoryQ[StringReplace[dirname<>"/Analysis","//"->"/"]]||
    safeDirectoryQ[StringReplace[dirname<>"/analysis","//"->"/"]]||
    safeDirectoryQ[StringReplace[dirname<>"/Psi4ModeDecomp","//"->"/"]]||
    safeDirectoryQ[StringReplace[dirname<>"/NinjaCleanPsi","//"->"/"]]
    
    
    HasModesFiles[dirname_?StringQ,pattern_?StringQ]:=Length@FileNames[pattern,dirname,2]>0
    
    
    ModesDirectory[dirname_?StringQ,pattern_?StringQ]:=Module[{sel},
    
    sel = Select[FileNames["*",dirname],FileType@# == Directory&];
    sel = Select[sel,HasModesFiles[#,pattern]&];
    
    If[Length@sel > 0, First@sel,""]
    ];
    
    
    BAMDataDirectories[runName_?StringQ,
       OptionsPattern[{
    "RunsRoot"        -> HomeDirectory,
    "ReducedRoot"     -> Global`BBHDataDir,
    "InitialDataRoot" -> HomeDirectory,
    "TraverseLevels" -> 4}]
    ]:=Module[{RunsRoot,ReducedRoot,TraverseLevels,runsFound,
    RunsDirectories,obsolete,reducedFound,ReducedDirectories,evolDirHasModes,reducedDirHasModes,
    evolutionModesDir,reducedModesDir,psidfile,InitialDataRoot,IDDir},
    
    RunsRoot        = OptionValue["RunsRoot"];
    ReducedRoot     = OptionValue["ReducedRoot"];
    InitialDataRoot = OptionValue["InitialDataRoot"];
    TraverseLevels  = OptionValue["TraverseLevels"];
    
    (* directory with the original run *)
    runsFound = FileNames[runName,RunsRoot,TraverseLevels];
    Print["Found runs directories:            ", RunsDirectories = Select[runsFound,FileType@#==Directory&]];
    Print[" - evolution directories:          ", RunsDirectories = Select[RunsDirectories,IsBAMEvolutionDirectory]];
    obsolete = Select[RunsDirectories,IsBAMObsoleteDirectory];
    Print[" - valid evolution directories:    ", RunsDirectories = Complement[RunsDirectories,obsolete]];
    
    If[Length@RunsDirectories == 1, 
      RunsDirectories = RunsDirectories[[1]];
      Print["Evolution directory has psi4 modes: ", evolDirHasModes = HasModesFiles[RunsDirectories,"psi3col*"]];, 
      Print["No unique runs directory found"];
      RunsDirectories = False;
      evolDirHasModes = False;
    ];
    
    
    (* directory with reduced data *)
    reducedFound = FileNames[RunDir,ReducedRoot,TraverseLevels];
    Print["Found reduced-data directories: ", ReducedDirectories=Select[reducedFound,FileType@#==Directory&]];
    
    If[Length@ReducedDirectories == 1, 
      ReducedDirectories = ReducedDirectories[[1]];
      Print["Reduced directory has psi4 modes:",  reducedDirHasModes = HasModesFiles[ReducedDirectories,"psi3col*"]];, 
      Print["No unique reduced-data directory found"];
      ReducedDirectories = False;
      reducedDirHasModes = False;
    ];
    
    If[evolDirHasModes,    evolutionModesDir = ModesDirectory[ReducedDirectories,"psi3*"], evolutionModesDir = False];
    If[reducedDirHasModes, reducedModesDir   = ModesDirectory[RunsDirectories,   "psi3*"],   reducedModesDir = False];
    
    
    If[StringQ@RunsDirectories,
     psidfile = BAMStringParameter[RunsDirectories,"punctures_ps_file"];
     psidfile = Last@StringSplit[psidfile,"/"];
     psidfile = StringReplace[psidfile," "->""];
     Print["psid-file determined from evolution parameter file as: ", psidfile];
     IDDir = LocateInitialDataDirectory[InitialDataRoot,psidfile];,
     IDDir = False;
    ];
    
    {"EvolutionDirectory" -> RunsDirectories, "ReducedDirectory" -> ReducedDirectories,
    "EvolutionDirectoryHasModes" -> evolDirHasModes,"ReducedDirectoryHasModes" -> reducedDirHasModes,
    "EvolutionModesDir"-> evolutionModesDir, "ReducedModesDir"-> reducedModesDir, "InitialDataDir" -> IDDir, 
    "PSIDFile"-> IDDir<> "/"<> psidfile}
    ]
    
    
    LocateInitialDataDirectory[rootDir_,psidfile_]:=Module[{file,psidFiles,levels=5,guess,result},
    file=StringReplace[psidfile," "->""];
    guess=rootDir<>"/"<>StringSplit[file,"."][[1]]<>"/"<>file;
    Print["Making guess for psid file:", guess];
    If[FileType@guess==File,
    Print["guessed id-file confirmed to be ",guess];
    result=rootDir<>"/"<>StringSplit[file,"."][[1]];
    Print["Taking as psid-result: ", result];
    Return[result]
    ];
    
    Print["Looking for psid file: ",file," in directory ",rootDir];
    psidFiles=Union@FileNames[file,rootDir,levels];
    Print["Found matching psid-files:",psidFiles];
    
    If[Length@psidFiles>1,Print["taking first psid-file"];
    ];
    
    If[psidFiles!={},
    result=DropLastDirectory[First@psidFiles];,Print["ERROR: LocateInitialDataDirectory could not locate ID-directory"];
    result="";
    ];
    
    ToString@result
    ];
    
    
    InitialDataParameters[psidfile_]:=Module[{file,m1,m2,madm,sep,abschi1,abschi2,s1x,s1y,s1z,s2x,s2y,s2z,
    x1,y1,z1,x2,y2,z2,px1,py1,pz1,px2,py2,pz2},
    
    file=psidfile;
    Print["InitialDataParameters: Processing psidfile ", file]; 
    m1=PSIDHashedNumberParameter[file,"M1"];
    m2=PSIDHashedNumberParameter[file,"M2"];
    madm=PSIDHashedNumberParameter[file,"Madm"];
    
    sep=PSIDHashedNumberParameter[file,"d"];
    abschi1=PSIDHashedNumberParameter[file,"S1/M1^2"];
    abschi2=PSIDHashedNumberParameter[file,"S2/M2^2"];
    
    s1x=PSIDNumberParameter[file,"bhsx1"];
    s1y=PSIDNumberParameter[file,"bhsy1"];
    s1z=PSIDNumberParameter[file,"bhsz1"];
    
    s2x=PSIDNumberParameter[file,"bhsx2"];
    s2y=PSIDNumberParameter[file,"bhsy2"];
    s2z=PSIDNumberParameter[file,"bhsz2"];
    
    x1=PSIDNumberParameter[file,"bhx1"];
    y1=PSIDNumberParameter[file,"bhy1"];
    z1=PSIDNumberParameter[file,"bhz1"];
    
    x2=PSIDNumberParameter[file,"bhx2"];
    y2=PSIDNumberParameter[file,"bhy2"];
    z2=PSIDNumberParameter[file,"bhz2"];
    
    px1=PSIDNumberParameter[file,"bhpx1"];
    py1=PSIDNumberParameter[file,"bhpy1"];
    pz1=PSIDNumberParameter[file,"bhpz1"];
    
    px2=PSIDNumberParameter[file,"bhpx2"];
    py2=PSIDNumberParameter[file,"bhpy2"];
    pz2=PSIDNumberParameter[file,"bhpz2"];
    
    {m1,m2,madm,sep,abschi1,abschi2,s1x,s1y,s1z,s2x,s2y,s2z,{x1,y1,z1},{x2,y2,z2},{px1,py1,pz1},{px2,py2,pz2}}
    ];
    
    
    InitialDataParameters[idDir_,psidfile_]:=Module[{file},
    
    file=idDir<>"/"<>psidfile;
    
    InitialDataParameters[file]
    ];
    
    
    ParfileToRules[filename_String]:=Module[{parlines,rules},
    	parlines=ReadList[filename,String];
    
    	rules=Flatten[DefinitionsFromString/@parlines];
    
    	(*modify selected rules*)
    	rules=ReplaceRuleInRuleList[rules,"amr_nxyz",ToExpression/@StringSplit["amr_nxyz"/.rules]];
    
    	rules
    ]
    
    
    BAMParfileToRules[filename_String]:=Module[{parlines,rules,par1,rules2,pos,rulesPn,var,value,list},
    	parlines=ReadList[filename,String];
    
        rules=Delete[parlines,Position[StringMatchQ[parlines,"*#*"],True]];
        rulesPn=StringTrim/@StringSplit[Flatten@StringSplit[Flatten[StringSplit[Select[parlines,StringMatchQ[#,"#$$*"]&],"#"]],"$$"],"="];
        (*pos=Flatten@Position[TakeColumn[rules2,1],StringMatchQ[" ",#]&,True];*)
      
    (*	
    	(*modify selected rules*)
    	rules=ReplaceRuleInRuleList[rules,"amr_nxyz",ToExpression/@StringSplit["amr_nxyz"/.rules]];
    *)
    	
        rules=StringTrim/@StringSplit[rules,"="];
        rulesPn=Flatten[ToExpression@StringReplace[ToString/@Select[StringSplit[ToString/@rulesPn," "],Length@#==2&],",,"->","],1];
        rules=Join[rules,rulesPn];
    
        var=ToString/@rules[[All,1]];
        value=ToExpression/@rules[[All,2]];
       
        
        list={};
        Do[
          If[IsFPNumberQ@value[[i]],
    	  value[[i]]=(StringToNumber@#)&/@value[[i]];
            ];
          list=AppendTo[list,{var[[i]]->value[[i]]}];
         ,
          {i,1,Length@value}];
       
         Flatten@list
    
    ]
    
    
    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
    
    ]
    
    
    RITMetaFilesToRules[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]]],"<"\[Rule]""];*)
    (*We do not have < in the MetaFiles of RIT*)
    
    (*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
    
    ]
    
    
    BAMMetaFilesToRules[filePath_]:=Module[{filepath,fileList,meta1,pos1,meta2,meta3,meta4,value,var,list,pos,rads,posdel},
    
    
    If[ListQ@filePath,filepath=filePath[[1]],filepath=filePath];
    If[Not@FileExistsQ@filepath,Print["File not found"];Return[]];
    
    (*Reading the file*);
    fileList=ToString/@ReadList[filepath,String];
    
    (*Delete comments*)
    meta1=Delete[fileList,Position[StringMatchQ[fileList,"*#*"],True]];
    
    (*Find = *)
    meta2=meta1[[Flatten@Position[StringMatchQ[meta1,"*=*"],True]]];
    meta3=StringSplit[meta2,"="];
    (*Select non-empty arrays*)
    meta4=Select[meta3,Length@#==2&];
    
    (*Delete spaces*)
    meta4=Transpose[{StringTrim/@meta4[[All,1]],StringTrim/@meta4[[All,2]]}];
    (*Delete empty elements*)
    meta4=Select[meta4, UnsameQ[#[[2]], {}] &];
    
    var=meta4[[All,1]];
    value=meta4[[All,2]];
    
    value=StringSplit[#,","]&/@value;
    
    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}];
    
    list=Flatten@list;
    
    (* Fix extraction radius repetition. Assume that the first entry is finite-radii *)
    pos=Position[Flatten@list,"extraction-radius"];
    rads=Select[Flatten@list[[pos[[All,1]],2]],NumberQ@#&];
    list[[pos[[1,1]],2]]=rads;
    posdel=Partition[TakeColumn[pos,1][[2;;-1]],1];
    list=Delete[list,posdel];
    
    (* Gather together the extraction radius per each mode *)
    pos=Select[DuplicatesPosition[list[[All,1]]],Length@#>1&];
    Do[list[[pos[[i,1]],2]]=Flatten@Join[{list[[pos[[i,1]],2]]},list[[pos[[i,2;;-1]],2]]],{i,Length@pos}];
    posdel=Partition[Flatten@Table[pos[[i,2;;-1]],{i,Length@pos}],1];
    
    list=Delete[list,posdel]
    
    ]
    
    
    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[{}]];
    
    (* Useful function to detect the duplicates *)
    positionDuplicates[listaux_]:=GatherBy[Range@Length[listaux],listaux[[#]]&];
    
    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]]]
    ]
    
    
    SXSParClassification2[sxsdir_?ListQ,ClassStr_?ListQ,OptionsPattern[{"\[Epsilon]"->0.001,"HighSpin"->0.8,"UnRepeated"->False,"Verbose"->False}]]:=Module[{metafiles,metadata,
    orbitStr="number-of-orbits",dStr="initial-separation",mass1Str="relaxed-mass1",mass2Str="relaxed-mass2",spin1Str="relaxed-spin1",spin2Str="relaxed-spin1",
    STTOV="relaxed-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,parmatrix,myindex},
    
    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"];
    
    
    
    (* Useful function to detect the duplicates *)
    positionDuplicates[listaux_]:=GatherBy[Range@Length[listaux],listaux[[#]]&];
    
    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)/(mass1*mass1);
    spin2Dim=((spin2Str/.#)&/@metadata)/(mass2*mass2);
    \[Chi]eff=massratio/(1.+massratio)*spin1Dim +1./(1+massratio)*spin2Dim;
    spinzDiff=Abs[(#[[2]]-#[[1]])&/@Transpose[{TakeColumn[spin1Dim,3],TakeColumn[spin2Dim,3]}]];
    A1=(1+3. massratio/(4.));
    A2=(1+3. /(4.*massratio) );
    eccentricity=Flatten[(eccStr/.#)&/@metadata];
    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;
    spin1Norm=Norm[#]&/@spin1Dim;
    spin2Norm=Norm[#]&/@spin2Dim;
    
    parmatrix=Transpose[{massratio,dist,orbit,precvalueNorm,spinzDiff,spin1Norm,spin2Norm,\[Chi]eff,spin1Dim,spin2Dim,eccentricity,sxsdirout}];
    
    condition=ClassStr[[All,1]];
    select=Table[If[Length@ClassStr[[i]]==2,ToExpression@(ClassStr[[i,2]]),"Null"],{i,1,Length@ClassStr}];
    
    Do[
    
    Which[condition[[i]]== "MassRatio",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=1;
    pos=Flatten@Position[parmatrix[[All,myindex]],_?(Evaluate[select[[i]]]&)];
    parmatrix=parmatrix[[pos]];,
    
    condition[[i]]== "Distance",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=2;
    parmatrix=Select[parmatrix,Evaluate[StringReplace[select[[i]],"#"->ToString@(#[[myindex]])&]]];,
    
    condition[[i]]== "Orbits",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=3;
    pos=Flatten@Position[parmatrix[[All,myindex]],_?(Evaluate[select[[i]]]&)];
    parmatrix=parmatrix[[pos]];,
    
    condition[[i]]== "Non-Precessing",
    
    myindex=4;
    
    parmatrix=Select[parmatrix,#[[myindex]]<\[Epsilon]&];,
    
    condition[[i]]== "Unequal",
    
    myindex=5;
    parmatrix=Select[parmatrix,#[[myindex]]>\[Epsilon]&];,
    
    condition[[i]]== "Precessing",
    
    myindex=6;
    parmatrix=Select[parmatrix,#[[myindex]]>=\[Epsilon]&];,
    
    condition[[i]]== "High-Spin",
    
    myindex=7;
    parmatrix=Select[parmatrix,#[[myindex]]>=highspin& ||#[[myindex+1]]>=highspin& ];,
    
    condition[[i]]== "\[Chi]eff",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=8;
    pos=Flatten@Position[parmatrix[[All,myindex]],_?(Evaluate[select[[i]]]&)];
    parmatrix=parmatrix[[pos]];,
    
    condition[[i]]== "\[Chi]1",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=9;
    pos=Flatten@Position[parmatrix[[All,myindex]],_?(Evaluate[select[[i]]]&)];
    parmatrix=parmatrix[[pos]];,
    
    condition[[i]]== "\[Chi]2",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=10;
    pos=Flatten@Position[parmatrix[[All,myindex]],_?(Evaluate[select[[i]]]&)];
    parmatrix=parmatrix[[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[#&/@parmatrix[[All,1]],0.1],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,9]],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,10]],parmatrix[[All,2]],parmatrix[[All,11]]}];
    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}];
    parmatrix=parmatrix[[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]];
    parmatrix=parmatrix[[posdistecc]];
    
    parmatrix=SortBy[parmatrix,#[[1]]&];
    sxsdiroutaux=Transpose[{#&/@parmatrix[[All,12]],Round[#&/@parmatrix[[All,1]],0.1],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,9]],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,10]],parmatrix[[All,2]],parmatrix[[All,11]]}];
    
    If[verbose, 
    Print[Prepend[Table[ToString@#&/@sxsdiroutaux[[i]],{i,1,Length@sxsdiroutaux}],{"Case","q","\[Chi]1","\[Chi]2","D","e"}]//TableForm]];
    
    parmatrix[[All,12]],
    
    sxsdiroutaux=Transpose[{#&/@parmatrix[[All,12]],Round[#&/@parmatrix[[All,1]],0.1],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,9]],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,10]],parmatrix[[All,2]],parmatrix[[All,11]]}];
    
    (*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]];
    
    parmatrix[[All,12]]]
    ]
    
    
    RITParClassification[ritdir_?ListQ,ClassStr_?ListQ,OptionsPattern[{"\[Epsilon]"->0.001,"HighSpin"->0.8,"UnRepeated"->False,"Verbose"->False}]]:=Module[{metafiles,metadata,
    orbitStr="number-of-cycles-22",dStr="initial-separation",mass1Str="initial-mass1",mass2Str="initial-mass2",spin1x="initial-bh-chi1x",spin1y="initial-bh-chi1y",spin1z="initial-bh-chi1z",
    spin2x="initial-bh-chi2x",spin2y="initial-bh-chi2y",spin2z="initial-bh-chi2z",eccStr="eccentricity",spin1Dim,spin2Dim,mass1,mass2,massratio,eccentricity,dist,orbit,select,pos,condition,A1,A2,precvalue,precvalueNorm,\[Epsilon],
    spin1Norm,spin2Norm,highspin,spintest,\[Chi]eff,ritdirout,spinz,spinzDiff,auxDist,posdup,posdupDist,posdistecc,unrepeated,verbose,ritdiroutaux,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 rit file names are not consistent with the metadata files"];
    Print["The spin definition is consitent with 'initial-spin' values and not relaxed ones"];
    
    (* Useful function to detect the duplicates *)
    positionDuplicates[listaux_]:=GatherBy[Range@Length[listaux],listaux[[#]]&];
    
    ritdirout=ritdir;
    
    \[Epsilon]=OptionValue["\[Epsilon]"];
    highspin=OptionValue["HighSpin"];
    unrepeated=OptionValue["UnRepeated"];
    verbose=OptionValue["Verbose"];
    
    
    metafiles=Flatten[FileNames["Metadata",#,2]&/@ritdirout,1];
    metadata=RITMetaFilesToRules[#]&/@metafiles;
    
    mass1=Flatten@((mass1Str/.#)&/@metadata);
    mass2=Flatten@((mass2Str/.#)&/@metadata);
    
    massratio=mass1/mass2;
    dist=Flatten@((dStr/.#)&/@metadata);
    orbit=Flatten@((orbitStr/.#)&/@metadata);
    spin1Dim=(({spin1x,spin1y,spin1z}/.#)&/@metadata)(*/(mass1*mass1)*); (* In RIT the spin is already adimensional *)
    spin2Dim=(({spin2x,spin2y,spin2z}/.#)&/@metadata)(*/(mass2*mass2)*);
    spin1Dim=Flatten[#,1]&/@spin1Dim;
    spin2Dim=Flatten[#,1]&/@spin2Dim;
    Do[  (* This is to set to 0 ths x and y spin components for aligned cases *)
     Do[
       If[Not@NumberQ[spin1Dim[[i,j]]],spin1Dim[[i,j]]=0];
       If[Not@NumberQ[spin2Dim[[i,j]]],spin2Dim[[i,j]]=0];
     ,{j,1,2}];
    ,{i,Length@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=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,1,Length@ritdirout}];
    
    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]];
    
    ritdirout=ritdirout[[pos]];,
    
    condition[[i]]== "Distance",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    
    pos=Flatten@Position[ritdirout=ritdirout[[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]];
    
    ritdirout=ritdirout[[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]];
    
    ritdirout=ritdirout[[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]];
    
    ritdirout=ritdirout[[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]];
    
    ritdirout=ritdirout[[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]];
    
    ritdirout=ritdirout[[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]];
    
    ritdirout=ritdirout[[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]];
    
    ritdirout=ritdirout[[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]];
    
    ritdirout=ritdirout[[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]];
    
    ritdirout=ritdirout[[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]];
    
    ritdirout=ritdirout[[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}];
    ritdirout=ritdirout[[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]];
    ritdirout=ritdirout[[posdistecc]];
    ritdiroutaux=SortBy[Table[Join[{ritdirout[[i]]},auxDist[[i]]],{i,1,Length@ritdirout}],#[[2]]&];
    
    If[verbose, 
    Print[Prepend[Table[ToString@#&/@ritdiroutaux[[i]],{i,1,Length@ritdiroutaux}],{"Case","q","\[Chi]1","\[Chi]2","D","e"}]//TableForm]];
    
    ritdiroutaux[[All,1]],
    
    auxDist=Transpose[{Round[#&/@massratio,0.1],Round[Chop[#,10^(-2)],0.01]&/@spin1Dim,Round[Chop[#,10^(-2)],0.01]&/@spin2Dim,dist,eccentricity}];
    (*ritdiroutaux=Table[Join[{ritdirout[[i]]},auxDist[[i]]],{i,1,Length@ritdirout}];*)
    ritdiroutaux=SortBy[Table[Join[{ritdirout[[i]]},auxDist[[i]]],{i,1,Length@ritdirout}],#[[2]]&];
    If[verbose, 
    Print[Prepend[Table[ToString@#&/@ritdiroutaux[[i]],{i,1,Length@ritdiroutaux}],{"Case","q","\[Chi]1","\[Chi]2","D","e"}]//TableForm]];
    
    ritdiroutaux[[All,1]]]
    ]
    
    
    ritParClassification2[ritdir_?ListQ,ClassStr_?ListQ,OptionsPattern[{"\[Epsilon]"->0.001,"HighSpin"->0.8,"UnRepeated"->False,"Verbose"->False}]]:=Module[{metafiles,metadata,
    orbitStr="number-of-orbits",dStr="initial-separation",mass1Str="relaxed-mass1",mass2Str="relaxed-mass2",spin1Str="relaxed-spin1",
    spin2Str="relaxed-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,ritdirout,spinz,spinzDiff,auxDist,posdup,posdupDist,posdistecc,unrepeated,verbose,ritdiroutaux,precvalue1,precvalue2,
    precvalueNorm1,precvalueNorm2,parmatrix,myindex},
    
    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 rit file names are not consistent with the metadata files"];
    
    (* Useful function to detect the duplicates *)
    positionDuplicates[listaux_]:=GatherBy[Range@Length[listaux],listaux[[#]]&];
    
    ritdirout=ritdir;
    
    \[Epsilon]=OptionValue["\[Epsilon]"];
    highspin=OptionValue["HighSpin"];
    unrepeated=OptionValue["UnRepeated"];
    verbose=OptionValue["Verbose"];
    
    
    metafiles=Flatten[FileNames["metadata.txt",#,4]&/@ritdirout,1];
    metadata=ritMetaFilesToRules[#]&/@metafiles;
    
    
    mass1=Flatten@((mass1Str/.#)&/@metadata);
    mass2=Flatten@((mass2Str/.#)&/@metadata);
    
    massratio=mass1/mass2;
    dist=Flatten@((dStr/.#)&/@metadata);
    orbit=Flatten@((orbitStr/.#)&/@metadata);
    spin1Dim=((spin1Str/.#)&/@metadata)/(mass1*mass1);
    spin2Dim=((spin2Str/.#)&/@metadata)/(mass2*mass2);
    \[Chi]eff=massratio/(1.+massratio)*spin1Dim +1./(1+massratio)*spin2Dim;
    spinzDiff=Abs[(#[[2]]-#[[1]])&/@Transpose[{TakeColumn[spin1Dim,3],TakeColumn[spin2Dim,3]}]];
    A1=(1+3. massratio/(4.));
    A2=(1+3. /(4.*massratio) );
    eccentricity=Flatten[(eccStr/.#)&/@metadata];
    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;
    spin1Norm=Norm[#]&/@spin1Dim;
    spin2Norm=Norm[#]&/@spin2Dim;
    
    parmatrix=Transpose[{massratio,dist,orbit,precvalueNorm,spinzDiff,spin1Norm,spin2Norm,\[Chi]eff,spin1Dim,spin2Dim,eccentricity,ritdirout}];
    
    condition=ClassStr[[All,1]];
    select=Table[If[Length@ClassStr[[i]]==2,ToExpression@(ClassStr[[i,2]]),"Null"],{i,1,Length@ClassStr}];
    
    Do[
    
    Which[condition[[i]]== "MassRatio",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=1;
    pos=Flatten@Position[parmatrix[[All,myindex]],_?(Evaluate[select[[i]]]&)];
    parmatrix=parmatrix[[pos]];,
    
    condition[[i]]== "Distance",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=2;
    parmatrix=Select[parmatrix,Evaluate[StringReplace[select[[i]],"#"->ToString@(#[[myindex]])&]]];,
    
    condition[[i]]== "Orbits",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=3;
    pos=Flatten@Position[parmatrix[[All,myindex]],_?(Evaluate[select[[i]]]&)];
    parmatrix=parmatrix[[pos]];,
    
    condition[[i]]== "Non-Precessing",
    
    myindex=4;
    
    parmatrix=Select[parmatrix,#[[myindex]]<\[Epsilon]&];,
    
    condition[[i]]== "Unequal",
    
    myindex=5;
    parmatrix=Select[parmatrix,#[[myindex]]>\[Epsilon]&];,
    
    condition[[i]]== "Precessing",
    
    myindex=6;
    parmatrix=Select[parmatrix,#[[myindex]]>=\[Epsilon]&];,
    
    condition[[i]]== "High-Spin",
    
    myindex=7;
    parmatrix=Select[parmatrix,#[[myindex]]>=highspin& ||#[[myindex+1]]>=highspin& ];,
    
    condition[[i]]== "\[Chi]eff",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=8;
    pos=Flatten@Position[parmatrix[[All,myindex]],_?(Evaluate[select[[i]]]&)];
    parmatrix=parmatrix[[pos]];,
    
    condition[[i]]== "\[Chi]1",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=9;
    pos=Flatten@Position[parmatrix[[All,myindex]],_?(Evaluate[select[[i]]]&)];
    parmatrix=parmatrix[[pos]];,
    
    condition[[i]]== "\[Chi]2",
    
    If[Length@ClassStr[[i]]!= 2,Print["Wrong input"];Break[]];
    myindex=10;
    pos=Flatten@Position[parmatrix[[All,myindex]],_?(Evaluate[select[[i]]]&)];
    parmatrix=parmatrix[[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[#&/@parmatrix[[All,1]],0.1],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,9]],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,10]],parmatrix[[All,2]],parmatrix[[All,11]]}];
    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}];
    parmatrix=parmatrix[[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]];
    parmatrix=parmatrix[[posdistecc]];
    
    parmatrix=SortBy[parmatrix,#[[1]]&];
    ritdiroutaux=Transpose[{#&/@parmatrix[[All,12]],Round[#&/@parmatrix[[All,1]],0.1],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,9]],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,10]],parmatrix[[All,2]],parmatrix[[All,11]]}];
    
    If[verbose, 
    Print[Prepend[Table[ToString@#&/@ritdiroutaux[[i]],{i,1,Length@ritdiroutaux}],{"Case","q","\[Chi]1","\[Chi]2","D","e"}]//TableForm]];
    
    parmatrix[[All,12]],
    
    ritdiroutaux=Transpose[{#&/@parmatrix[[All,12]],Round[#&/@parmatrix[[All,1]],0.1],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,9]],Round[Chop[#,10^(-2)],0.01]&/@parmatrix[[All,10]],parmatrix[[All,2]],parmatrix[[All,11]]}];
    
    (*ritdiroutaux=SortBy[Table[Join[{ritdirout[[i]]},auxDist[[i]]],{i,1,Length@ritdirout}],#[[2]]&];*)
    If[verbose, 
    Print[Prepend[Table[ToString@#&/@ritdiroutaux[[i]],{i,1,Length@ritdiroutaux}],{"Case","q","\[Chi]1","\[Chi]2","D","e"}]//TableForm]];
    
    parmatrix[[All,12]]]
    ]
    
    
    BAMStringParameter[directory_,parametername_]:=Module[{file,dirname,content},
    
    file = directory;
    
     If[TrueQ[FileType@directory == Directory],
     file=ZippedOrUnzipped@ParfileInDirectory[directory];
    ];
    
    Print["BAMStringParameter identified parameter file ",file];
    
    If[FileType@file==File,content=StringSplit[Import[file,"String"],EndOfLine],content={}];
    If[content=={},Print["parameter file does not exist"],content=Map[StringCases[#,StartOfLine~~parametername~~WhitespaceCharacter...~~"="~~WhitespaceCharacter...~~x__~~WhitespaceCharacter...~~EndOfLine->x]&,content];
    content=First@Flatten@content;
    content=First@StringSplit[content,EndOfLine];];
    content
    ];
    
    
    BAMNumberParameter[directory_,parametername_]:=Module[{file,dirname,content},
    
    file = directory;
    
    If[TrueQ[FileType@directory == Directory],
     file=ZippedOrUnzipped@ParfileInDirectory[directory];
    ];
    
    Print["BAMStringParameter identified parameter file ",file];
    
    If[FileType@file==File,content=StringSplit[Import[file,"String"],EndOfLine],content={}];
    
    Print["Read ",Length@content, " parameter file entries"];
    
    If[content== {},
    Print["parameter file does not exist"],
    content=Flatten@Map[StringCases[#,parametername~~WhitespaceCharacter...~~"="~~WhitespaceCharacter...~~x:NumberString-> x]&, content];
    content=Map[ToExpression,content];
    
    content=First@Select[content,NumberQ]
    ];
    content
    ];
    
    
    BAMNumberParameterInFile[file_,parametername_,default_]:=Module[{dirname,content},
    
    If[FileType@file==File,content=StringSplit[Import[file,"String"],EndOfLine],content={}];
    
    If[content== {},
    Print["parameter file does not exist"],
    content=Flatten@Map[StringCases[#,parametername~~WhitespaceCharacter...~~"="~~WhitespaceCharacter...~~x:NumberString-> x]&, content];
    content=Map[ToExpression,content];
    
    If[content == {},
     content = default,
     content=First@Select[content,NumberQ]
    ]
    ];
    content
    ];
    
    
    BAMNumberParameterInFile[file_,parametername_]:=BAMNumberParameterInFile[file,parametername,Indeterminate]
    
    
    BAMNumberParametersInFile[file_,parametername_]:=Module[{dirname,content,x},
    
    If[FileType@file==File,content=StringSplit[Import[file,"String"],EndOfLine],content={}];
    
    If[content== {},
    Print["parameter file does not exist"],
    content=Flatten@Map[StringCases[#,parametername~~WhitespaceCharacter...~~"="~~WhitespaceCharacter...~~x__..-> x]&, content];
    
    content=Map[StringSplit[#,"#"]&,content];
    content=Map[First,content];
    content=StringReplace[content,Whitespace-> ","];
    content=Map[StringReplace[#, ","~~EndOfString->""]&,content];
    content=Map["{"<>#<>"}"&,content];
    content=Map[ToExpression,content][[1]]
    ];
    content
    ];
    
    
    BAMNumberParameters[directory_,parametername_]:=Module[{file,dirname,content,x},
    
    file=ParfileInDirectory[directory];
    file=ZippedOrUnzipped@file;
    Print["Identified parameter file ",file];
    If[FileType@file==File,content=StringSplit[Import[file,"String"],EndOfLine],content={}];
    
    Print["Read ",Length@content, " parameter file entries"];
    
    If[content== {},
    Print["parameter file does not exist"],
    content=Flatten@Map[StringCases[#,parametername~~WhitespaceCharacter...~~"="~~WhitespaceCharacter...~~x__..-> x]&, content];
    
    content=Map[StringSplit[#,"#"]&,content];
    content=Map[First,content];
    content=StringReplace[content,Whitespace-> ","];
    content=Map[StringReplace[#, ","~~EndOfString->""]&,content];
    content=Map["{"<>#<>"}"&,content];
    content=Map[ToExpression,content][[1]]
    ];
    content
    ];
    
    
    PSIDHashedNumberParameter[file_,parametername_]:=Module[{content,x},
    
    	content = PSIDReadHeader[file];
    
    	If[content == {},
    		Print["psid file does not exist"],
    		content = Flatten@Map[StringCases[#,"# "~~parametername~~WhitespaceCharacter...~~"="~~WhitespaceCharacter...~~x:NumberString-> x]&, content];
    		content = Map[ToExpression,content];
    
    		content = First@Select[content,NumberQ]
    	];
    	content
    ];
    
    
    PSIDNumberParameter[file_,parametername_]:=Module[{content,x},
    
    	content = PSIDReadHeader[file];
    
    	If[content == {},
    		Print["psid file does not exist"],
    		content = Flatten@Map[StringCases[#,StartOfLine~~parametername~~WhitespaceCharacter...~~"="~~WhitespaceCharacter...~~x___-> x]&, content];
    		content = Map[StringToNumber,content];
    
    		content = First@Select[content,NumberQ]
    	];
    	content
    ];
    
    
    PSIDReadHeader[file_]:=Module[{content,psidStream,str},
      (* only read the header of the psid file until the ID section starting with "data xx xx xx" *)
    
    	If[FileType@file == File,
    		psidStream = OpenRead[file];
    		content = {}; str = {};
    		While[True,
    			str = Read[psidStream, String]; 
    			If[StringMatchQ[str, RegularExpression["data[\\s\\d+]+"]], 
    				Break[]
    			];
    			AppendTo[content, str];
    		];
    		Close[psidStream];
    		,
    		content = {};
    	];
    	content
    ];
    
    
    PSIDReadData[file_]:=Module[{header,data,str,dims},
    	(* helper function to partition flat data *)
    	unflatten[e_,{d__?((IntegerQ[#]&&Positive[#])&)}]:= Fold[Partition,e,Take[{d},{-1,2,-1}]] /;(Length[e]===Times[d]);
    
    	data = {};
    	If[FileType@file == File,
    		psidStream = OpenRead[file];
    		header = {}; data={}; str = {};
    
    		(* read header *)
    		While[True,
    			str = Read[psidStream, String]; 
    			If[StringMatchQ[str, RegularExpression["data[\\s\\d+]+"]], Break[]];
    			AppendTo[header, str];
    		];
    
    	(* get dimensions (nz,ny,nx) *)
    	dims = Reverse@ToExpression@StringCases[str, RegularExpression["\\d+"]];
    
    	(* read data *)
    	data = ReadList[psidStream,Real];
    
    	(* partition into array according to dimensions *)
    	data = unflatten[data, dims];
    	Close[psidStream];
    	];
    
    	data
    ];
    
    
    PSID2Rules[filename_?StringQ]:=Module[{M1,M2,x1,y1,z1,x2,y2,z2,px1,py1,pz1
    ,px2,py2,pz2,s1x,s1y,s1z,s2x,s2y,s2z,mtot,sep,sepInM,prel,xrel,Madm,m1,m2},
    
    M1=PSIDHashedNumberParameter[filename,"M1"];
    M2=PSIDHashedNumberParameter[filename,"M2"];
    
    Madm=PSIDHashedNumberParameter[filename,"Madm"];
    
    m1=PSIDNumberParameter[filename,"bhmass1"];
    m2=PSIDNumberParameter[filename,"bhmass2"];
    
    x1=PSIDNumberParameter[filename,"bhx1"];
    y1=PSIDNumberParameter[filename,"bhy1"];
    z1=PSIDNumberParameter[filename,"bhz1"];
    
    x2=PSIDNumberParameter[filename,"bhx2"];
    y2=PSIDNumberParameter[filename,"bhy2"];
    z2=PSIDNumberParameter[filename,"bhz2"];
    
    px1=PSIDNumberParameter[filename,"bhpx1"];
    py1=PSIDNumberParameter[filename,"bhpy1"];
    pz1=PSIDNumberParameter[filename,"bhpz1"];
    
    px2=PSIDNumberParameter[filename,"bhpx2"];
    py2=PSIDNumberParameter[filename,"bhpy2"];
    pz2=PSIDNumberParameter[filename,"bhpz2"];
    
    s1x=PSIDNumberParameter[filename,"bhsx1"];
    s1y=PSIDNumberParameter[filename,"bhsy1"];
    s1z=PSIDNumberParameter[filename,"bhsz1"];
    
    s2x=PSIDNumberParameter[filename,"bhsx2"];
    s2y=PSIDNumberParameter[filename,"bhsy2"];
    s2z=PSIDNumberParameter[filename,"bhsz2"];
    
    mtot=M1+M2;
    
    sep={x1,y1,z1}-{x2,y2,z2};
    sep=Sqrt[sep.sep];
    sepInM = sep/mtot;
    
    prel={px1,py1,pz1}-{px2,py2,pz2};
    xrel={x1,y1,z1}-{x2,y2,z2};
    
    
    {"M1"->M1,"M2"->M2,"M"->mtot,
    "x1"->x1,"y1"->y1,"z1"->z1,"px1"->px1,"py1"->py1,"pz1"->pz1,"s1x"->s1x,"s1y"->s1y,"s1z"->s1z,
    "x2"->x2,"y2"->y2,"z2"->z2,"px2"->px2,"py2"->py2,"pz2"->pz2,"s2x"->s2x,"s2y"->s2y,"s2z"->s2z,
    "D[M]"->sepInM, "Madm"-> Madm, "BHMassParam1" -> m1, "BHMassParam2" -> m2 
    }
    ]
    
    
    NMovingLevels[parRules_]:="amr_lmax"-"amr_move_lcube"/.parRules
    
    
    CactusThornsAvailable[cactusdir_]:=Module[{arrDir,arrangements,thorns,allThorns},
    
    arrDir=FileNameJoin[{cactusdir,"arrangements"}];
    arrangements=FileNames["*",arrDir];
    arrangements=FileNameTake[#,-2]& /@Select[arrangements,FileType@#==Directory&];
    
    
    thorns=FileNames["*",arrDir,2];
    thorns=FileNameTake[#,-2]& /@ Select[thorns,FileType@#==Directory&];
    allThorns=Select[Complement[thorns,arrangements],Not@StringMatchQ[#,".*"]&];
    
    
    {"Arrangements"-> arrangements, "Thorns"-> allThorns}
    ];
    
    
    NormalizeCactusParfile[parfile_,cactusDir_,OptionsPattern[
    {"ThornListOutputFile"-> "","ParameterOutputFile"-> ""}]]:=Module[{temp,arrangements,allThorns,
    import,activeThorns,flatThornList,fullNameThorns,
    outThornList,outcontent,content,i,thorn,outparfile},
    
    outThornList = OptionValue["ThornListOutputFile"];
    outparfile = OptionValue["ParameterOutputFile"];
    
    temp=CactusThornsAvailable@cactusDir;
    arrangements = "Arrangements"/.temp;
    allThorns="Thorns"/.temp;
    
    import = Import[parfile,"String"];
    import=StringReplace[import,{
    "\n "-> " ",
    "\n"~~Whitespace...~~"\""->"\"",
     "\n"~~Whitespace...~~"="->"\""," "...~~"="->" ="
    }];
    
    
    import=Flatten@StringSplit[StringSplit[import,EndOfLine],"\n"];
    import=Map[StringReplace[#,WhitespaceCharacter...~~x__~~WhitespaceCharacter...~~"#"~~y___->x]&,import];
    import=Map[StringReplace[#,WhitespaceCharacter...~~x__->x]&,import];
    import=Map[StringReplace[#,WhitespaceCharacter...~~"#"~~x___->""]&,import];
    import=Map[StringReplace[#,"\t"->" "]&,import];
    import=Select[import,StringLength@#>1&];
    
    activeThorns=Select[import,StringMatchQ[#,"*ActiveThorns*"]&];
    activeThorns=Map[StringReplace[#,
      "ActiveThorns"~~WhitespaceCharacter...~~"="~~WhitespaceCharacter...~~x__->x]&,
       activeThorns];
    activeThorns=Map[StringSplit[#," "]&,activeThorns];
    
    
    flatThornList=Sort@Union@(StringReplace[#,"\""-> ""]&/@Flatten@activeThorns);
    flatThornList=Select[flatThornList,# != ""& ];
    
    fullNameThorns=Flatten@Table[
    Select[allThorns,StringMatchQ[#,"*/"<>flatThornList[[i]]]&],{i,1,Length@flatThornList}];
    
    outcontent=StringJoin[#<>"\n"&/@fullNameThorns];
    
    If[outThornList!= "",
      Print["Exporting thornlist to file ", outThornList];
      Export[outThornList,outcontent,"Text"];
    ];
    
    temp=Table[thorn=flatThornList[[i]]; {"\nActiveThorns = \""<>thorn<>"\"", 
      Sort@Union@Select[import,StringMatchQ[#,thorn<>"::*"]&]} ,{i,1,Length@flatThornList}];
    temp = Flatten@{"# thornlist created by NormalizeCactusParfile from file "<>parfile, temp};
    
    If[outThornList!= "",
      Print["Exporting normalized parameters to file ", outparfile];
      Export[outparfile,Flatten@temp,"Text"];
    ];
    
    
    {"AvailableArrangements"-> arrangements,"AvailableThorns"-> allThorns,
    "ActiveThorns"-> activeThorns,"FullNameThorns"-> fullNameThorns,"Parameters"-> import,
    "NormalizedParfile"-> temp}
    ];
    
    
    CompleteCactusParfile[parfile_,cactusDir_,OptionsPattern[
    {"ThornListOutputFile"-> "","ParameterOutputFile"-> ""}]]:=Module[{temp,arrangements,allThorns,
    import,activeThorns,flatThornList,fullNameThorns,
    outThornList,outcontent,content,i,thorn,outparfile},
    
    outThornList = OptionValue["ThornListOutputFile"];
    outparfile = OptionValue["ParameterOutputFile"];
    
    temp=CactusThornsAvailable@cactusDir;
    arrangements = "Arrangements"/.temp;
    allThorns="Thorns"/.temp;
    
    import = Import[parfile,"String"];
    import=StringReplace[import,{
    "\n "-> " ",
    "\n"~~Whitespace...~~"\""->"\"",
     "\n"~~Whitespace...~~"="->"\""," "...~~"="->" ="
    }];
    
    
    import=Flatten@StringSplit[StringSplit[import,EndOfLine],"\n"];
    import=Map[StringReplace[#,WhitespaceCharacter...~~x__~~WhitespaceCharacter...~~"#"~~y___->x]&,import];
    import=Map[StringReplace[#,WhitespaceCharacter...~~x__->x]&,import];
    import=Map[StringReplace[#,WhitespaceCharacter...~~"#"~~x___->""]&,import];
    import=Map[StringReplace[#,"\t"->" "]&,import];
    import=Select[import,StringLength@#>1&];
    
    (*  *)
    activeThorns=Select[import,StringMatchQ[#,"*ActiveThorns*"]&];
    activeThorns=Map[StringReplace[#,
      "ActiveThorns"~~WhitespaceCharacter...~~"="~~WhitespaceCharacter...~~x__->x]&,
       activeThorns];
    activeThorns=Map[StringSplit[#," "]&,activeThorns];
    
    
    flatThornList=Sort@Union@(StringReplace[#,"\""-> ""]&/@Flatten@activeThorns);
    flatThornList=Select[flatThornList,# != ""& ];
    
    fullNameThorns=Flatten@Table[
    Select[allThorns,StringMatchQ[#,"*/"<>flatThornList[[i]]]&],{i,1,Length@flatThornList}];
    
    outcontent=StringJoin[#<>"\n"&/@fullNameThorns];
    
    If[outThornList!= "",
      Print["Exporting thornlist to file ", outThornList];
      Export[outThornList,outcontent,"Text"];
    ];
    
    
    temp=Table[thorn=flatThornList[[i]]; {"\nActiveThorns = \""<>thorn<>"\"", 
      Sort@Union@Select[import,StringMatchQ[#,thorn<>"::*"]&]} ,{i,1,Length@flatThornList}];
    temp = Flatten@{"# thornlist created by NormalizeCactusParfile from file "<>parfile, temp};
    
    If[outThornList!= "",
      Print["Exporting normalized parameters to file ", outparfile];
      Export[outparfile,Flatten@temp,"Text"];
    ];
    
    
    {"AvailableArrangements"-> arrangements,"AvailableThorns"-> allThorns,
    "ActiveThorns"-> activeThorns,"FullNameThorns"-> fullNameThorns,"Parameters"-> import,
    "NormalizedParfile"-> temp}
    ];
    
    
    SetParfileEntryValue[text_,key_,value_]:=
    StringReplace[text,key~~ws1:WhitespaceCharacter... ~~\!\(\*
    TagBox[
    StyleBox["\"\<=\>\"",
    ShowSpecialCharacters->False,
    ShowStringCharacters->True,
    NumberMarks->True],
    FullForm]\)~~ws2:WhitespaceCharacter...~~x__~~EndOfLine:> key<>ws1<>"="<>ws2<>ToString@value]
    
    
    SetParfileValue[text_,key_,value_]:=Module[{entry,temp,len,new,x,out},
    
    entry=StringCases[text,Shortest[key~~WhitespaceCharacter...~~"="~~WhitespaceCharacter...~~x___~~EndOfLine]];
    (* Print@entry; *)
    
    out = text;
    
    If[Length@entry > 0,
      entry=Last@entry;
    
      new = SetParfileEntryValue[entry,key,value];
    
      If[StringQ@entry,
         out = StringReplace[text,entry :>  new];
      ];
    ];
    
    out
    ];
    
    
    SetParfileVectorValue[text_,key_,component_,value_]:=Module[{entry,temp,len,new,s,x,out},
    s = "["~~WhitespaceCharacter...~~ ToString@component ~~WhitespaceCharacter...~~ "]"; 
    
    entry = StringCases[text,Shortest[key~~WhitespaceCharacter...~~ s ~~WhitespaceCharacter...~~"= "~~x__~~EndOfLine]];
    (* Print[entry]; *)
    
    out = text;
    
    If[Length@entry > 0,
    entry=Last@entry;
    
    s = "[" <> ToString@component <> "]"; 
    new = key <> s <> " = " <> ToString@value;
    
    If[StringQ@entry,
      out = StringReplace[text,entry ->   new];
    ];
    ];
    
    out
    ];
    
    
    CreateDataReduceDirectory[dirname_]:=CreateDirectory[dirname<>"/DataReduce"]
    
    
    CreateDataReduceDirectory[reduceroot_,dirname_]:=CreateDirectory[reduceroot<>"/"<>LastInPath@dirname<>"/DataReduce"]
    
    
    LocateMode[modeDir_,lmode_,mmode_]:=Module[{searchString,modeFiles,file,radii,outerradius,level},
    
    searchString="hmod.r*.l*.l"<>ToString@lmode<>".m"<>ToString@mmode<>"*";
    Print["Searching for files matching ", searchString];
    modeFiles=FileNames[searchString,modeDir];
    
    Print["Found mode files: ", modeFiles];
    
    level=Union@Map[IntegerPart,Map[levelFun,modeFiles]];
    Print["Available refinement levels: ", level];
    
    radii=Sort@TakeColumn[level,1];
    level=Sort@TakeColumn[level,2];
    
    outerradius={ToString@Last@radii,ToString@First@level};
    
    file="hmod.r"<>outerradius[[1]]<>".l"<>outerradius[[2]]<>".l"<>ToString@lmode<>".m"<>ToString@mmode;
    file=ZippedOrUnzipped[file];
    If[FileType@file!= File,Print["cannot find hmod file"];file="";];
    
    file
    ];
    
    
    LocateModes[modeDir_,lmode_,mmode_]:=Module[{searchString,modeFiles,file,radii,outerradius,level,levelFun},
    
    searchString="hmod.r*.l*.l"<>ToString@lmode<>".m"<>ToString@mmode<>"*";
    Print["Searching for files matching ", searchString];
    FileNames[searchString,modeDir]
    ];
    
    
    CopyL2mode[configStr_,modesdir_,reducedir_]:=Module[{command,modeFile,source,fname},
    modeFile=LocateMode[modesdir,2,2];
    Print["Called CopyL2mode with mode file ", modeFile];
    
    source=modesdir<>"/"<>modeFile;
    source=ZippedOrUnzipped@source;
    If[FileType@source==File,
    fname=Last@StringSplit[source,"/"];
    command="cp "<>source<>"  "<>reducedir<>"/"<>configStr<>"_"<>fname;
    Print["Running command ",command];
    Run[Evaluate@command];,
    Print["ERROR: Could not find hmod modes file"];
    ];
    
    source=ZippedOrUnzipped@StringReplace[source,"hmod"-> "psi3col"];
    If[FileType@source==File,
    fname=Last@StringSplit[source,"/"];
    command="cp "<>source<>"  "<>reducedir<>"/"<>configStr<>"_"<>fname;
    Print["Running command ",command];
    Run[Evaluate@command];,
    Print["ERROR: Could not find psi4 modes file"];
    ];
    ];
    
    
    CopyL2modes[configStr_,modesdir_,reducedir_]:=Module[{command,modeFiles,modeFile,source,fname,i},
    modeFiles=LocateModes[modesdir,2,2];
    Print["Called CopyL2modes with mode files ", modeFiles];
    
    Do[
    source=modeFiles[[i]];
    
    If[FileType@source==File,
    fname=Last@StringSplit[source,"/"];
    command="cp "<>source<>"  "<>reducedir<>"/"<>configStr<>"_"<>fname;
    Print["CopyL2modes running command ",command];
    Run[Evaluate@command];,
    Print["ERROR: Could not find hmod modes file"];
    ];
    source=ZippedOrUnzipped@StringReplace[source,"hmod"-> "psi3col"];
    If[FileType@source==File,
    fname=Last@StringSplit[source,"/"];
    command="cp "<>source<>"  "<>reducedir<>"/"<>configStr<>"_"<>fname;
    Print["CopyL2modes running command ",command];
    Run[Evaluate@command];,
    Print["ERROR: Could not find psi4 modes file"];
    ];,{i,1,Length@modeFiles}];
    ];
    
    
    CopyModes[configStr_,modesdir_,reducedir_,Lmode_,Mmode_]:=Module[{command,modeFiles,modeFile,source,fname,i},
    modeFiles=LocateModes[modesdir,Lmode,Mmode];
    Print["Called CopyModes with mode files ", modeFiles];
    
    Do[
    source=modeFiles[[i]];
    
    If[FileType@source==File,
    fname=Last@StringSplit[source,"/"];
    command="cp "<>source<>"  "<>reducedir<>"/"<>configStr<>"_"<>fname;
    Print["CopyModes running command ",command];
    Run[Evaluate@command];,
    Print["ERROR: Could not find hmod modes file"];
    ];
    source=ZippedOrUnzipped@StringReplace[source,"hmod"-> "psi3col"];
    If[FileType@source==File,
    fname=Last@StringSplit[source,"/"];
    command="cp "<>source<>"  "<>reducedir<>"/"<>configStr<>"_"<>fname;
    Print["CopyModes running command ",command];
    Run[Evaluate@command];,
    Print["ERROR: Could not find psi4 modes file"];
    ];,{i,1,Length@modeFiles}];
    ];
    
    
    CopyFiles[configStr_,modesdir_,reducedir_,patterns_]:=Module[{command,files,modeFile,source,fname,i,prefix},
    
    files=Map[FileNames[#,modesdir]&,patterns];
    
    (* Print["CopyFiles found files ", files]; *)
    
    files=Flatten@files;
    
    If[configStr == "", 
    prefix = "",
    prefix = configStr<>"_";
    ];
    
    Do[
    source=files[[i]];
    
    If[FileType@source==File || FileType@source==Directory,
    fname=Last@StringSplit[source,"/"];
    command="cp -r "<>source<>"  "<>reducedir<>"/"<>prefix<>fname;
    (* Print["CopyFiles running command ",command];*)
    Run[Evaluate@command];,
    Print["ERROR: CopyFiles could not find file"];
    ];,{i,1,Length@files}];
    
    
    ];
    
    
    FormatPunctureData[mp_,string_,M_]:=Module[{xl,yl,zl,vxl,vyl,vzl,rl,wl,\[Omega]l,times},xl=TakeColumn[Position/.mp,1];
    yl=TakeColumn[Position/.mp,2];
    zl=TakeColumn[Position/.mp,3];
    vxl=TakeColumn[Speed/.mp,1];
    vyl=TakeColumn[Speed/.mp,2];
    vzl=TakeColumn[Speed/.mp,3];
    times=Times/.mp;
    times=times/M;
    xl=xl/M;
    yl=yl/M;
    zl=zl/M;
    Clear[Evaluate["times"<>string],Evaluate["xlist"<>string],Evaluate["ylist"<>string],Evaluate["zlist"<>string],Evaluate["vxlist"<>string],Evaluate["vylist"<>string],Evaluate["vzlist"<>string],Evaluate["xf"<>string],Evaluate["yf"<>string],Evaluate["zf"<>string],Evaluate["vxf"<>string],Evaluate["vyf"<>string],Evaluate["vzf"<>string],Evaluate["rlist"<>string],Evaluate["rf"<>string],Evaluate["wlist"<>string],Evaluate["wf"<>string],Evaluate["\[Omega]list"<>string],Evaluate["\[Omega]f"<>string]];
    Print["Defining the quantities: \n",Evaluate["times"<>string]," , ",Evaluate["[xyz]list"<>string]," , ",Evaluate["v[xyz]list"<>string]," , ",Evaluate["v[xyz]f"<>string]," , ",Evaluate["rlist"<>string]," , ",Evaluate["rf"<>string]," , ",Evaluate["wlist"<>string]," , ",Evaluate["wf"<>string]];
    Evaluate@ToExpression["times"<>string]=times;
    Evaluate@ToExpression["xf"<>string]=Interpolation[Union@CombineColumns[times,xl]];
    Evaluate@ToExpression["yf"<>string]=Interpolation[Union@CombineColumns[times,yl]];
    Evaluate@ToExpression["zf"<>string]=Interpolation[Union@CombineColumns[times,zl]];
    Evaluate@ToExpression["vxf"<>string]=Interpolation[Union@CombineColumns[times,vxl]];
    Evaluate@ToExpression["vyf"<>string]=Interpolation[Union@CombineColumns[times,vyl]];
    Evaluate@ToExpression["vzf"<>string]=Interpolation[Union@CombineColumns[times,vzl]];
    Evaluate@ToExpression["xlist"<>string]=xl;
    Evaluate@ToExpression["ylist"<>string]=yl;
    Evaluate@ToExpression["zlist"<>string]=zl;
    Evaluate@ToExpression["vxlist"<>string]=vxl;
    Evaluate@ToExpression["vylist"<>string]=vyl;
    Evaluate@ToExpression["vzlist"<>string]=vzl;
    rl=Sqrt[xl^2+yl^2+zl^2];
    Evaluate@ToExpression["rf"<>string]=Interpolation[Union@CombineColumns[times,rl]];
    wl=Sqrt[vxl^2+vyl^2+vzl^2]/Sqrt[xl^2+yl^2+zl^2];
    Evaluate@ToExpression["wf"<>string]=Interpolation[Union@CombineColumns[times,wl]];
    \[Omega]l=(vyl xl-vxl yl)/(xl^2+yl^2+zl^2);
    Evaluate@ToExpression["\[Omega]f"<>string]=Interpolation[Union@CombineColumns[times,\[Omega]l]];
    Evaluate@ToExpression["rlist"<>string]=rl;
    Evaluate@ToExpression["wlist"<>string]=wl;
    Evaluate@ToExpression["\[Omega]list"<>string]=\[Omega]l;];
    
    
    
    SafeFormatPunctureData2[mp1_,mp2_,string_,m1_,m2_]:=Module[{M,xl1,yl1,zl1,xl2,yl2,zl2,xlc,ylc,zlc,dxl,dyl,dzl,
    vxl1,vyl1,vzl1,vxl2,vyl2,vzl2,vxlc,vylc,vzlc,rl,wl,\[Omega]l,vrxl,vryl,vrzl,wxl,wyl,wzl,vrl,vtl,times,timesSorted,timesSorted1,timesSorted2,pos},M=m1+m2;
    xl1=TakeColumn[Position/.mp1,1];
    yl1=TakeColumn[Position/.mp1,2];
    zl1=TakeColumn[Position/.mp1,3];
    vxl1=TakeColumn[Speed/.mp1,1];
    vyl1=TakeColumn[Speed/.mp1,2];
    vzl1=TakeColumn[Speed/.mp1,3];
    xl2=TakeColumn[Position/.mp2,1];
    yl2=TakeColumn[Position/.mp2,2];
    zl2=TakeColumn[Position/.mp2,3];
    vxl2=TakeColumn[Speed/.mp2,1];
    vyl2=TakeColumn[Speed/.mp2,2];
    vzl2=TakeColumn[Speed/.mp2,3];
    times=Times/.mp1;
    timesSorted1=Union@times;
    
    times=Times/.mp2;
    timesSorted2=Union@times;
    
    If[Length@timesSorted2 >= Length@timesSorted1, timesSorted = timesSorted1,timesSorted = timesSorted2]
    
    Clear[timesSorted1,timesSorted2];
    
    pos=Table[Position[times,timesSorted[[i]]],{i,1,Length@timesSorted}];
    pos=Flatten@Map[First,pos];
    times=Table[times[[pos[[i]]]],{i,1,Length@pos}];
    xl1=Table[xl1[[pos[[i]]]],{i,1,Length@pos}];
    yl1=Table[yl1[[pos[[i]]]],{i,1,Length@pos}];
    zl1=Table[zl1[[pos[[i]]]],{i,1,Length@pos}];
    xl2=Table[xl2[[pos[[i]]]],{i,1,Length@pos}];
    yl2=Table[yl2[[pos[[i]]]],{i,1,Length@pos}];
    zl2=Table[zl2[[pos[[i]]]],{i,1,Length@pos}];
    vxl1=Table[vxl1[[pos[[i]]]],{i,1,Length@pos}];
    vyl1=Table[vyl1[[pos[[i]]]],{i,1,Length@pos}];
    vzl1=Table[vzl1[[pos[[i]]]],{i,1,Length@pos}];
    vxl2=Table[vxl2[[pos[[i]]]],{i,1,Length@pos}];
    vyl2=Table[vyl2[[pos[[i]]]],{i,1,Length@pos}];
    vzl2=Table[vzl2[[pos[[i]]]],{i,1,Length@pos}];
    (*scale to units of M*)times=times/M;
    xl1=xl1/M;
    yl1=yl1/M;
    zl1=zl1/M;
    xl2=xl2/M;
    yl2=yl2/M;
    zl2=zl2/M;
    (*center-of-mass vector*)xlc=(m1 xl1+m2 xl2)/M;
    ylc=(m1 yl1+m2 yl2)/M;
    zlc=(m1 zl1+m2 zl2)/M;
    (*relative coordinates*)xl1=xl1-xlc;
    yl1=yl1-ylc;
    zl1=zl1-zlc;
    xl2=xl2-xlc;
    yl2=yl2-ylc;
    zl2=zl2-zlc;
    Clear[Evaluate["times"<>string],Evaluate["xlista"<>string],Evaluate["ylista"<>string],Evaluate["zlista"<>string],Evaluate["vxlista"<>string],Evaluate["vylista"<>string],Evaluate["vzlista"<>string],Evaluate["xfa"<>string],Evaluate["yfa"<>string],Evaluate["zfa"<>string],Evaluate["vxfa"<>string],Evaluate["vyfa"<>string],Evaluate["vzfa"<>string],Evaluate["rlist"<>string],Evaluate["rf"<>string],Evaluate["wlist"<>string],Evaluate["wf"<>string],Evaluate["\[Omega]list"<>string],Evaluate["\[Omega]f"<>string],Evaluate["xlistb"<>string],Evaluate["ylistb"<>string],Evaluate["zlistb"<>string],Evaluate["vxlistb"<>string],Evaluate["vylistb"<>string],Evaluate["vzlistb"<>string],Evaluate["xfb"<>string],Evaluate["yfb"<>string],Evaluate["zfb"<>string],Evaluate["vxfb"<>string],Evaluate["vyfb"<>string],Evaluate["vzfb"<>string]];
    Print["Defining the quantities: \n",Evaluate["times"<>string]," , ",Evaluate["[xyz]list"<>string]," , ",Evaluate["v[xyz]list"<>string]," , ",Evaluate["v[xyz]f"<>string]," , ",Evaluate["rlist"<>string]," , ",Evaluate["rf"<>string]," , ",Evaluate["wlist"<>string]," , ",Evaluate["wf"<>string]];
    Evaluate@ToExpression["times"<>string]=times;
    Evaluate@ToExpression["xfa"<>string]=Interpolation[CombineColumns[times,xl1]];
    Evaluate@ToExpression["yfa"<>string]=Interpolation[CombineColumns[times,yl1]];
    Evaluate@ToExpression["zfa"<>string]=Interpolation[CombineColumns[times,zl1]];
    Evaluate@ToExpression["vxfa"<>string]=Interpolation[CombineColumns[times,vxl1]];
    Evaluate@ToExpression["vyfa"<>string]=Interpolation[CombineColumns[times,vyl1]];
    Evaluate@ToExpression["vzfa"<>string]=Interpolation[CombineColumns[times,vzl1]];
    Evaluate@ToExpression["xlista"<>string]=xl1;
    Evaluate@ToExpression["ylista"<>string]=yl1;
    Evaluate@ToExpression["zlista"<>string]=zl1;
    Evaluate@ToExpression["vxlista"<>string]=vxl1;
    Evaluate@ToExpression["vylista"<>string]=vyl1;
    Evaluate@ToExpression["vzlista"<>string]=vzl1;
    Evaluate@ToExpression["xfb"<>string]=Interpolation[CombineColumns[times,xl2]];
    Evaluate@ToExpression["yfb"<>string]=Interpolation[CombineColumns[times,yl2]];
    Evaluate@ToExpression["zfb"<>string]=Interpolation[CombineColumns[times,zl2]];
    Evaluate@ToExpression["vxfb"<>string]=Interpolation[CombineColumns[times,vxl2]];
    Evaluate@ToExpression["vyfb"<>string]=Interpolation[CombineColumns[times,vyl2]];
    Evaluate@ToExpression["vzfb"<>string]=Interpolation[CombineColumns[times,vzl2]];
    Evaluate@ToExpression["xlistb"<>string]=xl2;
    Evaluate@ToExpression["ylistb"<>string]=yl2;
    Evaluate@ToExpression["zlistb"<>string]=zl2;
    Evaluate@ToExpression["vxlistb"<>string]=vxl2;
    Evaluate@ToExpression["vylistb"<>string]=vyl2;
    Evaluate@ToExpression["vzlistb"<>string]=vzl2;
    rl=Sqrt[(xl1-xl2)^2+(yl1-yl2)^2+(zl1-zl2)^2];
    rl=Table[Max[rl[[i]],$MachineEpsilon],{i,1,Length@rl}];
    Evaluate@ToExpression["rf"<>string]=Interpolation[CombineColumns[times,rl]];
    wl=Sqrt[(vxl1-vxl2)^2+(vyl1-vyl2)^2+(vzl1-vzl2)^2]/rl;
    Evaluate@ToExpression["wf"<>string]=Interpolation[CombineColumns[times,wl]];
    dxl=(xl2-xl1)/rl;
    dyl=(yl2-yl1)/rl;
    dzl=(zl2-zl1)/rl;
    (*w=drl\[Cross]\[CapitalDelta]yl*)wxl=dyl (vzl2-vzl1)-dzl (vyl2-vyl1);
    wyl=dzl (vxl2-vxl1)-dxl (vzl2-vzl1);
    wzl=dxl (vyl2-vyl1)-dyl (vxl2-vxl1);
    (*\[Omega]=(|drl\[Cross]\[CapitalDelta]yl|)/rl*)\[Omega]l=Sqrt[wxl^2+wyl^2+wzl^2]/rl;
    Evaluate@ToExpression["\[Omega]f"<>string]=Interpolation[CombineColumns[times,\[Omega]l]];
    Evaluate@ToExpression["rlist"<>string]=rl;
    Evaluate@ToExpression["wlist"<>string]=wl;
    Evaluate@ToExpression["\[Omega]list"<>string]=\[Omega]l;];
    
    
    SafeFormatPunctureDataCactus[mp_,string_,m1_,m2_]:=Module[{M,xl1,yl1,zl1,xl2,yl2,zl2,xlc,ylc,zlc,dxl,dyl,dzl,vxl1,vyl1,vzl1,vxl2,vyl2,vzl2,vxlc,vylc,vzlc,rl,wl,\[Omega]l,vrxl,vryl,vrzl,wxl,wyl,wzl,vrl,vtl,times,timesSorted,pos},M=m1+m2;
    xl1=TakeColumn[mp,23];
    yl1=TakeColumn[mp,33];
    zl1=TakeColumn[mp,43];
    xl2=TakeColumn[mp,24];
    yl2=TakeColumn[mp,34];
    zl2=TakeColumn[mp,44];
    times=TakeColumn[mp,9];
    timesSorted=Union@times;
    pos=Table[Position[times,timesSorted[[i]]],{i,1,Length@timesSorted}];
    pos=Flatten@Map[First,pos];
    times=Table[times[[pos[[i]]]],{i,1,Length@pos}];
    xl1=Table[xl1[[pos[[i]]]],{i,1,Length@pos}];
    yl1=Table[yl1[[pos[[i]]]],{i,1,Length@pos}];
    zl1=Table[zl1[[pos[[i]]]],{i,1,Length@pos}];
    xl2=Table[xl2[[pos[[i]]]],{i,1,Length@pos}];
    yl2=Table[yl2[[pos[[i]]]],{i,1,Length@pos}];
    zl2=Table[zl2[[pos[[i]]]],{i,1,Length@pos}];
    times=times/M;
    vx1=D[Interpolation[CombineColumns[times,xl1]][t],t]/.ff1_[tt_]->ff1;
    vy1=D[Interpolation[CombineColumns[times,yl1]][t],t]/.ff2_[tt_]->ff2;
    vz1=D[Interpolation[CombineColumns[times,zl1]][t],t]/.ff3_[tt_]->ff3;
    vx2=D[Interpolation[CombineColumns[times,xl2]][t],t]/.ff4_[tt_]->ff4;
    vy2=D[Interpolation[CombineColumns[times,yl2]][t],t]/.ff5_[tt_]->ff5;
    vz2=D[Interpolation[CombineColumns[times,zl2]][t],t]/.ff6_[tt_]->ff6;
    vxl1=Map[vx1,times];
    vyl1=Map[vy1,times];
    vzl1=Map[vz1,times];
    vxl2=Map[vx2,times];
    vyl2=Map[vy2,times];
    vzl2=Map[vz2,times];
    xl1=xl1/M;
    yl1=yl1/M;
    zl1=zl1/M;
    xl2=xl2/M;
    yl2=yl2/M;
    zl2=zl2/M;
    xlc=m1 xl1+m2 xl2;
    ylc=m1 yl1+m2 yl2;
    zlc=m1 zl1+m2 zl2;
    xl1=xl1-xlc;
    yl1=yl1-ylc;
    zl1=zl1-zlc;
    xl2=xl2-xlc;
    yl2=yl2-ylc;
    zl2=zl2-zlc;
    Clear[Evaluate["times"<>string],Evaluate["xlista"<>string],Evaluate["ylista"<>string],Evaluate["zlista"<>string],Evaluate["vxlista"<>string],Evaluate["vylista"<>string],Evaluate["vzlista"<>string],Evaluate["xfa"<>string],Evaluate["yfa"<>string],Evaluate["zfa"<>string],Evaluate["vxfa"<>string],Evaluate["vyfa"<>string],Evaluate["vzfa"<>string],Evaluate["rlist"<>string],Evaluate["rf"<>string],Evaluate["wlist"<>string],Evaluate["wf"<>string],Evaluate["\[Omega]list"<>string],Evaluate["\[Omega]f"<>string],Evaluate["xlistb"<>string],Evaluate["ylistb"<>string],Evaluate["zlistb"<>string],Evaluate["vxlistb"<>string],Evaluate["vylistb"<>string],Evaluate["vzlistb"<>string],Evaluate["xfb"<>string],Evaluate["yfb"<>string],Evaluate["zfb"<>string],Evaluate["vxfb"<>string],Evaluate["vyfb"<>string],Evaluate["vzfb"<>string]];
    Print["Defining the quantities: \n",Evaluate["times"<>string]," , ",Evaluate["[xyz]list"<>string]," , ",Evaluate["rlist"<>string]," , ",Evaluate["rf"<>string]," , ",Evaluate["wlist"<>string]," , ",Evaluate["wf"<>string]];
    Evaluate@ToExpression["times"<>string]=times;
    Evaluate@ToExpression["xfa"<>string]=Interpolation[CombineColumns[times,xl1]];
    Evaluate@ToExpression["yfa"<>string]=Interpolation[CombineColumns[times,yl1]];
    Evaluate@ToExpression["zfa"<>string]=Interpolation[CombineColumns[times,zl1]];
    Evaluate@ToExpression["xfb"<>string]=Interpolation[CombineColumns[times,xl2]];
    Evaluate@ToExpression["yfb"<>string]=Interpolation[CombineColumns[times,yl2]];
    Evaluate@ToExpression["zfb"<>string]=Interpolation[CombineColumns[times,zl2]];
    Evaluate@ToExpression["xlista"<>string]=xl1;
    Evaluate@ToExpression["ylista"<>string]=yl1;
    Evaluate@ToExpression["zlista"<>string]=zl1;
    Evaluate@ToExpression["xlistb"<>string]=xl2;
    Evaluate@ToExpression["ylistb"<>string]=yl2;
    Evaluate@ToExpression["zlistb"<>string]=zl2;
    rl=Sqrt[(xl1-xl2)^2+(yl1-yl2)^2+(zl1-zl2)^2];
    rl=Table[Max[rl[[i]],$MachineEpsilon],{i,1,Length@rl}];
    Evaluate@ToExpression["rf"<>string]=Interpolation[CombineColumns[times,rl]];
    wl=Sqrt[(vxl1-vxl2)^2+(vyl1-vyl2)^2+(vzl1-vzl2)^2]/rl;
    Evaluate@ToExpression["wf"<>string]=Interpolation[CombineColumns[times,wl]];
    dxl=(xl2-xl1)/rl;
    dyl=(yl2-yl1)/rl;
    dzl=(zl2-zl1)/rl;
    wxl=dyl (vzl2-vzl1)-dzl(vyl2-vyl1);
    wyl=dzl (vxl2-vxl1)-dxl(vzl2-vzl1);
    wzl=dxl (vyl2-vyl1)-dyl(vxl2-vxl1);
    \[Omega]l=Sqrt[wxl^2+wyl^2+wzl^2]/rl;
    Evaluate@ToExpression["\[Omega]f"<>string]=Interpolation[CombineColumns[times,\[Omega]l]];
    Evaluate@ToExpression["rlist"<>string]=rl;
    Evaluate@ToExpression["wlist"<>string]=wl;
    Evaluate@ToExpression["\[Omega]list"<>string]=\[Omega]l;];
    
    
    BAMModesFilesTo3Col[modesFile_,OptionsPattern[{"DeleteSourceFiles" -> False}]]:=Module[{rootString,path,noPath,reData,imData,rad,lev,modeString,timesRe,timesIm,
    col3Data,file,col3data,deleteSourceFiles,lm,commonTimes},
    deleteSourceFiles = OptionValue["DeleteSourceFiles"];
    
    path = FileNameTake[modesFile,FileNameDepth[modesFile]-1];
    noPath=FileNameTake[modesFile,-1];
    rootString = StringCases[modesFile,{"r","i"}~~"psi4mode"~~postfix__:>  {"rpsi4mode"<>postfix,"ipsi4mode"<>postfix}];
    
    If[rootString != {},
    
      {rad,lev}=First@StringCases[modesFile,"r"~~rad:NumberString~~".t"~~lev:NumberString~~___ :>  {rad,lev}];
    
      modeString = First@StringCases[modesFile,{"r","i"}~~"psi4mode"~~mode__~~"_r"~~postfix__:> mode];
      lm = Switch[modeString,
       "20",  {2, 0},
       "2m1", {2,-1},
       "21",  {2, 1},
       "2m2", {2,-2},
       "22",  {2, 2}
      ];
    
      file = FileNameJoin[{path,rootString[[1,1]]}];
      If[FileType@file !=  File,
        Print["BAMModesFilesTo3Col:: File not found: ", file]; Return[];,
    
        reData = TimeUnion@NRFiles`ReadColData[file,2,1];
        If[deleteSourceFiles,DeleteFile[file];];
        timesRe = NRLists`TakeColumn[reData,1];
    
        If[lm[[2]]!=0,
          file = FileNameJoin[{path,rootString[[1,2]]}];
    
          If[FileType@file !=  File,
            Print["BAMModesFilesTo3Col:: File not found: ", file];
            imData = {};,
            timesIm = {};
            imData = TimeUnion@NRFiles`ReadColData[file,2,1];
             If[deleteSourceFiles,DeleteFile[file];];
            timesIm = NRLists`TakeColumn[imData,1];
           ];,
           timesIm = timesRe;
           imData = Transpose@{timesIm,0*timesIm};
        ];
    
       If[timesIm==timesRe,
         col3data = Transpose@{timesRe, NRLists`TakeColumn[reData,2], NRLists`TakeColumn[imData,2]};
          file = path<>"/psi3col_bam_r"<>ToString@rad<>".l"<>ToString@lev<>".l"<>ToString@lm[[1]]<>".m"<>ToString@lm[[2]]<>".gz";
          Print["Exporting file ", file];
          Export[file,col3data,{"GZIP","Table"}],
          Print["BAMModesFilesTo3Col:: fixing inconsistent lengths of real and imaginary part"];
          Print["Length@Re@wave = ", Length@timesRe, ", ", "Length@Im@wave = ", Length@timesIm];
          commonTimes = Intersection[timesIm,timesRe];
          reData = Select[reData,MemberQ[commonTimes,#[[1]]]&];
          imData = Select[imData,MemberQ[commonTimes,#[[1]]]&];
    
          col3data = Transpose@{commonTimes, NRLists`TakeColumn[reData,2], NRLists`TakeColumn[imData,2]};
          file = path<>"/psi3col_bam_r"<>ToString@rad<>".l"<>ToString@lev<>".l"<>ToString@lm[[1]]<>".m"<>ToString@lm[[2]]<>".gz";
          Print["Exporting file ", file];
         Export[file,col3data,{"GZIP","Table"}]
       ];
       ];
    ];
    ];
    
    
    BAMTrajectoryFileTo4Col[inFile_,OptionsPattern[{"DeleteSourceFiles" -> False}]]:=Module[{path,noPath,data,
    lev,bh,modeString,timesRe,times,
    file,col4data,deleteSourceFiles},
    
    deleteSourceFiles = OptionValue["DeleteSourceFiles"];
    
    If[FileType@inFile == File,
    
    path = path = FileNameTake[inFile,FileNameDepth[inFile]-1];
    noPath = FileNameTake[inFile,-1];
    
    {bh,lev} = First@StringCases[noPath,
         "moving_puncture_integrate"~~bh:NumberString~~".txyz"~~lev:NumberString~~X___ :>  {bh,lev}];
    
    lev = StringReplace[lev, "."->""];
    
    data = Import@inFile;
    If[deleteSourceFiles,DeleteFile[inFile];];
    col4data = TimeUnion@NRLists`TakeColumn[data,{7,1,2,3}];
    
    file = path<>"/traj_"<>ToString@bh<>".l"<> ToString@lev<>".gz";
    Print["Exporting file ", file];
    Export[file,col4data,{"GZIP","Table"}];
    ];
    ];
    
    
    BAMHorizonFilesToNRARFiles[inFile_,OptionsPattern[{"DeleteSourceFiles" -> False}]]:=Module[{path,noPath,data,lev,bh,modeString,timesRe,times,
    file,col4data,deleteSourceFiles,
    hmass,hspin,totalSpin,htraj,
    t,m},
    
    deleteSourceFiles = OptionValue["DeleteSourceFiles"];
    
    If[FileType@inFile == File,
    
    path = path = FileNameTake[inFile,FileNameDepth[inFile]-1];
    noPath = FileNameTake[inFile,-1];
    Print["BAMHorizonFilesToNRARFiles: filename: ", noPath];
    
    bh = First@StringCases[noPath,"horizon_"~~bh:NumberString :>  bh];
    Print["BAMHorizonFilesToNRARFiles:  bh: ", bh];
    bh = ToExpression@bh + 1;
    
    data = TimeUnion@ReadColData[inFile,9,1];
    If[deleteSourceFiles, DeleteFile[inFile];];
    
    times = TakeColumn[data,1];
    hmass = TakeColumn[data,5];
    hspin = TakeColumn[data,{1,6,7,8}];
    
    
    totalSpin = TakeColumn[data,9];
    
    hmass = Sqrt[hmass^2+1/4 totalSpin^2/hmass^2];
    hmass = CombineColumns[times,hmass];
    
    
    htraj = TakeColumn[data,{1,2,3,4}];
    
    
    file = path<>"/hmass_"<>ToString@bh<>".gz";
    Print["Exporting file ", file];
    Export[file,hmass,{"GZIP","Table"}];
    
    file = path<>"/htraj_"<>ToString@bh<>".gz";
    Print["Exporting file ", file];
    Export[file,htraj,{"GZIP","Table"}];
    
    file = path<>"/hspin_"<>ToString@bh<>".gz";
    Print["Exporting file ", file];
    Export[file,hspin,{"GZIP","Table"}];
    ];
    ];
    
    
    WaveExtractionRadii[rootDir_,modesDir_]:=Module[{admFiles,adm,i,admradii,outerradius,level,massScale=1,len},
    
    admFiles=FileNames["hmod.r*.l*.l*.m*",modesDir];
    
    If[admFiles == {},
    Print["ERROR: could not find any hmod data in directory ", LastInPath@modesDir];
    Return[];
    ];
    
    level=Union@Map[IntegerPart,Map[levelFun,admFiles]];
    Print["WaveExtractionRadii: Available refinement levels: ", level];
    
    outerradius=TakeColumn[level,1];
    Print["Wave extraction indices: ", outerradius];
    outerradius=Max@Map[ToExpression,outerradius];
    Print["Maximal wave extraction radius: ", outerradius];
    
    admradii=BAMNumberParameters[rootDir,"invariants_modes_r"];
    Print["Wave extraction radii: ", admradii];
    
    admradii
    ];
    
    
    ADMReduce[rootDir_,reduceDir_]:=Module[{admFiles,adm,i,admradii,outerradius,level,levelFun,massScale=1,len},
    
    admFiles=FileNames["*ADM_mass_r*",rootDir];
    
    Print["ADM data files: ", admFiles];
    
    If[admFiles == {},
    Print["ERROR: ADMReduce could not find any ADM data in directory ", LastInPath@rootDir];Return[];,
    Print["Starting ADM analysis"];
    ];
    
    levelFun[str_]:=ToExpression@First@StringCases[str,"ADM_mass_r"~~NumberString..~~".t"~~x:NumberString-> x];
    
    level=Map[IntegerPart,Map[levelFun,admFiles]];
    Print["Available refinement levels: ", level];
    
    outerradius=Map[Last@StringSplit[#,"/"]&,admFiles];
    outerradius=Union@Map[StringSplit[#,"_"][[3]]&,outerradius];
    outerradius=Union@Map[StringSplit[#,"."][[1]]&,outerradius];
    outerradius=Union@StringReplace[outerradius,"r"-> ""];
    
    Print["ADM-style integrals extraction indices: ", outerradius];
    outerradius=Max@Map[ToExpression,outerradius];
    Print["Maximal ADM-stlye integral extraction radius: ", outerradius];
    
    admradii=BAMNumberParameters[rootDir,"ADM_mass_r"];
    Print["ADM-style integrals extraction radii: ", admradii];
    
    len=Length@admradii;
    level=First@Union@Flatten@Select[myGather[level],Length@#==len&];
    Print["All extraction radii seem present at level ", level];
    
    
    admFiles=FileNames["*ADM_mass_r*.t"<>ToString@level<>"*",rootDir];
    Print["Reading data from files ", admFiles, " at level ", level];
    adm=Table[{admFiles[[i]],admradii[[i]]} , {i,1,Length@admFiles}];
    
    ADMFit1Res[adm,massScale];
    ];
    
    
    BBHDataReduce[modesdir_,IDroot_,ReduceRoot_]:=Module[{modesroot,reducedir,outerradius,psidfile,configStr,
    idDir,m1,m2,madm,sep,abschi1,abschi2,s1x,s1y,s1z,s2x,s2y,s2z,a1x,a1y,a1z,a2x,a2y,a2z,q,eta,mtot,qString,
    dString,chiString,myround,tmpString,modeDecompConfig,bitant,lmax=6,cleanModeConfig,extractionRadii,
    extractionRadiiString,modesDecompDir,levels,radii,stdout},
    
    tmpString = Map[ToString,{m1,m2,madm,sep,abschi1,abschi2,q}];
    
    modesroot = ParentDirectory@modesdir;
    Print["========================================================"];
    Print["Start case modesroot = ", modesroot];
    
    
    If[ReduceRoot==".",
      reducedir=modesroot<>"/DataReduce",
      reducedir=ReduceRoot<>"/"<>Last@StringSplit[modesroot,"/"]<>"/DataReduce"
    ];
    
    Print["Data reduction directory is ", reducedir];
    CreateDirectory[reducedir];
    
    modesDecompDir = reducedir<>"/ModeDecomp";
    Print["Mode decomposition directory is ", modesDecompDir];
    CreateDirectory[modesDecompDir];
    
    outerradius = FileNames["*hmod.r*.l*",modesdir];
    
    outerradius = Map[Last@StringSplit[#,"/"]&,outerradius];
    levels = Map[IntegerPart,Union@Map[levelFun,outerradius]];
    Print["Levels ", levels];
    radii  = TakeColumn[levels,1];
    levels = TakeColumn[levels,2];
    outerradius = Union@Map[StringSplit[#,"."][[2]]&,outerradius];
    outerradius = Union@StringReplace[outerradius,"r"-> ""];
    
    extractionRadii = outerradius;
    Print["Radiation extraction radii: ", extractionRadii];
    outerradius = Max@Map[ToExpression,outerradius];
    Print["Maximal radiation extraction radius: ", outerradius];
    
    psidfile = BAMStringParameter[modesroot,"punctures_ps_file"];
    psidfile = Last@StringSplit[psidfile,"/"];
    psidfile = StringReplace[psidfile," "->""];
    Print["psid-file determined from evolution parameter file as: ", psidfile];
    idDir = LocateInitialDataDirectory[IDroot,psidfile];
    {m1,m2,madm,sep,abschi1,abschi2,s1x,s1y,s1z,s2x,s2y,s2z} = InitialDataParameters[idDir,psidfile];
    mtot = m1+m2;
    q = Max[m1,m2]/Min[m1,m2];
    {a1x,a1y,a1z}={s1x,s1y,s1z}/m1^2;
    {a2x,a2y,a2z}={s2x,s2y,s2z}/m2^2;
    
    Map[ValPrint,tmpString];
    Print["a1 =", {a1x,a1y,a1z}];
    Print["a2 =", {a2x,a2y,a2z}];
    
    myround[x_]:=N[Round[10000x]/10000];
    
    qString = ToString@myround[q];
    dString = ToString@myround[sep/mtot];
    
    bitant = False;
    If[s1x==s1y==s2x==s2y==0,
    Print["Aligned spins!"];
    bitant = True;
    chiString = ToString@myround@a1z<>"_"<>ToString@myround@a2z;,
    Print["Precessing spins!"];
    chiString = "Precess"<>ToString@myround@a1x<>"_"<>ToString@myround@a1y<>ToString@myround@a1z<>"_"<>ToString@myround@a2x<>ToString@myround@a1y<>"_"<>ToString@myround@a2z;
    ];
    
    configStr = "D"<>dString<>"_q"<> qString<>"_"<>chiString;
    Print["Generated configuration string ", configStr];
    Print["Finished Initialization, starting analysis"];
    
    modeDecompConfig="dirnames = {\""<>modesroot<>"\"};\n\n"<>
    "mass1 = {"<>ToString@m1<>"};\n"<>
    "mass2 = {"<>ToString@m2<>"};\n\n"<>
    "minradius = "<>ToString@Min@radii<>";\n"<>
    "maxradius = "<>ToString@Max@radii<>";\n"<>
    "level = "<>ToString@Max@levels<>";\n\n"<> (* This is not really correct, want a list here. *)
    "lmin = 2;\n"<>
    "lmax = "<>ToString@lmax<>";\n\n"<>
    "BITANT = "<>ToString@bitant<>";";
    
    Export[modesDecompDir<>"/allmodesH.in.m",modeDecompConfig,"Text"];
    
    extractionRadii=WaveExtractionRadii[modesroot,modesdir];
    extractionRadiiString=ToString@extractionRadii;
    
    cleanModeConfig="dirnames = {\""<>modesDecompDir<>"\"};\n\n"<>
    "RadiusValue = "<>ToString@extractionRadiiString<>";\n"<>
    "winoffset = -1;\n\n"<>
    "minradius = "<>ToString@Min@radii<>";\n"<>
    "maxradius = "<>ToString@Max@radii<>";\n"<>
    "level = "<>ToString@Max@levels<>";\n\n"<> (* This is not really correct, want a list here. *)
    "lmin = 2;\n"<>
    "lmax = "<>ToString@lmax<>";\n\n"<>
    "dtResample = -1;\n"<>
    "mintime = 150\n"<>
    "maxtime = 5000";
    
    Export[modesDecompDir<>"/cleanPsi.in.m",cleanModeConfig,"Text"];
    
    ADMReduce[modesroot,reducedir];
    
    (* CopyL2mode[configStr,modesdir,reducedir];
    CopyL2modes[configStr,modesdir,reducedir];
    CopyModes[configStr,modesdir,reducedir,2,1];
    CopyModes[configStr,modesdir,reducedir,3,2];
    CopyModes[configStr,modesdir,reducedir,3,3];
    CopyModes[configStr,modesdir,reducedir,4,4];*)
    
    CopyFiles[configStr,modesdir,reducedir,
    {
    "hmod.r*",
    "psi3col.r*"
    }];
    
    stdout = Last@FileNameSplit@First@FileNames["stdout.*",modesroot];
    Print["Found example stdout file: ", stdout];
    
    (* overwrite modes files if we find some in the main directory *)
    CopyFiles[configStr,modesroot,reducedir,
    {
    "d*dt_r*.t*","ADM_mass_r*.t*",
    "moving_puncture_integrate1.txyz*",
    "moving_puncture_integrate2.txyz*",
    "hmod.r*",
    "psi3col.r*",
    "system.log*",
    stdout
    }];
    
    Print["========================================================\n\n"];
    ];
    
    
    SymmetriesFromParfile[parfile_]:=Module[{grid,bitant, quadrant},
    
    grid = BAMStringParameter[FileNameDrop[parfile,-1], "grid"];
    
    If[StringMatchQ[grid, "*quadrant*"], quadrant = True, quadrant = False];
    If[StringMatchQ[grid, "*bitant*"],   bitant   = True, bitant   = False];
    
    {bitant, quadrant}
    ];
    
    
    WriteModeDecompConfigFile[DirectoryRules_?ListQ,
     OptionsPattern[{"LMin" -> 2,"LMax" -> 6, "RunDecomp"-> False, "ModesTargetDirectory" -> "","JustPretendCalculation"-> False}]]:=Module[{modeDecompConfig,
    evDir,targetDir,evDirHasModes,redDirHasModes,modesroot,parfile,
    evModesDir,idDir,exportFile,
    psidFile,psidRules,
    M1,M2,radii,bitant,quadrant,
    modesDecompDir,levels,lmin,lmax,runDecomp,command,justPretendCalculation},
    
    Print["======= entering WriteModeDecompConfigFile ======="];
    
    lmin      = OptionValue["LMin"];
    lmax      = OptionValue["LMax"];
    runDecomp = OptionValue["RunDecomp"];
    targetDir = OptionValue["ModesTargetDirectory"];
    justPretendCalculation = OptionValue["JustPretendCalculation"];
    
    evDir     = "EvolutionDirectory" /. DirectoryRules;
    If[TrueQ[targetDir == ""], targetDir = evDir];
    
    evDirHasModes  = "EvolutionDirectoryHasModes"/. DirectoryRules;
    (* targetDirHasModes = "ModesTargetDirectoryHasModes"  /. DirectoryRules;*)
    
    evModesDir  = "EvolutionModesDir" /. DirectoryRules;
    
    idDir    = "InitialDataDir" /. DirectoryRules;
    psidFile = "PSIDFile"       /. DirectoryRules;
    
    psidRules = PSID2Rules[psidFile];
    
    modesDecompDir = FileNameJoin[{targetDir,"Psi4ModeDecomp"}];
    
    If[FileType@modesDecompDir != None, modesDecompDir == modesDecompDir <> "_new"];
    If[FileType@modesDecompDir != None, Print["modes directory already exists!"]; Return[]];
    
    CreateDirectory[modesDecompDir];
    
    M1 = "M1" /. psidRules;
    M2 = "M2" /. psidRules;
    
    radii = Range@Length@BAMExtractionRadii[evDir];
    
    levels = Union@Flatten@TakeColumn[BAMExtractionRadii[evDir],2];
    
    parfile = ParfileInDirectory[evDir];
    {bitant, quadrant} = SymmetriesFromParfile[parfile];
    
    modesroot = evDir;
    modeDecompConfig =
    "(* run as cd " <> modesDecompDir <> "; math < " <> Global`BAMToolsDir <> "/ModeDecomp/allmodes.m > out.m 2> err.m & *)\n\n" <>
    "dirnames = {\""<>modesroot<>"\"};\n\n"<>
    "mass1 = {"<> ToString@M1 <>"};\n"<>
    "mass2 = {"<> ToString@M2 <>"};\n\n"<>
    "minradius = "<>ToString@Min@radii<>";\n"<>
    "maxradius = "<> ToString@Max@radii<>";\n"<>
    "level = "    <> ToString@Max@levels<>";\n\n"<> (* This is not really correct, want a list here. *)
    "lmin = "     <> ToString@lmin <>";\n"<>
    "lmax = "     <> ToString@lmax<>";\n\n"<>
    "PISYM = "    <> ToString@quadrant<>";\n" <>
    "BITANT = "   <> ToString@bitant<>";\n\n" <>
    "ExportRPsi4IPsi4 = False;";
    
    If[justPretendCalculation,
       modeDecompConfig = modeDecompConfig <> "\n\nJustPretendCalculation = True;";
    ];
    
    (* export mode decomp config file *)
    exportFile = modesDecompDir<>"/allmodes.in.m";
    Print["Exporting file: ", exportFile];
    Export[exportFile,modeDecompConfig,"Text"];
    
    
    CopyFile[Global`BAMToolsDir <> "/ModeDecomp/allmodes.m", modesDecompDir <>"/allmodes.m"];  (* copy the mode decomp file *)
    
    (* now write a shell script *)
    modeDecompConfig =
    "#!/bin/bash\n"                                                       <>
    "rm -f ModeDecomp_Done\n\n"                                           <>
    "cd " <> modesDecompDir <> "; math < allmodes.m > out.m 2> err.m\n\n"  <>
    "rm -f *psi*m-0.gz";
    
    If[Not@justPretendCalculation,
      modeDecompConfig = modeDecompConfig <> "\ntouch ModeDecomp_Done";  
    ];
    
    
    exportFile = modesDecompDir<>"/runModeDecomp.sh";
    Print["Exporting file: ", exportFile];
    Export[exportFile,modeDecompConfig,"Text"];
    
    (* change file permissions to execute *)
    command = "cd " <> modesDecompDir <> "; chmod u+x runModeDecomp.sh";
    Print["Running command ", command];
    Run[command];
    
    If[runDecomp,
    Print["Running mode decomposition in batch mode"];
    If[justPretendCalculation,
     command = "cd " <> modesDecompDir <> "; ./runModeDecomp.sh > out 2> err";,   (*        wait for result *)
     command = "cd " <> modesDecompDir <> "; ./runModeDecomp.sh > out 2> err &";  (* do not wait for result *)
    ];
    
    Print["Running command ", command];
    Run[command];
    ];
    
    Print["======= exiting WriteModeDecompConfigFile ======="];
    ];
    
    
    BAMExtractionRadii[evolutiondir_]:=Module[{outerradius,levels,radii,extractionRadii,levelfun,
    evolutionParfile,collectLevels,r,l,i},
    
    levelfun[str_]:=ToExpression@First@StringCases[str,"rpsi4.r"~~r:NumberString..~~".l"~~l:NumberString-> {r,l}];
    
    collectLevels[level_,levels_] := Map[Last,Select[levels,#[[1]]==level&]];
    
    outerradius = FileNames["rpsi4.r*.l*",evolutiondir];
    
    outerradius = Map[Last@StringSplit[#,"/"]&,outerradius];
    
    levels = Map[IntegerPart,Union@Map[levelfun,outerradius]];
    
    outerradius = Union@Map[StringSplit[#,"."][[2]]&,outerradius];
    outerradius = Union@StringReplace[outerradius,"r"-> ""];
    
    extractionRadii = outerradius;
    Print["Radiation extraction radii from psi4 file names: ", extractionRadii];
    
    evolutionParfile = ParfileInDirectory[evolutiondir];
    
    radii = BAMNumberParametersInFile[evolutionParfile, "invariants_modes_r"];
    
    If[Length@radii >  Length@extractionRadii,
     Print["Found inconsistent extraction radii, assume outermost radius did not fit on extraction levels"];
     Print["BAMExtractionRadii: radii = ", radii];
     Print["BAMExtractionRadii: radii = ", extractionRadii];
    ];
    
    If[Length@radii <  Length@extractionRadii,
     Print["Found inconsistent extraction radii!"];
     Print["radii : extractionRadii ", Length@radii, " : ", Length@extractionRadii];
     Print["BAMExtractionRadii: radii = ", radii];
     Print["BAMExtractionRadii: radii = ", extractionRadii];
    Return[{False,False}];
    ];
    
    Table[{radii[[i]],collectLevels[i,levels]},{i,1,Length@radii}]
    ];
    
    
    CurateBBHData[DirectoryRules_?ListQ,
     OptionsPattern[{"DestinationDirectory" -> ".", "SubmitterEmail" -> "Sascha Husa <sascha.husa@gmail.zebra.elephant.com>",
     "ModeDecomp"-> False,"VivekModesSource"-> {} }]]:=Module[{reducedir,outerradius,psidfile,configString,filePrefix,
    idDir,m1,m2,madm,sep,abschi1,abschi2,s1x,s1y,s1z,s2x,s2y,s2z,a1x,a1y,a1z,a2x,a2y,a2z,q,eta,mtot,qString,
    dString,chiString,myround,tmpString,modeDecompConfig,bitant,lmax,extractionRadii,
    extractionRadiiString,modesDecompDir,levels,radii,stdout,evDir,redDir,evDirHasModes,redDirHasModes,modesroot,parfile,
    evModesDir,redModesDir,exportFile,
    psidFile,psidRules,modesdir,modesString,timer,destinationDir,myDestinationDir,X1,X2,P1,P2,
    modesFiles,selectFiles,str,
    rad,lev,l,m,tmp,rules,
    modesFilesAndRadii,trajectoriesSection,metadataFile,baseSection,psi4BAMModesSections,psi4ModesSections,metadataContent,physicsSection,
    evParfileRules,resolution,L,eccentricityData,metaFile,metaFile2,
    radius,startTime,nCycles,startFreq,cycleInfo,
    afterJunkTime=140,
    tt,sx,sy,sz,
    pickHighestLevel,
    BAMLogDir,EnergyMomentumDir,Psi4ModesDir,Psi4BAMModesDir,
    submitterEmail,modeDecomp,canCreateModeFiles,VivekModesSource},
    
    Print["Entering Function CurateBBHData"];
    
    (* determine where we can find data, and where we will put data *)
    destinationDir = OptionValue["DestinationDirectory"];
    submitterEmail = OptionValue["SubmitterEmail"];
    modeDecomp     = OptionValue["ModeDecomp"];
    
    Print["DirectoryRules: ", DirectoryRules];
    
    evDir  = "EvolutionDirectory"/. DirectoryRules;
    redDir = "ReducedDirectory"  /. DirectoryRules;
    
    VivekModesSource = OptionValue["VivekModesSource"];
    Print["Value of VivekModesSource: ", VivekModesSource];
    
    configString = FileNameTake[evDir,-1];
    Print["Configuration Name = ", configString];
    
    evDirHasModes  = "EvolutionDirectoryHasModes"/. DirectoryRules;
    redDirHasModes = "ReducedDirectoryHasModes"  /. DirectoryRules;
    
    evModesDir  = "EvolutionModesDir" /. DirectoryRules;
    redModesDir = "ReducedModesDir"   /. DirectoryRules;
    
    idDir    = "InitialDataDir" /. DirectoryRules;
    psidFile = "PSIDFile"       /. DirectoryRules;
    
    parfile = ParfileInDirectory[evDir];
    
    (* create uuid and basic.bbh files in evolution directory, if they do not alrady exist *)
    tmp = FileNames["uuid",evDir];
    If[Length@tmp == 0,
      Run["cd "<> evDir <> "; uuidgen > uuid"];
    ];
    
    tmp = FileNames["basic.bbh",evDir];
    If[Length@tmp == 0,
      Run["cd "<> evDir <> "; cp " <> FileNameJoin[{Global`BBHReduceDir, "bam_metadata_template.bbh"}] <> " basic.bbh"];
    ];
    
    
    (* determine from where we will take the psi4-modes*)
    modesString = "psi3col*";
    If[StringQ@redModesDir,
       Print["Using modes from reduced-data directory: ", redModesDir]; 
       modesdir = redModesDir,
      If[StringQ@evModesDir,
       Print["Using modes from reduced-data directory: ", evModesDir]; 
       modesdir = evModesDir,
       Print["Using l=2 modes from BAM on-the-fly integration in directory: ", evDir]; 
       modesString = "*psi4*mode*";
    ]; 
    ]; 
    
    
    myDestinationDir = FileNameJoin[{destinationDir,configString}];
    Print["file destination directory: ", myDestinationDir];
    
    If[FileType@destinationDir == Directory,
      Print["Destination directory already exists!"];,
      CreateDirectory[destinationDir];
    ];
    
    If[FileType@myDestinationDir == Directory,
      Print["Configuration-specific destination directory already exists"];,
      CreateDirectory[myDestinationDir];
      Print[myDestinationDir]
    ];
    
    CreateDirectory[myDestinationDir];
    
    
    (* determine simulation parameters *)
    psidRules = PSID2Rules[psidFile];
    evParfileRules = ParfileToRules[parfile];
    
    tmpString = Map[ToString,{m1,m2,madm,sep,abschi1,abschi2,q}];
    {m1,m2,madm,sep,abschi1,abschi2,s1x,s1y,s1z,s2x,s2y,s2z,X1,X2,P1,P2} = InitialDataParameters[psidFile];
    
    mtot = m1 + m2;
    
    q = Max[m1,m2]/Min[m1,m2];
    
    {a1x,a1y,a1z}={s1x,s1y,s1z}/m1^2;
    {a2x,a2y,a2z}={s2x,s2y,s2z}/m2^2;
    Print["a1 =", {a1x,a1y,a1z}];
    Print["a2 =", {a2x,a2y,a2z}];
    
    resolution = "UNDEFINED";
    
    If[IsFPNumberQ@ToString["nxyz" /. evParfileRules],
     resolution = "nxyz" /. evParfileRules;
    ];
    
    If[IsFPNumberQ@ToString["amr_move_nxyz" /. evParfileRules],
     resolution = "amr_move_nxyz" /. evParfileRules;
    ];
    
    tmp = "amr_nxyz" /. evParfileRules;
    If[ListQ@tmp,
      If[IsFPNumberQ@ToString@Last@tmp,
        resolution = Last["amr_nxyz" /. evParfileRules];
      ];
    ];
    
    L = (X1-X2)\[Cross]P1;
    
    physicsSection = {
    "submitter-email"              ->  submitterEmail,
    "simulation-name"              ->  configString,
    "resolution"                   ->  resolution,
    "initial-ADM-energy"           ->  "Madm" /. psidRules,
    "initial-angular-momentumx"    ->  L[[1]],
    "initial-angular-momentumy"    ->  L[[2]],
    "initial-angular-momentumz"    ->  L[[3]],
    "initial-separation"           ->  sep,
    (*     *)
    "eccentricity"                 ->  "",
    (*     *)
    "freq-start-22"                ->  "",
    "number-of-cycles-22"          ->  "",
    (*     *)
    "phase-error"                  -> "",
    "amplitude-error-relative"     -> "",
    (*     *)
    "after-junkradiation-time"     ->  afterJunkTime,
    (*     *)
    "mass1"                        ->  m1,
    "mass2"                        ->  m2,
    (*     *)
    "initial-bh-position1x" ->   X1[[1]],
    "initial-bh-position1y" ->   X1[[2]],
    "initial-bh-position1z" ->   X1[[3]],
    (*     *)
    "initial-bh-position2x" ->   X2[[1]],
    "initial-bh-position2y" ->   X2[[2]],
    "initial-bh-position2z" ->   X2[[3]],
    (*     *)
    "initial-bh-momentum1x"  ->   P1[[1]],
    "initial-bh-momentum1y"  ->   P1[[2]],
    "initial-bh-momentum1z"  ->   P1[[3]],
    (*     *)
    "initial-bh-momentum2x"  ->   P2[[1]],
    "initial-bh-momentum2y"  ->   P2[[2]],
    "initial-bh-momentum2z"  ->   P2[[3]],
    (*     *)
    "initial-bh-spin1x" ->   s1x,
    "initial-bh-spin1y" ->   s1y,
    "initial-bh-spin1z" ->   s1z,
    (*     *)
    "initial-bh-spin2x" ->   s2x,
    "initial-bh-spin2y" ->   s2y,
    "initial-bh-spin2z" ->   s2z,
    (*     *)
    "after-junkradiation-spin1x"   ->  "",
    "after-junkradiation-spin1y"   ->  "",
    "after-junkradiation-spin1z"   ->  "",
    (*     *)
    "after-junkradiation-spin2x"   ->  "",
    "after-junkradiation-spin2y"   ->  "",
    "after-junkradiation-spin2z"   ->  ""
    };
    physicsSection = Drop[ComposeConfigSection["dummy",physicsSection],1];
    
    Map[ValPrint,tmpString];
    
    myround[x_]:=N[Round[10000x]/10000];
    
    qString = ToString@myround[q];
    dString = ToString@myround[sep/mtot];
    
    bitant = False;
    If[s1x==s1y==s2x==s2y==0,
    Print["Aligned spins!"];,  (* TODO:  check consistency with settings in evolution-parfile *)
    Print["Precessing spins!"];
    ];
    
    (* copy and create files *)
    filePrefix = ""; (* do not prefix files *)
    
    (* stdout and timer files: only copy 2 files each *)
    stdout = FileNames["stdout.*",evDir];
    If[Length@stdout >= 1,
      stdout = {Last@FileNameSplit@First@stdout, Last@FileNameSplit@Last@stdout}
    ];
    Print["Found example stdout file: ", stdout];
    
    timer = FileNames["timer.*",evDir];
    If[Length@timer >= 1,
      timer = {Last@FileNameSplit@First@FileNames["timer.*",evDir], Last@FileNameSplit@Last@FileNames["timer.*",evDir]};
    ];
    Print["Found example timer file: ", timer];
    
    BAMLogDir = FileNameJoin[{myDestinationDir,"BAMLogs"}];
    CreateDirectory[BAMLogDir];
    
    If[Length@stdout > 0, CopyFiles[filePrefix, evDir, BAMLogDir, stdout]];
    If[Length@timer  > 0, CopyFiles[filePrefix, evDir, BAMLogDir, timer]];
    
    (* GW modes *)
    Psi4BAMModesDir = FileNameJoin[{myDestinationDir,"Psi4BAMModes"}];
    CreateDirectory[Psi4BAMModesDir];
    
    (* Psi4ModesDir = FileNameJoin[{myDestinationDir,"Psi4Modes"}];*)
    Psi4ModesDir = myDestinationDir;
    (*CreateDirectory[Psi4ModesDir];*)
    
    
    
    Print["Copying GW modes of the form ", modesString];
    CopyFiles[filePrefix,evDir,Psi4BAMModesDir, {modesString}];
    
    modesFiles = FileNames[modesString,myDestinationDir,2];
    
    selectFiles= Select[modesFiles,StringMatchQ[#,"*/rpsi4mode*"]&];
    
    Map[BAMModesFilesTo3Col[#,"DeleteSourceFiles"-> True]&,selectFiles];
    modesString = "psi3col*";
    
    (* now that we have a unified naming convention for psi4-files, compose modes metadata sections *)
    modesFiles = FileNames[modesString,myDestinationDir,2];
    
    Print["Calling BAMExtractionRadii"];
    radii = TakeColumn[BAMExtractionRadii[evDir],1];
    Print["BAMExtractionRadii found radii = ", radii];
    
    modesFilesAndRadii = Table[
    str = modesFiles[[i]];
    {StringCases[str,"r"~~rad:NumberString~~".l"~~lev:NumberString~~".l"~~l:NumberString~~".m"~~m:NumberString~~___ :> 
            {rad,lev}],str},{i,1,Length@modesFiles}];
    
    pickHighestLevel[xxx_?ListQ]:=\[NonBreakingSpace]Last@Gather[xxx, First@#1 ==First@#2 &];
    
    tmp = Gather[modesFilesAndRadii, First@First@First@#1 == First@First@First@#2 &];
    modesFilesAndRadii = Map[pickHighestLevel,tmp];
    
    
    psi4BAMModesSections = Table[
      modesFiles = TakeColumn[modesFilesAndRadii[[i]],2];
      {rad,lev}  = TakeColumn[modesFilesAndRadii[[i]], 1][[1, 1]];
    
      str = radii[[ToExpression@rad]];
    
      tmp=ComposeStrainModesSection[First@modesFiles,
         "Verbose"-> False,"SectionHeader"-> "psi4t-data",
         "Modes"-> "All","Path" -> "Psi4BAMModes/"];
      tmp[[2]] = {tmp[[2]], "extraction-radius = "<>  ToString@str <>"\n"};
      AppendTo[tmp, "\n\n"];
      Flatten@tmp,
    {i,1,Length@modesFilesAndRadii}];
    ExportText[FileNameJoin[{myDestinationDir,"psi4BAMmodes.bbh"}], Flatten@psi4BAMModesSections,"Overwrite"->True];
    
    (* NOW THE WHOLE THING AGAIN FOR THE FULL MODES *)
    
    (* First, check if Vivek-style modes are available *)
    VivekModesFiles = FileNames[modesString,VivekModesSource,2];
    Print["Found Vivek-style modes files: ", VivekModesFiles];
    
    If[Length@VivekModesFiles > 0,
    modesFiles=VivekModesFiles;
    modesDecompDir = FileNameJoin[{Psi4ModesDir,"Psi4ModeDecomp"}];
    (*CreateDirectory[modesDecompDir];*)
    
    CopyDirectory[VivekModesSource<>"/data",modesDecompDir];
    
    canCreateModeFiles = False;,
    
    (* check for psi4 source files and whether mode decomp can be run *)
    WriteModeDecompConfigFile[DirectoryRules,"RunDecomp"-> True,
         "ModesTargetDirectory"-> Psi4ModesDir,"JustPretendCalculation"-> True];
    
    modesFiles = FileNames[modesString,Psi4ModesDir<>"/Psi4ModeDecomp",2];
    If[Length@modesFiles == 0,
     Print["Mode decomposition yields no results - no mode files can be created."];
     canCreateModeFiles = False;,
     canCreateModeFiles = True;
    ];
    ];
    
    modesFilesAndRadii = Table[
    str = modesFiles[[i]];
    {StringCases[str,"r"~~rad:NumberString~~".l"~~lev:NumberString~~".l"~~l:NumberString~~".m"~~m:NumberString~~___ :> 
            {rad,lev}],str},{i,1,Length@modesFiles}];
    
    pickHighestLevel[xxx_?ListQ]:=\[NonBreakingSpace]Last@Gather[xxx, First@#1 ==First@#2 &];
    
    tmp = Gather[modesFilesAndRadii, First@First@First@#1 == First@First@First@#2 &];
    modesFilesAndRadii = Map[pickHighestLevel,tmp];
    
    psi4ModesSections = Table[
      modesFiles = TakeColumn[modesFilesAndRadii[[i]],2];
      {rad,lev}  = TakeColumn[modesFilesAndRadii[[i]], 1][[1, 1]];
    
      str = radii[[ToExpression@rad]];
    
      tmp=ComposeStrainModesSection[First@modesFiles,
         "Verbose"-> False,"SectionHeader"-> "psi4t-data",
         "Modes"-> "All","Path" -> "Psi4ModeDecomp/"];
      tmp[[2]] = {tmp[[2]], "extraction-radius = "<>  ToString@str <>"\n"};
      AppendTo[tmp, "\n\n"];
      Flatten@tmp,
    {i,1,Length@modesFilesAndRadii}];
    
    Export[FileNameJoin[{myDestinationDir,"psi4modes.bbh"}], Flatten@psi4ModesSections,"Text"];
    
    
    (* overwrite modes files if we find some in the main directory *)
    CopyFiles[filePrefix,evDir,myDestinationDir,
    {
    "moving_puncture_integrate1.txyz*",
    "moving_puncture_integrate2.txyz*",
    "psi3col.r*",
    "*.par*",
    "horizon_*",
    "ah.xy*",
    "basic.bbh",
    "uuid"
    }];
    
    
    CopyFiles[filePrefix,evDir,BAMLogDir,
    {"system.log*"}];
    
    EnergyMomentumDir = FileNameJoin[{myDestinationDir,"EnergyMomentum"}];
    CreateDirectory[EnergyMomentumDir];
    
    CopyFiles[filePrefix,evDir,EnergyMomentumDir,
    {"d*dt_r*.t*","ADM_mass_r*.t*"}];
    
    (* compress some files *)
    Run["cd "<> myDestinationDir <> "; " <>
        "gzip moving_puncture_integrate* system.log* ah.xy*/*"];
    
    Run["cd "<> EnergyMomentumDir <> "; " <>
        "gzip d*dt_r*.t* ADM_mass_r*.t*"];
    
    Run["cd "<> BAMLogDir <> "; " <>
        "gzip *"];
    
    (* convert horizon and puncture track files to NRAR format and add entries to metadata file *)
    rules={};
    
    tmp = FileNames["horizon_*",myDestinationDir,2];
    If[Length@tmp == 0,
    Print["Could not find horizon files!"];,
    
    Map[BAMHorizonFilesToNRARFiles[#,"DeleteSourceFiles"-> False]&,tmp];
    
    tmp = LastInPath@First@FileNames["hmass_1.*",myDestinationDir,2];
    AppendTo[rules, "horizon-mass1" -> tmp];
    
    tmp = LastInPath@First@FileNames["hmass_2.*",myDestinationDir,2];
    AppendTo[rules, "horizon-mass2" -> tmp];
    
    tmp = LastInPath@First@FileNames["hspin_1.*",myDestinationDir,2];
    AppendTo[rules, "spin1" -> tmp];
    
    tmp = LastInPath@First@FileNames["hspin_2.*",myDestinationDir,2];
    AppendTo[rules, "spin2" -> tmp];
    
    tmp = LastInPath@First@FileNames["htraj_1.*",myDestinationDir,2];
    AppendTo[rules, "horizon-center1" -> tmp];
    
    tmp = LastInPath@First@FileNames["htraj_2.*",myDestinationDir,2];
    AppendTo[rules, "horizon-center2" -> tmp];
    ];
    
    tmp = Join[FileNames["moving_puncture_integrate1.txyz*",myDestinationDir,2],
               FileNames["moving_puncture_integrate2.txyz*",myDestinationDir,2]];
    
    If[Length@tmp == 0,
    Print["Could not find puncture track files!"];,
    Print["Found puncture track files: ", tmp];
    Map[BAMTrajectoryFileTo4Col[#,"DeleteSourceFiles"-> False]&,tmp];
    
    (* below, Last is used to get data at the finest available level *)
    tmp = Last@FileNames["traj_1.l*",myDestinationDir,2];
    AppendTo[rules,"trajectory1" -> FileNameTake[tmp,-1]];
    
    tmp = Last@FileNames["traj_2.l*",myDestinationDir,2];
    AppendTo[rules, "trajectory2" -> FileNameTake[tmp,-1]];
    ];
    
    
    trajectoriesSection = ComposeConfigSection["body-data", Flatten@rules, "Verbose" -> True];
    
    tmp=FileNames["basic.bbh",myDestinationDir,2];
    If[Length@tmp == 0,
      Print["Missing metadata file basic.bbh!"];
      baseSection = {"[metadata]","# empty"};,
      metadataFile = First@tmp;
      If[FileType@metadataFile == File,
       tmp = ConfigFileToRules@metadataFile;
       baseSection = Flatten@StringSplit[StringSplit[Import[metadataFile,"String"],EndOfLine],"\n"];
       DeleteFile@metadataFile;
    ];
    ];
    
    
    metadataContent = Flatten@Join[baseSection,physicsSection,trajectoriesSection,psi4BAMModesSections];
    metaFile = FileNameJoin[{myDestinationDir, configString<>".raw.bbh"}];
    ExportText[metaFile,metadataContent,"Overwrite"-> True];
    
    metadataContent = Flatten@Join[baseSection,physicsSection,trajectoriesSection,psi4ModesSections];
    metaFile2 = FileNameJoin[{myDestinationDir, configString<>".bbh"}];
    ExportText[metaFile2,metadataContent,"Overwrite"-> True];
    
    
    (* compute eccentricity and modify metadata *)
    
    (* use the maximal finite extraction radius *)
    eccentricityData = NinjaBase`NinjaEccentricity[metaFile,"SectionInstance" -> "MaxRad","Verbose"-> True];
    ChangeConfigEntry[metaFile,metaFile, "eccentricity", ToString@eccentricityData[[1]]];
    ChangeConfigEntry[metaFile2,metaFile2, "eccentricity", ToString@eccentricityData[[1]]];
    
    Print["select max from metaFile ", metaFile];
    radius=ReadNinjaModes[FileNameTake[metaFile,FileNameDepth@metaFile-1],FileNameTake[metaFile,-1],
            "MModes"-> {2,2}, "DataSection" -> "psi4t-data","SectionInstance" -> "MaxRad", "Tag" -> "mode22","Verbose" -> True,
            "OnlyComputeExtractionInfo" -> True];
    
    Print[];
    Print["Found extraction radii = ", FullForm@radius];
    Print["Found max radius = ", radius=Max[StringToNumber /@ radius]];
    
    
    (* compute start frequency and number of cycles, and modify metadata *)
    startTime = afterJunkTime + radius;
    Print["Calling NinjaCyclesStartFreq with startTime = ", startTime];
    cycleInfo = NinjaCyclesStartFreq[metaFile,startTime];
    
    startFreq = "freq-start-22"    /.cycleInfo;
    nCycles = "number-of-cycles-22" /. cycleInfo;
    
    Print["found start freq = ", startFreq, ", ", "#\[NonBreakingSpace]of cycles = ", nCycles ];
    
    (* tmp = ListPlot[{Abs@ninjlm["mode22",2,2],Re@ninjlm["mode22",2,2],
      {{afterJunkTime + radius,0},{afterJunkTime + radius,Max[TakeColumn[Abs@ninjlm["mode22",2,2],2]]}}},
    PlotRange\[Rule] All,Joined\[Rule] True];
    Export[FileNameJoin[{myDestinationDir,"timedomain_psi4_22.pdf"}],tmp]; *)
    
    ChangeConfigEntry[metaFile,metaFile, "freq-start-22", ToString@startFreq];
    ChangeConfigEntry[metaFile,metaFile, "number-of-cycles-22", ToString@nCycles];
    
    ChangeConfigEntry[metaFile2,metaFile2, "freq-start-22", ToString@startFreq];
    ChangeConfigEntry[metaFile2,metaFile2, "number-of-cycles-22", ToString@nCycles];
    
    
    Print["changed metadata entries for freq-start-22 and number-of-cycles-22"];
    
    (* compute spins after junk-radiation, and modify metadata *)
    Clear[tt,sx,sy,sz];
    tmp = ReadNinjaData[metaFile,"body-data", "spin1"];
    If[ListQ@tmp,
      tmp = tmp /. {tt_?NumberQ,sx_?NumberQ,sy_?NumberQ,sz_?NumberQ} -> {tt,{sx,sy,sz}};
      tmp = Interpolation@tmp;
      tmp = Chop/@ tmp@afterJunkTime;
      Print["spin1 after junk = ", tmp];
      ChangeConfigEntry[metaFile,metaFile, "after-junkradiation-spin1x", ToString@tmp[[1]]];
      ChangeConfigEntry[metaFile,metaFile, "after-junkradiation-spin1y", ToString@tmp[[2]]];
      ChangeConfigEntry[metaFile,metaFile, "after-junkradiation-spin1z", ToString@tmp[[3]]];
    
      ChangeConfigEntry[metaFile2,metaFile2, "after-junkradiation-spin1x", ToString@tmp[[1]]];
      ChangeConfigEntry[metaFile2,metaFile2, "after-junkradiation-spin1y", ToString@tmp[[2]]];
      ChangeConfigEntry[metaFile2,metaFile2, "after-junkradiation-spin1z", ToString@tmp[[3]]];
    ];
    
    
    tmp = ReadNinjaData[metaFile,"body-data", "spin2"]  /. {tt_,sx_,sy_,sz_} -> {tt,{sx,sy,sz}};
    If[ListQ@tmp,
      tmp = Interpolation@tmp;
      tmp = Chop /@ tmp@afterJunkTime;
      Print["spin2 after junk = ", tmp];
      ChangeConfigEntry[metaFile,metaFile, "after-junkradiation-spin2x", ToString@tmp[[1]]];
      ChangeConfigEntry[metaFile,metaFile, "after-junkradiation-spin2y", ToString@tmp[[2]]];
      ChangeConfigEntry[metaFile,metaFile, "after-junkradiation-spin2z", ToString@tmp[[3]]];
    
      ChangeConfigEntry[metaFile2,metaFile2, "after-junkradiation-spin2x", ToString@tmp[[1]]];
      ChangeConfigEntry[metaFile2,metaFile2, "after-junkradiation-spin2y", ToString@tmp[[2]]];
      ChangeConfigEntry[metaFile2,metaFile2, "after-junkradiation-spin2z", ToString@tmp[[3]]];
    ];
    
    Print["Deleting file ", FileNameJoin[{myDestinationDir,"psi4BAMmodes.bbh"}]];
    DeleteFile[FileNameJoin[{myDestinationDir,"psi4BAMmodes.bbh"}]];
    
    If[canCreateModeFiles,
     WriteModeDecompConfigFile[DirectoryRules,"RunDecomp"-> modeDecomp,"ModesTargetDirectory"-> Psi4ModesDir];,
     CopyFile[metaFile, metaFile2];
    ];
    
    Print["========================================================\n\n"];
    ];
    
    
    NRARPsi4ToStrain[metaFile_,LMAX_,\[Omega]GWStart_]:=Module[{newMetaFile,dir,newMeta,tmp,psi4tInstances,readOut,filesRead,originalFile,newFile,times,data,col3data,file},
    
    newMetaFile=StringReplace[metaFile,".bbh"->".hFFI.bbh"];
    dir=FileNameTake[metaFile,FileNameDepth[metaFile]-1];
    
    myDestinationDir=FileNameJoin[{FileNameDrop[metaFile,-1],"FFIStrainModes"}];
    
    If[FileType@myDestinationDir==None,
      CreateDirectory[myDestinationDir],
      RenameDirectory[myDestinationDir,myDestinationDir<>"_bak"];
    ];
    
    Run["cp "<>metaFile<>" "<>newMetaFile];
    
    newMeta=Import[newMetaFile,"Text"];
    tmp=StringReplace[newMeta,"[psi4t-data]"->"[ht-data]"];
    tmp=StringReplace[tmp,"Psi4ModeDecomp/psi"->"FFIStrainModes/h"];
    Export[newMetaFile,tmp,"Text"];
    
    psi4tInstances=Length@Cases[ConfigSectionHeaders[metaFile,"Union"-> False],"psi4t-data"];
    
    Do[
    readOut=ReadNinjaModes[dir,FileNameTake[metaFile,-1],"DataSection"->"psi4t-data","LMax"->LMAX,"SectionInstance"->i,"Tag"->"dummy"];
    filesRead=readOut[[2]];
    
    Print@readOut;
    
    Do[
    
    If[m!=0,strain=hFromPsi4FFI[ninjlm["dummy",l,m],(Abs@m/2)\[Omega]GWStart/(2 \[Pi])];,strain=hFromPsi4FFI[ninjlm["dummy",l,m],(1/2)\[Omega]GWStart/(2 \[Pi])];];
    
    originalFile=Last@First@Select[filesRead,#[[1]]=={l,m}&];
    newFile=StringReplace[originalFile,"Psi4ModeDecomp/psi"->"FFIStrainModes/h"];
    
    times=TakeColumn[strain,1];
    data=TakeColumn[strain,2];
    
    col3data=Transpose@{times,Re@data,Im@data};
    file=FileNameJoin[{FileNameDrop[myDestinationDir,-1],newFile}];
    Print["Exporting file ",file];
    Export[file,col3data,{"GZIP","Table"}],{l,2,LMAX},{m,-l,l}]
    ,{i,psi4tInstances}]
    ]
    
    
    Options[SXSLuminosityFromMetaFiles]={"PrecessingTolerance"->0.001`,"OutputTag"->"SXSLuminosities_","OutputDir"->".","ExtrapOrder"->{2},"Verbose"->"False"}
    
    
    SXSLuminosityFromMetaFiles[metadata_?ListQ,modes_?ListQ,OptionsPattern[]]:=Module[{tags,metarules,\[CapitalOmega]startMetadata,M1,M2,\[Chi]1vect,\[Chi]2vect,\[Chi]1,\[Chi]2,M,\[Eta],q,A1,A2,Lx,Ly,Lz,Lvect,
    LvectNorm,precvalue1,precvalue2,precvalueNorm1,precvalueNorm2,precvalueNorm,precessing,\[Chi]1p,\[Chi]2p,Initial\[Chi]Plane,\[Chi]1vect\[Chi]2vect,fnameBase,myFile,mymodes,h5file,dir,
    ta,test\[Psi],pos,test\[Psi]22,domain,times,f0,mm,lum,sumModes,sumModesDominant,sumModesLum,sumModesLumDom,tmaxSum,tmaxSumDom,modesMaxima,data,int,verbose,lmax,ls,lsms,
    lsmsm1,dommodes,testdommodes,posdom,domsneg,precessingtolerance,outputtag,outputdir,xx,yy,extraporder},
    
    mymodes=modes;
    
    verbose=OptionValue["Verbose"];
    precessingtolerance=OptionValue["PrecessingTolerance"];
    outputtag=OptionValue["OutputTag"];
    outputdir=OptionValue["OutputDir"];
    extraporder=OptionValue["ExtrapOrder"];
    
    (* read metadata and convert to physical quantities *)
    tags=(FileNameTake[#1,{-3}]&)/@metadata;
    metarules=SXSMetaFilesToRules/@metadata;
    \[CapitalOmega]startMetadata[tags]="initial-orbital-frequency"/. metarules;
    M1=Flatten["initial-mass1"/. metarules,1];
    M2=Flatten["initial-mass2"/. metarules,1];
    M=M1+M2;\[Eta]=(M1 M2)^2/M^2;
    q=(Max[#1,1]&)/@(M1/M2);
    
    \[Chi]1vect=("initial-spin1"/. metarules)/M1^2;
    \[Chi]2vect=("initial-spin2"/. metarules)/M2^2;
    \[Chi]1vect\[Chi]2vect=Transpose[{\[Chi]1vect,\[Chi]2vect}];
    \[Chi]1=Norm/@\[Chi]1vect;
    \[Chi]2=Norm/@\[Chi]2vect;
    
    Lvect="initial-ADM-angular-momentum"/. metarules;
    LvectNorm=Lvect/(Norm/@Lvect);
    Initial\[Chi]Plane=Table[{\[Chi]1vect\[Chi]2vect[[i,1]]-Lvect[[i]] Lvect[[i]].\[Chi]1vect\[Chi]2vect[[i,2]],\[Chi]1vect\[Chi]2vect[[i,2]]-Lvect[[i]] Lvect[[i]].\[Chi]1vect\[Chi]2vect[[i,2]]},{i,Length[Lvect]}];
    \[Chi]1p=Norm/@Initial\[Chi]Plane[[All,1]];
    \[Chi]2p=Norm/@Initial\[Chi]Plane[[All,2]];
    
    
    (* Is it precessing? *)
    A1=1+3/4. q;
    A2=1+3/(4. q);
    precvalue1=Table[(\[Chi]1vect[[i]] A1[[i]])\[Cross]LvectNorm[[i]],{i,Length[tags]}];
    precvalue2=Table[(\[Chi]2vect[[i]] A2[[i]])\[Cross]LvectNorm[[i]],{i,Length[tags]}];
    precvalueNorm1=Norm/@precvalue1;
    precvalueNorm2=Norm/@precvalue2;
    precvalueNorm=precvalueNorm1^2+precvalueNorm2^2;
    Do[If[precvalueNorm[[i]]<0.001,precessing[tags[[i]]]=False,precessing[tags[[i]]]=True],{i,Length[tags]}];
    
    (* Compute Luminosity *)
    Table[
        (* Find rpsi4 files *)
        fnameBase=ToString[outputtag];
        Print["Computing Luminosity from: ",tags[[i]]];
        myFile=ToString[outputdir]<>"/"<>fnameBase<>ToString[tags[[i]]]<>".dat";
    
        dir=FileNameDrop[metadata[[i]],-2];
        h5file=First[FileNames["rMPsi4_Asymptotic_GeometricUnits.h5",dir,2]];
    
        If[FileExistsQ[h5file],Print["Found ",h5file," . Continue"];,Print["Not found ",h5file," . Bye"];Return[]];
        If[precessing[tags[[i]]],Print["Precessing run. Non simetry. Taking all modes"];,Print["Non-Precessing run. Assuming simetry. Taking only m>0 modes"];mymodes=Select[mymodes,#1[[2]]>=0&];];
        Print["Modes selection : ",mymodes];
    
        (* How many dominant modes are present? eg. {{2,2},{2,1},{2,-2}}? *)
        lmax=Max[mymodes[[All,1]]];
        ls=DeleteDuplicates[mymodes[[All,1]]];
        lsms=CombineColumns[ls,ls];
        lsmsm1=CombineColumns[ls,ls-1];
        dommodes=Sort[Join[lsms,lsmsm1]];
        domsneg=dommodes/. {xx_,yy_}->{xx,-yy};
        dommodes=Sort[Join[dommodes,domsneg]];
        testdommodes=(MemberQ[mymodes,#1]&)/@dommodes;
        posdom=Flatten[Position[testdommodes,True]];
        dommodes=dommodes[[posdom]];Print["Dominant modes : ",dommodes];
    
        ta=Table[
                (* Load and process data *)
                Print["Extrapolation order : ",j];
                test\[Psi]=GetAsymptoticMultiMode[h5file,j,mymodes];
                Print["Modes Loaded"];
                
                pos=First[First[Position[mymodes,{2,2}]]];
                test\[Psi]22=test\[Psi][[pos]];
                domain={First[First[test\[Psi]22]],First[Last[test\[Psi]22]]};
                times=Range[First[domain],Last[domain],0.5];
                test\[Psi]=Interpolation/@test\[Psi];
                test\[Psi]=Table[CombineColumns[times,test\[Psi][[i]][times]],{i,Length[test\[Psi]]}];
    
                (* 22 mode guess frequency *)
                f0=Abs[First[GuessFFIf0[test\[Psi][[pos]],"MinTime"->First[First[test\[Psi][[pos]]]]+250]]];
    
                (* compute Luminosoties *)
    
                Do[If[mymodes[[l,2]]!=0,mm=mymodes[[l,2]]/2,mm=(mymodes[[l,2]]+1)/2]; 
                   lum[mymodes[[l,1]],mymodes[[l,2]]]=Abs[TakeColumn[NewsFromPsi4FFI[test\[Psi][[l]],0.75f0 Abs[mm]],2]]^2;Print[mymodes[[l]]];,{l,Length[mymodes]}];
               
                If[precessing[tags[[i]]],
                             sumModes=1/(16 \[Pi])Sum[lum[mymodes[[l,1]],mymodes[[l,2]]],{l,Length@mymodes}];
                             sumModesDominant=1/(16 \[Pi])Sum[lum[dommodes[[l,1]],dommodes[[l,2]]],{l,Length@dommodes}];
                             ,
                             sumModes=1/(8 \[Pi])Sum[lum[mymodes[[l,1]],mymodes[[l,2]]],{l,Length@mymodes}];
                             sumModesDominant=1/(8 \[Pi])Sum[lum[dommodes[[l,1]],dommodes[[l,2]]],{l,Length@dommodes}];];
    
    
                 sumModesLum=CombineColumns[times,sumModes];
                 sumModesLumDom=CombineColumns[times,sumModesDominant];
                 tmaxSum=TimeOfMaximum[sumModesLum];
                 tmaxSumDom=TimeOfMaximum[sumModesLumDom];
                 modesMaxima=Table[data=CombineColumns[times,lum[mymodes[[l,1]],mymodes[[l,2]]]];
                                   int=Interpolation[data];
                                   {ValueOfMaximum[data],int[tmaxSum],int[tmaxSumDom]},{l,Length[mymodes]}];
    
    If[verbose,Print[ListPlot[{sumModesLum,sumModesLumDom},Joined->True,PlotRange->All]];];
    Flatten[{tags[[i]],\[Eta][[i]],\[Chi]1[[i]],\[Chi]2[[i]],\[Chi]1p[[i]],\[Chi]2p[[i]],Initial\[Chi]Plane[[i]],j,ValueOfMaximum[sumModesLum],tmaxSum,ValueOfMaximum[sumModesLumDom],tmaxSumDom,Flatten[modesMaxima]}],{j,extraporder}];
    
    Print["Exporting results to : ",myFile];Export[myFile,ta];
    ,{i,Length[tags]}];
    
    ta]
    
    
    Options[BAMLuminosityFromMetaFiles]={"PrecessingTolerance"->0.001,"OutputTag"->"BAMLuminosities_","OutputDir"->".","ExtractionRadius"->{},"Verbose"->"False"};
    
    
    BAMLuminosityFromMetaFiles[metadata_?ListQ,modes_?ListQ,OptionsPattern[]]:=Module[{tags,metarules,\[CapitalOmega]startMetadata,M1,M2,\[Chi]1vect,\[Chi]2vect,\[Chi]1,\[Chi]2,M,\[Eta],q,A1,A2,Lx,Ly,Lz,Lvect,
    LvectNorm,precvalue1,precvalue2,precvalueNorm1,precvalueNorm2,precvalueNorm,precessing,\[Chi]1p,\[Chi]2p,Initial\[Chi]Plane,\[Chi]1vect\[Chi]2vect,fnameBase,myFile,mymodes,h5file,dir,
    ta,test\[Psi],pos,test\[Psi]22,domain,times,f0,mm,lum,sumModes,sumModesDominant,sumModesLum,sumModesLumDom,tmaxSum,tmaxSumDom,modesMaxima,data,int,verbose,lmax,ls,lsms,
    lsmsm1,dommodes,testdommodes,posdom,domsneg,precessingtolerance,outputtag,outputdir,xx,yy,extractionradius,p1,p2,r1,r2,InitialOrbitalAngularMomentum,InitialAngularMomentum,
    extraddef,myRad,myrads},
    
    mymodes=modes;
    
    verbose=OptionValue["Verbose"];
    precessingtolerance=OptionValue["PrecessingTolerance"];
    outputtag=OptionValue["OutputTag"];
    outputdir=OptionValue["OutputDir"];
    extractionradius=OptionValue["ExtractionRadius"];
    
    (* read metadata and convert to physical quantities *)
    tags=FileNameTake[#,{-3}]&/@metadata;
    metarules=BAMMetaFilesToRules/@metadata;
    
    M1=Flatten["mass1"/. metarules,1];
    M2=Flatten["mass2"/. metarules,1];
    M=M1+M2;
    \[Eta]=(M1 M2)^2/M^2;
    q=(Max[#1,1]&)/@(M2/M1);
    
    \[Chi]1vect=Transpose@{(("initial-bh-spin1x"/. metarules)/M1^2)[[All,1]],(("initial-bh-spin1y"/. metarules)/M1^2)[[All,1]],(("initial-bh-spin1z"/. metarules)/M1^2)[[All,1]]};
    
    \[Chi]2vect=Transpose@{(("initial-bh-spin2x"/. metarules)/M2^2)[[All,1]],(("initial-bh-spin2y"/. metarules)/M2^2)[[All,1]],(("initial-bh-spin2z"/. metarules)/M2^2)[[All,1]]};
    \[Chi]1vect\[Chi]2vect=Transpose[{\[Chi]1vect,\[Chi]2vect}];
    \[Chi]1=Norm/@\[Chi]1vect;
    \[Chi]2=Norm/@\[Chi]2vect;
    
    
    r1=Transpose[{("initial-bh-position1x"/.metarules)[[All,1]],("initial-bh-position1y"/.metarules)[[All,1]],("initial-bh-position1z"/.metarules)[[All,1]]}];
    r2=Transpose[{("initial-bh-position2x"/.metarules)[[All,1]],("initial-bh-position2y"/.metarules)[[All,1]],("initial-bh-position2z"/.metarules)[[All,1]]}];
    p1=Transpose[{("initial-bh-momentum1x"/.metarules)[[All,1]],("initial-bh-momentum1y"/.metarules)[[All,1]],("initial-bh-momentum1z"/.metarules)[[All,1]]}];
    p2=Transpose[{("initial-bh-momentum2x"/.metarules)[[All,1]],("initial-bh-momentum2y"/.metarules)[[All,1]],("initial-bh-momentum2z"/.metarules)[[All,1]]}];
    
    Lvect= Table[p1[[i]]\[Cross]r1[[i]] + p2[[i]]\[Cross]r2[[i]],{i,Length@metadata}];
    LvectNorm=Lvect/(Norm/@Lvect);
    Initial\[Chi]Plane=Table[{\[Chi]1vect\[Chi]2vect[[i,1]]-Lvect[[i]] Lvect[[i]].\[Chi]1vect\[Chi]2vect[[i,2]],\[Chi]1vect\[Chi]2vect[[i,2]]-Lvect[[i]] Lvect[[i]].\[Chi]1vect\[Chi]2vect[[i,2]]},{i,Length[Lvect]}];
    
    \[Chi]1p=Norm/@Initial\[Chi]Plane[[All,1]];
    \[Chi]2p=Norm/@Initial\[Chi]Plane[[All,2]];
    
    
    (* Is it precessing? *)
    A1=1+3/4. q;
    A2=1+3/(4. q);
    precvalue1=Table[(\[Chi]1vect[[i]] A1[[i]])\[Cross]LvectNorm[[i]],{i,Length[tags]}];
    precvalue2=Table[(\[Chi]2vect[[i]] A2[[i]])\[Cross]LvectNorm[[i]],{i,Length[tags]}];
    precvalueNorm1=Norm/@precvalue1;
    precvalueNorm2=Norm/@precvalue2;
    precvalueNorm=precvalueNorm1^2+precvalueNorm2^2;
    Do[If[precvalueNorm[[i]]<0.001,precessing[tags[[i]]]=False,precessing[tags[[i]]]=True],{i,Length[tags]}];
    
    extraddef=Table["extraction-radius"/.metarules[[i]],{i,Length@metarules}];
    
    If[Length@extractionradius!=0,extraddef==TakeColumn[extraddef,extractionradius]];
    
    (* Compute Luminosity *)
    Table[
        (* Find rpsi4 files *)
        fnameBase=ToString[outputtag];
        Print["Computing Luminosity from: ",tags[[i]]];
        myFile=ToString[outputdir]<>"/"<>fnameBase<>ToString[tags[[i]]]<>".dat";
    
        dir=FileNameTake[metadata[[i]],FileNameDepth[metadata[[i]]]-1]; 
        myrads=extraddef[[i]];
      
        If[FileExistsQ[dir],Print["Found ",dir," . Continue"];,Print["Not found ",dir," . Bye"];Return[]];
        If[precessing[tags[[i]]],Print["Precessing run. Non simetry. Taking all modes"];,Print["Non-Precessing run. Assuming simetry. Taking only m>0 modes"];mymodes=Select[mymodes,#1[[2]]>=0&];];
        Print["Modes selection : ",mymodes];
     
        (* How many dominant modes are present? eg. {{2,2},{2,1},{2,-2}}? *)
        lmax=Max[mymodes[[All,1]]];
        ls=DeleteDuplicates[mymodes[[All,1]]];
        lsms=CombineColumns[ls,ls];
        lsmsm1=CombineColumns[ls,ls-1];
        dommodes=Sort[Join[lsms,lsmsm1]];
        domsneg=dommodes/. {xx_,yy_}->{xx,-yy};
        dommodes=Sort[Join[dommodes,domsneg]];
        testdommodes=(MemberQ[mymodes,#1]&)/@dommodes;
        posdom=Flatten[Position[testdommodes,True]];
        dommodes=dommodes[[posdom]];Print["Dominant modes : ",dommodes];
    
        ta=Table[
                (* Load and process data *)
                Clear[ninjlm];            
                Print["Extraction radius : ",myrads[[j]]];
                
                Do[myReadNinjaModes[dir,FileNameTake[metadata[[i]],-1],"DataSection"-> "psi4t-data","LMax" -> 2,"MModes" -> {mymodes[[l,1]],mymodes[[l,2]]}, "SectionInstance" -> j,
                   "Tag"-> "dummy"];,{l,Length@mymodes}];
                Print["Modes Loaded"];
      
                (* 22 mode guess frequency *)
                test\[Psi]=TimeUnion@ninjlm["dummy",2,2];
                f0=Abs@First@GuessFFIf0[test\[Psi], "MinTime"-> First@First@test\[Psi]+250];
                times=TakeColumn[test\[Psi],1];
    
                (* compute Luminosoties *)
                Do[If[mymodes[[l,2]]!=0,mm=mymodes[[l,2]]/2,mm=(mymodes[[l,2]]+1)/2];
                   test\[Psi]=TimeUnion@ninjlm["dummy",mymodes[[l,1]],mymodes[[l,2]]]; 
                   lum[mymodes[[l,1]],mymodes[[l,2]]]=Abs[TakeColumn[NewsFromPsi4FFI[test\[Psi],0.75f0 Abs[mm]],2]]^2;Print[mymodes[[l]]];,{l,Length[mymodes]}];
               
                If[precessing[tags[[i]]],
                             sumModes=1/(16 \[Pi])Sum[lum[mymodes[[l,1]],mymodes[[l,2]]],{l,Length@mymodes}];
                             sumModesDominant=1/(16 \[Pi])Sum[lum[dommodes[[l,1]],dommodes[[l,2]]],{l,Length@dommodes}];
                             ,
                             sumModes=1/(8 \[Pi])Sum[lum[mymodes[[l,1]],mymodes[[l,2]]],{l,Length@mymodes}];
                             sumModesDominant=1/(8 \[Pi])Sum[lum[dommodes[[l,1]],dommodes[[l,2]]],{l,Length@dommodes}];];
    
    
                 sumModesLum=CombineColumns[times,sumModes];
                 sumModesLumDom=CombineColumns[times,sumModesDominant];
                 tmaxSum=TimeOfMaximum[sumModesLum];
                 tmaxSumDom=TimeOfMaximum[sumModesLumDom];
                 modesMaxima=Table[data=CombineColumns[times,lum[mymodes[[l,1]],mymodes[[l,2]]]];
                                   int=Interpolation[data];
                                   {ValueOfMaximum[data],int[tmaxSum],int[tmaxSumDom]},{l,Length[mymodes]}];
    
    If[verbose,Print[ListPlot[{sumModesLum,sumModesLumDom},Joined->True,PlotRange->All]];];
    Flatten[{tags[[i]],\[Eta][[i]],\[Chi]1[[i]],\[Chi]2[[i]],\[Chi]1p[[i]],\[Chi]2p[[i]],Initial\[Chi]Plane[[i]],myrads[[j]],ValueOfMaximum[sumModesLum],tmaxSum,ValueOfMaximum[sumModesLumDom],tmaxSumDom,Flatten[modesMaxima]}],{j,Length@myrads}];
    
    Print["Exporting results to : ",myFile];
    Export[myFile,ta];
    
    ,{i,Length[tags]}];
    
    ta
    ]
    
    
    AHBAMsmall[q_,s_]:=Sqrt[0.7599216762029171` -0.6443660299643744` s^2]/(1+q)
    AHBAMbig[q_,s_]:=(q Sqrt[0.9466225527668473` -0.8578738205662785` s^2])/(1+q)
    
    
    AHBAMsmall2017[q_,s_]:=(0.9097788740970632` Sqrt[1-0.8630715321085808` s^2]-0.013097939946145838` q Sqrt[1-0.8630715321085808` s^2])1/(1+q)
    AHBAMbig2017[q_,s_]:=(0.9092695479821113` Sqrt[1-0.9004695426988402` s^2]+0.025102624533684333` q Sqrt[1-0.9004695426988402` s^2])q/(1+q)
    
    
    End[];
    EndPackage[];