diff --git a/watpy/utils/coreh5.py b/watpy/utils/coreh5.py
index 60f33486115c3bfff962221cffcbc75ce19507f7..459142f282fb92c05406a01f286915967b8b6cbe 100644
--- a/watpy/utils/coreh5.py
+++ b/watpy/utils/coreh5.py
@@ -12,9 +12,18 @@ import os.path
import re
import numpy as np
-from ..wave.wave import wfile_parse_name
+from ..wave.wave import wfile_parse_name, rinf_float_to_str, rinf_str_to_float, rInf, write_headstr
from .viz import wplot
+def write_keyh5(l, m, r):
+ """
+ Writes key string 'l#_m#_r#'
+ """
+ if(r==rInf):
+ key = 'l{}_m{}_rInf'.format(l,m)
+ else:
+ key = 'l{}_m{}_r{:05d}'.format(l,m,int(r))
+ return key
class CoRe_h5():
"""
@@ -140,24 +149,25 @@ class CoRe_h5():
rad = self.dset_radii(fn,group=group, det=det)
if group.startswith('rh_'):
l,m = self.lm_from_group(group)
- if(rad==1e6):
- dset = fn[group]['Rh_l{}_m{}_rInf.txt'.format(l,m)][()]
- else:
- dset = fn[group]['Rh_l{}_m{}_r{:05d}.txt'.format(l,m,int(rad))][()]
+ key = write_keyh5(l,m,rad)
+ filename = 'Rh_'+key+'.txt'
+ dset = fn[group][filename][()]
+ #dset = fn[group]['Rh_l{}_m{}_r{:05d}.txt'.format(l,m,int(rad))][()]
elif group.startswith('rpsi4_'):
l,m = self.lm_from_group(group)
- if(rad==1e6):
- dset = fn[group]['Rh_l{}_m{}_rInf.txt'.format(l,m)][()]
- else:
- dset = fn[group]['Rpsi4_l{}_m{}_r{:05d}.txt'.format(l,m,int(rad))][()]
+ key = write_keyh5(l,m,rad)
+ filename = 'Rpsi4_'+key+'.txt'
+ dset = fn[group][filename][()]
+ #dset = fn[group]['Rpsi4_l{}_m{}_r{:05d}.txt'.format(l,m,int(rad))][()]
elif group.startswith('EJ_'):
- if(rad==1e6):
- dset = fn[group]['EJ__rInf.txt'][()]
- else:
- dset = fn[group]['EJ__r{:05d}.txt'.format(int(rad))][()]
+ rad_str = rinf_float_to_str(rad)
+ filename = 'EJ_r'+rad_str+'.txt'
+ dset = fn[group][filename][()]
+ #dset = fn[group]['EJ__r{:05d}.txt'.format(int(rad)][()]
else:
raise ValueError("Unknown group {}".format(group))
return np.array(dset)
+
def dump(self):
"""
@@ -182,11 +192,9 @@ class CoRe_h5():
"""
radii = []
for ds in fp[group].keys():
- rad_str = ds[-8:-4]
- if(rad_str=='rInf'):
- radii = np.append(radii,1e6)
- else:
- radii = np.append(radii,float(ds[-8:-4]))
+ rad = rinf_str_to_float(ds[-8:-4])
+ radii = np.append(radii,rad)
+ #radii = np.append(radii,float(ds[-8:-4]))
if det in radii:
return det
else:
@@ -203,12 +211,10 @@ class CoRe_h5():
group = 'rh_{}{}'.format(l,m)
if group not in fn.keys(): continue
for f in fn[group]:
- rad_str = f[-8:-4]
- if(rad_str=='rInf'):
- rad = 1e6 #for simulations extrapolated at infinity
- else:
- rad = float(f[-8:-4])
- headstr = "r=%e\nM=%e\n " % (rad, mass)
+ #rad = float(f[-8:-4])
+ rad = rinf_str_to_float(f[-8:-4])
+ #headstr = "r=%e\nM=%e\n " % (rad, mass)
+ headstr = write_headstr(rad,mass)
headstr += "u/M:0 Reh/M:1 Imh/M:2 Redh/M:3 Imdh/M:4 Momega:5 A/M:6 phi:7 t:8"
dset = fn[group][f]
data = np.c_[dset[:,0],dset[:,1],dset[:,2],
@@ -225,24 +231,22 @@ class CoRe_h5():
mass = float(self.mdata.data['id_mass'])
with h5py.File(os.path.join(self.path,self.dfile), 'r') as fn:
for l,m in lm:
- group = 'rh_{}{}'.format(l,m)
+ group = 'rpsi4_{}{}'.format(l,m)
if group not in fn.keys(): continue
for f in fn[group]:
- rad_str = f[-8:-4]
- if(rad_str=='rInf'):
- rad = 1e6 #for simulations extrapolated at infinity
- else:
- rad = float(f[-8:-4])
- headstr = "r=%e\nM=%e\n " % (rad, mass)
+ #rad = float(f[-8:-4])
+ #headstr = "r=%e\nM=%e\n " % (rad, mass)
+ rad = rinf_str_to_float(f[-8:-4])
+ headstr = write_headstr(rad,mass)
dset = fn[group][f]
- try:
- headstr += "u/M:0 RePsi4/M:1 ImPsi4/M:2 Momega:3 A/M:4 phi:5 t:6"
+ try:
data = np.c_[dset[:,0],dset[:,1],dset[:,2],
- dset[:,3],dset[:,4],dset[:,5],dset[:,6]]
+ dset[:,3],dset[:,4],dset[:,5],dset[:,6]]
+ headstr += "u/M:0 RePsi4/M:1 ImPsi4/M:2 Momega:3 A/M:4 phi:5 t:6"
except:
- headstr += "u/M:0 RePsi4/M:1 ImPsi4/M:2 t:4"
data = np.c_[dset[:,0],dset[:,1],dset[:,2],
- dset[:,3]]
+ dset[:,3]]
+ headstr += "u/M:0 RePsi4/M:1 ImPsi4/M:2 t:4"
np.savetxt(os.path.join(self.path,f),
data, header=headstr)
@@ -258,21 +262,19 @@ class CoRe_h5():
print("No group {}".format(group))
return
for f in fn[group]:
- rad_str = f[-8:-4]
- if(rad_str=='rInf'):
- rad = 1e6 #for simulations extrapolated at infinity
- else:
- rad = float(f[-8:-4])
- headstr = "r=%e\nM=%e\n " % (rad, mass)
+ #rad = float(f[-8:-4])
+ #headstr = "r=%e\nM=%e\n " % (rad, mass)
+ rad = rinf_str_to_float(f[-8:-4])
+ headstr = write_headstr(rad,mass)
dset = fn['energy'][f]
try:
- headstr += "J_orb:0 E_b:1 u/M:2 E_rad:3 J_rad:4 t:5"
data = np.c_[dset[:,0],dset[:,1],dset[:,2],
dset[:,3],dset[:,4],dset[:,5]]
+ headstr += "J_orb:0 E_b:1 u/M:2 E_rad:3 J_rad:4 t:5"
except:
- headstr += "J_orb:0 E_b:1 u/M:2 E_rad:3 J_rad:4"
data = np.c_[dset[:,0],dset[:,1],dset[:,2],
dset[:,3],dset[:,4]]
+ headstr += "J_orb:0 E_b:1 u/M:2 E_rad:3 J_rad:4"
np.savetxt(os.path.join(self.path,f),
data, header=headstr)
@@ -300,16 +302,16 @@ class CoRe_h5():
rad = self.dset_radii(fn,group=group, det=det)
if group.startswith('rh_'):
l,m = self.lm_from_group(group)
- if(rad==1e6):
- dset = fn[group]['Rh_l{}_m{}_rInf.txt'.format(l,m)]
- else:
- dset = fn[group]['Rh_l{}_m{}_r{:05d}.txt'.format(l,m,int(rad))]
+ subkey = write_keyh5(l,m,rad)
+ key = 'Rh_'+subkey+'.txt'
+ dset = fn[group][key]
+ #dset = fn[group]['Rh_l{}_m{}_r{:05d}.txt'.format(l,m,int(rad))]
elif group.startswith('rpsi4_'):
l,m = self.lm_from_group(group)
- if(rad==1e6):
- dset = fn[group]['Rpsi4_l{}_m{}_rInf.txt'.format(l,m)]
- else:
- dset = fn[group]['Rpsi4_l{}_m{}_r{:05d}.txt'.format(l,m,int(rad))]
+ subkey = write_keyh5(l,m,rad)
+ key = 'Rpsi4_'+subkey+'.txt'
+ dset = fn[group][key]
+ #dset = fn[group]['Rpsi4_l{}_m{}_r{:05d}.txt'.format(l,m,int(rad))]
else:
raise ValueError("Group {} is now a waveform".format(group))
x = dset[:,0]
diff --git a/watpy/wave/gwutils.py b/watpy/wave/gwutils.py
index b433f282dbfb7192b1ebc8b7a90646972cd3aa48..0cdc98fd6c2122f3debd21acd44ad58d2206a6a0 100644
--- a/watpy/wave/gwutils.py
+++ b/watpy/wave/gwutils.py
@@ -404,7 +404,10 @@ def ret_time(t,r,M=1.):
"""
Retarded time on Schwarzschild
"""
- rs = r + 2.*M*np.log(r/(2.*M) -1.)
+ if(r==1.0 or r==-1.0):
+ rs = 0 # For the case when r=-1 (extrapolated at infinity)
+ else:
+ rs = r + 2.*M*np.log(r/(2.*M) -1.)
return t - rs
diff --git a/watpy/wave/wave.py b/watpy/wave/wave.py
index 4f4f21413e275cdea3cc1f3418d0832ac08cff59..d52ee70f8862b726cadd87f2fc008303c2f875f0 100644
--- a/watpy/wave/wave.py
+++ b/watpy/wave/wave.py
@@ -2,12 +2,55 @@ from ..utils.ioutils import *
from ..utils.num import diff1, diffo
from ..utils.viz import wplot
from .gwutils import fixed_freq_int_2, waveform2energetics, ret_time
+import numpy as np
# ------------------------------------------------------------------
# Routines for waveform files
# ------------------------------------------------------------------
+# Value for extraction radius at infinity
+rInf = -1.
+
+def write_headstr(radius,mass):
+ if(radius==rInf):
+ headstr = "r=Infinity\nM=%e\n" % (mass)
+ else:
+ headstr = "r=%e\nM=%e\n" % (radius, mass)
+ return headstr
+
+def write_key(l, m, r):
+ """
+ Writes key string 'l#_m#_r#'
+ """
+ if(r==rInf):
+ key = 'l%d_m%d_rInf' % (l, m)
+ else:
+ key = 'l%d_m%d_r%05d' % (l, m, r)
+ return key
+
+def rinf_float_to_str(radius):
+ """
+ Converts detector radius (float) to string
+ """
+ if(radius==rInf):
+ rad_str = 'Inf'
+ else:
+ rad_str = '%05d' % int(radius)
+ return rad_str
+
+
+def rinf_str_to_float(rad_str):
+ """
+ Converts detector radius (string) to float
+ """
+ opt = ['rInf','Inf','Infinity']
+ if(rad_str in opt):
+ radius = rInf
+ else:
+ radius = float(rad_str)
+ return radius
+
def wave_prop_default():
"""
@@ -33,10 +76,11 @@ def wfile_parse_name(fname):
fname : Name of the file to parse for information
"""
#FIXME: negative modes?!
- t = ['bam','cactus','core','core-energy']
+ t = ['bam','cactus','core','core','core-energy']
s = [r'R(\w+)mode(\d)(\d)_r(\d+).l(\d+)',
r'mp_(\w+)_l(\d)_m(\d)_r(\d+\.\d\d).asc',
r'R(\w+)_l(\d+)_m(\d+)_r(\d+).txt',
+ r'R(\w+)_l(\d+)_m(\d+)_r(\w+).txt',
r'EJ_r(\d+).txt']
vlmr = None
for tp, sm in zip(t,s):
@@ -44,14 +88,16 @@ def wfile_parse_name(fname):
if name is not None:
if tp == 'core-energy':
v = 'EJ'
- r = float(name.group(1))
+ #r = float(name.group(1))
+ r = rinf_str_to_float(name.group(1))
vlmr = (v,None,None,r,tp)
return vlmr
else:
v = name.group(1)
l = int(name.group(2))
m = int(name.group(3))
- r = float(name.group(4))
+ #r = float(name.group(4))
+ r = rinf_str_to_float(name.group(4))
vlmr = (v,l,m,r,tp)
return vlmr
return vlmr
@@ -69,7 +115,13 @@ def wfile_get_detrad(fname):
fname : Name of the file to parse for information
"""
s = extract_comments(fname, '#')
- return float(re.findall("\d+\.\d+",s[0])[0])
+ #return float(re.findall("\d+\.\d+",s[0])[0])
+ try:
+ #FIXME: Doesn't identify the exponential, returns only the 1st digit
+ rad_str = re.findall("\d+\.\d+",s[0])[1]
+ except:
+ rad_str = re.findall("\w+",s[0])[1]
+ return rinf_str_to_float(rad_str)
def wfile_get_mass(fname):
@@ -86,7 +138,6 @@ def wfile_get_mass(fname):
# BAM specials
-
def wfile_get_detrad_bam(fname):
"""
Get radius from wf BAM file header
@@ -97,7 +148,12 @@ def wfile_get_detrad_bam(fname):
fname : Name of the file to parse for information
"""
s = extract_comments(fname, '"')
- return float(re.findall("\d+\.\d+",s[0])[0])
+ #return float(re.findall("\d+\.\d+",s[0])[0])
+ try:
+ rad_str = re.findall("\d+\.\d+",s[0])[2]
+ except:
+ rad_str = re.findall("\w+",s[0])[2]
+ return rinf_str_to_float(rad_str)
# Cactus/THC specials
@@ -121,7 +177,9 @@ def cactus_to_core(path, prop):
var = np.array([])
for seg in os.listdir(path):
if re.match(seg_tmpl, seg):
- fpath = '/data/mp_Psi4_l%d_m%d_r%.2f.asc' % (l, m, det)
+ key = write_key(l,m,det)
+ fpath = '/data/mp_Psi4_'+key+'.asc'
+ #fpath = '/data/mp_Psi4_l%d_m%d_r%.2f.asc' % (l, m, det)
raw = np.loadtxt(os.path.join(path,seg+fpath))
t = np.append(t, raw[:,0])
var = np.append(var, raw[:,1]+raw[:,2]*1.0j)
@@ -211,7 +269,9 @@ class mwaves(object):
self.radii.add(r)
#self.dtype.add(tp)
- key = "%s_l%d_m%d_r%.2f" % (var, l, m, r)
+ subkey = write_key(l,m,r)
+ key = var+"_"+subkey
+ #key = "%s_l%d_m%d_r%.2f" % (var, l, m, r)
if key in self.data:
self.data[key].append(fname)
else:
@@ -254,7 +314,9 @@ class mwaves(object):
if m not in self.mmode:
raise ValueError("Unknown m-index {}".format(m))
- key = "%s_l%d_m%d_r%.2f" % (var, l, m, r)
+ #key = "%s_l%d_m%d_r%.2f" % (var, l, m, r)
+ subkey = write_key(l,m,r)
+ key = var+"_"+subkey
return wave(path = self.path, code = self.code, filename = self.data[key][0],
mass = self.mass, f0 = self.f0)
@@ -284,11 +346,13 @@ class mwaves(object):
self.jorb = (jadm - self.j) / (m1*m2)
if path_out:
- headstr = "r=%e\nM=%e\n" % (rad, self.mass)
+ #headstr = "r=%e\nM=%e\n" % (rad, self.mass)
+ headstr = write_headstr(rad,self.mass)
headstr += "J_orb:0 E_b:1 u/M:2 E_rad:3 J_rad:4 t:5"
data = np.c_[self.jorb, self.eb, u[(2,2)]/self.mass,
self.e, self.j, t]
- fname = "EJ_r%05d.txt" % int(rad)
+ rad_str = rinf_float_to_str(rad)
+ fname = "EJ_r"+rad_str+".txt"
np.savetxt('{}/{}'.format(path_out,fname), data, header=headstr)
@@ -390,13 +454,16 @@ class wave(object):
"""
M = self.prop["mass"]
R = self.prop['detector.radius']
- headstr = "r=%e\n" %(self.prop["detector.radius"])
- headstr += "M=%e\n" %(M)
+ #headstr = "r=%e\n" %(self.prop["detector.radius"])
+ #headstr += "M=%e\n" %(M)
+ headstr = write_headstr(R,M)
if var == 'Psi4':
headstr += "u/M:0 RePsi4/M:1 ImPsi4/M:2 Momega:3 A/M:4 phi:5 t:6"
data = np.c_[self.time_ret()/M, self.p4.real*R/M, self.p4.imag*R/M, M*self.phase_diff1(var),
self.amplitude(var)*R/M, self.phase(var), self.time]
- fname = 'Rpsi4_l%d_m%d_r%05d.txt' % (self.prop['lmode'], self.prop['mmode'], R)
+ key = write_key(self.prop['lmode'], self.prop['mmode'], R)
+ fname = 'Rpsi4_'+key+'.txt'
+ #fname = 'Rpsi4_l%d_m%d_r%05d.txt' % (self.prop['lmode'], self.prop['mmode'], R)
else:
headstr += "u/M:0 Reh/M:1 Imh/M:2 Momega:3 A/M:4 phi:5 t:6"
data = np.c_[self.time_ret()/M, self.h.real*R/M, self.h.imag*R/M, M*self.phase_diff1(var),
@@ -444,8 +511,8 @@ class wave(object):
"""
Print the properties
"""
- for key, vak in self.prop.items():
- print(key+":"+val)
+ for key, val in self.prop.items():
+ print(key+":"+str(val))
def prop_set(self, key, val):
"""
@@ -536,7 +603,7 @@ class wave(object):
"""
Retarded time based on tortoise Schwarzschild coordinate
"""
- return ret_time(self.time,self.prop["detector.radius"], self.prop["mass"])
+ return ret_time(self.time,np.abs(self.prop["detector.radius"]), self.prop["mass"])
def data_interp1(self, timei, useu=0, kind='linear'):
"""