Select Git revision
-
Oliver Bock authoredOliver Bock authored
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[];