Load a BAM hdf5 file, and given an inclination value, determine the angle $\theta_{\mathrm{JN}}$

The goal of this script is to take in INCLINATION and then output $\theta_{\mathrm{JN}}$



In [46]:

    
# Import useful standard things
from os import system, remove, makedirs, path
from os.path import dirname, basename, isdir, realpath
from numpy import array,ones,pi,dot,cos,sin,savetxt
from numpy import arccos as acos
from numpy.linalg import norm
from matplotlib.pyplot import *
from os.path import expanduser
# Import nrutils (See: https://github.com/llondon6/nrutils_dev)
from nrutils import alert,nr2h5,scsearch,gwylm,alert
from nrutils.tools.unit.conversion import *

Use nrutils to locate the simulation of interest.



In [20]:

    
Y = scsearch(keyword='q1.2_dc2dcp2',verbose=True,unique=True)[0]









    



[scsearch]>> Found keyword (='q1.2_dc2dcp2') keyword.
[scsearch]>> Found unique (=True) keyword.
[scsearch]>> Found verbose (=True) keyword.
## Found 1 unique simulations:
[0001] q1.2_dc2dcp2: p-q1.20

Define function for $\theta_{\mathrm{JN}}$ calculation



In [27]:

    
def iota2thetajn(A,                # scnetry object from nrutils
                 IOTA,             # The angle between L and the line of sight 
                 verbose=False):   # be verbose toggle
    
    # Import useful things
    from numpy import array,ones,pi,dot,cos,sin,vstack,hstack
    from numpy import arccos as acos
    from numpy.linalg import norm
    
    # Extract L, S, and J (NOTE: these are after-junk quantities)
    L = A.L1 + A.L2
    S = A.S1 + A.S2 
    J = L + S

    # Calculate the direction of J and L using the loaded information
    info = A.raw_metadata
    L_bbh = array( [ info.initial_angular_momentumx, info.initial_angular_momentumy, info.initial_angular_momentumz ] )
    S1 = array( [ info.after_junkradiation_spin1x, info.after_junkradiation_spin1y, info.after_junkradiation_spin1z ] )
    S2 = array( [ info.after_junkradiation_spin2x, info.after_junkradiation_spin2y, info.after_junkradiation_spin2z ] )
    S = S1+S2
    J_bbh = L_bbh+S

    # The magnitude of For now, apply norm(L_bbh)

    if False: # verbose:
        print 'L from the bbh files is %s' % L_bbh
        print 'L after junk from nrutils is %s' % L
        print 40*'--'

    # Find unit vectors
    L_hat = L / norm(L)
    J_hat = J / norm(J)

    # Extract oribital separation (NOTE: this is also referenced after-junk)
    n_hat = ( A.R1 - A.R2 ) / norm( A.R1 - A.R2 )

    #
    theta = IOTA

    # per LAL convention
    u_hat = [1,0,0]

    #
    phi = acos( dot( n_hat,u_hat ) )

    #
    N_hat = array( [  cos(phi)*sin(theta), sin(phi)*sin(theta), cos(theta)  ] )

    #
    theta_JN = acos( dot(N_hat,J_hat) )
    theta_LN = acos( dot(N_hat,L_hat) )
    theta_JL = acos( dot(J_hat,L_hat) )

    #
    if verbose:
        print '(iota,theta_JN) = (%1.4f,%1.4f)' % (IOTA*180/pi,theta_JN*180/pi)
    if False: # verbose:
        print 'phi = %1.4f' % phi
        print 'n_hat = %s' % n_hat
        print 'N_hat = %s' % N_hat
        print 'L_hat = %s' % L_hat
        print 'J_hat = %s' % J_hat
        print 'theta_JN = %1.4f degrees' % (theta_JN*180/pi)
        print 'theta_LN = %1.4f degrees' % (theta_LN*180/pi)
        print 'theta_JL = %1.4f degrees' % (theta_JL*180/pi)
        print 'theta_JZ = %1.4f degrees' % ( acos(dot(J_hat,[0,0,1])) )
        
    #
    return theta_JN

Define inclinations of interest.



In [30]:

    
NUM = 9 # Number of points in iota to use 
iota = array(  [ (k*pi/(NUM-1.0)) for k in range(0,NUM) ]  )
print iota









    



[ 0.          0.39269908  0.78539816  1.17809725  1.57079633  1.96349541
  2.35619449  2.74889357  3.14159265]

Evaluate above function for the given list of runs (scentry objects).



In [40]:

    
theta_JN = array(  [ iota2thetajn(Y,k,verbose=True) for k in iota ]  )









    



(iota,theta_JN) = (0.0000,5.3149)
(iota,theta_JN) = (22.5000,24.9079)
(iota,theta_JN) = (45.0000,47.1415)
(iota,theta_JN) = (67.5000,69.5223)
(iota,theta_JN) = (90.0000,91.9350)
(iota,theta_JN) = (112.5000,114.3452)
(iota,theta_JN) = (135.0000,136.7142)
(iota,theta_JN) = (157.5000,158.8802)
(iota,theta_JN) = (180.0000,174.6851)

Make array of iota and $\theta_{\mathrm{JN}}$ values



In [44]:

    
output_data = array( [iota,theta_JN] ).T
print output_data









    



[[ 0.          0.09276337]
 [ 0.39269908  0.43472526]
 [ 0.78539816  0.82277365]
 [ 1.17809725  1.2133932 ]
 [ 1.57079633  1.60456847]
 [ 1.96349541  1.99570061]
 [ 2.35619449  2.38611276]
 [ 2.74889357  2.77298228]
 [ 3.14159265  3.04882929]]

Save the array to a txt file.



In [49]:

    
output_file_path = '/Users/book/GARREG/REPOS/cbc_svn_nr_systematics/nr_systematics/scripts/thetaJN_'+Y.setname+'.txt'
savetxt( output_file_path, output_data, fmt='%1.8f', delimiter='\t' )



In [ ]: