Manually load and inspect Ylm_depict data

In [1]:

    

# Low-level import 
from numpy import array,loadtxt,linspace,zeros,exp,ones,unwrap,angle,pi

# Setup ipython environment
%load_ext autoreload
%autoreload 2
%matplotlib inline

# Import useful things from kerr
from kerr.formula.ksm2_cw import CW as cwfit
from kerr.formula.ksm2_sc import SC as scfit
from kerr.pttools import leaver_workfunction as lvrwork
from kerr import leaver,rgb
from kerr.models import mmrdns 

#
from nrutils import scsearch,gwylm,gwf

# Setup plotting backend
import matplotlib as mpl
mpl.rcParams['lines.linewidth'] = 0.8
mpl.rcParams['font.family'] = 'serif'
mpl.rcParams['font.size'] = 12
mpl.rcParams['axes.labelsize'] = 20
from matplotlib.pyplot import *


    

    




kerr## Found formula module "ksm2_cw"
kerr## Found formula module "ksm2_sc"
kerr## Found formula module "ksm2_slm_norm"
kerr## Found formula module "mmrdns_amplitudes"
kerr## Found formula module "mmrdns_Mfjf"
The highest level init for nrutils is located at: /Users/book/JOKI/Libs/KOALA/nrutils_dev/nrutils/__init__.pyc

>> Initiating nrutils ...

>> Sub-Packages to be imported:
   -> core
   -> generate
   -> manipulate
   -> tools
>> Please note style conventions:                  
   * lower case function/method/variable names                  
   * no underscore in names unless there are repeated letters, or counfounded syllables                  
   * information is implicitely in time domain unless explicitely stated.                  
   * frequency domain information will start with "fd".

nrutils:

  .core: 
      .basics*
      .basics
      .nrsc
  .generate: 
  .manipulate: 
      .bundlers
      .rotate
  .tools: 
    .unit: 
      .conversion

In [7]:

    

T0 = 10
ll,mm = 3,2
q = 2
data_file_string = '/Users/book/GARREG/Spectroscopy/Ylm_Depictions/NonPrecessing/MULTI_DATA_6/T0_%i_nmax2_Mmin97ll_Mmin_r75_qref1.50__p1_17-Mar-2014gnx_/Data_Sets/HRq-series/D9_q%1.1f_a0.0_m160/DEPICTION_INFO::NODD_INPUT_ll%i_mm%i_r75.asc'%(T0,q,ll,mm)
data = loadtxt(data_file_string)


    

In [8]:

    

# Collect raw fit data for later convenience
rfdata = {}
for k,row in enumerate(data):
    #
    ll,mm,q,m1,m2,x1,x2,jf,Mf,qid,rew,imw,rewfit,imwfit,reA,imA,reAmean,imAmean,minA,maxA,T1,dT,match,rmse,reB,imB,reBmean,imBmean,minB,maxB = row
    ll,mm = int(ll),int(mm)
    A = reA+1j*imA
    cw = rew + 1j*imw
    try:
        l,m,n,p = mmrdns.calc_z(qid)
    except:
        l,m,n,p,l2,m2,n2,p2 = mmrdns.calc_z(qid)
    rfdata[(l,m,n,p)] = {}
    rfdata[(l,m,n,p)]['ll'],rfdata[(l,m,n,p)]['mm'],rfdata[(l,m,n,p)]['A'],rfdata[(l,m,n,p)]['cw'] = ll,mm,A,cw


    

In [9]:

    

print angle( rfdata[(mm,mm,0,1)]['A'] * rfdata[(ll,mm,0,1)]['A'].conj() )
print angle( mmrdns.Afit(mm,mm,0,mmrdns.q2eta(q))  * mmrdns.Afit(ll,mm,0,mmrdns.q2eta(q)).conj() )


    

    




-0.758233993633
0.727670507856

In [10]:

    

#
def rawfit(t):
    y = zeros( t.shape, dtype=complex )
    for k,row in enumerate(data):
        #
        ll,mm,q,m1,m2,x1,x2,jf,Mf,qid,rew,imw,rewfit,imwfit,reA,imA,reAmean,imAmean,minA,maxA,T1,dT,match,rmse,reB,imB,reBmean,imBmean,minB,maxB = row
        ll,mm = int(ll),int(mm)
        A = reA+1j*imA
        cw = rew + 1j*imw
        try:
            l,m,n,p = mmrdns.calc_z(qid)
        except:
            l,m,n,p,l2,m2,n2,p2 = mmrdns.calc_z(qid)
        # NOTE that the amplitudes are for Psi4 here
        if True: # (l,m,n,p) in [ (2,2,0,1) ,(2,2,1,1) ]  :
            y += A*exp( 1j*cw*(t-T0) )
            
    #
    a = gwf( array( [t,y.real,-y.imag] ).T )
    
    #
    return a,q

_,q = rawfit( linspace(T0,50) )

#
A = scsearch( keyword='sxs', q=q, nonspinning=True,verbose=True )[0]

#
imrnr = gwylm( A, lm=([ll,mm],[2,2]), verbose=True, dt=0.5 )
nr = imrnr.ringdown(T0=T0)

y,_ = rawfit( nr.lm[(ll,mm)]['psi4'].t )

#
eta = mmrdns.q2eta(q)
h = mmrdns.meval_spherical_mode(ll,mm,eta,kind='psi4',gwfout=True)(nr.ylm[0].t)

h.align(nr.lm[(ll,mm)]['psi4'],method='average-phase',mask=nr.ylm[0].t<60)
y.align(nr.lm[(ll,mm)]['psi4'],method='average-phase',mask=nr.ylm[0].t<60)

nr.lm[(ll,mm)]['psi4'].plot()
y.plot()
h.plot()

fig = figure( figsize=2*array([5,3]) )
gca().set_yscale("log", nonposy='clip')

plot( nr.ylm[0].t, nr.lm[(ll,mm)]['psi4'].amp, color=0.5*ones((3,)), label=None )
plot( nr.ylm[0].t, y.amp, '--k', label=None )
plot( nr.ylm[0].t, h.amp, 'k', alpha=0.2, linewidth=6, label=None )

plot( nr.ylm[0].t, nr.lm[(ll,mm)]['psi4'].plus, color=0.5*ones((3,)), label='NR' )
plot( nr.ylm[0].t, y.plus, '--k', label='RAW-FIT' )
plot( nr.ylm[0].t, h.plus, 'k', alpha=0.2, linewidth=6,label='MMRDNS' )

# plot( nr.ylm[0].t, nr.lm[(ll,mm)]['psi4'].cross, color=0.5*ones((3,)), label='NR', alpha=0.8 )
# plot( nr.ylm[0].t, y.cross, '--k', label='RAW-FIT', alpha=0.8 )
# plot( nr.ylm[0].t, h.cross, 'k', alpha=0.1, linewidth=6, label='MMRDNS' )

ylim( [max(nr.lm[(ll,mm)]['psi4'].amp)*1e-5,1.2*max(nr.lm[(ll,mm)]['psi4'].amp)] )
xlim( [T0,150] )

xlabel(r'$(t-{t^{\mathrm{Peak}}}_{\mathrm{Lum.}})/M$')
ylabel(r'${rM}\psi_{%i%i}$'%(ll,mm))

legend(frameon=False)
title( nr.label )

savefig('mmrdns_psi4_comparison_%s_ll%imm%i.pdf'%(nr.label.replace('-','_'),ll,mm))

# gca().set_yscale("log", nonposy='clip')


    

    




[scsearch]>> Found keyword (='sxs') keyword.
[scsearch]>> Found nonspinning (=True) keyword.
[scsearch]>> Found q (=2.000003) keyword.
[scsearch]>> Found validate_remnant (=False) keyword.
[scsearch]>> Found verbose (=True) keyword.
(scsearch)>> List of keywords or string keyword found: ALL scentry objects matching will be passed. To pass ANY entries matching the keywords, input the keywords using an iterable of not of type list.
## Found 2 possibly degenerate simulations:
[0001][sxs] SXS0184: qc-ns-q2.00	(SXS0184)
[0002][sxs] SXS0169: qc-ns-q2.00	(SXS0169)

(gwylm)>> Found clean (=False) keyword.
(gwylm)>> Found dt (=0.5) keyword.
(gwylm)>> Found lm (=([3, 2], [2, 2])) keyword.
(gwylm)>> Found load (=True) keyword.
(gwylm)>> Found lowpass (=False) keyword.
(gwylm)>> Found scentry_obj (=<nrutils.core.nrsc.scentry instance at 0x109bddf38>) keyword.
(gwylm)>> Found verbose (=True) keyword.
(load)>> Loading: rMPsi4_Y_l3_m2.asc
(gwylm.setfields)>> Interpolating data to dt=0.500000
(gwylm.setfields)>> Interpolating data to dt=0.500000
(load)>> Loading: rMPsi4_Y_l2_m2.asc
(gwylm.setfields)>> Interpolating data to dt=0.500000
(gwylm.setfields)>> Interpolating data to dt=0.500000
(gwylm)>> Using w22 from a PN estimate to calculate strain multipoles [see pnw0 in basics.py, and/or arxiv:1310.1528v4].
* w0(w22) = 0.031781 (this is the lower frequency used for FFI method [arxiv:1006.1632v3])
(gwylm.calchlm)>> The user should note that there is no minus sign used in front of the double time integral for strain (i.e. Eq 4 of arxiv:1006.1632). This differs from Eq 3.4 of arxiv:0707.4654v3. The net effect is a rotation of the overall polarization of pi degrees. The user should also note that there is no minus sign applied to h_cross meaning that the user must be mindful to write h_pluss-1j*h_cross when appropriate.
* w0(w22) = 0.031781 (this is the lower frequency used for FFI method [arxiv:1006.1632v3])
(gwylm.calchlm)>> The user should note that there is no minus sign used in front of the double time integral for strain (i.e. Eq 4 of arxiv:1006.1632). This differs from Eq 3.4 of arxiv:0707.4654v3. The net effect is a rotation of the overall polarization of pi degrees. The user should also note that there is no minus sign applied to h_cross meaning that the user must be mindful to write h_pluss-1j*h_cross when appropriate.
* w0(w22) = 0.031781 (this is the lower frequency used for FFI method [arxiv:1006.1632v3])
* w0(w22) = 0.031781 (this is the lower frequency used for FFI method [arxiv:1006.1632v3])






    




/Users/book/JOKI/Libs/KOALA/nrutils_dev/nrutils/core/basics.py:1028: RuntimeWarning: invalid value encountered in subtract
  iNeq = where(  ( abs(yTemp[1:]-yTemp[:-1])>1e-12 )  *  ( yFinite[:-1]+yFinite[1:] )  )






    













    













    













    

In [ ]:

    

nr.simdir


    

In [ ]: