#!/usr/bin/python

from DiscreteFunction import *
from WaveFunction import *
from DataAnalysis import *

# Maximum l-mode to process
lmax = 8
# FFI cutoff frequency. This must be choosen
# smaller than any physically expected frequency.
f0 = 0.03/(2*pi)

# Initialize a mode array for psi4
WF = InitModeArray(lmax)
# Initialize a mode array for h
WFint = InitModeArray(lmax)

# Load each psi4-mode into a WaveFunction object and store it in mode array
for l in range(2,lmax+1):
    for m in range(-l, l+1):
	print "(l,m) = ", l, m
	WF[l][m] = WaveFunction([], []) 
	WF[l][m].Load("rPsi4_l"+str(l)+"m"+str(m)+".dat")


# Integrate 2,2 mode using FFI
WFint[2][2] = GethFromPsi4(WF[2][2], f0)
# Get time of merger from maximum of 2,2-amplitude
tmerger = WFint[2][2].Amplitude().FindAbsMaxInterpolated(1e-4)

# Integrate all remaining modes, shift them according to tmerger, and
# store them in a file.
print "tmerger =", tmerger
for l in range(2,lmax+1):
    for m in range(-l, l+1):
	print "(l,m) = ", l, m
	WFint[l][m] = GethFromPsi4(WF[l][m], f0*m*0.5)
	WFint[l][m].x -= tmerger
	fout = open("h_from_rPsi4_l"+str(l)+"m"+str(m)+".dat", 'w')
	fout.write("# 1:time   2:h+   3:hx \n# FFI cut-off: omega = "+str(f0*pi*m)+"\n")
	for ii in range(0, WFint[l][m].Length()):
	    fout.write(str(WFint[l][m].x[ii])+" "+str(WFint[l][m].f[ii].real)+" "+str(WFint[l][m].f[ii].imag)+"\n")
	fout.close()