Skip to content
Snippets Groups Projects
Commit 41e41181 authored by Alejandra Gonzalez's avatar Alejandra Gonzalez
Browse files

Added and adapted functions to deal with infinite radius, plus some other minor corrections

parent 27de8773
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
watpy/utils/coreh5.py
+ 54
− 52
View file @ 41e41181
......@@ -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]
......
watpy/wave/gwutils.py
+ 4
− 1
View file @ 41e41181
......@@ -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
......
watpy/wave/wave.py
+ 84
− 17
View file @ 41e41181
......@@ -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'):
"""
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment