import CactusTool

name = 'BBH000'
Level = ['Lev3', 'Lev4', 'Lev5']
n = len(Level)

mp = []
bbh = []
carpet = []
for l in Level:
    sim = CactusTool.load(name+l, '/Volumes/simulations/BBH_Catalog/')
    mp.append(sim.ThornOutput('multipole'))
    bbh.append(sim.ThornOutput('TwoPunctures'))
    carpet.append(sim.ThornOutput('carpet').dx)

Q = 6
h = [dx[len(dx)-1][0] for dx in carpet]
factor = round((h[0]**Q - h[1]**Q) / (h[1]**Q - h[2]**Q), 2)
factor

2.02

import matplotlib.pyplot as plt
import numpy as np

rPsi4 = []
for i in range(n):
    rPsi4.append(mp[i].rPsi4(bbh[i].ADMMass, (2,2), 500))

tstart = max([v.tstart for v in rPsi4])
tend = min([v.tend for v in rPsi4])
dt = min([v.dt for v in rPsi4])
t = np.arange(tstart, tend, dt)
rpsi4 = [v.resample(t).real.y for v in rPsi4]

fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)
plt.xlabel(r'$(t-r*)~[M]$')
plt.ylabel(r'$Re[r\Psi_{4}^{(2,2)}]$')

for i in range(n-1):
    plt.plot(t, np.abs(rpsi4[i+1]-rpsi4[i]), alpha = 0.7, label='|{}-{}|'.format(Level[i+1], Level[i]))

plt.plot(t, factor * np.abs(rpsi4[2]-rpsi4[1]), ':', alpha = 0.7,label='{}|{}-{}|'.format(factor, Level[2], Level[1]))

ax.set_yscale('log')
# plt.ylim(1e-6,1e-2)
# plt.xlim(0,1500)
plt.legend()
plt.title(name)

plt.show()

import matplotlib.pyplot as plt
import numpy as np

Strain = []
for i in range(n):
    Strain.append(mp[i].Strain(bbh[i].ADMMass, bbh[i].CutoffFrequency, (2,2), -1))

tstart = max([v.tstart for v in Strain])
tend = min([v.tend for v in Strain])
dt = min([v.dt for v in Strain])
t = np.arange(tstart, tend, dt)

phase = [v.phase.resample(t).y for v in Strain]
amp = [v.abs().resample(t).y for v in Strain]

fig = plt.figure(figsize=(10, 12))
ax = fig.add_subplot(211)
ax.set_xlabel(r'$(t-r*)~[M]$')
ax.set_ylabel(r'$\Delta \phi^{(2,2)}$')
ax.set_yscale('log')

for i in range(n-1):
    plt.plot(t, np.abs(phase[i+1]-phase[i]), alpha = 0.7, label='|{}-{}|'.format(Level[i+1], Level[i]))

plt.plot(t, factor * np.abs(phase[2]-phase[1]), ':', alpha = 0.7,label='{}|{}-{}|'.format(factor, Level[2], Level[1]))

plt.legend()
plt.title(name)

ax = fig.add_subplot(212)
ax.set_xlabel(r'$(t-r*)~[M]$')
ax.set_ylabel(r'$\Delta |h|^{(2,2)}$')
ax.set_yscale('log')

for i in range(n-1):
    plt.plot(t, np.abs(amp[i+1]-amp[i]), alpha = 0.7, label='|{}-{}|'.format(Level[i+1], Level[i]))

plt.plot(t, factor * np.abs(amp[2]-amp[1]), ':', alpha = 0.7,label='{}|{}-{}|'.format(factor, Level[2], Level[1]))

plt.legend()
plt.title(name)

plt.show()

import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)
ax.set_xlabel('t [M]')
ax.set_ylabel(r'$h^{(2,2)}_{+}$')

for i in range(n):
    Strain = mp[i].Strain(bbh[i].ADMMass, bbh[i].CutoffFrequency, (2,2), -1).real
    plt.plot(Strain.t, Strain.y, alpha = 0.7, label='{}'.format(Level[i]))

# plt.xlim(900, 1000)
plt.legend()