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: