[Commits] [svn:einsteintoolkit] pyGWAnalysis/trunk/ (Rev. 15)
reisswig at tapir.caltech.edu
reisswig at tapir.caltech.edu
Mon Nov 12 05:00:27 CST 2012
User: reisswig
Date: 2012/11/12 05:00 AM
Modified:
/trunk/DataAnalysis/
DataAnalysis.py, RadiatedQuantities.py
/trunk/DiscreteFunction/
DiscreteFunction.py, WaveFunction.py
Log:
1) Radiated energy flux from RWZM.
2) Some operators for DiscreteFunction
3) Load mode array by specifying file format.
File Changes:
Directory: /trunk/DataAnalysis/
===============================
File [modified]: DataAnalysis.py
Delta lines: +21 -5
===================================================================
--- trunk/DataAnalysis/DataAnalysis.py 2012-05-11 00:23:43 UTC (rev 14)
+++ trunk/DataAnalysis/DataAnalysis.py 2012-11-12 11:00:27 UTC (rev 15)
@@ -427,20 +427,36 @@
res[-1].append([])
return res
-def LoadModeArray(basename, lmax, col_time=1, col_re=2, col_im=3):
+def LoadModeArray(basename, lmax, col_time=1, col_re=2, col_im=3, filenameformat="cce_conversionscript"):
"""Loads an array of modes up to lmax from disk.
We assume that files begin with 'basename' and end with _l*m*.dat"""
res = InitModeArray(lmax)
for ll in range(2, lmax+1):
for mm in range(-ll, ll+1):
- res[ll][mm] = WaveFunction([], [])
- res[ll][mm].Load(basename+"_l"+str(ll)+"m"+str(mm)+".dat", col_time, col_re, col_im)
-
+ if (filenameformat == "cce_conversion_script"):
+ res[ll][mm] = WaveFunction([], [])
+ res[ll][mm].Load(basename+"_l"+str(ll)+"m"+str(mm)+".dat", col_time, col_re, col_im)
+ continue
+ if (filenameformat == "pittnull"):
+ res[ll][-mm] = WaveFunction([], [])
+ # Get spaces from Fortran outpout filenames right
+ tmp1 = "%02d" % ll
+ if mm < 0:
+ tmp2 = "m"
+ else:
+ tmp2 = "p"
+ tmp2 += "%02d" % abs(mm)
+ # Take into account that nullcode has a weird convention for psi4!
+ res[ll][-mm].Load(basename+".L"+tmp1+"M"+tmp2+".asc", col_time, col_re, col_im)
+ res[ll][-mm].f *= 2.0*(-1.0)**(mm+1)
+ res[ll][-mm].f = map(lambda x: x.conjugate(), res[ll][-mm].f)
+ continue
+ print("Error loading mode array: filenameformat not recognized!")
+
return res
-
def ReconFromsYlm(arrayWF, theta, phi, lmax):
"""Given an array of WF modes h_lm reconstruct h at angle theta, phi up to l=lmax"""
res = WaveFunction(arrayWF[2][0].x, arrayWF[2][0].f*0.0)
File [modified]: RadiatedQuantities.py
Delta lines: +24 -0
===================================================================
--- trunk/DataAnalysis/RadiatedQuantities.py 2012-05-11 00:23:43 UTC (rev 14)
+++ trunk/DataAnalysis/RadiatedQuantities.py 2012-11-12 11:00:27 UTC (rev 15)
@@ -39,6 +39,30 @@
+def RadiatedEnergyFluxRWZM(Qeven, Qodd, lmax):
+ """ Given the time integral of rPsi4 (l,m) modes (a DiscreteFunction mode array),
+ this function computes the radiated energy flux
+ that is emitted as a function of time in units of [M]."""
+
+ npoints = Qeven[2][0].Length()
+ Edot = DiscreteFunction(Qeven[2][0].x, zeros(npoints, float))
+
+ dtQeven = InitModeArray(2)
+ for l in range(2, lmax+1):
+ for m in range(-l, l+1):
+ dtQeven[l][m] = Qeven[l][m].DifferentiateFunction()
+
+ for ii in range(Edot.Length()):
+ e = 0.0
+ for l in range(2, lmax+1):
+ for m in range(-l, l+1):
+ e = e + abs(dtQeven[l][m].f[ii])**2 + abs(Qodd[l][m].f[ii])**2
+ Edot.f[ii] = e/(32.*pi)
+
+ return Edot
+
+
+
def RadiatedEnergyFlux(int_rPsi4_ModeArray, lmax):
""" Given the time integral of rPsi4 (l,m) modes (a DiscreteFunction mode array),
this function computes the radiated energy flux
Directory: /trunk/DiscreteFunction/
===================================
File [modified]: DiscreteFunction.py
Delta lines: +68 -2
===================================================================
--- trunk/DiscreteFunction/DiscreteFunction.py 2012-05-11 00:23:43 UTC (rev 14)
+++ trunk/DiscreteFunction/DiscreteFunction.py 2012-11-12 11:00:27 UTC (rev 15)
@@ -32,10 +32,13 @@
return False
- def __init__(self, x_, f_):
+ def __init__(self, x_, f_, g_=None):
"""Initialise a 1D transformation function"""
self.x = array(x_)
- self.f = array(f_)
+ if g_ != None:
+ self.f = array(f_) + 1j*array(g_)
+ else:
+ self.f = array(f_)
if self.HasUniformSpacing() == True:
self.dx = self.x[1]-self.x[0]
else:
@@ -383,12 +386,31 @@
return xroot
+ def DifferentiateFunction(self):
+ """ Differentiate this function and return a function"""
+
+ npoints = self.Length()
+ res = self.__class__(self.x, self.f*0.0)
+
+ # TODO: fail if not uniform spacing
+ self.dx = self.x[1]-self.x[0]
+ res.dx = self.dx
+
+ # Simple second order centred differentiation
+
+ for ii in range (1, self.Length()-1):
+ res.f[ii] = (self.f[ii+1] - self.f[ii-1])/(2*res.dx)
+
+ return res
+
+
def IntegrateFunction(self):
""" Integrate this function and return a function"""
npoints = self.Length()
res = self.__class__(self.x, self.f*0.0)
+ # TODO: fail if not uniform spacing
self.dx = self.x[1]-self.x[0]
res.dx = self.dx
@@ -572,6 +594,12 @@
return res
+ def __add__(self, rval):
+ res = self.Union(rval)
+ for ii in range(0, res.Length()):
+ res.f[ii] = self.Interpolate(res.x[ii]) + rval.Interpolate(res.x[ii])
+ return res
+
def __sub__(self, rval):
res = self.Union(rval)
for ii in range(0, res.Length()):
@@ -591,7 +619,45 @@
res.f[ii] = self.Interpolate(res.x[ii]) / rval.Interpolate(res.x[ii])
return res
+ def __abs__(self):
+ npoints = len(self.x)
+ res = self.__class__(self.x, zeros(npoints, float))
+ for ii in range(0, Length()):
+ res.f[ii] = abs(self.f[ii])
+ return res
+ def __pos__(self):
+ npoints = len(self.x)
+ res = self.__class__(self.x, self.f)
+ return res
+
+ def __neg__(self):
+ npoints = len(self.x)
+ res = self.__class__(self.x, -self.f)
+ return res
+
+ def real(self):
+ npoints = len(self.x)
+ res = self.__class__(self.x, zeros(npoints, float))
+ for ii in range(0, self.Length()):
+ res.f[ii] = self.f[ii].real
+ return res
+
+ def imag(self):
+ npoints = len(self.x)
+ res = self.__class__(self.x, zeros(npoints, float))
+ for ii in range(0, self.Length()):
+ res.f[ii] = self.f[ii].imag
+ return res
+
+ def conjugate(self):
+ npoints = len(self.x)
+ res = self.__class__(self.x, zeros(npoints, type(self.f[0])))
+ for ii in range(0, self.Length()):
+ res.f[ii] = self.f[ii].conjugate()
+ return res
+
+
# implement iterator
def __iter__(self):
self.index = 0
File [modified]: WaveFunction.py
Delta lines: +5 -2
===================================================================
--- trunk/DiscreteFunction/WaveFunction.py 2012-05-11 00:23:43 UTC (rev 14)
+++ trunk/DiscreteFunction/WaveFunction.py 2012-11-12 11:00:27 UTC (rev 15)
@@ -19,10 +19,13 @@
class WaveFunction(DiscreteFunction):
"""Class describing discrete 1d wave functions"""
- def __init__(self, x_, f_):
+ def __init__(self, x_, f_, g_=None):
"""Initialise a 1D wave function"""
self.x = array(x_)
- self.f = array(f_)
+ if g_ != None:
+ self.f = array(f_) + 1j*array(g_)
+ else:
+ self.f = array(f_)
if self.HasUniformSpacing() == True:
self.dx = self.x[1]-self.x[0]
else:
More information about the Commits
mailing list