pulsarpputils.py

Go to the documentation of this file.

# -*- coding: utf-8 -*-
#
#       pulsarpputils.py
#
#       Copyright 2012
#       Matthew Pitkin <matthew.pitkin@ligo.org>
#
#
#       This program is free software; you can redistribute it and/or modify
#       it under the terms of the GNU General Public License as published by
#       the Free Software Foundation; either version 2 of the License, or
#       (at your option) any later version.
#
#       This program is distributed in the hope that it will be useful,
#       but WITHOUT ANY WARRANTY; without even the implied warranty of
#       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#       GNU General Public License for more details.
#
#       You should have received a copy of the GNU General Public License
#       along with this program; if not, write to the Free Software
#       Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
#       MA 02110-1301, USA.
 
# known pulsar analysis post-processing utilities
 
# Many functions in this a taken from, or derived from equivalents available in
# the PRESTO pulsar software package http://www.cv.nrao.edu/~sransom/presto/
 
from __future__ import print_function, division
 
import sys
import math
import cmath
import os
import numpy as np
import struct
import re
import h5py
 
try:
    from scipy.integrate import trapezoid, cumulative_trapezoid
except ImportError:
    # FIXME: Remove this once we require scipy >=1.6.0.
    from scipy.integrate import trapz as trapezoid, cumtrapz as cumulative_trapezoid
 
from scipy.interpolate import interp1d
 
try:
    from scipy.special import logsumexp
except ImportError:  # scipy < 0.19.0
    from scipy.misc import logsumexp
 
# some common constants taken from psr_constants.py in PRESTO
ARCSECTORAD = float("4.8481368110953599358991410235794797595635330237270e-6")
RADTOARCSEC = float("206264.80624709635515647335733077861319665970087963")
SECTORAD = float("7.2722052166430399038487115353692196393452995355905e-5")
RADTOSEC = float("13750.987083139757010431557155385240879777313391975")
RADTODEG = float("57.295779513082320876798154814105170332405472466564")
DEGTORAD = float("1.7453292519943295769236907684886127134428718885417e-2")
RADTOHRS = float("3.8197186342054880584532103209403446888270314977710")
HRSTORAD = float("2.6179938779914943653855361527329190701643078328126e-1")
PI = float("3.1415926535897932384626433832795028841971693993751")
TWOPI = float("6.2831853071795864769252867665590057683943387987502")
PIBYTWO = float("1.5707963267948966192313216916397514420985846996876")
SECPERDAY = float("86400.0")
SECPERJULYR = float("31557600.0")
KMPERPC = float("3.0856776e13")
KMPERKPC = float("3.0856776e16")
Tsun = float("4.925490947e-6")  # sec
Msun = float("1.9891e30")  # kg
Mjup = float("1.8987e27")  # kg
Rsun = float("6.9551e8")  # m
Rearth = float("6.378e6")  # m
SOL = float("299792458.0")  # m/s
MSUN = float("1.989e+30")  # kg
G = float("6.673e-11")  # m^3/s^2/kg
C = SOL
KPC = float("3.0856776e19")  # kiloparsec in metres
I38 = float("1e38")  # moment of inertia kg m^2
 
 
# some angle conversion functions taken from psr_utils.py in PRESTO
def rad_to_dms(rad):
    """
    rad_to_dms(rad):
       Convert radians to degrees, minutes, and seconds of arc.
    """
    if rad < 0.0:
        sign = -1
    else:
        sign = 1
    arc = RADTODEG * np.fmod(np.fabs(rad), math.pi)
    d = int(arc)
    arc = (arc - d) * 60.0
    m = int(arc)
    s = (arc - m) * 60.0
    if sign == -1 and d == 0:
        return (sign * d, sign * m, sign * s)
    else:
        return (sign * d, m, s)
 
 
def dms_to_rad(deg, mins, sec):
    """
    dms_to_rad(deg, min, sec):
       Convert degrees, minutes, and seconds of arc to radians.
    """
    if deg < 0.0:
        sign = -1
    elif deg == 0.0 and (mins < 0.0 or sec < 0.0):
        sign = -1
    else:
        sign = 1
    return (
        sign
        * ARCSECTORAD
        * (60.0 * (60.0 * np.fabs(deg) + np.fabs(mins)) + np.fabs(sec))
    )
 
 
def dms_to_deg(deg, mins, sec):
    """
    dms_to_deg(deg, min, sec):
       Convert degrees, minutes, and seconds of arc to degrees.
    """
    return RADTODEG * dms_to_rad(deg, mins, sec)
 
 
def rad_to_hms(rad):
    """
    rad_to_hms(rad):
       Convert radians to hours, minutes, and seconds of arc.
    """
    rad = np.fmod(rad, 2.0 * math.pi)
    if rad < 0.0:
        rad = rad + 2.0 * math.pi
    arc = RADTOHRS * rad
    h = int(arc)
    arc = (arc - h) * 60.0
    m = int(arc)
    s = (arc - m) * 60.0
    return (h, m, s)
 
 
def hms_to_rad(hour, mins, sec):
    """
    hms_to_rad(hour, min, sec):
       Convert hours, minutes, and seconds of arc to radians
    """
    if hour < 0.0:
        sign = -1
    else:
        sign = 1
    return (
        sign * SECTORAD * (60.0 * (60.0 * np.fabs(hour) + np.fabs(mins)) + np.fabs(sec))
    )
 
 
def coord_to_string(h_or_d, m, s):
    """
    coord_to_string(h_or_d, m, s):
       Return a formatted string of RA or DEC values as
       'hh:mm:ss.ssss' if RA, or 'dd:mm:ss.ssss' if DEC.
    """
    retstr = ""
    if h_or_d < 0:
        retstr = "-"
    elif abs(h_or_d) == 0:
        if (m < 0.0) or (s < 0.0):
            retstr = "-"
 
    h_or_d, m, s = abs(h_or_d), abs(m), abs(s)
    if s >= 9.9995:
        return retstr + "%.2d:%.2d:%.4f" % (h_or_d, m, s)
    else:
        return retstr + "%.2d:%.2d:0%.4f" % (h_or_d, m, s)
 
 
def rad_to_string(rad, ra_or_dec):
    """
    rad_to_string(rad, ra_or_dec):
       Convert an angle in radians to hours/degrees, minutes seconds and output
       it as a string in the format 'hh:mm:ss.ssss' if RA, or 'dd:mm:ss.ssss' if DEC.
       Whether to use hours or degrees is set by whether ra_or_dec is 'RA' or 'DEC'
    """
    if ra_or_dec.upper() == "RA":
        v, m, s = rad_to_hms(rad)
    elif ra_or_dec.upper() == "DEC":
        v, m, s = rad_to_dms(rad)
    else:
        raise ValueError(
            "Unrecognised option: Expected 'ra_or_dec' to be 'RA' or 'DEC'"
        )
 
    return coord_to_string(v, m, s)
 
 
def ra_to_rad(ra_string):
    """
    ra_to_rad(ar_string):
       Given a string containing RA information as
       'hh:mm:ss.ssss', return the equivalent decimal
       radians. Also deal with cases where input
       string is just hh:mm, or hh.
    """
    hms = ra_string.split(":")
    if len(hms) == 3:
        return hms_to_rad(int(hms[0]), int(hms[1]), float(hms[2]))
    elif len(hms) == 2:
        return hms_to_rad(int(hms[0]), int(hms[1]), 0.0)
    elif len(hms) == 1:
        return hms_to_rad(float(hms[0]), 0.0, 0.0)
    else:
        print("Problem parsing RA string %s" % ra_string, file=sys.stderr)
        sys.exit(1)
 
 
def dec_to_rad(dec_string):
    """
    dec_to_rad(dec_string):
       Given a string containing DEC information as
       'dd:mm:ss.ssss', return the equivalent decimal
       radians. Also deal with cases where input string
       is just dd:mm or dd
    """
    dms = dec_string.split(":")
    if "-" in dms[0] and float(dms[0]) == 0.0:
        m = "-"
    else:
        m = ""
 
    if len(dms) == 3:
        return dms_to_rad(int(dms[0]), int(m + dms[1]), float(m + dms[2]))
    elif len(dms) == 2:
        return dms_to_rad(int(dms[0]), int(m + dms[1]), 0.0)
    elif len(dms) == 1:
        return dms_to_rad(float(dms[0]), 0.0, 0.0)
    else:
        print("Problem parsing DEC string %s" % dec_string, file=sys.stderr)
        sys.exit(1)
 
 
def p_to_f(p, pd, pdd=None):
    """
    p_to_f(p, pd, pdd=None):
      Convert period, period derivative and period second
      derivative to the equivalent frequency counterparts.
      Will also convert from f to p.
    """
    f = 1.0 / p
    fd = -pd / (p * p)
    if pdd == None:
        return [f, fd]
    else:
        if pdd == 0.0:
            fdd = 0.0
        else:
            fdd = 2.0 * pd * pd / (p**3.0) - pdd / (p * p)
 
        return [f, fd, fdd]
 
 
def pferrs(porf, porferr, pdorfd=None, pdorfderr=None):
    """
    pferrs(porf, porferr, pdorfd=None, pdorfderr=None):
       Calculate the period or frequency errors and
       the pdot or fdot errors from the opposite one.
    """
    if pdorfd == None:
        return [1.0 / porf, porferr / porf**2.0]
    else:
        forperr = porferr / porf**2.0
        fdorpderr = np.sqrt(
            (4.0 * pdorfd**2.0 * porferr**2.0) / porf**6.0
            + pdorfderr**2.0 / porf**4.0
        )
        [forp, fdorpd] = p_to_f(porf, pdorfd)
 
        return [forp, forperr, fdorpd, fdorpderr]
 
 
# class to read in a pulsar par file - this is heavily based on the function
# in parfile.py in PRESTO
float_keys = [
    "F",
    "F0",
    "F1",
    "F2",
    "F3",
    "F4",
    "F5",
    "F6",
    "F7",
    "F8",
    "F9",
    "F10",
    "F11",
    "F12",
    "PEPOCH",
    "POSEPOCH",
    "DM",
    "START",
    "FINISH",
    "NTOA",
    "TRES",
    "TZRMJD",
    "TZRFRQ",
    "TZRSITE",
    "NITS",
    "A1",
    "XDOT",
    "E",
    "ECC",
    "EDOT",
    "T0",
    "PB",
    "PBDOT",
    "OM",
    "OMDOT",
    "EPS1",
    "EPS2",
    "EPS1DOT",
    "EPS2DOT",
    "TASC",
    "LAMBDAPIN",
    "BETA",
    "RA_RAD",
    "DEC_RAD",
    "GAMMA",
    "SINI",
    "M2",
    "MTOT",
    "FB0",
    "FB1",
    "FB2",
    "ELAT",
    "ELONG",
    "PMRA",
    "PMDEC",
    "DIST",
    "PB_2",
    "PB_3",
    "T0_2",
    "T0_3",
    "A1_2",
    "A1_3",
    "OM_2",
    "OM_3",
    "ECC_2",
    "ECC_3",
    "PX",
    "KIN",
    "KOM",
    "A0",
    "B0",
    "D_AOP",
    # GW PARAMETERS
    "H0",
    "COSIOTA",
    "PSI",
    "PHI0",
    "THETA",
    "I21",
    "I31",
    "C22",
    "HPLUS",
    "HCROSS",
    "C21",
    "PHI22",
    "PHI21",
    "SNR",
    "COSTHETA",
    "IOTA",
    "Q22",
]
str_keys = [
    "FILE",
    "PSR",
    "PSRJ",
    "NAME",
    "RAJ",
    "DECJ",
    "RA",
    "DEC",
    "EPHEM",
    "CLK",
    "BINARY",
    "UNITS",
]
 
 
class psr_par:
    def __init__(self, parfilenm):
        """
        This class parses a TEMPO(2)-style pulsar parameter file. If possible all parameters
        are converted into SI units and angles are in radians. Epochs will be converted from
        MJD values into GPS times.
        """
        self.FILE = parfilenm
        pf = open(parfilenm)
        for line in pf.readlines():
            # ignore empty lines (i.e. containing only whitespace)
            if not line.strip():
                continue
 
            # Convert any 'D-' or 'D+' to 'E-' or 'E+'
            line = line.replace("D-", "E-")
            line = line.replace("D+", "E+")
            # also check for lower case
            line = line.replace("d-", "e-")
            line = line.replace("d+", "e+")
            splitline = line.split()
 
            # get all upper case version in case lower case in par file
            key = splitline[0].upper()
 
            if key in str_keys:
                setattr(self, key, splitline[1])
            elif key in float_keys:
                try:
                    setattr(self, key, float(splitline[1]))
                except:
                    continue
 
            if (
                len(splitline) == 3
            ):  # Some parfiles don't have flags, but do have errors
                if splitline[2] not in ["0", "1"]:
                    setattr(self, key + "_ERR", float(splitline[2]))
                    setattr(self, key + "_FIT", 0)  # parameter was not fit
 
            if len(splitline) == 4:
                if splitline[2] == "1":  # parameter was fit
                    setattr(self, key + "_FIT", 1)
                else:
                    setattr(self, key + "_FIT", 0)
                setattr(self, key + "_ERR", float(splitline[3]))
 
        # sky position
        if hasattr(self, "RAJ"):
            setattr(self, "RA_RAD", ra_to_rad(self.RAJ))
 
            # set RA error in rads (rather than secs)
            if hasattr(self, "RAJ_ERR"):
                setattr(self, "RA_RAD_ERR", hms_to_rad(0, 0, self.RAJ_ERR))
 
                if hasattr(self, "RAJ_FIT"):
                    setattr(self, "RA_RAD_FIT", self["RAJ_FIT"])
        if hasattr(self, "DECJ"):
            setattr(self, "DEC_RAD", dec_to_rad(self.DECJ))
 
            # set DEC error in rads (rather than arcsecs)
            if hasattr(self, "DECJ_ERR"):
                setattr(self, "DEC_RAD_ERR", dms_to_rad(0, 0, self.DECJ_ERR))
 
                if hasattr(self, "DECJ_FIT"):
                    setattr(self, "DEC_RAD_FIT", self["DECJ_FIT"])
 
        # convert proper motions to rads/sec from mas/year
        for pv in ["RA", "DEC"]:
            if hasattr(self, "PM" + pv):
                pmv = self["PM" + pv]
                setattr(self, "PM" + pv + "_ORIGINAL", pmv)  # save original value
                setattr(
                    self, "PM" + pv, pmv * np.pi / (180.0 * 3600.0e3 * 365.25 * 86400.0)
                )
 
                if hasattr(self, "PM" + pv + "_ERR"):
                    pmve = self["PM" + pv + "_ERR"]
                    setattr(
                        self, "PM" + pv + "_ERR_ORIGINAL", pmv
                    )  # save original value
                    setattr(
                        self,
                        "PM" + pv + "_ERR",
                        pmve * np.pi / (180.0 * 3600.0e3 * 365.25 * 86400.0),
                    )
 
        # periods and frequencies
        if hasattr(self, "P"):
            setattr(self, "P0", self.P)
        if hasattr(self, "P0"):
            setattr(self, "F0", 1.0 / self.P0)
        if hasattr(self, "F0"):
            setattr(self, "P0", 1.0 / self.F0)
        if hasattr(self, "FB0"):
            setattr(self, "PB", (1.0 / self.FB0) / 86400.0)
        if hasattr(self, "P0_ERR"):
            if hasattr(self, "P1_ERR"):
                f, ferr, fd, fderr = pferrs(self.P0, self.P0_ERR, self.P1, self.P1_ERR)
                setattr(self, "F0_ERR", ferr)
                setattr(self, "F1", fd)
                setattr(self, "F1_ERR", fderr)
            else:
                if hasattr(self, "P1"):
                    (
                        f,
                        fd,
                    ) = p_to_f(self.P0, self.P1)
                    setattr(self, "F0_ERR", self.P0_ERR / (self.P0 * self.P0))
                    setattr(self, "F1", fd)
        if hasattr(self, "F0_ERR"):
            if hasattr(self, "F1_ERR"):
                p, perr, pd, pderr = pferrs(self.F0, self.F0_ERR, self.F1, self.F1_ERR)
                setattr(self, "P0_ERR", perr)
                setattr(self, "P1", pd)
                setattr(self, "P1_ERR", pderr)
            else:
                if hasattr(self, "F1"):
                    (
                        p,
                        pd,
                    ) = p_to_f(self.F0, self.F1)
                    setattr(self, "P0_ERR", self.F0_ERR / (self.F0 * self.F0))
                    setattr(self, "P1", pd)
 
        # convert epochs (including binary epochs) to GPS if possible
        try:
            from lalpulsar import TTMJDtoGPS
 
            for epoch in [
                "PEPOCH",
                "POSEPOCH",
                "DMEPOCH",
                "T0",
                "TASC",
                "T0_2",
                "T0_3",
            ]:
                if hasattr(self, epoch):
                    setattr(
                        self, epoch + "_ORIGINAL", self[epoch]
                    )  # save original value
                    setattr(self, epoch, TTMJDtoGPS(self[epoch]))
 
                    if hasattr(
                        self, epoch + "_ERR"
                    ):  # convert errors from days to seconds
                        setattr(
                            self, epoch + "_ERR_ORIGINAL", self[epoch + "_ERR"]
                        )  # save original value
                        setattr(self, epoch + "_ERR", self[epoch + "_ERR"] * SECPERDAY)
        except:
            print(
                "Could not convert epochs to GPS times. They are all still MJD values.",
                file=sys.stderr,
            )
 
        # distance and parallax (distance: kpc -> metres, parallax: mas -> rads)
        convfacs = {"DIST": KPC, "PX": 1e-3 * ARCSECTORAD}
        for item in convfacs:
            if hasattr(self, item):  # convert kpc to metres
                setattr(self, item + "_ORIGINAL", self[item])  # save original value
                setattr(self, item, self[item] * convfacs[item])
 
                if hasattr(self, item + "_ERR"):
                    setattr(
                        self, item + "_ERR_ORIGINAL", self[item + "_ERR"]
                    )  # save original value
                    setattr(self, item + "_ERR", self[item + "_ERR"] * convfacs[item])
 
        # binary parameters
        # omega (angle of periastron) parameters (or others requiring conversion from degs to rads)
        for om in ["OM", "OM_2", "OM_3", "KIN", "KOM"]:
            if hasattr(self, om):  # convert OM from degs to rads
                setattr(self, om, self[om] / RADTODEG)
 
                if hasattr(self, om + "_ERR"):
                    setattr(self, om + "_ERR", self[om + "_ERR"] / RADTODEG)
 
        # period
        for pb in ["PB", "PB_2", "PB_3"]:
            if hasattr(self, pb):  # convert PB from days to seconds
                setattr(self, pb + "_ORIGINAL", self[pb])  # save original value
                setattr(self, pb, self[pb] * SECPERDAY)
 
                if hasattr(self, pb + "_ERR"):
                    setattr(
                        self, pb + "_ERR_ORIGINAL", self[pb + "_ERR"]
                    )  # save original value
                    setattr(self, pb + "_ERR", self[pb + "_ERR"] * SECPERDAY)
 
        # OMDOT
        if hasattr(self, "OMDOT"):  # convert from deg/year to rad/sec
            setattr(self, "OMDOT_ORIGINAL", self["OMDOT"])  # save original value
            setattr(self, "OMDOT", self["OMDOT"] / (RADTODEG * 365.25 * SECPERDAY))
 
            if hasattr(self, "OMDOT_ERR"):
                setattr(
                    self, "OMDOT_ERR_ORIGINAL", self["OMDOT_ERR"]
                )  # save original value
                setattr(
                    self,
                    "OMDOT_ERR",
                    self["OMDOT_ERR"] / (RADTODEG * 365.25 * SECPERDAY),
                )
 
        if hasattr(self, "EPS1") and hasattr(self, "EPS2"):
            ecc = math.sqrt(self.EPS1 * self.EPS1 + self.EPS2 * self.EPS2)
            omega = math.atan2(self.EPS1, self.EPS2)
            setattr(self, "E", ecc)
            setattr(self, "OM", omega)  # omega in rads
            setattr(self, "T0", self.TASC + self.PB * omega / TWOPI)
        if (
            hasattr(self, "PB")
            and hasattr(self, "A1")
            and not (hasattr(self, "E") or hasattr(self, "ECC"))
        ):
            setattr(self, "E", 0.0)
        if (
            hasattr(self, "T0")
            and not hasattr(self, "TASC")
            and hasattr(self, "OM")
            and hasattr(self, "PB")
        ):
            setattr(self, "TASC", self.T0 - self.PB * self.OM / TWOPI)
 
        # binary unit conversion for small numbers (TEMPO2 checks if these are > 1e-7 and if so then the units are in 1e-12) - errors are not converted
        for binu in ["XDOT", "PBDOT", "EDOT", "EPS1DOT", "EPS2DOT", "XPBDOT"]:
            if hasattr(self, binu):
                setattr(self, binu + "_ORIGINAL", self[binu])  # save original value
                if np.abs(self[binu]) > 1e-7:
                    setattr(self, binu, self[binu] * 1.0e-12)
 
                    # check value is not extremely large due to ATNF catalogue error
                    if self[binu] > 10000.0:  # set to zero
                        setattr(self, binu, 0.0)
 
        # masses
        for mass in ["M2", "MTOT"]:
            if hasattr(self, mass):  # convert solar masses to kg
                setattr(self, mass + "_ORIGINAL", self[mass])  # save original value
                setattr(self, mass, self[mass] * MSUN)
 
                if hasattr(self, mass + "_ERR"):
                    setattr(
                        self, mass + "_ERR_ORIGINAL", self[mass + "_ERR"]
                    )  # save original value
                    setattr(self, mass + "_ERR", self[mass + "_ERR"] * MSUN)
 
        # D_AOP
        if hasattr(self, "D_AOP"):  # convert from inverse arcsec to inverse radians
            setattr(self, "D_AOP_ORIGINAL", self["D_AOP"])  # save original value
            setattr(self, "D_AOP", self["D_AOP"] * RADTODEG * 3600.0)
 
        pf.close()
 
    def __getitem__(self, key):
        try:
            par = getattr(self, key)
        except:
            par = None
 
        return par
 
    def __str__(self):
        out = ""
        for k, v in self.__dict__.items():
            if k[:2] != "__":
                if isinstance(self.__dict__[k], string_type):
                    out += "%10s = '%s'\n" % (k, v)
                else:
                    out += "%10s = %-20.15g\n" % (k, v)
 
        return out
 
 
# class to read in a nested sampling prior file
class psr_prior:
    def __init__(self, priorfilenm):
        self.FILE = priorfilenm
        pf = open(priorfilenm)
        for line in pf.readlines():
            splitline = line.split()
 
            # get all upper case version in case lower case in par file
            key = splitline[0].upper()
 
            if key in str_keys:
                # everything in a prior files should be numeric
                setattr(self, key, [float(splitline[2]), float(splitline[3])])
            elif key in float_keys:
                setattr(self, key, [float(splitline[2]), float(splitline[3])])
 
        # get sky positions in rads as strings 'dd/hh:mm:ss.s'
        if hasattr(self, "RA"):
            hl, ml, sl = rad_to_hms(self.RA[0])
            rastrl = coord_to_string(hl, ml, sl)
            hu, mu, su = rad_to_hms(self.RA[1])
            rastru = coord_to_string(hu, mu, su)
            setattr(self, "RA_STR", [rastrl, rastru])
 
        if hasattr(self, "DEC"):
            dl, ml, sl = rad_to_dms(self.DEC[0])
            decstrl = coord_to_string(dl, ml, sl)
            du, mu, su = rad_to_dms(self.DEC[1])
            decstru = coord_to_string(du, mu, su)
            setattr(self, "DEC_STR", [decstrl, decstru])
 
        pf.close()
 
    def __getitem__(self, key):
        try:
            atr = getattr(self, key)
        except:
            atr = None
 
        return atr
 
    def __str__(self):
        out = ""
        for k, v in self.__dict__.items():
            if k[:2] != "__":
                out += "%10s = %-20.15g, %-20.15g\n" % (k, float(v[0]), float(v[1]))
 
        return out
 
 
# Function to return a pulsar's strain spin-down limit given its spin frequency
# (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 = np.sqrt((5.0 / 2.0) * (G / C**3) * I38 * np.fabs(fdot) / freq) / (
        dist * KPC
    )
 
    return hsd
 
 
# 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.0 * dist * KPC / (16.0 * 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 = (
        np.sqrt(15.0 / (8.0 * np.pi))
        * h0
        * C**4.0
        * dist
        * KPC
        / (16.0 * 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.0 * np.pi) / 15.0)
        * 16.0
        * np.pi**2
        * G
        * freq**2
        / (C**4.0 * 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):
    chainlen = len(phipchain)
 
    phichain = []
    psichain = []
 
    theta = math.atan2(1, 2)
    ct = math.cos(theta)
    st = math.sin(theta)
 
    for i in range(0, chainlen):
        phi0 = (1 / (2 * st)) * phipchain[i] - (1 / (2 * st)) * psipchain[i]
        psi = (1 / (2 * ct)) * phipchain[i] + (1 / (2 * ct)) * psipchain[i]
 
        # put psi between +/-pi/4
        if math.fabs(psi) > math.pi / 4.0:
            # shift phi0 by pi
            phi0 = phi0 + math.pi
 
            # wrap around psi
            if psi > math.pi / 4.0:
                psi = -(math.pi / 4.0) + math.fmod(psi + (math.pi / 4.0), math.pi / 2.0)
            else:
                psi = (math.pi / 4.0) - math.fmod((math.pi / 4.0) - psi, math.pi / 2.0)
 
        # get phi0 into 0 -> 2pi range
        if phi0 > 2.0 * math.pi:
            phi0 = math.fmod(phi0, 2.0 * math.pi)
        else:
            phi0 = 2.0 * math.pi - math.fmod(2.0 * math.pi - phi0, 2.0 * math.pi)
 
        phichain.append(phi0)
        psichain.append(psi)
 
    return phichain, psichain
 
 
# function to create histogram plot of the 1D posterior (potentially for
# multiple IFOs) for a parameter (param). If an upper limit is given then
# that will be output
def plot_posterior_hist(
    poslist,
    param,
    ifos,
    parambounds=[float("-inf"), float("inf")],
    nbins=50,
    upperlimit=0,
    overplot=False,
    parfile=None,
    mplparams=False,
):
    import matplotlib
    from matplotlib import pyplot as plt
    from lalpulsar.pulsarhtmlutils import paramlatexdict
 
    # create list of figures
    myfigs = []
 
    # create a list of upper limits
    ulvals = []
 
    # set some matplotlib defaults for hist
    if not mplparams:
        mplparams = {
            "backend": "Agg",
            "text.usetex": True,  # use LaTeX for all text
            "axes.linewidth": 0.5,  # set axes linewidths to 0.5
            "axes.grid": True,  # add a grid
            "grid.linewidth": 0.5,
            "font.family": "serif",
            "font.size": 12,
        }
 
    matplotlib.rcParams.update(mplparams)
 
    # ifos line colour specs
    coldict = {"H1": "r", "H2": "c", "L1": "g", "V1": "b", "G1": "m", "Joint": "k"}
 
    # param name for axis label
    try:
        paraxis = paramlatexdict[param.upper()]
    except:
        paraxis = param
 
    ymax = []
 
    # if a par file object is given expect that we have an injection file
    # containing the injected value to overplot on the histgram
    parval = None
    if parfile:
        parval = parfile[param.upper()]
 
    if ifos == None:
        # default to just output colour for H1
        ifos = ["H1"]
 
    # loop over ifos
    for idx, ifo in enumerate(ifos):
        # check whether to plot all figures on top of each other
        if overplot and idx == 0:
            myfig = plt.figure(figsize=(4, 4), dpi=200)
            plt.hold(True)
        elif not overplot:
            myfig = plt.figure(figsize=(4, 4), dpi=200)
 
        pos = poslist[idx]
 
        pos_samps = pos[param].samples
 
        # get a normalised histogram for each
        n, bins = hist_norm_bounds(
            pos_samps, int(nbins), parambounds[0], parambounds[1]
        )
 
        # plot histogram
        plt.step(bins, n, color=coldict[ifo])
 
        if "h0" not in param:
            plt.xlim(parambounds[0], parambounds[1])
 
        plt.xlabel(r"" + paraxis, fontsize=14, fontweight=100)
        plt.ylabel(r"Probability Density", fontsize=14, fontweight=100)
        myfig.subplots_adjust(left=0.18, bottom=0.15)  # adjust size
 
        if not overplot:
            plt.ylim(0, n.max() + 0.1 * n.max())
            # plt.legend(ifo)
            # set background colour of axes
            ax = plt.gca()
            ax.set_axis_bgcolor("#F2F1F0")
            myfigs.append(myfig)
        else:
            ymax.append(n.max() + 0.1 * n.max())
 
        # if upper limit is needed then integrate posterior using trapezium rule
        if upperlimit != 0:
            ct = cumulative_trapezoid(n, bins)
 
            # prepend a zero to ct
            ct = np.insert(ct, 0, 0)
 
            # use linear interpolation to find the value at 'upper limit'
            ctu, ui = np.unique(ct, return_index=True)
            intf = interp1d(ctu, bins[ui], kind="linear")
            ulvals.append(intf(float(upperlimit)))
 
    # plot parameter values
    if parval:
        if not overplot:
            plt.hold(True)
            plt.axvline(parval, color="k", ls="--", linewidth=1.5)
        else:
            plt.axvline(parval, color="k", ls="--", linewidth=1.5)
 
    if overplot:
        plt.ylim(0, max(ymax))
        # plt.legend(ifos)
        ax = plt.gca()
        ax.set_axis_bgcolor("#F2F1F0")
        plt.hold(False)
        myfigs.append(myfig)
 
    return myfigs, ulvals
 
 
# function to return an upper limit from a posteriors: pos is an array of posteriors samples for a
# particular parameter
def upper_limit(
    pos, upperlimit=0.95, parambounds=[float("-inf"), float("inf")], nbins=50
):
    ulval = 0
 
    # get a normalised histogram of posterior samples
    n, bins = hist_norm_bounds(pos, int(nbins), parambounds[0], parambounds[1])
 
    # if upper limit is needed then integrate posterior using trapezium rule
    if upperlimit != 0:
        ct = cumulative_trapezoid(n, bins)
 
        # prepend a zero to ct
        ct = np.insert(ct, 0, 0)
 
        # use linear interpolation to find the value at 'upper limit'
        ctu, ui = np.unique(ct, return_index=True)
        intf = interp1d(ctu, bins[ui], kind="linear")
        ulval = intf(float(upperlimit))
 
    return ulval
 
 
def upper_limit_greedy(pos, upperlimit=0.95, nbins=100):
    n, binedges = np.histogram(pos, bins=nbins)
    dbins = binedges[1] - binedges[0]  # width of a histogram bin
 
    frac = 0.0
    j = 0
    for nv in n:
        prevfrac = frac
        frac += float(nv) / len(pos)
        j += 1
        if frac > upperlimit:
            break
 
    # linearly interpolate to get upper limit
    ul = binedges[j - 1] + (upperlimit - prevfrac) * (dbins / (frac - prevfrac))
 
    return ul
 
 
# function to plot a posterior chain (be it MCMC chains or nested samples)
# the input should be a list of posteriors for each IFO, and the parameter
# required, the list of IFO. grr is a list of dictionaries giving
# the Gelman-Rubins statistics for the given parameter for each IFO.
# If withhist is set then it will also output a histgram, with withhist number
# of bins
def plot_posterior_chain(poslist, param, ifos, grr=None, withhist=0, mplparams=False):
    import matplotlib
    from matplotlib import pyplot as plt
    from lalpulsar.pulsarhtmlutils import paramlatexdict
 
    try:
        from matplotlib import gridspec
    except:
        return None
 
    if not mplparams:
        mplparams = {
            "backend": "Agg",
            "text.usetex": True,  # use LaTeX for all text
            "axes.linewidth": 0.5,  # set axes linewidths to 0.5
            "axes.grid": True,  # add a grid
            "grid.linewidth": 0.5,
            "font.family": "sans-serif",
            "font.sans-serif": "Avant Garde, Helvetica, Computer Modern Sans serif",
            "font.size": 15,
        }
 
    matplotlib.rcParams.update(mplparams)
 
    coldict = {"H1": "r", "H2": "c", "L1": "g", "V1": "b", "G1": "m", "Joint": "k"}
 
    # param name for axis label
    try:
        if param == "iota":
            p = "cosiota"
        else:
            p = param
 
        paryaxis = paramlatexdict[p.upper()]
    except:
        paryaxis = param
 
    if grr:
        legendvals = []
 
    maxiter = 0
    maxn = 0
    minsamp = float("inf")
    maxsamp = -float("inf")
 
    for idx, ifo in enumerate(ifos):
        if idx == 0:
            myfig = plt.figure(figsize=(12, 4), dpi=200)
            myfig.subplots_adjust(bottom=0.15)
 
            if withhist:
                gs = gridspec.GridSpec(1, 4, wspace=0)
                ax1 = plt.subplot(gs[:-1])
                ax2 = plt.subplot(gs[-1])
 
        pos = list(poslist)[idx]
 
        # check for cosiota
        if "iota" == param:
            pos_samps = np.cos(pos["iota"].samples)
        else:
            pos_samps = pos[param].samples
 
        if np.min(pos_samps) < minsamp:
            minsamp = np.min(pos_samps)
        if np.max(pos_samps) > maxsamp:
            maxsamp = np.max(pos_samps)
 
        if withhist:
            ax1.plot(pos_samps, ".", color=coldict[ifo], markersize=1)
 
            n, binedges = np.histogram(pos_samps, withhist, density=True)
            n = np.append(n, 0)
            ax2.step(n, binedges, color=coldict[ifo])
 
            if np.max(n) > maxn:
                maxn = np.max(n)
        else:
            plt.plot(pos_samps, ".", color=coldict[ifo], markersize=1)
 
        if grr:
            try:
                legendvals.append(r"$R = %.2f$" % grr[idx][param])
            except:
                legendvals = []
        else:
            legendvals = None
 
        if len(pos_samps) > maxiter:
            maxiter = len(pos_samps)
 
    if not withhist:
        ax1 = plt.gca()
 
    bounds = [minsamp, maxsamp]
 
    ax1.set_ylabel(r"" + paryaxis, fontsize=16, fontweight=100)
    ax1.set_xlabel(r"Iterations", fontsize=16, fontweight=100)
 
    ax1.set_xlim(0, maxiter)
    ax1.set_ylim(bounds[0], bounds[1])
 
    if withhist:
        ax2.set_ylim(bounds[0], bounds[1])
        ax2.set_xlim(0, maxn + 0.1 * maxn)
        ax2.set_yticklabels([])
        ax2.set_axis_bgcolor("#F2F1F0")
 
    # add gelman-rubins stat data
    if legendvals:
        ax1.legend(legendvals, title="Gelman-Rubins test")
 
    return myfig
 
 
# function to read in and plot a 2D histogram from a binary file, where the files
# structure is of the form:
#   a header of six doubles with:
#     - minimum of N-dim
#     - the step size in N-dim
#     - N - number of N-dim values
#     - minimum of M-dim
#     - the step size in M-dim
#     - M - number of M-dim values
#   followed by an NxM double array of values.
# The information returned are two lists of the x and y-axis bin centres, along
# with a numpy array containing the histogram values.
def read_hist_from_file(histfile):
    # read in 2D binary file
    try:
        fp = open(histfile, "rb")
    except:
        print("Could not open prior file %s" % histfile, file=sys.stderr)
        return None, None, None
 
    try:
        pd = fp.read()  # read in all the data
    except:
        print("Could not read in data from prior file %s" % histfile, file=sys.stderr)
        return None, None, None
 
    fp.close()
 
    # read in the header (6 doubles)
    # try:
    header = struct.unpack("d" * 6, pd[: 6 * 8])
    # except:
    #  print >> sys.stderr, "Could not read in header data from prior file %s" % histfile
    #  return None, None, None
 
    # read in histogram grid (NxM doubles)
    # try:
    grid = struct.unpack("d" * int(header[2]) * int(header[5]), pd[6 * 8 :])
    # except:
    #  print >> sys.stderr, "Could not read in header data from prior file %s" % histfile
    #  return None, None, None
 
    header = list(header)
 
    # convert grid into numpy array
    g = list(grid)  # convert to list
    histarr = np.array([g[: int(header[5])]], ndmin=2)
    for i in range(int(header[2]) - 1):
        histarr = np.append(
            histarr, [g[(i + 1) * int(header[5]) : (i + 2) * int(header[5])]], axis=0
        )
 
    xbins = np.linspace(
        header[0], header[0] + header[1] * (header[2] - 1), int(header[2])
    )
    ybins = np.linspace(
        header[3], header[3] + header[4] * (header[5] - 1), int(header[5])
    )
 
    return xbins, ybins, histarr
 
 
# Function to plot a 2D histogram of from a binary file containing it.
# Also supply the label names for the n-dimension (x-axis) and m-dimension
# (y-axis). If margpars is true marginalised plots of both dimensions will
# also be returned.
def plot_2Dhist_from_file(
    histfile, ndimlabel, mdimlabel, margpars=True, mplparams=False
):
    import matplotlib
    from matplotlib import pyplot as plt
    from lalpulsar.pulsarhtmlutils import paramlatexdict
 
    # read in 2D h0 vs cos(iota) binary prior file
    xbins, ybins, histarr = read_hist_from_file(histfile)
 
    if not xbins.any():
        print("Could not read binary histogram file", file=sys.stderr)
        return None
 
    figs = []
 
    # set some matplotlib defaults for amplitude spectral density
    if not mplparams:
        mplparams = {
            "backend": "Agg",
            "text.usetex": True,  # use LaTeX for all text
            "axes.linewidth": 0.5,  # set axes linewidths to 0.5
            "axes.grid": True,  # add a grid
            "grid.linewidth": 0.5,
            "font.family": "serif",
            "font.size": 12,
        }
 
    matplotlib.rcParams.update(mplparams)
 
    fig = plt.figure(figsize=(4, 4), dpi=200)
 
    # param name for axis label
    try:
        parxaxis = paramlatexdict[ndimlabel.upper()]
    except:
        parxaxis = ndimlabel
 
    try:
        paryaxis = paramlatexdict[mdimlabel.upper()]
    except:
        paryaxis = mdimlabel
 
    plt.xlabel(r"" + parxaxis, fontsize=14, fontweight=100)
    plt.ylabel(r"" + paryaxis, fontsize=14, fontweight=100, rotation=270)
 
    dx = xbins[1] - xbins[0]
    dy = ybins[1] - ybins[0]
 
    extent = [
        xbins[0] - dx / 2.0,
        xbins[-1] + dx / 2.0,
        ybins[-1] + dy / 2.0,
        ybins[0] - dy / 2.0,
    ]
 
    plt.imshow(
        np.transpose(histarr),
        aspect="auto",
        extent=extent,
        interpolation="bicubic",
        cmap="gray_r",
    )
 
    fig.subplots_adjust(left=0.18, bottom=0.15)  # adjust size
    fax = (fig.get_axes())[0].axis()
 
    figs.append(fig)
 
    if margpars:
        # marginalise over y-axes and produce plots
        xmarg = []
        for i in range(len(xbins)):
            xmarg.append(trapezoid(histarr[:][i], x=ybins))
 
        # normalise
        xarea = trapezoid(xmarg, x=xbins)
        xmarg = map(lambda x: x / xarea, xmarg)
 
        ymarg = []
        for i in range(len(ybins)):
            ymarg.append(trapezoid(np.transpose(histarr)[:][i], x=xbins))
 
        # normalise
        yarea = trapezoid(ymarg, x=ybins)
        ymarg = map(lambda x: x / yarea, ymarg)
 
        # plot x histogram
        figx = plt.figure(figsize=(4, 4), dpi=200)
        plt.step(xbins, xmarg, color="k")
        plt.ylim(0, max(xmarg) + 0.1 * max(xmarg))
        plt.xlim(fax[0], fax[1])
        plt.xlabel(r"" + parxaxis, fontsize=14, fontweight=100)
        plt.ylabel(r"Probability Density", fontsize=14, fontweight=100)
        figx.subplots_adjust(left=0.18, bottom=0.15)  # adjust size
        ax = plt.gca()
        ax.set_axis_bgcolor("#F2F1F0")
 
        # plot y histogram
        figy = plt.figure(figsize=(4, 4), dpi=200)
        plt.step(ybins, ymarg, color="k")
        plt.ylim(0, max(ymarg) + 0.1 * max(ymarg))
        plt.xlim(fax[3], fax[2])
        plt.xlabel(r"" + paryaxis, fontsize=14, fontweight=100)
        plt.ylabel(r"Probability Density", fontsize=14, fontweight=100)
        figy.subplots_adjust(left=0.18, bottom=0.15)  # adjust size
        ax = plt.gca()
        ax.set_axis_bgcolor("#F2F1F0")
 
        figs.append(figx)
        figs.append(figy)
 
    return figs
 
 
# using a binary file containing a histogram oh h0 vs cos(iota) calculate the
# upper limit on h0
def h0ul_from_prior_file(priorfile, ulval=0.95):
    # read in 2D h0 vs cos(iota) binary prior file
    h0bins, cibins, histarr = read_hist_from_file(priorfile)
 
    if not h0bins.any():
        print("Could not read binary histogram file", file=sys.stderr)
        return None
 
    # marginalise over cos(iota)
    h0marg = []
    for i in range(len(h0bins)):
        h0marg.append(trapezoid(histarr[:][i], x=cibins))
 
    # normalise h0 posterior
    h0bins = h0bins - (h0bins[1] - h0bins[0]) / 2
    h0area = trapezoid(h0marg, x=h0bins)
    h0margnorm = map(lambda x: x / h0area, h0marg)
 
    # get cumulative probability
    ct = cumulative_trapezoid(h0margnorm, h0bins)
 
    # prepend a zero to ct
    ct = np.insert(ct, 0, 0)
 
    # use spline interpolation to find the value at 'upper limit'
    ctu, ui = np.unique(ct, return_index=True)
    intf = interp1d(ctu, h0bins[ui], kind="linear")
    return intf(float(ulval))
 
 
# function to create a histogram plot of the 2D posterior
def plot_posterior_hist2D(
    poslist,
    params,
    ifos,
    bounds=None,
    nbins=[50, 50],
    parfile=None,
    overplot=False,
    mplparams=False,
):
    import matplotlib
    from matplotlib import pyplot as plt
    from lalpulsar.pulsarhtmlutils import paramlatexdict
 
    if len(params) != 2:
        print("Require 2 parameters", file=sys.stderr)
        sys.exit(1)
 
    # set some matplotlib defaults for amplitude spectral density
    if not mplparams:
        mplparams = {
            "backend": "Agg",
            "text.usetex": True,  # use LaTeX for all text
            "axes.linewidth": 0.5,  # set axes linewidths to 0.5
            "axes.grid": True,  # add a grid
            "grid.linewidth": 0.5,
            "font.family": "serif",
            "font.size": 12,
        }
 
    matplotlib.rcParams.update(mplparams)
 
    if overplot:
        coldict = {
            "H1": "Reds",
            "H2": "GnBu",
            "L1": "Greens",
            "V1": "Blues",
            "G1": "BuPu",
            "Joint": "Greys",
        }
 
    myfigs = []
 
    # param name for axis label
    try:
        parxaxis = paramlatexdict[params[0].upper()]
    except:
        parxaxis = params[0]
 
    try:
        paryaxis = paramlatexdict[params[1].upper()]
    except:
        paryaxis = params[1]
 
    parval1 = None
    parval2 = None
 
    if parfile:
        parval1 = parfile[params[0].upper()]
        parval2 = parfile[params[1].upper()]
 
    if ifos == None:
        ifos = ["H1"]
 
    for idx, ifo in enumerate(ifos):
        if overplot:
            if idx == 0:
                # if overplotting just have one figure
                myfig = plt.figure(figsize=(4, 4), dpi=200)
 
            if len(ifos) == 1:
                alpha = 1.0
                cmap = plt.cm.get_cmap("Greys")
            else:
                alpha = 0.5  # set transparency
                cmap = plt.cm.get_cmap(coldict[ifo])  # set colormap
                # cmap.set_bad(alpha=0.) # set 'bad' masked values to be transparent
        else:
            myfig = plt.figure(figsize=(4, 4), dpi=200)
            alpha = 1.0  # no transparency
            cmap = plt.cm.get_cmap("Greys")
 
        posterior = poslist[idx]
 
        a = np.squeeze(posterior[params[0]].samples)
        b = np.squeeze(posterior[params[1]].samples)
 
        # Create 2D bin array
        par1pos_min = a.min()
        par2pos_min = b.min()
 
        par1pos_max = a.max()
        par2pos_max = b.max()
 
        plt.xlabel(r"" + parxaxis, fontsize=14, fontweight=100)
        plt.ylabel(r"" + paryaxis, fontsize=14, fontweight=100, rotation=270)
 
        H, xedges, yedges = np.histogram2d(a, b, nbins, normed=True)
 
        # if overplot and len(ifos) > 1:
        #  # mask values that are zero, so that they are
        #  Hm = np.ma.masked_where(np.transpose(H) == 0, np.transpose(H))
        # else:
        #  Hm = np.transpose(H)
        Hm = np.transpose(H)
 
        extent = [xedges[0], xedges[-1], yedges[-1], yedges[0]]
        plt.imshow(
            Hm,
            aspect="auto",
            extent=extent,
            interpolation="bilinear",
            cmap=cmap,
            alpha=alpha,
        )
        # plt.colorbar()
 
        if bounds:
            plt.xlim(bounds[0][0], bounds[0][1])
            plt.ylim(bounds[1][0], bounds[1][1])
 
        # plot injection values if given
        if parval1 and parval2:
            if overplot and idx == len(ifos) - 1:
                # only plot after final posterior has been plotted
                plt.hold(True)
                plt.plot(parval1, parval2, "rx", markersize=8, mew=2)
            else:
                plt.hold(True)
                plt.plot(parval1, parval2, "rx", markersize=8, mew=2)
 
        if overplot:
            plt.hold(True)
        if not overplot:
            myfig.subplots_adjust(left=0.18, bottom=0.15)  # adjust size
            myfigs.append(myfig)
 
    if overplot:
        myfig.subplots_adjust(left=0.18, bottom=0.15)  # adjust size
        myfigs.append(myfig)
 
    return myfigs
 
 
# a function that creates and normalises a histograms of samples, with nbins
# between an upper and lower bound a upper and lower bound. The values at the
# bin points (with length nbins+2) will be returned as numpy arrays
def hist_norm_bounds(samples, nbins, low=float("-inf"), high=float("inf")):
    # get histogram
    n, binedges = np.histogram(samples, nbins)
 
    # get bin width
    binwidth = binedges[1] - binedges[0]
 
    # create bin centres
    bincentres = np.array([])
    for i in range(0, len(binedges) - 1):
        bincentres = np.append(bincentres, binedges[i] + binwidth / 2)
 
    # if histogram points are not close to boundaries (i.e. within a bin of the
    # boundaries) then add zeros to histrogram edges
    if bincentres[0] - binwidth > low:
        # prepend a zero to n
        n = np.insert(n, 0, 0)
 
        # prepend a new bin centre at bincentres[0] - binwidth
        bincentres = np.insert(bincentres, 0, bincentres[0] - binwidth)
    else:
        # we're  closer to the boundary edge than the bin width then, so set a new
        # bin on the boundary with a value linearly extrapolated from the
        # gradiant of the adjacent points
        dx = bincentres[0] - low
 
        # prepend low value to bins
        bincentres = np.insert(bincentres, 0, low)
 
        dn = n[1] - n[0]
 
        nbound = n[0] - (dn / binwidth) * dx
 
        n = n.astype(float)  # convert to floats
 
        # prepend to n
        n = np.insert(n, 0, nbound)
 
    # now the other end!
    if bincentres[-1] + binwidth < high:
        # append a zero to n
        n = np.append(n, 0)
 
        # append a new bin centre at bincentres[end] + binwidth
        bincentres = np.append(bincentres, bincentres[-1] + binwidth)
    else:
        dx = high - bincentres[-1]
 
        # prepend low value to bins
        bincentres = np.append(bincentres, high)
 
        dn = n[-1] - n[-2]
 
        nbound = n[-1] + (dn / binwidth) * dx
 
        n = n.astype(float)  # convert to floats
 
        # prepend to n
        n = np.append(n, nbound)
 
    # now calculate area and normalise
    area = trapezoid(n, x=bincentres)
 
    ns = np.array([])
    for i in range(0, len(bincentres)):
        ns = np.append(ns, float(n[i]) / area)
 
    return ns, bincentres
 
 
# create a Tukey window of length N
def tukey_window(N, alpha=0.5):
    # if alpha >= 1 just return a Hanning window
    if alpha >= 1:
        return np.hanning(N)
 
    # get x values at which to calculate window
    x = np.linspace(0, 1, N)
 
    # initial square window
    win = np.ones(x.shape)
 
    # get the left-hand side of the window  0 <= x < alpha/2
    lhs = x < alpha / 2
    win[lhs] = 0.5 * (1 + np.cos(2 * np.pi / alpha * (x[lhs] - alpha / 2)))
 
    # get right hand side condition 1 - alpha / 2 <= x <= 1
    rhs = x >= (1 - alpha / 2)
    win[rhs] = 0.5 * (1 + np.cos(2 * np.pi / alpha * (x[rhs] - 1 + alpha / 2)))
 
    return win
 
 
# create a function for plotting the absolute value of Bk data (read in from
# data files) and an averaged 1 day "two-sided" amplitude spectral density
# spectrogram for each IFO. Bkdata should be a dictionary contains the files
# keyed to the IFO
def plot_Bks_ASDs(
    Bkdata,
    delt=86400,
    plotpsds=True,
    plotfscan=False,
    removeoutlier=None,
    mplparams=False,
):
    import matplotlib
    from matplotlib.mlab import specgram
    from matplotlib import colors
    from matplotlib import pyplot as plt
 
    # create list of figures
    Bkfigs = []
    psdfigs = []
    fscanfigs = []
    asdlist = []
    sampledt = []
 
    # set some matplotlib defaults for amplitude spectral density
    if not mplparams:
        mplparams = {
            "backend": "Agg",
            "text.usetex": True,  # use LaTeX for all text
            "axes.linewidth": 0.5,  # set axes linewidths to 0.5
            "axes.grid": True,  # add a grid
            "grid.linewidth": 0.5,
            "font.family": "sans-serif",
            "font.sans-serif": "Avant Garde, Helvetica, Computer Modern Sans serif",
            "font.size": 15,
        }
 
    matplotlib.rcParams.update(mplparams)
    matplotlib.rcParams["text.latex.preamble"] = r"\usepackage{xfrac}"
 
    # ifos line colour specs
    coldict = {"H1": "r", "H2": "c", "L1": "g", "V1": "b", "G1": "m"}
    colmapdic = {
        "H1": "Reds",
        "H2": "PuBu",
        "L1": "Greens",
        "V1": "Blues",
        "G1": "PuRd",
    }
 
    # there should be data for each ifo
    for ifo in Bkdata:
        # get data for given ifo
        try:
            Bk = np.loadtxt(Bkdata[ifo], comments=["#", "%"])
        except:
            print("Could not open file %s" % Bkdata[ifo], file=sys.stderr)
            sys.exit(-1)
 
        # sort the array in time
        Bk = Bk[Bk[:, 0].argsort()]
 
        # make sure times are unique
        uargs = np.unique(Bk[:, 0], return_index=True)[1]
        Bk = Bk[uargs, :]
 
        # should be three lines in file
        gpstime = []
 
        # remove outliers at Xsigma by working out sigma from the peak of the
        # distribution of the log absolute value
        if removeoutlier:
            n, binedges = np.histogram(
                np.log(np.fabs(np.concatenate((Bk[:, 1], Bk[:, 2])))), 50
            )
            j = n.argmax(0)
 
            # standard devaition estimate
            stdest = math.exp((binedges[j] + binedges[j + 1]) / 2)
 
            # get values within +/-8 sigma
            # test real parts
            vals = np.where(np.fabs(Bk[:, 1]) < removeoutlier * stdest)
            Bknew = Bk[vals, :][-1]
            # test imag parts
            vals = np.where(np.fabs(Bknew[:, 2]) < removeoutlier * stdest)
            Bk = Bknew[vals, :][-1]
 
        gpstime = Bk[:, 0]
 
        # minimum time step between points (should generally be 60 seconds)
        mindt = min(np.diff(gpstime))
        sampledt.append(mindt)
 
        Bkabs = np.sqrt(Bk[:, 1] ** 2 + Bk[:, 2] ** 2)
 
        # plot the time series of the data
        Bkfig = plt.figure(figsize=(11, 3.5), dpi=200)
        Bkfig.subplots_adjust(bottom=0.15, left=0.09, right=0.94)
 
        tms = list(map(lambda x: x - gpstime[0], gpstime))
 
        plt.plot(tms, Bkabs, ".", color=coldict[ifo], markersize=1)
        plt.xlabel(r"GPS - %d" % int(gpstime[0]), fontsize=14, fontweight=100)
        plt.ylabel(r"$|B_k|$", fontsize=14, fontweight=100)
        # plt.title(r'$B_k$s for ' + ifo.upper(), fontsize=14)
        plt.xlim(tms[0], tms[-1])
 
        Bkfigs.append(Bkfig)
 
        if plotpsds or plotfscan:
            # create PSD by splitting data into days, padding with zeros to give a
            # sample a second, getting the PSD for each day and combining them
            totlen = gpstime[-1] - gpstime[0]  # total data length
 
            # check mindt is an integer and greater than 1
            if math.fmod(mindt, 1) != 0.0 or mindt < 1:
                print(
                    "Error... time steps between data points must be integers",
                    file=sys.stderr,
                )
                sys.exit(-1)
 
            count = 0
 
            # zero pad the data and bin each point in the nearest bin
            datazeropad = np.zeros(int(math.ceil(totlen / mindt)) + 1, dtype=complex)
 
            idx = list(map(lambda x: int(math.floor((x / mindt) + 0.5)), tms))
            for i in range(0, len(idx)):
                datazeropad[idx[i]] = complex(Bk[i, 1], Bk[i, 2])
 
            win = tukey_window(math.floor(delt / mindt), alpha=0.1)
 
            Fs = 1.0 / mindt  # sample rate in Hz
 
            fscan, freqs, t = specgram(
                datazeropad,
                NFFT=int(math.floor(delt / mindt)),
                Fs=Fs,
                window=win,
                noverlap=int(math.floor(delt / (2.0 * mindt))),
            )
 
            if plotpsds:
                fshape = fscan.shape
 
                totalasd = np.zeros(fshape[0])
 
                for i in range(0, fshape[0]):
                    scanasd = np.sqrt(fscan[i, :])
                    nonzasd = np.nonzero(scanasd)
                    scanasd = scanasd[nonzasd]
 
                    # median amplitude spectral density
                    # totalasd[i] = hmean(scanasd)
                    totalasd[i] = np.median(scanasd)
 
                # average amplitude spectral density
                asdlist.append(totalasd)
 
                # plot PSD
                psdfig = plt.figure(figsize=(4, 3.5), dpi=200)
                psdfig.subplots_adjust(left=0.18, bottom=0.15)
 
                plt.plot(freqs, totalasd, color=coldict[ifo])
                plt.xlim(freqs[0], freqs[-1])
                plt.xlabel(r"Frequency (Hz)", fontsize=14, fontweight=100)
                plt.ylabel(r"$h/\sqrt{\rm Hz}$", fontsize=14, fontweight=100)
 
                # convert frequency labels to fractions
                ax = plt.gca()
                xt = [-Fs / 2.0, -Fs / 4.0, 0.0, Fs / 2.0, Fs / 4.0]
                ax.set_xticks(xt)
                xl = []
                for item in xt:
                    if item == 0:
                        xl.append("0")
                    else:
                        if item < 0:
                            xl.append(r"$-\sfrac{1}{%d}$" % (-1.0 / item))
                        else:
                            xl.append(r"$\sfrac{1}{%d}$" % (1.0 / item))
                ax.set_xticklabels(xl)
                # plt.setp(ax.get_xticklabels(), fontsize=16)  # increase font size
                plt.tick_params(axis="x", which="major", labelsize=14)
 
                psdfigs.append(psdfig)
 
            if plotfscan:
                fscanfig = plt.figure(figsize=(11, 3.5), dpi=200)
                fscanfig.subplots_adjust(bottom=0.15, left=0.09, right=0.94)
 
                extent = [tms[0], tms[-1], freqs[0], freqs[-1]]
                plt.imshow(
                    np.sqrt(np.flipud(fscan)),
                    aspect="auto",
                    extent=extent,
                    interpolation=None,
                    cmap=colmapdic[ifo],
                    norm=colors.Normalize(),
                )
                plt.ylabel(r"Frequency (Hz)", fontsize=14, fontweight=100)
                plt.xlabel(r"GPS - %d" % int(gpstime[0]), fontsize=14, fontweight=100)
 
                # convert frequency labels to fractions
                ax = plt.gca()
                yt = [-Fs / 2.0, -Fs / 4.0, 0.0, Fs / 2.0, Fs / 4.0]
                ax.set_yticks(yt)
                yl = []
                for item in yt:
                    if item == 0:
                        yl.append("0")
                    else:
                        if item < 0:
                            yl.append(r"$-\sfrac{1}{%d}$" % (-1.0 / item))
                        else:
                            yl.append(r"$\sfrac{1}{%d}$" % (1.0 / item))
                ax.set_yticklabels(yl)
                plt.tick_params(axis="y", which="major", labelsize=14)
 
                fscanfigs.append(fscanfig)
 
    return Bkfigs, psdfigs, fscanfigs, asdlist, sampledt
 
 
# a function to create a histogram of log10(results) for a list of a given parameter
# (if a list with previous value is given they will be plotted as well)
#  - lims is a dictionary of lists for each IFO
def plot_limits_hist(
    lims, param, ifos, prevlims=None, bins=20, overplot=False, mplparams=False
):
    import matplotlib
    from matplotlib import pyplot as plt
    from lalpulsar.pulsarhtmlutils import paramlatexdict
 
    if not mplparams:
        mplparams = {
            "backend": "Agg",
            "text.usetex": True,  # use LaTeX for all text
            "axes.linewidth": 0.5,  # set axes linewidths to 0.5
            "axes.grid": True,  # add a grid
            "grid.linewidth": 0.5,
            "font.family": "serif",
            "font.size": 12,
        }
 
    matplotlib.rcParams.update(mplparams)
 
    myfigs = []
 
    # ifos line colour specs
    coldict = {"H1": "r", "H2": "c", "L1": "g", "V1": "b", "G1": "m", "Joint": "k"}
 
    try:
        parxaxis = paramlatexdict[param.upper()]
    except:
        parxaxis = param
 
    paryaxis = "Number of Pulsars"
 
    # get log10 of previous results
    logprevlims = None
    if prevlims is not None:
        logprevlims = np.log10(prevlims)
 
    # get the limits of the histogram range
    hrange = None
    stacked = False
    if overplot:
        highbins = []
        lowbins = []
 
        # combine dataset to plot as stacked
        stackdata = []
        stacked = True
 
        for j, ifo in enumerate(ifos):
            theselims = lims[ifo]
 
            # remove any None's
            for i, val in enumerate(theselims):
                if val == None:
                    del theselims[i]
 
            loglims = np.log10(theselims)
 
            stackdata.append(loglims)
 
            highbins.append(np.max(loglims))
            lowbins.append(np.min(loglims))
 
        if logprevlims is not None:
            highbins.append(max(logprevlims))
            lowbins.append(min(logprevlims))
 
        hrange = (min(lowbins), max(highbins))
 
    for j, ifo in enumerate(ifos):
        # get log10 of results
        theselims = lims[ifo]
 
        # remove any None's
        for i, val in enumerate(theselims):
            if val == None:
                del theselims[i]
 
        loglims = np.log10(theselims)
 
        if not overplot:
            stackdata = loglims
 
        # if not overplot or (overplot and j == 0):
        myfig = plt.figure(figsize=(4, 4), dpi=200)
        maxlims = []
        minlims = []
 
        if overplot:
            edgecolor = []
            facecolor = []
            for ifoname in ifos:
                edgecolor.append(coldict[ifoname])
                facecolor.append(coldict[ifoname])
        else:
            edgecolor = [coldict[ifo]]
            facecolor = [coldict[ifo]]
 
        # plt.hist(loglims, bins, range=hrange, histtype='step', fill=True, edgecolor=coldict[ifo], facecolor=coldict[ifo], alpha=0.6)
        n, bins, patches = plt.hist(
            stackdata,
            bins,
            range=hrange,
            histtype="step",
            stacked=stacked,
            fill=True,
            color=edgecolor,
            alpha=0.6,
        )
        for i, patch in enumerate(patches):
            plt.setp(patch, "facecolor", facecolor[i])
 
        if not overplot:
            maxlims.append(np.max(loglims))
            minlims.append(np.min(loglims))
 
        if logprevlims is not None:
            if not overplot:
                maxlims.append(max(logprevlims))
                minlims.append(min(logprevlims))
 
                maxlim = max(maxlims)
                minlim = min(minlims)
            else:
                maxlim = hrange[1]
                minlim = hrange[0]
 
            plt.hold(True)
            plt.hist(
                logprevlims,
                bins,
                range=hrange,
                edgecolor="k",
                lw=2,
                histtype="step",
                fill=False,
            )
            plt.hold(False)
        else:
            if not overplot:
                maxlim = max(maxlims)
                minlim = min(minlims)
            else:
                maxlim = hrange[1]
                minlim = hrange[0]
 
        # set xlabels to 10^x
        ax = plt.gca()
 
        # work out how many ticks to set
        tickvals = range(int(math.ceil(maxlim) - math.floor(minlim)))
        tickvals = map(lambda x: x + int(math.floor(minlim)), tickvals)
        ax.set_xticks(tickvals)
        tls = map(lambda x: "$10^{%d}$" % x, tickvals)
        ax.set_xticklabels(tls)
 
        plt.ylabel(r"" + paryaxis, fontsize=14, fontweight=100)
        plt.xlabel(r"" + parxaxis, fontsize=14, fontweight=100)
 
        myfig.subplots_adjust(left=0.18, bottom=0.15)  # adjust size
 
        myfigs.append(myfig)
 
        if overplot:
            break
 
    return myfigs
 
 
# a function plot upper limits verses GW frequency in log-log space. If upper and
# lower estimates for the upper limit are supplied these will also be plotted as a
# band. If previous limits are supplied these will be plotted as well.
# h0lims is a dictionary of lists of h0 upper limits for each of the IFOs given in
# the list of ifos
# ulesttop and ulestbot are the upper and lower ranges of the estimates limit given
# as a dictionary if list for each ifo.
# prevlim is a list of previous upper limits at the frequencies prevlimf0gw
# xlims is the frequency range for the plot
# overplot - set to true to plot different IFOs on same plot
def plot_h0_lims(
    h0lims,
    f0gw,
    ifos,
    xlims=[10, 1500],
    ulesttop=None,
    ulestbot=None,
    prevlim=None,
    prevlimf0gw=None,
    overplot=False,
    mplparams=False,
):
    import matplotlib
    from matplotlib import pyplot as plt
 
    if not mplparams:
        mplparams = {
            "backend": "Agg",
            "text.usetex": True,  # use LaTeX for all text
            "axes.linewidth": 0.5,  # set axes linewidths to 0.5
            "axes.grid": True,  # add a grid
            "grid.linewidth": 0.5,
            "font.family": "serif",
            "font.size": 12,
        }
 
    matplotlib.rcParams.update(mplparams)
 
    myfigs = []
 
    # ifos line colour specs
    coldict = {"H1": "r", "H2": "c", "L1": "g", "V1": "b", "G1": "m", "Joint": "k"}
 
    parxaxis = "Frequency (Hz)"
    paryaxis = "$h_0$"
 
    for j, ifo in enumerate(ifos):
        h0lim = h0lims[ifo]
 
        if len(ifos) > 1:
            f0s = f0gw[ifo]
        else:
            f0s = f0gw
 
        if len(h0lim) != len(f0s):
            return None  # exit if lengths aren't correct
 
        if not overplot or (overplot and j == 0):
            myfig = plt.figure(figsize=(7, 5.5), dpi=200)
 
        plt.hold(True)
        if ulesttop is not None and ulestbot is not None:
            ult = ulesttop[ifo]
            ulb = ulestbot[ifo]
 
            if len(ult) == len(ulb) and len(ult) == len(f0s):
                for i in range(len(ult)):
                    plt.loglog(
                        [f0s[i], f0s[i]], [ulb[i], ult[i]], ls="-", lw=3, c="lightgrey"
                    )
 
        if (
            prevlim is not None
            and prevlimf0gw is not None
            and len(prevlim) == len(prevlimf0gw)
        ):
            if not overplot or (overplot and j == len(ifos) - 1):
                plt.loglog(
                    prevlimf0gw,
                    prevlim,
                    marker="*",
                    ms=10,
                    alpha=0.7,
                    mfc="None",
                    mec="k",
                    ls="None",
                )
 
        # plot current limits
        plt.loglog(
            f0s, h0lim, marker="*", ms=10, mfc=coldict[ifo], mec=coldict[ifo], ls="None"
        )
 
        plt.ylabel(r"" + paryaxis, fontsize=14, fontweight=100)
        plt.xlabel(r"" + parxaxis, fontsize=14, fontweight=100)
 
        plt.xlim(xlims[0], xlims[1])
 
        if not overplot or (overplot and j == len(ifos) - 1):
            plt.hold(False)
            myfig.subplots_adjust(left=0.12, bottom=0.10)  # adjust size
            myfigs.append(myfig)
 
    return myfigs
 
 
# a function to create the signal model for a heterodyned triaxial pulsar given
# a signal GPS start time, signal duration (in seconds), sample interval
# (seconds), detector, and the parameters h0, cos(iota), psi (rads), initial
# phase phi0 (rads), right ascension (rads) and declination (rads) in a
# dictionary. The list of time stamps, and the real and imaginary parts of the
# signal are returned
def heterodyned_triaxial_pulsar(starttime, duration, dt, detector, pardict):
    # create a list of times stamps
    ts = []
    tmpts = starttime
 
    # create real and imaginary parts of the signal
    s = np.array([], dtype=complex)
 
    sphi = np.sin(pardict["phi0"])
    cphi = np.cos(pardict["phi0"])
 
    Xplus = 0.25 * (1.0 + pardict["cosiota"] * pardict["cosiota"]) * pardict["h0"]
    Xcross = 0.5 * pardict["cosiota"] * pardict["h0"]
    Xpsinphi = Xplus * sphi
    Xcsinphi = Xcross * sphi
    Xpcosphi = Xplus * cphi
    Xccosphi = Xcross * cphi
 
    i = 0
    while tmpts < starttime + duration:
        ts.append(starttime + (dt * i))
 
        # get the antenna response
        fp, fc = antenna_response(
            ts[i], pardict["ra"], pardict["dec"], pardict["psi"], detector
        )
 
        # create real part of signal
        s = np.append(
            s, (fp * Xpcosphi + fc * Xcsinphi) + 1j * (fp * Xpsinphi - fc * Xccosphi)
        )
 
        tmpts = ts[i] + dt
 
        i = i + 1
 
    return ts, s
 
 
# a function to create a heterodyned signal model for a neutron star using the complex
# amplitude parameters C22, phi22, C21 and phi21 and the orientation parameters cos(iota)
# and psi. If both C22 and C21 are non-zero then a signal at both the rotation frequency
# and twice the rotation frequency will be generated.
#
# Note that for the purely triaxial signal the waveform, as defined in Jones (2015)
# http://arxiv.org/abs/1501.05832 (and subsequently Pitkin et al (2015) http://arxiv.org/abs/1508.00416),
# has the opposite sign to the convention in e.g. JKS (which is used for signal generation in hardware
# injections). Hence, when converting h0, from the conventional model, to C22 a minus sign needs to be
# introduced. Also, phi0 here is defined as the rotational phase of the source, so when converting to phi22
# a factor of 2 is needed.
def heterodyned_pulsar_signal(
    pardict, detector, starttime=900000000.0, duration=86400.0, dt=60.0, datatimes=None
):
    if "cosiota" in pardict:
        cosiota = pardict["cosiota"]
        iota = math.acos(cosiota)
        siniota = math.sin(iota)
    else:
        raise KeyError("cos(iota) not defined!")
 
    if "psi" in pardict:
        psi = pardict["psi"]
    else:
        raise KeyError("psi not defined!")
 
    if "C22" in pardict:
        C22 = pardict["C22"]
    else:
        if "h0" in pardict:
            C22 = -pardict["h0"] / 2.0
        else:
            C22 = 0.0
 
    if "phi22" in pardict:
        phi22 = pardict["phi22"]
    else:
        if "phi0" in pardict:
            phi22 = 2.0 * pardict["phi0"]
        else:
            phi22 = 0.0
 
    ePhi22 = cmath.exp(phi22 * 1j)
 
    if "C21" in pardict:
        C21 = pardict["C21"]
    else:
        C21 = 0.0
 
    if "phi21" in pardict:
        phi21 = pardict["phi21"]
    else:
        phi21 = 0.0
 
    ePhi21 = cmath.exp(phi21 * 1j)
 
    s = []  # signal
    ts = []  # times
 
    if "C21" in pardict and "C22" in pardict and "h0" not in pardict:
        freqs = [1.0, 2.0]
    else:
        freqs = [2.0]
 
    for f in freqs:
        i = 0
        if datatimes is None:
            datatimes = np.arange(starttime, starttime + duration, dt)
 
        # get the antenna response
        fp, fc = antenna_response(
            datatimes, pardict["ra"], pardict["dec"], pardict["psi"], detector
        )
 
        if f == 1.0:
            sf = (
                -(C21 / 4.0) * ePhi21 * siniota * cosiota * fp
                + 1j * (C21 / 4.0) * ePhi21 * siniota * fc
            )
        elif f == 2.0:
            sf = (
                -(C22 / 2.0) * ePhi22 * (1.0 + cosiota**2.0) * fp
                + 1j * (C22) * ePhi22 * cosiota * fc
            )
 
        ts.append(datatimes)
        s.append(sf)
 
    return ts, s
 
 
# a function to convert the pinned superfluid model parameters to the complex
# amplitude parameters using the Eqns 76-79 defined in Jones 2012 LIGO DCC T1200265-v3
def convert_model_parameters(pardict):
    costheta = pardict["costheta"]
    theta = np.arccos(costheta)
    sintheta = math.sin(theta)
    sin2theta = math.sin(2.0 * theta)
    sinlambda = math.sin(pardict["lambdapin"])
    coslambda = math.cos(pardict["lambdapin"])
    sin2lambda = math.sin(2.0 * pardict["lambdapin"])
 
    phi0 = pardict["phi0"]
 
    I21 = pardict["I21"]
    I31 = pardict["I31"]
 
    A22 = I21 * (sinlambda**2 - (coslambda * costheta) ** 2) - I31 * sintheta**2
    B22 = I21 * sin2lambda * costheta
 
    C22 = 2.0 * math.sqrt(A22**2 + B22**2)
 
    A21 = I21 * sin2lambda * sintheta
    B21 = (I21 * coslambda**2 - I31) * sin2theta
 
    C21 = 2.0 * math.sqrt(A21**2 + B21**2)
 
    phi22 = np.fmod(2.0 * phi0 - math.atan2(B22, A22), 2.0 * np.pi)
    phi21 = np.fmod(phi0 - math.atan2(B21, A21), 2.0 * np.pi)
 
    outvals = {"C22": C22, "C21": C21, "phi22": phi22, "phi21": phi21}
 
    return outvals
 
 
# a function to create the signal model for a heterodyned pinned superfluid
# model pulsar a signal GPS start time, signal duration (in seconds),
# sample interval (seconds), detector, and the parameters I21, I31,
# cos(theta), lambda, cos(iota), psi (rads), initial phase phi0 (for l=m=2 mode in rads), right
# ascension (rads), declination (rads), distance (kpc) and frequency (f0) in a
# dictionary. The list of time stamps, and the real and imaginary parts of the
# 1f and 2f signals are returned in an array
def heterodyned_pinsf_pulsar(starttime, duration, dt, detector, pardict):
    iota = np.arccos(pardict["cosiota"])
    theta = np.arccos(pardict["costheta"])
    siniota = math.sin(iota)
    sintheta = math.sin(theta)
    sin2theta = math.sin(2.0 * theta)
    sinlambda = math.sin(pardict["lambdapin"])
    coslambda = math.cos(pardict["lambdapin"])
    sin2lambda = math.sin(2.0 * pardict["lambdapin"])
 
    ePhi = cmath.exp(0.5 * pardict["phi0"] * 1j)
    e2Phi = cmath.exp(pardict["phi0"] * 1j)
 
    f2_r = pardict["f0"] ** 2 / pardict["dist"]
 
    """
      This model is a complex heterodyned time series for a pinned superfluid
      neutron star emitting at its roation frequency and twice its rotation
      frequency (as defined in Eqns 35-38 of Jones 2012 LIGO DCC T1200265-v3)
    """
 
    Xplusf = -(f2_r / 2.0) * siniota * pardict["cosiota"]
    Xcrossf = (f2_r / 2.0) * siniota
 
    Xplus2f = -f2_r * (1.0 + pardict["cosiota"] ** 2)
    Xcross2f = 2.0 * f2_r * pardict["cosiota"]
 
    A21 = pardict["I21"] * sin2lambda * sintheta
    B21 = (pardict["I21"] * coslambda**2 - pardict["I31"]) * sin2theta
 
    A22 = (
        pardict["I21"] * (sinlambda**2 - coslambda**2 * pardict["costheta"] ** 2)
        - pardict["I31"] * sintheta**2
    )
    B22 = pardict["I21"] * sin2lambda * pardict["costheta"]
 
    # create a list of times stamps
    ts1 = []
    ts2 = []
    tmpts = starttime
 
    # create real and imaginary parts of the 1f signal
    s1 = np.array([], dtype=complex)
    s2 = np.array([], dtype=complex)
 
    i = 0
    while tmpts < starttime + duration:
        ts1.append(starttime + (dt * i))
        ts2.append(starttime + (dt * i))
 
        # get the antenna response
        fp, fc = antenna_response(
            ts1[i], pardict["ra"], pardict["dec"], pardict["psi"], detector
        )
 
        # create the complex signal amplitude model at 1f
        s1 = np.append(
            s1,
            (fp * Xplusf * ePhi * (A21 - 1j * B21))
            + (fc * Xcrossf * ePhi * (B21 + 1j * A21)),
        )
 
        # create the complex signal amplitude model at 2f
        s2 = np.append(
            s2,
            (fp * Xplus2f * e2Phi * (A22 - 1j * B22))
            + (fc * Xcross2f * e2Phi * (B22 + 1j * A22)),
        )
 
        tmpts = ts1[i] + dt
 
        i = i + 1
 
    # combine data into 1 array
    ts = np.vstack([ts1, ts2])
    s = np.vstack([s1, s2])
 
    return ts, s
 
 
# function to get the antenna response for a given detector. It takes in a GPS time or list (or numpy array)
# of GPS times, right  ascension (rads), declination (rads), polarisation angle (rads) and a detector name
# e.g. H1, L1, V1. The plus and cross polarisations are returned.
def antenna_response(gpsTime, ra, dec, psi, det):
    import lal
 
    # create detector-name map
    detMap = {
        "H1": lal.LALDetectorIndexLHODIFF,
        "H2": lal.LALDetectorIndexLHODIFF,
        "L1": lal.LALDetectorIndexLLODIFF,
        "G1": lal.LALDetectorIndexGEO600DIFF,
        "V1": lal.LALDetectorIndexVIRGODIFF,
        "T1": lal.LALDetectorIndexTAMA300DIFF,
        "AL1": lal.LALDetectorIndexLLODIFF,
        "AH1": lal.LALDetectorIndexLHODIFF,
        "AV1": lal.LALDetectorIndexVIRGODIFF,
        "E1": lal.LALDetectorIndexE1DIFF,
        "E2": lal.LALDetectorIndexE2DIFF,
        "E3": lal.LALDetectorIndexE3DIFF,
    }
 
    try:
        detector = detMap[det]
    except KeyError:
        raise KeyError("ERROR. Key {} is not a valid detector name.".format(det))
 
    # get detector
    detval = lal.CachedDetectors[detector]
    response = detval.response
 
    # check if ra and dec are floats or strings (if strings in hh/dd:mm:ss.s format then convert to rads)
    if not isinstance(ra, float):
        if isinstance(ra, str):
            try:
                ra = ra_to_rad(ra)
            except:
                raise ValueError(
                    "Right ascension string '{}' not formatted properly".format(ra)
                )
        else:
            raise ValueError(
                "Right ascension must be a 'float' in radians or a string of the format 'hh:mm:ss.s'"
            )
 
    if not isinstance(dec, float):
        if isinstance(dec, str):
            try:
                dec = dec_to_rad(dec)
            except:
                raise ValueError(
                    "Declination string '{}' not formatted properly".format(ra)
                )
        else:
            raise ValueError(
                "Declination must be a 'float' in radians or a string of the format 'dd:mm:ss.s'"
            )
 
    # check if gpsTime is just a float or int, and if so convert into an array
    if isinstance(gpsTime, float) or isinstance(gpsTime, int):
        gpsTime = np.array([gpsTime])
    else:  # make sure it's a numpy array
        gpsTime = np.copy(gpsTime)
 
    # make sure times are floats
    gpsTime = gpsTime.astype("float64")
 
    # if gpsTime is a list of regularly spaced values then use ComputeDetAMResponseSeries
    if len(gpsTime) == 1 or np.unique(np.diff(gpsTime)).size == 1:
        gpsStart = lal.LIGOTimeGPS(gpsTime[0])
        dt = 0.0
        if len(gpsTime) > 1:
            dt = gpsTime[1] - gpsTime[0]
        fp, fc = lal.ComputeDetAMResponseSeries(
            response, ra, dec, psi, gpsStart, dt, len(gpsTime)
        )
 
        # return elements from Time Series
        return fp.data.data, fc.data.data
    else:  # we'll have to work out the value at each point in the time series
        fp = np.zeros(len(gpsTime))
        fc = np.zeros(len(gpsTime))
 
        for i in range(len(gpsTime)):
            gps = lal.LIGOTimeGPS(gpsTime[i])
            gmst_rad = lal.GreenwichMeanSiderealTime(gps)
 
            # actual computation of antenna factors
            fp[i], fc[i] = lal.ComputeDetAMResponse(response, ra, dec, psi, gmst_rad)
 
        return fp, fc
 
 
# a function to inject a heterodyned pulsar signal into noise of a
# given level for detectors. it will return the signal + noise and the optimal
# SNR. If an snrscale value is passed to the function the signal will be scaled
# to that SNR. nsigs
def inject_pulsar_signal(
    starttime,
    duration,
    dt,
    detectors,
    pardict,
    freqfac=[2.0],
    npsds=None,
    snrscale=None,
):
    # if detectors is just a string (i.e. one detector) then make it a list
    if isinstance(detectors, str):
        detectors = [detectors]
 
    # if not noise sigma's are given then generate a noise level from the given
    # detector noise curves
    if npsds is None:
        npsds = []
 
        for det in detectors:
            for frf in freqfac:
                psd = detector_noise(det, frf * pardict["f0"])
 
                # convert to time domain standard devaition shared between real and
                # imaginary signal
                ns = np.sqrt((psd / 2.0) / (2.0 * dt))
 
                npsds.append(ns)
    else:
        # convert input psds into time domain noise standard deviation
        tmpnpsds = []
 
        count = 0
        for j, det in enumerate(detectors):
            for frf in freqfac:
                if len(npsds) == 1:
                    tmpnpsds.append(np.sqrt((npsds[0] / 2.0) / (2.0 * dt)))
                else:
                    tmpnpsds.append(np.sqrt((npsds[count] / 2.0) / (2.0 * dt)))
 
                count = count + 1
 
        npsds = tmpnpsds
 
    if len(freqfac) == 1 and len(detectors) != len(npsds):
        raise ValueError(
            "Number of detectors {} not the same as number of noises {}".format(
                len(detectors), len(npsds)
            )
        )
 
    if len(freqfac) == 2 and 2 * len(detectors) != len(npsds):
        raise ValueError(
            "Number of detectors {} not half the number of noises {}".format(
                len(detectors), len(npsds)
            )
        )
 
    tss = np.array([])
    ss = np.array([])
 
    snrtot = 0
    for j, det in enumerate(detectors):
        # create the pulsar signal
        if len(freqfac) == 1 and freqfac[0] == 2.0:
            ts, s = heterodyned_pulsar_signal(pardict, det, starttime, duration, dt)
 
            if j == 0:
                tss = np.append(tss, ts)
                ss = np.append(ss, s)
            else:
                tss = np.vstack([tss, ts])
                ss = np.vstack([ss, s])
 
            # get SNR
            snrtmp = get_optimal_snr(s[0], npsds[j])
        elif len(freqfac) == 2:
            ts, s = heterodyned_pulsar_signal(pardict, det, starttime, duration, dt)
 
            snrtmp = 0
            for k, frf in enumerate(freqfac):
                if j == 0 and k == 0:
                    tss = np.append(tss, ts[k][:])
                    ss = np.append(ss, s[k][:])
                else:
                    tss = np.vstack([tss, ts[k][:]])
                    ss = np.vstack([ss, s[k][:]])
 
                snrtmp2 = get_optimal_snr(s[k][:], npsds[2 * j + k])
                snrtmp = snrtmp + snrtmp2 * snrtmp2
 
            snrtmp = np.sqrt(snrtmp)
 
        snrtot = snrtot + snrtmp * snrtmp
 
    # total multidetector/data stream snr
    snrtot = np.sqrt(snrtot)
 
    # add noise and rescale signals if necessary
    if snrscale is not None:
        if snrscale != 0:
            snrscale = snrscale / snrtot
        # print snrscale
    else:
        snrscale = 1
 
    signalonly = ss.copy()
 
    i = 0
    for det in detectors:
        # for triaxial model
        if len(freqfac) == 1:
            # generate random numbers
            rs = np.random.randn(len(ts[0]), 2)
 
            for j, t in enumerate(ts[0]):
                if len(tss.shape) == 1:
                    signalonly[j] = (snrscale * ss[j].real) + 1j * (
                        snrscale * ss[j].imag
                    )
                    ss[j] = (snrscale * ss[j].real + npsds[i] * rs[j][0]) + 1j * (
                        snrscale * ss[j].imag + npsds[i] * rs[j][1]
                    )
                else:
                    signalonly[i][j] = (snrscale * ss[i][j].real) + 1j * (
                        snrscale * ss[i][j].imag
                    )
                    ss[i][j] = (snrscale * ss[i][j].real + npsds[i] * rs[j][0]) + 1j * (
                        snrscale * ss[i][j].imag + npsds[i] * rs[j][1]
                    )
 
            i = i + 1
        elif len(freqfac) == 2:
            # generate random numbers
            rs = np.random.randn(len(ts[0][:]), 4)
 
            for j, t in enumerate(ts[0][:]):
                signalonly[i][j] = (snrscale * ss[i][j].real) + 1j * (
                    snrscale * ss[i][j].imag
                )
                ss[i][j] = (snrscale * ss[i][j].real + npsds[i] * rs[j][0]) + 1j * (
                    snrscale * ss[i][j].imag + npsds[i] * rs[j][1]
                )
                signalonly[i + 1][j] = (snrscale * ss[i + 1][j].real) + 1j * (
                    snrscale * ss[i + 1][j].imag
                )
                ss[i + 1][j] = (
                    snrscale * ss[i + 1][j].real + npsds[i + 1] * rs[j][2]
                ) + 1j * (snrscale * ss[i + 1][j].imag + npsds[i + 1] * rs[j][3])
 
            i = i + 2
        else:
            print("Something wrong with injection", file=sys.stderr)
            sys.exit(1)
 
    snrtot = snrtot * snrscale
 
    return tss, ss, signalonly, snrtot, snrscale, npsds
 
 
# function to create a time domain PSD from theoretical
# detector noise curves. It takes in the detector name and the frequency at
# which to generate the noise.
#
# The noise models are taken from those in lalsimulation/lib/LALSimNoisePSD.c
def detector_noise(det, f):
    import lalsimulation
 
    if det == "AV1":  # Advanced Virgo
        return lalsimulation.SimNoisePSDAdvVirgo(f)
    elif det == "H1" or det == "L1":  # iLIGO SRD
        return lalsimulation.SimNoisePSDiLIGOSRD(f)
    elif det == "H2":
        return lalsimulation.SimNoisePSDiLIGOSRD(f) * 2.0
    elif det == "G1":  # GEO_600
        return lalsimulation.SimNoisePSDGEO(f)
    elif det == "V1":  # initial Virgo
        return lalsimulation.SimNoisePSDVirgo(f)
    elif det == "T1":  # TAMA
        return lalsimulation.SimNoisePSDTAMA(f)
    elif det == "K1":  # KAGRA
        return lalsimulation.SimNoisePSDKAGRA(f)
    elif det == "AL1" or det == "AH1":
        return lalsimulation.SimNoisePSDaLIGOZeroDetHighPower(f)
    else:
        raise ValueError("{} is not a recognised detector".format(det))
 
 
# function to calculate the optimal SNR of a heterodyned pulsar signal - it
# takes in a complex signal model and noise standard deviation
def get_optimal_snr(s, sig):
    ss = 0
    # sum square of signal
    for val in s:
        ss = ss + val.real**2 + val.imag**2
 
    return np.sqrt(ss / sig**2)
 
 
# use the Gelman-Rubins convergence test for MCMC chains, where chains is a list
# of MCMC numpy chain arrays - this copies the gelman_rubins function in
# bayespputils
def gelman_rubins(chains):
    chainMeans = [np.mean(data) for data in chains]
    chainVars = [np.var(data) for data in chains]
 
    BoverN = np.var(chainMeans)
    W = np.mean(chainVars)
 
    sigmaHat2 = W + BoverN
 
    m = len(chains)
 
    VHat = sigmaHat2 + BoverN / m
 
    R = VHat / W
 
    return R
 
 
# function to convert MCMC files output from pulsar_parameter_estimation into
# a posterior class object - also outputting:
#  - the mean effective sample size of each parameter chain
#  - the Gelman-Rubins statistic for each parameter (a dictionary)
#  - the original length of each chain
# Th input is a list of MCMC chain files
def pulsar_mcmc_to_posterior(chainfiles):
    from lalinference import bayespputils as bppu
 
    cl = []
    neffs = []
    grr = {}
 
    mcmc = []
 
    for cfile in chainfiles:
        if os.path.isfile(cfile):
            # load MCMC chain
            mcmcChain = read_pulsar_mcmc_file(cfile)
 
            # if no chain found then exit with None's
            if mcmcChain is None:
                return None, None, None, None
 
            # find number of effective samples for the chain
            neffstmp = []
            for j in range(1, mcmcChain.shape[1]):
                neff, acl, acf = bppu.effectiveSampleSize(mcmcChain[:, j])
                neffstmp.append(neff)
 
            # get the minimum effective sample size
            # neffs.append(min(neffstmp))
            # get the mean effective sample size
            neffs.append(math.floor(np.mean(neffstmp)))
 
            # nskip = math.ceil(mcmcChain.shape[0]/min(neffstmp))
            nskip = int(math.ceil(mcmcChain.shape[0] / np.mean(neffstmp)))
 
            # output every nskip (independent) value
            mcmc.append(mcmcChain[::nskip, :])
            cl.append(mcmcChain.shape[0])
        else:
            print("File %s does not exist!" % cfile, file=sys.stderr)
            return None, None, None, None
 
    # output data to common results format
    # get first line of MCMC chain file for header names
    cf = open(chainfiles[0], "r")
    headers = cf.readline()
    headers = cf.readline()  # column names are on the second line
    # remove % from start
    headers = re.sub("%", "", headers)
    # remove rads
    headers = re.sub("rads", "", headers)
    # remove other brackets e.g. around (iota)
    headers = re.sub("[()]", "", headers)
    cf.close()
 
    # get Gelman-Rubins stat for each parameter
    for idx, parv in enumerate(headers.split()):
        lgr = []
        if parv != "logL":
            for j in range(0, len(mcmc)):
                achain = mcmc[j]
                singlechain = achain[:, idx]
                lgr.append(singlechain)
            grr[parv.lower()] = gelman_rubins(lgr)
 
    # logL in chain is actually log posterior, so also output the posterior
    # values (can be used to estimate the evidence)
    headers = headers.replace("\n", "\tpost\n")
 
    # output full data to common format
    comfile = chainfiles[0] + "_common_tmp.dat"
    try:
        cf = open(comfile, "w")
    except:
        print("Can't open common posterior file!", file=sys.stderr)
        sys.exit(0)
 
    cf.write(headers)
    for narr in mcmc:
        for j in range(0, narr.shape[0]):
            mline = narr[j, :]
            # add on posterior
            mline = np.append(mline, np.exp(mline[0]))
 
            strmline = " ".join(str(x) for x in mline) + "\n"
            cf.write(strmline)
    cf.close()
 
    # read in as common object
    peparser = bppu.PEOutputParser("common")
    cf = open(comfile, "r")
    commonResultsObj = peparser.parse(cf)
    cf.close()
 
    # remove temporary file
    os.remove(comfile)
 
    # create posterior class
    pos = bppu.Posterior(commonResultsObj, SimInspiralTableEntry=None, votfile=None)
 
    # convert iota back to cos(iota)
    # create 1D posterior class of cos(iota) values
    cipos = None
    cipos = bppu.PosteriorOneDPDF("cosiota", np.cos(pos["iota"].samples))
 
    # add it back to posterior
    pos.append(cipos)
 
    # remove iota samples
    pos.pop("iota")
 
    return pos, neffs, grr, cl
 
 
# a function that attempt to load an pulsar MCMC chain file: first it tries using
# numpy.loadtxt; if it fails it tries reading line by line and checking
# for consistent line numbers, skipping lines that are inconsistent; if this
# fails it returns None
def read_pulsar_mcmc_file(cf):
    cfdata = None
 
    # first try reading in with loadtxt (skipping lines starting with %)
    try:
        cfdata = np.loadtxt(cf, comments="%")
    except:
        try:
            fc = open(cf, "r")
 
            # read in header lines and count how many values each line should have
            headers = fc.readline()
            headers = fc.readline()  # column names are on the second line
            fc.close()
            # remove % from start
            headers = re.sub("%", "", headers)
            # remove rads
            headers = re.sub("rads", "", headers)
            # remove other brackets e.g. around (iota)
            headers = re.sub("[()]", "", headers)
 
            lh = len(headers.split())
 
            cfdata = np.array([])
 
            lines = cf.readlines()
 
            for i, line in enumerate(lines):
                if "%" in line:  # skip lines containing %
                    continue
 
                lvals = line.split()
 
                # skip line if number of values isn't consistent with header
                if len(lvals) != lh:
                    continue
 
                # convert values to floats
                try:
                    lvalsf = map(float, lvals)
                except:
                    continue
 
                # add values to array
                if i == 0:
                    cfdata = np.array(lvalsf)
                else:
                    cfdata = np.vstack((cfdata, lvalsf))
 
            if cfdata.size == 0:
                cfdata = None
        except:
            cfdata = None
 
    return cfdata
 
 
def pulsar_nest_to_posterior(postfile, nestedsamples=False, removeuntrig=True):
    """
    This function will import a posterior sample file created by `lalapps_nest2pos` (or a nested
    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
        The file name of the posterior or nested sample file. In general this should be a HDF5 file with the extension
        '.hdf' or '.h5', although older format ascii files are still supported at the moment.
    nestedsamples : bool, default: False
        If the file being input contains nested samples, rather than posterior samples, then this flag must be set to
        True
    removeuntrig : bool, default: True
        If this is True then any parameters that are sines or cosines of a value will have the value removed if present
        e.g. if cosiota and iota exist then iota will be removed.
    """
 
    from lalinference import bayespputils as bppu
    from lalinference.bayespputils import replace_column
 
    fe = os.path.splitext(postfile)[-1].lower()  # file extension
 
    # combine multiple nested sample files for an IFO into a single posterior (copied from lalapps_nest2pos)
    if fe == ".h5" or fe == ".hdf":  # HDF5 file
        if nestedsamples:
            # use functions from lalapps_nest2pos to read values from nested sample files i.e. not a posterior file created by lalapps_nest2pos
            try:
                from lalinference.io import read_samples
                from lalinference import LALInferenceHDF5NestedSamplesDatasetName
 
                samples = read_samples(
                    postfile, tablename=LALInferenceHDF5NestedSamplesDatasetName
                )
                params = samples.colnames
 
                # make everything a float, since that's what's excected of a CommonResultsObj
                for param in params:
                    replace_column(samples, param, samples[param].astype(float))
 
                nsResultsObject = (
                    samples.colnames,
                    samples.as_array().view(float).reshape(-1, len(params)),
                )
            except:
                raise IOError("Could not import nested samples")
        else:  # assume posterior file has been created with lalapps_nest2pos
            peparser = bppu.PEOutputParser("hdf5")
            nsResultsObject = peparser.parse(postfile)
    elif fe == ".gz":  # gzipped file
        import gzip
 
        peparser = bppu.PEOutputParser("common")
        nsResultsObject = peparser.parse(gzip.open(postfile, "r"))
    else:  # assume an ascii text file
        peparser = bppu.PEOutputParser("common")
        nsResultsObject = peparser.parse(open(postfile, "r"))
 
    pos = bppu.Posterior(nsResultsObject, SimInspiralTableEntry=None)
 
    # remove any unchanging variables and randomly shuffle the rest
    pnames = pos.names
    nsamps = len(pos[pnames[0]].samples)
    permarr = np.arange(nsamps)
    np.random.shuffle(permarr)
    posdist = None
    posfreqs = None
    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
        ):
            if pname == "f0_fixed":
                # try getting a fixed f0 value (for calculating h0 from Q22)
                posfreqs = pos[pname].samples[0]
            elif pname == "dist":
                # try getting a fixed distance value (for calculating h0 from Q22 if required)
                posdist = pos[pname].samples[0] / KPC
 
            pos.pop(pname)
        else:
            # shuffle
            shufpos = None
            shufpos = bppu.PosteriorOneDPDF(pname, pos[pname].samples[permarr])
            pos.pop(pname)
            pos.append(shufpos)
 
    # check whether iota has been used
    posIota = None
    if "iota" in pos.names:
        posIota = pos["iota"].samples
 
    if posIota is not None:
        cipos = None
        cipos = bppu.PosteriorOneDPDF("cosiota", np.cos(posIota))
        pos.append(cipos)
        if removeuntrig:
            pos.pop("iota")
 
    # check whether sin(i) binary parameter has been used
    posI = None
    if "i" in pos.names:
        posI = pos["i"].samples
 
    if posI is not None:
        sinipos = None
        sinipos = bppu.PosteriorOneDPDF("sini", np.sin(posI))
        pos.append(sinipos)
        if removeuntrig:
            pos.pop("i")
 
    # convert C22 back into h0, and phi22 back into phi0 if required
    posC21 = None
    if "c21" in pos.names:
        posC21 = pos["c21"].samples
 
    posC22 = None
    if "c22" in pos.names:
        posC22 = pos["c22"].samples
 
    posphi22 = None
    if "phi22" in pos.names:
        posphi22 = pos["phi22"].samples
 
    # convert C22 back into h0, and phi22 back into phi0 if required
    if posC22 is not None and posC21 is None:
        h0pos = None
        h0pos = bppu.PosteriorOneDPDF("h0", 2.0 * pos["c22"].samples)
 
        pos.append(h0pos)
        pos.pop("c22")
 
        if posphi22 is not None:
            phi0pos = None
            phi0pos = bppu.PosteriorOneDPDF(
                "phi0", np.fmod(pos["phi22"].samples + math.pi, 2.0 * math.pi)
            )
 
            pos.append(phi0pos)
            pos.pop("phi22")
 
    # convert Q22 to h0 is required
    posQ22 = None
    if "q22" in pos.names:
        posQ22 = pos["q22"].samples
 
    if "dist" in pos.names:
        posdist = (
            pos["dist"].samples / KPC
        )  # distance in kpc (for use in converting Q22 to h0)
 
        # convert distance samples to kpc
        pos.pop("dist")
        distpos = bppu.PosteriorOneDPDF("dist", posdist)
        pos.append(distpos)
 
    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
        hdf = h5py.File(postfile, "r")
        a = hdf["lalinference"]["lalinference_nest"]
        sigev = a.attrs["log_evidence"]
        noiseev = a.attrs["log_noise_evidence"]
        hdf.close()
    else:
        B = np.loadtxt(postfile.replace(".gz", "") + "_B.txt")
        sigev = B[1]
        noiseev = B[2]
 
    # return posterior samples, signal evidence (B[1]) and noise evidence (B[2])
    return pos, sigev, noiseev
 
 
# function to add two exponentiated log values and return the log of the result
def logplus(x, y):
    return np.logaddexp(x, y)
 
 
def logtrapz(lnf, dx):
    """
    Perform trapezium rule integration for the logarithm of a function on a regular grid.
 
    Inputs
    ------
    lnf - a numpy array of values that are the natural logarithm of a function
    dx  - a float giving the step size between values in the function
 
    Returns
    -------
    The natural logarithm of the area under the function.
    """
 
    return np.log(dx / 2.0) + logsumexp([logsumexp(lnf[:-1]), logsumexp(lnf[1:])])
 
 
# function to marginalise over a parameter
def marginalise(like, pname, pnames, ranges):
    # like - a copy of the loglikelihood array
    # pname - the parameter to marginalise over
    # pnames - a copy of the list of parameters that haven't been marginalised over
    # ranges - a copy of the dictionary of ranges that haven't been marginalised over
 
    places = ranges[pname]
    axis = pnames.index(pname)
 
    # perform marginalisation
    if len(places) > 1:
        x = np.apply_along_axis(logtrapz, axis, like, places[1] - places[0])
    elif len(places) == 1:
        # no marginalisation required just remove the specific singleton dimension via reshaping
        z = like.shape
        q = np.arange(0, len(z)).astype(int) != axis
        newshape = tuple((np.array(list(z)))[q])
        x = np.reshape(like, newshape)
 
    # remove name from pnames and ranges
    del ranges[pname]
    pnames.remove(pname)
 
    return x
 
 
# return the log marginal posterior on a given parameter (if 'all' marginalise over all parameters)
def marginal(lnlike, pname, pnames, ranges, multflatprior=False):
    from copy import deepcopy
 
    # make copies of everything
    lnliketmp = deepcopy(lnlike)
    pnamestmp = deepcopy(pnames)
    rangestmp = deepcopy(ranges)
 
    for name in pnames:
        if name != pname:
            if multflatprior:  # multiply by flat prior range
                lnliketmp = marginalise(
                    lnliketmp
                    + np.log(1.0 / (rangestmp[name][-1] - rangestmp[name][0])),
                    name,
                    pnamestmp,
                    rangestmp,
                )
            else:
                lnliketmp = marginalise(lnliketmp, name, pnamestmp, rangestmp)
 
    return lnliketmp
 
 
# a function to chop the data time series, ts, up into contigous segments with a maximum length, chunkMax
def get_chunk_lengths(ts, chunkMax):
    i = j = count = 0
 
    length = len(ts)
 
    chunkLengths = []
 
    dt = np.min(np.diff(ts))  # the time steps in the data
 
    # create vector of data segment length
    while True:
        count += 1  # counter
 
        # break clause
        if i > length - 2:
            # set final value of chunkLength
            chunkLengths.append(count)
            break
 
        i += 1
 
        t1 = ts[i - 1]
        t2 = ts[i]
 
        # if consecutive points are within two sample times of each other count as in the same chunk
        if t2 - t1 > 2.0 * dt or count == chunkMax:
            chunkLengths.append(count)
            count = 0  # reset counter
 
            j += 1
 
    return chunkLengths
 
 
def pulsar_estimate_snr(source, det, tstart, duration, Sn, dt=1800):
    """
    A function to estimate the signal-to-noise ratio of a pulsar signal (for a triaxial neutron
    star emitting at twice the rotation frequency from the l=m=2 quadrupole mode) in a given
    detector over a given time range, for a given one-sided power spectral density.
 
    Inputs
    ------
        source - a dictionary containing 'h0', 'cosiota', 'psi', 'ra' (rads or string), 'dec' (rads or string)
           det - the detector, e.g. 'H1'
        tstart - a GPS start time
      duration - a signal duration (seconds)
            Sn - a one-sided power spectral density
            dt - timestep for antenna response calculation [default: 1800s]
 
    Returns
    -------
    rho        - an estimate of the optimal matched filter SNR
    """
 
    # if duration is greater than a sidereal day then just calculate the antenna pattern for 1 sidereal day plus the remainder
    sidday = 86164.0905  # one sidereal day in seconds
    Ndays = duration / sidday
    Nsamps = np.floor(sidday / dt)  # number of time samples in 1 sidereal day
    tend = tstart + duration
 
    sFp = 0
    sFc = 0
 
    # get antenna pattern and sum over the square of it
    if Ndays > 1:
        tts = np.arange(tstart, tstart + sidday, dt)
 
        # get antenna pattern over one sidereal day
        Fp, Fc = antenna_response(tts, source["ra"], source["dec"], source["psi"], det)
 
        sFp = np.floor(Ndays) * np.sum(Fp**2)
        sFc = np.floor(Ndays) * np.sum(Fc**2)
 
    # add on extra fractional part of a day
    if np.floor(Ndays) * sidday > dt:
        tts = np.arange(tstart + sidday * np.floor(Ndays), tend, dt)
 
        Fp, Fc = antenna_response(tts, source["ra"], source["dec"], source["psi"], det)
 
        sFp += np.sum(Fp**2)
        sFc += np.sum(Fc**2)
 
    # get snr squared
    snr2 = (
        (source["h0"] ** 2 / Sn)
        * dt
        * (
            sFp * (0.5 * (1 + source["cosiota"] ** 2)) ** 2
            + sFc * source["cosiota"] ** 2
        )
    )
 
    return np.sqrt(snr2)
 
 
def pulsar_estimate_h0_from_snr(snr, source, det, tstart, duration, Sn, dt=600):
    """
    A function to estimate the signal amplitude of a pulsar signal (for a triaxial neutron star emitting at twice
    the rotation frequency from the l=m=2 quadrupole mode) for a given SNR in a particular detector over
    a given time range, for a given one-sided power spectral density.
 
    Inputs
    ------
           snr - the optimal matched filter signal-to-noise ratio
        source - a dictionary containing 'cosiota', 'psi', 'ra' (rads or string), 'dec' (rads or string)
           det - the detector, e.g. 'H1'
        tstart - a GPS start time for the signal
      duration - a signal duration (seconds)
            Sn - a one-sided power spectral density
            dt - timestep for antenna response calculation [default: 600s]
 
    Returns
    -------
    h0         - an estimate of signal amplitude required to give the input SNR
    """
 
    # if duration is greater than a sidereal day then just calculate the antenna pattern for 1 sidereal day plus the remainder
    sidday = 86164.0905  # one sidereal day in seconds
    Ndays = duration / sidday
    Nsamps = np.floor(sidday / dt)  # number of time samples in 1 sidereal day
    tend = tstart + duration
 
    sFp = 0
    sFc = 0
 
    # get antenna pattern and sum over the square of it
    if Ndays > 1:
        tts = np.arange(tstart, tstart + sidday, dt)
 
        # get antenna pattern over one sidereal day
        Fp, Fc = antenna_response(tts, source["ra"], source["dec"], source["psi"], det)
 
        sFp = np.floor(Ndays) * np.sum(Fp**2)
        sFc = np.floor(Ndays) * np.sum(Fc**2)
 
    # add on extra fractional part of a day
    if np.floor(Ndays) * sidday > dt:
        tts = np.arange(tstart + sidday * np.floor(Ndays), tend, dt)
 
        Fp, Fc = antenna_response(tts, source["ra"], source["dec"], source["psi"], det)
 
        sFp += np.sum(Fp**2)
        sFc += np.sum(Fc**2)
 
    # get h0 squared
    h02 = (snr**2 * Sn) / (
        dt
        * (
            sFp * (0.5 * (1 + source["cosiota"] ** 2)) ** 2
            + sFc * source["cosiota"] ** 2
        )
    )
 
    return np.sqrt(h02)
 
 
def pulsar_posterior_grid(
    dets,
    ts,
    data,
    ra,
    dec,
    sigmas=None,
    paramranges={},
    datachunks=30,
    chunkmin=5,
    ngrid=50,
    outputlike=False,
):
    """
    A function to calculate the 4-parameter posterior probability density for a continuous wave signal
    given a set of processed data from a set of detectors.
 
    Inputs
    ------
           dets - a list of strings containing the detectors being used in the likelihood calculation
             ts - a dictionary of 1d numpy arrays containing the GPS times of the data for each detector
           data - a dictionary of 1d numpy arrays containing the complex data values for each detector
             ra - the right ascension (in rads) of the source
            dec - the declination (in rads) of the source
         sigmas - a dictionary of 1d numpy arrays containing the times series of standard deviation
                  values for each detector. If this is not given a Student's t likelihood will be used,
                  but if it is given a Gaussian likelihood will be used (default: None)
    paramranges - a dictionary of tuples for each parameter ('h0', 'phi0', 'psi' and 'cosiota') giving
                  the lower and upper ranges of the parameter grid and the number of grid points. If not
                  given then defaults will be used.
     datachunks - in the calculation split the data into stationary chunks with a maximum length given
                  by this value. If set to 0 or inf then the data will not be split. (default: 30)
       chunkmin - this is the shortest chunk length allowed to be included in the likelihood calculation
                  (default: 5)
          ngrid - the number of grid points to use for each dimension of the likelihood calculation. This
                  is used if the values are not specified in the paramranges argument (default: 50)
     outputlike - output the log likelihoods rather than posteriors (default: False)
 
    Returns
    -------
    L           - The 4d posterior (or likelihood) over all parameters
    h0pdf       - The 1d marginal posterior for h0
    phi0pdf     - The 1d marginal posterior for phi0 (the rotation frequency, not GW frequency)
    psipdf      - The 1d marginal posterior for psi
    cosiotapdf  - The 1d marginal posterior for cosiota
    lingrids    - A dictionary of the grid points for each parameter
    sigev       - The log evidence for the signal model
    noiseev     - The log evidence for the noise model
 
    An example would be:
    # set the detectors
    dets = ['H1', 'L1']
 
    # set the time series and data
    ts = {}
    data = {}
    for det in dets:
      ts[det] = np.arange(900000000., 921000843., 60.)
      data[det] = np.random.randn(len(ts[det])) + 1j*np.random.randn(len(ts[det]))
 
    # set the parameter ranges
    ra = 0.2
    dec = -0.8
    paramranges = {}
    paramranges['h0'] = (0., 2., 50)
    paramranges['psi'] = (0., np.pi/2., 50)
    paramranges['phi0'] = (0., np.pi, 50)
    paramranges['cosiota'] = (-1., 1., 50)
 
    L, h0pdf, phi0pdf, psipdf, cosiotapdf, grid, sigev, noiseev = pulsar_posterior_grid(dets, ts, data, ra, dec,
                                                                                        paramranges=paramranges)
    """
 
    # import numpy
    import numpy as np
    import sys
    from scipy.special import gammaln
 
    # set the likelihood to either Student's or Gaussian
    if sigmas == None:
        liketype = "studentst"
    else:
        liketype = "gaussian"
 
    # if dets is just a single string
    if not isinstance(dets, list):
        if isinstance(dets, str):
            dets = [dets]  # make into list
        else:
            print("Detector not, or incorrectly, set", file=sys.stderr)
            return
 
    # allowed list of detectors
    alloweddets = ["H1", "L1", "V1", "G1", "T1", "H2"]
 
    # check consistency of arguments
    for det in dets:
        if det not in alloweddets:
            print(
                "Detector not in list of allowed detectors ("
                + ",".join(alloweddets)
                + ")",
                file=sys.stderr,
            )
            return
 
        # checks on time stamps
        if det not in ts:
            print("No time stamps given for detector %s" % det, file=sys.stderr)
            return
 
        # checks on data
        if det not in data:
            print("No data time series given for detector %s" % det, file=sys.stderr)
            return
 
        # checks on sigmas
        if sigmas != None:
            if det not in sigmas:
                print(
                    "No sigma time series given for detector %s" % det, file=sys.stderr
                )
                return
 
        # check length consistency
        if len(ts[det]) != len(data[det]):
            print(
                "Length of times stamps array and data array are inconsistent for %s"
                % det,
                file=sys.stderr,
            )
 
        if sigmas != None:
            if len(ts[det]) != len(sigmas[det]):
                print(
                    "Length of times stamps array and sigma array are inconsistent for %s"
                    % det,
                    file=sys.stderr,
                )
 
    # setup grid on parameter space
    params = ["h0", "phi0", "psi", "cosiota"]
    defaultranges = {}  # default ranges for use if not specified
    defaultranges["phi0"] = (0.0, np.pi, ngrid)  # 0 to pi for phi0
    defaultranges["psi"] = (0.0, np.pi / 2.0, ngrid)  # 0 to pi/2 for psi
    defaultranges["cosiota"] = (-1.0, 1.0, ngrid)  # -1 to 1 for cos(iota)
    # 0 to 2 times the max standard deviation of the data (scaled to the data length)
    defaultranges["h0"] = (
        0.0,
        2.0
        * np.max([np.std(data[dv]) for dv in data])
        * np.sqrt(1440.0 / np.max([len(ts[dtx]) for dtx in ts])),
        ngrid,
    )
    lingrids = {}
    shapetuple = ()
    for param in params:
        if param in paramranges:
            if len(paramranges[param]) != 3:
                paramranges[param] = defaultranges[param]  # set to default range
            else:
                if (
                    paramranges[param][1] < paramranges[param][0]
                    or paramranges[param][2] < 1
                ):
                    print(
                        "Parameter ranges wrong for %s, reverting to defaults" % param,
                        file=sys.stderr,
                    )
                    paramranges[param] = defaultranges[param]
        else:  # use defaults
            paramranges[param] = defaultranges[param]
 
        lingrids[param] = np.linspace(
            paramranges[param][0], paramranges[param][1], paramranges[param][2]
        )
        shapetuple += (paramranges[param][2],)
 
    # set up meshgrid on parameter space
    [H0S, PHI0S, PSIS, COSIS] = np.meshgrid(
        lingrids["h0"],
        lingrids["phi0"],
        lingrids["psi"],
        lingrids["cosiota"],
        indexing="ij",
    )
 
    APLUS = 0.5 * (1.0 + COSIS**2)
    ACROSS = COSIS
    SINPHI = np.sin(2.0 * PHI0S)  # multiply phi0 by two to get GW frequency
    COSPHI = np.cos(2.0 * PHI0S)
    SIN2PSI = np.sin(2.0 * PSIS)
    COS2PSI = np.cos(2.0 * PSIS)
 
    Apsinphi_2 = 0.5 * APLUS * SINPHI
    Acsinphi_2 = 0.5 * ACROSS * SINPHI
    Apcosphi_2 = 0.5 * APLUS * COSPHI
    Accosphi_2 = 0.5 * ACROSS * COSPHI
    Acpsphicphi_2 = 0.5 * APLUS * ACROSS * COSPHI * SINPHI
 
    AP2 = APLUS**2
    AC2 = ACROSS**2
    TWOAPACSPCP = 2.0 * APLUS * ACROSS * SINPHI * COSPHI
    SP2 = SINPHI**2
    CP2 = COSPHI**2
 
    C2PS2P = COS2PSI * SIN2PSI
    C2P2 = COS2PSI**2
    S2P2 = SIN2PSI**2
 
    H02 = H0S**2
 
    # initialise loglikelihood
    like = np.zeros(shapetuple)
    noiselike = 0.0
 
    # get flat prior
    logprior = 0.0
    for p in params:
        logprior -= np.log(paramranges[p][1] - paramranges[p][0])
 
    # loop over detectors and calculate the log likelihood
    for det in dets:
        if liketype == "gaussian":
            nstd = sigmas[det]
        else:
            nstd = np.ones(len(ts[det]))
 
        # get the antenna pattern
        [As, Bs] = antenna_response(ts[det], ra, dec, 0.0, det)
 
        # split the data into chunks if required
        if np.isfinite(datachunks) and datachunks > 0:
            chunklengths = get_chunk_lengths(ts[det], datachunks)
        else:
            chunklengths = [len(ts[det])]
 
        startidx = 0
        for cl in chunklengths:
            endidx = startidx + cl
 
            # check for shortest allowed chunks
            if cl < chunkmin:
                continue
 
            thisdata = data[det][startidx:endidx] / nstd[startidx:endidx]
            ast = As[startidx:endidx] / nstd[startidx:endidx]
            bst = Bs[startidx:endidx] / nstd[startidx:endidx]
 
            A = np.sum(ast**2)
            B = np.sum(bst**2)
            AB = np.dot(ast, bst)
 
            dA1real = np.dot(thisdata.real, ast)
            dA1imag = np.dot(thisdata.imag, ast)
 
            dB1real = np.dot(thisdata.real, bst)
            dB1imag = np.dot(thisdata.imag, bst)
 
            dd1real = np.sum(thisdata.real**2)
            dd1imag = np.sum(thisdata.imag**2)
 
            P1 = A * C2P2 + B * S2P2 + 2.0 * AB * C2PS2P
            P2 = B * C2P2 + A * S2P2 - 2.0 * AB * C2PS2P
            P3 = AB * (C2P2 - S2P2) - A * C2PS2P + B * C2PS2P
 
            hr2 = (Apcosphi_2**2) * P1 + (Acsinphi_2**2) * P2 + Acpsphicphi_2 * P3
            hi2 = (Apsinphi_2**2) * P1 + (Accosphi_2**2) * P2 - Acpsphicphi_2 * P3
 
            d1hr = Apcosphi_2 * (dA1real * COS2PSI + dB1real * SIN2PSI) + Acsinphi_2 * (
                dB1real * COS2PSI - dA1real * SIN2PSI
            )
            d1hi = Apsinphi_2 * (dA1imag * COS2PSI + dB1imag * SIN2PSI) - Accosphi_2 * (
                dB1imag * COS2PSI - dA1imag * SIN2PSI
            )
 
            chiSq = dd1real + dd1imag + H02 * (hr2 + hi2) - 2.0 * H0S * (d1hr + d1hi)
 
            if liketype == "gaussian":
                like -= 0.5 * (chiSq)
                like -= (cl * np.log(2.0 * np.pi)) + 2.0 * np.sum(
                    np.log(nstd[startidx:endidx])
                )
                noiselike -= 0.5 * (dd1real + dd1imag)
                noiselike -= (cl * np.log(2.0 * np.pi)) + 2.0 * np.sum(
                    np.log(nstd[startidx:endidx])
                )
            else:
                like -= float(cl) * np.log(chiSq)
                like += gammaln(cl) - np.log(2.0) - cl * np.log(np.pi)
                noiselike -= float(cl) * np.log(dd1real + dd1imag)
                noiselike += gammaln(cl) - np.log(2.0) - cl * np.log(np.pi)
 
            startidx += cl  # updated start index
 
    # convert to posterior
    like += logprior
 
    # get signal evidence
    sigev = marginal(like, "all", params, lingrids)
 
    # normalise posterior
    if not outputlike:
        like -= sigev
    else:
        like -= logprior  # remove prior if outputting likelihood
 
    # get marginal posteriors for each parameter
    posts = {}
    for p in params:
        if not outputlike:
            posts[p] = np.exp(marginal(like, p, params, lingrids))
        else:
            posts[p] = marginal(
                like, p, params, lingrids, multflatprior=True
            )  # output individual log likelihoods
 
    # odds ratio for signal versus noise
    evrat = sigev - noiselike
 
    return (
        like,
        posts["h0"],
        posts["phi0"],
        posts["psi"],
        posts["cosiota"],
        lingrids,
        sigev,
        noiselike,
    )
 
 
# current version of the ATNF pulsar catalogue
ATNF_VERSION = "1.58"
 
 
def get_atnf_info(psr):
    """
    Get the pulsar (psr) distance (DIST in kpc), proper motion corrected period derivative (P1_I) and any association
    (ASSOC e.g. GC) from the ATNF catalogue.
    """
 
    import requests
 
    psrname = re.sub(r"\+", "%2B", psr)  # switch '+' for unicode character
 
    atnfurl = (
        "http://www.atnf.csiro.au/people/pulsar/psrcat/proc_form.php?version="
        + ATNF_VERSION
    )
    atnfurl += "&Dist=Dist&Assoc=Assoc&P1_i=P1_i"  # set parameters to get
    atnfurl += (
        "&startUserDefined=true&c1_val=&c2_val=&c3_val=&c4_val=&sort_attr=jname&sort_order=asc&condition=&pulsar_names="
        + psrname
    )
    atnfurl += "&ephemeris=selected&submit_ephemeris=Get+Ephemeris&coords_unit=raj%2Fdecj&radius=&coords_1=&coords_2="
    atnfurl += "&style=Long+with+last+digit+error&no_value=*&fsize=3&x_axis=&x_scale=linear&y_axis=&y_scale=linear&state=query"
 
    try:
        urldat = requests.get(atnfurl)  # read ATNF url
        predat = re.search(
            r"<pre[^>]*>([^<]+)</pre>", str(urldat.content)
        )  # to extract data within pre environment (without using BeautifulSoup) see e.g. http://stackoverflow.com/a/3369000/1862861 and http://stackoverflow.com/a/20046030/1862861
        pdat = (
            predat.group(1).strip().split(r"\n")
        )  # remove preceeding and trailing new lines and split lines
    except:
        print(
            "Warning... could not get information from ATNF pulsar catalogue.",
            file=sys.stderr,
        )
        return None
 
    # check whether information could be found and get distance, age and association from data
    dist = None
    p1_I = None
    assoc = None
    for line in [p for p in pdat if len(p) != 0]:
        if "WARNING" in line or "not in catalogue" in line:
            return None
        vals = line.split()
        if "DIST" in vals[0]:
            dist = float(vals[1])
        if "P1_I" in vals[0]:
            age = float(vals[1])
        if "ASSOC" in vals[0]:
            assoc = vals[1]
 
    return (dist, p1_I, assoc, atnfurl)
 
 
def cov_to_cor(cov):
    """
    Convert a covariance matrix to a correlation coefficient matrix.
    Return the correlation coefficient matrix and the standard deviations
    from the covariance matrix.
    """
 
    # get the standard deviations from the covariance matrix
    sigmas = np.sqrt(np.diag(cov))
 
    # convert these into a diagonal matrix
    D = sigmas * np.identity(cov.shape[0])
 
    # invert D
    Dinv = np.linalg.inv(D)
 
    # get correlation coefficient matrix
    Ccor = np.dot(np.dot(Dinv, cov), Dinv)
 
    return Ccor, sigmas