Commit 745ba2cf authored by Matthew David Pitkin's avatar Matthew David Pitkin
Browse files

pulsarpputils.py: convert Q22 to h0 in posterior file reading function

parent c30bde72
......@@ -543,7 +543,7 @@ class psr_prior:
#(Hz), spin-down (Hz/s) and distance (kpc). The canonical value of moment of
# inertia of 1e38 kg m^2 is used
def spin_down_limit(freq, fdot, dist):
hsd = math.sqrt((5./2.)*(G/C**3)*I38*math.fabs(fdot)/freq)/(dist*KPC)
hsd = np.sqrt((5./2.)*(G/C**3)*I38*np.fabs(fdot)/freq)/(dist*KPC)
return hsd
......@@ -551,18 +551,25 @@ def spin_down_limit(freq, fdot, dist):
# Function to convert a pulsar stain into ellipticity assuming the canonical
# moment of inertia
def h0_to_ellipticity(h0, freq, dist):
ell = h0*C**4.*dist*KPC/(16.*math.pi**2*G*I38*freq**2)
ell = h0*C**4.*dist*KPC/(16.*np.pi**2*G*I38*freq**2)
return ell
# Function to convert a pulsar strain into a mass quadrupole moment
def h0_to_quadrupole(h0, freq, dist):
q22 = math.sqrt(15./(8.*math.pi))*h0*C**4.*dist*KPC/(16.*math.pi**2*G*freq**2)
q22 = np.sqrt(15./(8.*np.pi))*h0*C**4.*dist*KPC/(16.*np.pi**2*G*freq**2)
return q22
# Function to conver quadrupole moment to strain
def quadrupole_to_h0(q22, freq, dist):
h0 = q22*np.sqrt((8.*np.pi)/15.)*16.*np.pi**2*G*freq**2/(C**4.*dist*KPC)
return h0
# function to convert the psi' and phi0' coordinates used in nested sampling
# into the standard psi and phi0 coordinates (using vectors of those parameters
def phipsiconvert(phipchain, psipchain):
......@@ -2417,6 +2424,9 @@ def pulsar_nest_to_posterior(postfile, nestedsamples=False, removeuntrig=True):
sample file). It will be returned as a Posterior class object from `bayespputils`. The
signal evidence and noise evidence are also returned.
If there are samples in 'Q22' and not 'H0', and also a distance, or distance samples, then
the Q22 samples will be converted into equivalent H0 values.
Parameters
----------
postfile : str, required
......@@ -2474,7 +2484,13 @@ def pulsar_nest_to_posterior(postfile, nestedsamples=False, removeuntrig=True):
for pname in pnames:
# check if all samples are the same
if pos[pname].samples.tolist().count(pos[pname].samples[0]) == len(pos[pname].samples):
pos.pop(pname)
if pname == 'f0_fixed':
# try getting a fixed f0 value (for calculating h0 from Q22)
posfreqs = pos[pname].samples[0]
else:
# if distance is in the file then don't remove it (for calculating h0 from Q22 if required)
if pname != 'dist':
pos.pop(pname)
else:
# shuffle
shufpos = None
......@@ -2483,10 +2499,9 @@ def pulsar_nest_to_posterior(postfile, nestedsamples=False, removeuntrig=True):
pos.append(shufpos)
# check whether iota has been used
try:
posIota = None
if 'iota' in pos.names:
posIota = pos['iota'].samples
except:
posIota = None
if posIota is not None:
cipos = None
......@@ -2496,10 +2511,9 @@ def pulsar_nest_to_posterior(postfile, nestedsamples=False, removeuntrig=True):
pos.pop('iota')
# check whether sin(i) binary parameter has been used
try:
posI = None
if 'i' in pos.names:
posI = pos['i'].samples
except:
posI = None
if posI is not None:
sinipos = None
......@@ -2509,20 +2523,17 @@ def pulsar_nest_to_posterior(postfile, nestedsamples=False, removeuntrig=True):
pos.pop('i')
# convert C22 back into h0, and phi22 back into phi0 if required
try:
posC21 = None
if 'c21' in pos.names:
posC21 = pos['c21'].samples
except:
posC21 = None
try:
posC22 = None
if 'c22' in pos.names:
posC22 = pos['c22'].samples
except:
posC22 = None
try:
posphi22 = None
if 'phi22' in pos.names:
posphi22 = pos['phi22'].samples
except:
posphi22 = None
# convert C22 back into h0, and phi22 back into phi0 if required
if posC22 is not None and posC21 is None:
......@@ -2539,6 +2550,25 @@ def pulsar_nest_to_posterior(postfile, nestedsamples=False, removeuntrig=True):
pos.append(phi0pos)
pos.pop('phi22')
# convert Q22 to h0 is required
posQ22 = None
if 'q22' in pos.names:
posQ22 = pos['q22'].samples
posdist = None
if 'dist' in pos.names:
posdist = pos['dist'].samples
if 'f0' in pos.names:
posfreqs = pos['f0'].samples
# if not found this will have been set from the f0_fixed value above
if posQ22 is not None and posdist is not None and 'h0' not in pos.names:
h0pos = None
h0pos = bppu.PosteriorOneDPDF('h0', quadrupole_to_h0(posQ22, posfreqs, posdist))
pos.append(h0pos)
# get evidence (as in lalapps_nest2pos)
if fe == '.h5' or fe == '.hdf':
# read evidence from HDF5 file
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment