Computes integrated limb darkening profiles for ck2004 atmospheres.
        These are used for intensity-to-flux transformations. The evaluated
        integral is:

        ldint = 2 \pi \int_0^1 Imu mu dmu
339008
def compute_ck2004_ldints(self):
        

        if 'ck2004_all' not in self.content:
            print('Castelli & Kurucz (2004) intensities are not computed yet. Please compute those first.')
            return None

        ldaxes = self._ck2004_intensity_axes
        ldtable = self._ck2004_Imu_energy_grid
        pldtable = self._ck2004_Imu_photon_grid

        self._ck2004_ldint_energy_grid = np.nan*np.ones((len(ldaxes[0]), len(ldaxes[1]), len(ldaxes[2]), 1))
        self._ck2004_ldint_photon_grid = np.nan*np.ones((len(ldaxes[0]), len(ldaxes[1]), len(ldaxes[2]), 1))

        mu = ldaxes[3]
        Imu = 10**ldtable[:,:,:,:]/10**ldtable[:,:,:,-1:]
        pImu = 10**pldtable[:,:,:,:]/10**pldtable[:,:,:,-1:]

        # To compute the fluxes, we need to evaluate \int_0^1 2pi Imu mu dmu.

        for a in range(len(ldaxes[0])):
            for b in range(len(ldaxes[1])):
                for c in range(len(ldaxes[2])):

                    ldint = 0.0
                    pldint = 0.0
                    for i in range(len(mu)-1):
                        ki = (Imu[a,b,c,i+1]-Imu[a,b,c,i])/(mu[i+1]-mu[i])
                        ni = Imu[a,b,c,i]-ki*mu[i]
                        ldint += ki/3*(mu[i+1]**3-mu[i]**3) + ni/2*(mu[i+1]**2-mu[i]**2)

                        pki = (pImu[a,b,c,i+1]-pImu[a,b,c,i])/(mu[i+1]-mu[i])
                        pni = pImu[a,b,c,i]-pki*mu[i]
                        pldint += pki/3*(mu[i+1]**3-mu[i]**3) + pni/2*(mu[i+1]**2-mu[i]**2)

                    self._ck2004_ldint_energy_grid[a,b,c] = 2*ldint
                    self._ck2004_ldint_photon_grid[a,b,c] = 2*pldint

        self.content.append('ck2004_ldint')
1037164
def _calc_specpp_hiE(self, Egamma):
        
        # Fixed quad with n=40 is about 15 times faster and is always within
        # 0.5% of the result of adaptive quad for Egamma>0.1
        # WARNING: It also produces artifacts for steep distributions (e.g.
        # Maxwellian) at ~500 GeV. Reverting to adaptative quadrature
        # from scipy.integrate import fixed_quad
        # result=c*fixed_quad(self._photon_integrand, 0., 1., args = [Egamma,
        # ], n = 40)[0]
        from scipy.integrate import quad

        Egamma = Egamma.to("TeV").value
        specpp = (
            c.cgs.value
            * quad(
                self._photon_integrand,
                0.0,
                1.0,
                args=Egamma,
                epsrel=1e-3,
                epsabs=0,
            )[0]
        )

        return specpp * u.Unit("1/(s TeV)")
71867
def _compute_intensity(ccube, bexpcube):
        
        bexp_data = np.sqrt(bexpcube.data[0:-1, 0:] * bexpcube.data[1:, 0:])
        intensity_data = ccube.data / bexp_data
        intensity_map = HpxMap(intensity_data, ccube.hpx)
        return intensity_map
179920
def read_spectral_data(hdu):
    
    ebins = read_energy_bounds(hdu)
    fluxes = np.ndarray((len(ebins)))
    try:
        fluxes[0:-1] = hdu.data.field("E_MIN_FL")
        fluxes[-1] = hdu.data.field("E_MAX_FL")[-1]
        npreds = hdu.data.field("NPRED")
    except:
        fluxes = np.ones((len(ebins)))
        npreds = np.ones((len(ebins)))
    return ebins, fluxes, npreds
103646
def set_pd_mag_from_flux_density(photodict,
                                 fd='',
                                 efd='',
                                 lefd='',
                                 uefd='',
                                 sig=DEFAULT_UL_SIGMA):
    
    with localcontext() as ctx:
        if lefd == '' or uefd == '':
            lefd = efd
            uefd = efd
        prec = max(
            get_sig_digits(str(fd), strip_zeroes=False),
            get_sig_digits(str(lefd), strip_zeroes=False),
            get_sig_digits(str(uefd), strip_zeroes=False)) + 1
        ctx.prec = prec
        dlefd = Decimal(str(lefd))
        duefd = Decimal(str(uefd))
        if fd != '':
            dfd = Decimal(str(fd))
        dsig = Decimal(str(sig))
        if fd == '' or float(fd) < DEFAULT_UL_SIGMA * float(uefd):
            photodict[PHOTOMETRY.UPPER_LIMIT] = True
            photodict[PHOTOMETRY.UPPER_LIMIT_SIGMA] = str(sig)
            photodict[PHOTOMETRY.MAGNITUDE] = str(Decimal('23.9') - D25 * (
                dsig * duefd).log10())
            if fd:
                photodict[PHOTOMETRY.E_UPPER_MAGNITUDE] = str(D25 * (
                    (dfd + duefd).log10() - dfd.log10()))
        else:
            photodict[PHOTOMETRY.MAGNITUDE] = str(Decimal('23.9') - D25 *
                                                  dfd.log10())
            photodict[PHOTOMETRY.E_UPPER_MAGNITUDE] = str(D25 * (
                (dfd + duefd).log10() - dfd.log10()))
            photodict[PHOTOMETRY.E_LOWER_MAGNITUDE] = str(D25 * (
                dfd.log10() - (dfd - dlefd).log10()))
1030025
def do_diagnostics(self):
        ''''''
        self.OLR = self.subprocess['LW'].flux_to_space
        self.LW_down_sfc = self.subprocess['LW'].flux_to_sfc
        self.LW_up_sfc = self.subprocess['LW'].flux_from_sfc
        self.LW_absorbed_sfc = self.LW_down_sfc - self.LW_up_sfc
        self.LW_absorbed_atm = self.subprocess['LW'].absorbed
        self.LW_emission = self.subprocess['LW'].emission
            #  contributions to OLR from surface and atm. levels
            #self.diagnostics['OLR_sfc'] = self.flux['sfc2space']
            #self.diagnostics['OLR_atm'] = self.flux['atm2space']
        self.ASR = (self.subprocess['SW'].flux_from_space -
                    self.subprocess['SW'].flux_to_space)
        #self.SW_absorbed_sfc = (self.subprocess['surface'].SW_from_atm -
        #                            self.subprocess['surface'].SW_to_atm)
        self.SW_absorbed_atm = self.subprocess['SW'].absorbed
        self.SW_down_sfc = self.subprocess['SW'].flux_to_sfc
        self.SW_up_sfc = self.subprocess['SW'].flux_from_sfc
        self.SW_absorbed_sfc = self.SW_down_sfc - self.SW_up_sfc
        self.SW_up_TOA = self.subprocess['SW'].flux_to_space
        self.SW_down_TOA = self.subprocess['SW'].flux_from_space
        self.planetary_albedo = (self.subprocess['SW'].flux_to_space /
                                 self.subprocess['SW'].flux_from_space)
921375
def Fierz_to_EOS_lep(C, ddll, parameters):
    """"""
    p = parameters
    V = ckmutil.ckm.ckm_tree(p["Vus"], p["Vub"], p["Vcb"], p["delta"])
    Vtb = V[2,2]
    Vts = V[2,1]
    ind = ddll.replace('l_','').replace('nu_','')
    ind2 = ddll.replace('l_','').replace('nu_','')[2::]
    e = sqrt(4* pi * parameters['alpha_e'])
    dic = {
        'b->s' + ind2 + '::c9' : (16 * pi**2) / e**2 * C['F' + ind + '9'],
        'b->s' + ind2 + "::c9'" : (16 * pi**2) / e**2 * C['F' + ind + '9p'],
        'b->s' + ind2 + "::c10" : (16 * pi**2) / e**2 * C['F' + ind + '10'],
        'b->s' + ind2 + "::c10'" : (16 * pi**2) / e**2 * C['F' + ind + '10p'],
        'b->s' + ind2 + "::cS" : (16 * pi**2) / e**2  * C['F' + ind + 'S'],
        'b->s' + ind2 + "::cS'" : (16 * pi**2) / e**2 * C['F' + ind + 'Sp'],
        'b->s' + ind2 + "::cP" : (16 * pi**2) / e**2  * C['F' + ind + 'P'],
        'b->s' + ind2 + "::cP'" : (16 * pi**2) / e**2  * C['F' + ind + 'Pp'],
        'b->s' + ind2 + "::cT"  : (16 * pi**2) / e**2  * C['F' + ind + 'T'],
        'b->s' + ind2 + "::cT5"  : (16 * pi**2) / e**2  * C['F' + ind + 'T5']
        }
    prefactor = sqrt(2)/p['GF']/Vtb/Vts.conj()/4
    return {k: prefactor * v for k,v in dic.items()}
989047
def expected_magx(RAW_IMU, ATTITUDE, inclination, declination):
    ''''''
    v = expected_mag(RAW_IMU, ATTITUDE, inclination, declination)
    return v.x
501714
def expected_magy(RAW_IMU, ATTITUDE, inclination, declination):
    ''''''
    v = expected_mag(RAW_IMU, ATTITUDE, inclination, declination)
    return v.y
719271
def _compute_heating_rates(self):
        ''''''
        (ncol, nlay, icld, permuteseed, irng, idrv, cp,
                play, plev, tlay, tlev, tsfc,
                h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
                cfc11vmr, cfc12vmr, cfc22vmr, ccl4vmr, emis,
                inflglw, iceflglw, liqflglw,
                cldfrac, ciwp, clwp, reic, relq, tauc, tauaer,) = self._prepare_lw_arguments()
        if icld == 0:  # clear-sky only
            cldfmcl = np.zeros((ngptlw,ncol,nlay))
            ciwpmcl = np.zeros((ngptlw,ncol,nlay))
            clwpmcl = np.zeros((ngptlw,ncol,nlay))
            reicmcl = np.zeros((ncol,nlay))
            relqmcl = np.zeros((ncol,nlay))
            taucmcl = np.zeros((ngptlw,ncol,nlay))
        else:
            #  Call the Monte Carlo Independent Column Approximation (McICA, Pincus et al., JC, 2003)
            (cldfmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl, taucmcl) = \
                _rrtmg_lw.climlab_mcica_subcol_lw(
                            ncol, nlay, icld,
                            permuteseed, irng, play,
                            cldfrac, ciwp, clwp, reic, relq, tauc)
            #  Call the RRTMG_LW driver to compute radiative fluxes
        (uflx, dflx, hr, uflxc, dflxc, hrc, duflx_dt, duflxc_dt) = \
            _rrtmg_lw.climlab_rrtmg_lw(ncol, nlay, icld, idrv,
                 play, plev, tlay, tlev, tsfc,
                 h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
                 cfc11vmr, cfc12vmr, cfc22vmr, ccl4vmr, emis,
                 inflglw, iceflglw, liqflglw, cldfmcl,
                 taucmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl,
                 tauaer)
        #  Output is all (ncol,nlay+1) or (ncol,nlay)
        self.LW_flux_up = _rrtm_to_climlab(uflx) + 0.*self.LW_flux_up
        self.LW_flux_down = _rrtm_to_climlab(dflx) + 0.*self.LW_flux_down
        self.LW_flux_up_clr = _rrtm_to_climlab(uflxc) + 0.*self.LW_flux_up_clr
        self.LW_flux_down_clr = _rrtm_to_climlab(dflxc) + 0.*self.LW_flux_down_clr
        #  Compute quantities derived from fluxes, including OLR
        self._compute_LW_flux_diagnostics()
        #  calculate heating rates from flux divergence
        LWheating_Wm2 = np.array(np.diff(self.LW_flux_net, axis=-1)) + 0.*self.Tatm
        LWheating_clr_Wm2 = np.array(np.diff(self.LW_flux_net_clr, axis=-1)) + 0.*self.Tatm
        self.heating_rate['Ts'] = np.array(-self.LW_flux_net[..., -1, np.newaxis]) + 0.*self.Ts
        self.heating_rate['Tatm'] = LWheating_Wm2
        #  Convert to K / day
        Catm = self.Tatm.domain.heat_capacity
        self.TdotLW = LWheating_Wm2 / Catm * const.seconds_per_day
        self.TdotLW_clr = LWheating_clr_Wm2 / Catm * const.seconds_per_day
216573
def f_aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
    
    return pvlib.irradiance.aoi(surface_tilt, surface_azimuth,
                                solar_zenith, solar_azimuth)
596623
def eclipseprob(self):
        
        #TODO: incorporate eccentricity/omega for exact calculation?
        s = self.stars
        return ((s['radius_1'] + s['radius_2'])*RSUN /
                (semimajor(s['P'],s['mass_1'] + s['mass_2'])*AU))
1103362
def get_magicc7_to_openscm_variable_mapping(inverse=False):
    

    def get_openscm_replacement(in_var):
        if in_var.endswith("_INVERSE_EMIS"):
            prefix = "Inverse Emissions"
        elif in_var.endswith("_EMIS"):
            prefix = "Emissions"
        elif in_var.endswith("_CONC"):
            prefix = "Atmospheric Concentrations"
        elif in_var.endswith("_RF"):
            prefix = "Radiative Forcing"
        elif in_var.endswith("_OT"):
            prefix = "Optical Thickness"
        else:
            raise ValueError("This shouldn't happen")

        variable = in_var.split("_")[0]
        # I hate edge cases
        if variable.endswith("EQ"):
            variable = variable.replace("EQ", " Equivalent")

        if "GHG" in variable:
            variable = variable.replace("GHG", "Greenhouse Gases")

        if "BIOMASSAER" in variable:
            variable = variable.replace("BIOMASSAER", "Aerosols|MAGICC AFOLU")

        if "CO2CH4N2O" in variable:
            variable = variable.replace("CO2CH4N2O", "CO2, CH4 and N2O")

        aggregate_indicators = {
            "KYOTO": "Kyoto Gases",
            "FGASSUM": "F Gases",
            "MHALOSUM": "Montreal Protocol Halogen Gases",
        }
        for agg_indicator, long_name in aggregate_indicators.items():
            if variable.startswith(agg_indicator):
                stripped_var = variable.replace(agg_indicator, "")
                if stripped_var:
                    variable = DATA_HIERARCHY_SEPARATOR.join([stripped_var, long_name])
                else:
                    variable = long_name

        edge_case_B = variable.upper() in ("HCFC141B", "HCFC142B")
        if variable.endswith("I"):
            variable = DATA_HIERARCHY_SEPARATOR.join(
                [variable[:-1], "MAGICC Fossil and Industrial"]
            )
        elif variable.endswith("B") and not edge_case_B:
            variable = DATA_HIERARCHY_SEPARATOR.join([variable[:-1], "MAGICC AFOLU"])

        case_adjustments = {
            "SOX": "SOx",
            "NOX": "NOx",
            "HFC134A": "HFC134a",
            "HFC143A": "HFC143a",
            "HFC152A": "HFC152a",
            "HFC227EA": "HFC227ea",
            "HFC236FA": "HFC236fa",
            "HFC245FA": "HFC245fa",
            "HFC365MFC": "HFC365mfc",
            "HCFC141B": "HCFC141b",
            "HCFC142B": "HCFC142b",
            "CH3CCL3": "CH3CCl3",
            "CCL4": "CCl4",
            "CH3CL": "CH3Cl",
            "CH2CL2": "CH2Cl2",
            "CHCL3": "CHCl3",
            "CH3BR": "CH3Br",
            "HALON1211": "Halon1211",
            "HALON1301": "Halon1301",
            "HALON2402": "Halon2402",
            "HALON1202": "Halon1202",
            "SOLAR": "Solar",
            "VOLCANIC": "Volcanic",
            "EXTRA": "Extra",
        }
        variable = apply_string_substitutions(variable, case_adjustments)

        return DATA_HIERARCHY_SEPARATOR.join([prefix, variable])

    magicc7_suffixes = ["_EMIS", "_CONC", "_RF", "_OT", "_INVERSE_EMIS"]
    magicc7_base_vars = MAGICC7_EMISSIONS_UNITS.magicc_variable.tolist() + [
        "SOLAR",
        "VOLCANIC",
        "CO2EQ",
        "KYOTOCO2EQ",
        "FGASSUMHFC134AEQ",
        "MHALOSUMCFC12EQ",
        "GHG",
        "KYOTOGHG",
        "FGASSUM",
        "MHALOSUM",
        "BIOMASSAER",
        "CO2CH4N2O",
        "EXTRA",
    ]
    magicc7_vars = [
        base_var + suffix
        for base_var in magicc7_base_vars
        for suffix in magicc7_suffixes
    ]

    replacements = {m7v: get_openscm_replacement(m7v) for m7v in magicc7_vars}

    replacements.update(
        {
            "SURFACE_TEMP": "Surface Temperature",
            "TOTAL_INCLVOLCANIC_RF": "Radiative Forcing",
            "VOLCANIC_ANNUAL_RF": "Radiative Forcing|Volcanic",
            "TOTAL_ANTHRO_RF": "Radiative Forcing|Anthropogenic",
            "TOTAER_DIR_RF": "Radiative Forcing|Aerosols|Direct Effect",
            "CLOUD_TOT_RF": "Radiative Forcing|Aerosols|Indirect Effect",
            "MINERALDUST_RF": "Radiative Forcing|Mineral Dust",
            "STRATOZ_RF": "Radiative Forcing|Stratospheric Ozone",
            "TROPOZ_RF": "Radiative Forcing|Tropospheric Ozone",
            "CH4OXSTRATH2O_RF": "Radiative Forcing|CH4 Oxidation Stratospheric H2O",  # what is this
            "LANDUSE_RF": "Radiative Forcing|Land-use Change",
            "BCSNOW_RF": "Radiative Forcing|Black Carbon on Snow",
            "CO2PF_EMIS": "Land to Air Flux|CO2|MAGICC Permafrost",
            # "CH4PF_EMIS": "Land to Air Flux|CH4|MAGICC Permafrost",  # TODO: test and then add when needed
        }
    )

    agg_ocean_heat_top = "Aggregated Ocean Heat Content"
    heat_content_aggreg_depths = {
        "HEATCONTENT_AGGREG_DEPTH{}".format(i): "{}{}Depth {}".format(
            agg_ocean_heat_top, DATA_HIERARCHY_SEPARATOR, i
        )
        for i in range(1, 4)
    }
    replacements.update(heat_content_aggreg_depths)
    replacements.update({"HEATCONTENT_AGGREG_TOTAL": agg_ocean_heat_top})

    ocean_temp_layer = {
        "OCEAN_TEMP_LAYER_{0:03d}".format(i): "Ocean Temperature{}Layer {}".format(
            DATA_HIERARCHY_SEPARATOR, i
        )
        for i in range(1, 999)
    }
    replacements.update(ocean_temp_layer)

    if inverse:
        return {v: k for k, v in replacements.items()}
    else:
        return replacements
406880
def calcSMAfromT(self, epsilon=0.7):
        

        return eq.MeanPlanetTemp(self.albedo(), self.star.T, self.star.R, epsilon, self.T).a
1184815
def imrphenomc_tmplt(**kwds):
    
    # Pull out the input arguments
    f_min = float128(kwds['f_lower'])
    f_max = float128(kwds['f_final'])
    delta_f = float128(kwds['delta_f'])
    distance = float128(kwds['distance'])
    mass1 = float128(kwds['mass1'])
    mass2 = float128(kwds['mass2'])
    spin1z = float128(kwds['spin1z'])
    spin2z = float128(kwds['spin2z'])

    if 'out' in kwds:
        out = kwds['out']
    else:
        out = None

    # Calculate binary parameters
    M = mass1 + mass2
    eta = mass1 * mass2 / (M * M)
    Xi = (mass1 * spin1z / M) + (mass2 * spin2z / M)
    Xisum = 2.*Xi
    Xiprod = Xi*Xi
    Xi2 = Xi*Xi

    m_sec = M * lal.MTSUN_SI;
    piM = lal.PI * m_sec;

    ## The units of distance given as input is taken to pe Mpc. Converting to SI
    distance *= (1.0e6 * lal.PC_SI / (2. * sqrt(5. / (64.*lal.PI)) * M * lal.MRSUN_SI * M * lal.MTSUN_SI))

    # Check if the value of f_max is correctly given, else replace with the fCut
    # used in the PhenomB code in lalsimulation. The various coefficients come
    # from Eq.(4.18) of http://arxiv.org/pdf/0710.2335 and
    # Table I of http://arxiv.org/pdf/0712.0343
    if not f_max:
        f_max = (1.7086 * eta * eta - 0.26592 * eta + 0.28236) / piM

    # Transform the eta, chi to Lambda parameters, using Eq 5.14, Table II of Main
    # paper.
    z101 = -2.417e-03
    z102 = -1.093e-03
    z111 = -1.917e-02
    z110 = 7.267e-02
    z120 = -2.504e-01

    z201 = 5.962e-01
    z202 = -5.600e-02
    z211 = 1.520e-01
    z210 = -2.970e+00
    z220 = 1.312e+01

    z301 = -3.283e+01
    z302 = 8.859e+00
    z311 = 2.931e+01
    z310 = 7.954e+01
    z320 = -4.349e+02

    z401 = 1.619e+02
    z402 = -4.702e+01
    z411 = -1.751e+02
    z410 = -3.225e+02
    z420 = 1.587e+03

    z501 = -6.320e+02
    z502 = 2.463e+02
    z511 = 1.048e+03
    z510 = 3.355e+02
    z520 = -5.115e+03

    z601 = -4.809e+01
    z602 = -3.643e+02
    z611 = -5.215e+02
    z610 = 1.870e+03
    z620 = 7.354e+02

    z701 = 4.149e+00
    z702 = -4.070e+00
    z711 = -8.752e+01
    z710 = -4.897e+01
    z720 = 6.665e+02

    z801 = -5.472e-02
    z802 = 2.094e-02
    z811 = 3.554e-01
    z810 = 1.151e-01
    z820 = 9.640e-01

    z901 = -1.235e+00
    z902 = 3.423e-01
    z911 = 6.062e+00
    z910 = 5.949e+00
    z920 = -1.069e+01

    eta2 = eta*eta
    Xi2 = Xiprod

    # Calculate alphas, gamma, deltas from Table II and Eq 5.14 of Main paper
    a1 = z101 * Xi + z102 * Xi2 + z111 * eta * Xi + z110 * eta + z120 * eta2
    a2 = z201 * Xi + z202 * Xi2 + z211 * eta * Xi + z210 * eta + z220 * eta2
    a3 = z301 * Xi + z302 * Xi2 + z311 * eta * Xi + z310 * eta + z320 * eta2
    a4 = z401 * Xi + z402 * Xi2 + z411 * eta * Xi + z410 * eta + z420 * eta2
    a5 = z501 * Xi + z502 * Xi2 + z511 * eta * Xi + z510 * eta + z520 * eta2
    a6 = z601 * Xi + z602 * Xi2 + z611 * eta * Xi + z610 * eta + z620 * eta2

    g1 = z701 * Xi + z702 * Xi2 + z711 * eta * Xi + z710 * eta + z720 * eta2

    del1 = z801 * Xi + z802 * Xi2 + z811 * eta * Xi + z810 * eta + z820 * eta2
    del2 = z901 * Xi + z902 * Xi2 + z911 * eta * Xi + z910 * eta + z920 * eta2

    # Get the spin of the final BH
    afin = FinalSpin( Xi, eta )
    Q = Qa( abs(afin) )

    # Get the fRD
    frd = fRD( abs(afin), M)
    Mfrd = frd * m_sec

    # Define the frequencies where SPA->PM->RD
    f1 = 0.1 * frd
    Mf1 = m_sec * f1
    f2 = frd
    Mf2 = m_sec * f2
    d1 = 0.005
    d2 = 0.005
    f0 = 0.98 * frd
    Mf0 = m_sec * f0
    d0 = 0.015

    # Now use this frequency for calculation of betas
    # calculate beta1 and beta2, that appear in Eq 5.7 in the main paper.
    b2 = ((-5./3.)* a1 * pow(Mfrd,(-8./3.)) - a2/(Mfrd*Mfrd) - \
      (a3/3.)*pow(Mfrd,(-4./3.)) + (2./3.)* a5 * pow(Mfrd,(-1./3.)) + a6)/eta

    psiPMrd = (a1 * pow(Mfrd,(-5./3.)) + a2/Mfrd + a3 * pow(Mfrd,(-1./3.)) + \
      a4 + a5 * pow(Mfrd,(2./3.)) + a6 * Mfrd)/eta
    b1 = psiPMrd - (b2 * Mfrd)

    ### Calculate the PN coefficients, Eq A3 - A5 of main paper ###
    pfaN = 3.0/(128.0 * eta)
    pfa2 = (3715./756.) + (55.*eta/9.0)
    pfa3 = -16.0*lal.PI + (113./3.)*Xi - 38.*eta*Xisum/3.
    pfa4 = (152.93365/5.08032) - 50.*Xi2 + eta*(271.45/5.04 + 1.25*Xiprod) + \
        3085.*eta2/72.
    pfa5 = lal.PI*(386.45/7.56 - 65.*eta/9.) - \
        Xi*(735.505/2.268 + 130.*eta/9.) + Xisum*(1285.0*eta/8.1 + 170.*eta2/9.) - \
        10.*Xi2*Xi/3. + 10.*eta*Xi*Xiprod
    pfa6 = 11583.231236531/4.694215680 - 640.0*lal.PI*lal.PI/3. - \
        6848.0*lal.GAMMA/21. - 684.8*log(64.)/6.3 + \
        eta*(2255.*lal.PI*lal.PI/12. - 15737.765635/3.048192) + \
        76.055*eta2/1.728 - (127.825*eta2*eta/1.296) + \
        2920.*lal.PI*Xi/3. - (175. - 1490.*eta)*Xi2/3. - \
        (1120.*lal.PI/3. - 1085.*Xi/3.)*eta*Xisum + \
        (269.45*eta/3.36 - 2365.*eta2/6.)*Xiprod

    pfa6log = -6848./63.

    pfa7 = lal.PI*(770.96675/2.54016 + 378.515*eta/1.512 - 740.45*eta2/7.56) - \
        Xi*(20373.952415/3.048192 + 1509.35*eta/2.24 - 5786.95*eta2/4.32) + \
        Xisum*(4862.041225*eta/1.524096 + 1189.775*eta2/1.008 - 717.05*eta2*eta/2.16 - 830.*eta*Xi2/3. + 35.*eta2*Xiprod/3.) - \
        560.*lal.PI*Xi2 + 20.*lal.PI*eta*Xiprod + \
        Xi2*Xi*(945.55/1.68 - 85.*eta) + Xi*Xiprod*(396.65*eta/1.68 + 255.*eta2)


    xdotaN = 64.*eta/5.
    xdota2 = -7.43/3.36 - 11.*eta/4.
    xdota3 = 4.*lal.PI - 11.3*Xi/1.2 + 19.*eta*Xisum/6.
    xdota4 = 3.4103/1.8144 + 5*Xi2 + eta*(13.661/2.016 - Xiprod/8.) + 5.9*eta2/1.8
    xdota5 = -lal.PI*(41.59/6.72 + 189.*eta/8.) - Xi*(31.571/1.008 - 116.5*eta/2.4) + \
          Xisum*(21.863*eta/1.008 - 79.*eta2/6.) - 3*Xi*Xi2/4. + \
          9.*eta*Xi*Xiprod/4.
    xdota6 = 164.47322263/1.39708800 - 17.12*lal.GAMMA/1.05 + \
          16.*lal.PI*lal.PI/3 - 8.56*log(16.)/1.05 + \
          eta*(45.1*lal.PI*lal.PI/4.8 - 561.98689/2.17728) + \
          5.41*eta2/8.96 - 5.605*eta*eta2/2.592 - 80.*lal.PI*Xi/3. + \
          eta*Xisum*(20.*lal.PI/3. - 113.5*Xi/3.6) + \
          Xi2*(64.153/1.008 - 45.7*eta/3.6) - \
          Xiprod*(7.87*eta/1.44 - 30.37*eta2/1.44)

    xdota6log = -856./105.

    xdota7 = -lal.PI*(4.415/4.032 - 358.675*eta/6.048 - 91.495*eta2/1.512) - \
          Xi*(252.9407/2.7216 - 845.827*eta/6.048 + 415.51*eta2/8.64) + \
          Xisum*(158.0239*eta/5.4432 - 451.597*eta2/6.048 + 20.45*eta2*eta/4.32 + 107.*eta*Xi2/6. - 5.*eta2*Xiprod/24.) + \
          12.*lal.PI*Xi2 - Xi2*Xi*(150.5/2.4 + eta/8.) + \
          Xi*Xiprod*(10.1*eta/2.4 + 3.*eta2/8.)


    AN = 8.*eta*sqrt(lal.PI/5.)
    A2 = (-107. + 55.*eta)/42.
    A3 = 2.*lal.PI - 4.*Xi/3. + 2.*eta*Xisum/3.
    A4 = -2.173/1.512 - eta*(10.69/2.16 - 2.*Xiprod) + 2.047*eta2/1.512
    A5 = -10.7*lal.PI/2.1 + eta*(3.4*lal.PI/2.1)

    A5imag = -24.*eta

    A6 = 270.27409/6.46800 - 8.56*lal.GAMMA/1.05 + \
      2.*lal.PI*lal.PI/3. + \
      eta*(4.1*lal.PI*lal.PI/9.6 - 27.8185/3.3264) - \
      20.261*eta2/2.772 + 11.4635*eta*eta2/9.9792 - \
      4.28*log(16.)/1.05

    A6log = -428./105.

    A6imag = 4.28*lal.PI/1.05

    ### Define other parameters needed by waveform generation ###
    kmin = int(f_min / delta_f)
    kmax = int(f_max / delta_f)
    n = kmax + 1;

    if not out:
        htilde = FrequencySeries(zeros(n,dtype=numpy.complex128), delta_f=delta_f, copy=False)
    else:
        if type(out) is not Array:
            raise TypeError("Output must be an instance of Array")
        if len(out) < kmax:
            raise TypeError("Output array is too small")
        if out.dtype != complex64:
            raise TypeError("Output array is the wrong dtype")
        htilde = FrequencySeries(out, delta_f=delta_f, copy=False)

    phenomC_kernel(htilde.data[kmin:kmax], kmin, delta_f, eta, Xi, distance,
                                       m_sec,  piM,  Mfrd,
                                       pfaN,  pfa2,  pfa3,  pfa4, pfa5,  pfa6,  pfa6log,  pfa7,
                                       a1,  a2,  a3,  a4, a5,  a6,  b1,  b2,
                                       Mf1,  Mf2,  Mf0, d1,  d2,  d0,
                                       xdota2,  xdota3,  xdota4, xdota5,  xdota6,  xdota6log,
                                       xdota7,  xdotaN,  AN, A2,  A3,  A4,  A5,
                                       A5imag,  A6,  A6log,  A6imag,
                                       g1,  del1,  del2,  Q )
    hp = htilde
    hc = htilde * 1j
    return hp, hc
1251204
def densm(alt, d0, xm, tz, mn3, zn3, tn3, tgn3, mn2, zn2, tn2, tgn2):
    '''

'''
    xs = [0.0]*10
    ys = [0.0]*10
    y2out = [0.0]*10
    rgas = 831.4
    #rgas = 831.44621    #maybe make this a global constant?
    densm_tmp=d0
    if (alt>zn2[0]): # pragma: no cover
        if(xm==0.0):
            return tz[0]
        else:
            return d0

    #/* STRATOSPHERE/MESOSPHERE TEMPERATURE */
    if (alt>zn2[mn2-1]):
        z=alt
    else:
        z=zn2[mn2-1];
    mn=mn2;
    z1=zn2[0];
    z2=zn2[mn-1];
    t1=tn2[0];
    t2=tn2[mn-1];
    zg = zeta(z, z1);
    zgdif = zeta(z2, z1);

    #/* set up spline nodes */
    for k in range(mn):
        xs[k]=zeta(zn2[k],z1)/zgdif;
        ys[k]=1.0 / tn2[k];
    yd1=-tgn2[0] / (t1*t1) * zgdif;
    yd2=-tgn2[1] / (t2*t2) * zgdif * (pow(((re[0]+z2)/(re[0]+z1)),2.0));

    #/* calculate spline coefficients */
    spline (xs, ys, mn, yd1, yd2, y2out);   #No need to change this
    x = zg/zgdif;
    y = [0.0]
    splint (xs, ys, y2out, mn, x, y);

    #/* temperature at altitude */
    tz[0] = 1.0 / y[0];
    if (xm!=0.0):
        #/* calaculate stratosphere / mesospehere density */
        glb = gsurf[0] / (pow((1.0 + z1/re[0]),2.0));
        gamm = xm * glb * zgdif / rgas;

        #/* Integrate temperature profile */
        yi = [0.0]
        splini(xs, ys, y2out, mn, x, yi);
        expl=gamm*yi[0];
        if (expl>50.0): # pragma: no cover
            expl=50.0

        #/* Density at altitude */
        densm_tmp = densm_tmp * (t1 / tz[0]) * exp(-expl);

    if (alt>zn3[0]):
        if (xm==0.0):
            return tz[0]
        else:
            return densm_tmp

    #/* troposhere / stratosphere temperature */
    z = alt;
    mn = mn3;
    z1=zn3[0];
    z2=zn3[mn-1];
    t1=tn3[0];
    t2=tn3[mn-1];
    zg=zeta(z,z1);
    zgdif=zeta(z2,z1);



    #/* set up spline nodes */
    for k in range(mn):
        xs[k] = zeta(zn3[k],z1) / zgdif;
        ys[k] = 1.0 / tn3[k];
    
    yd1=-tgn3[0] / (t1*t1) * zgdif;
    yd2=-tgn3[1] / (t2*t2) * zgdif * (pow(((re[0]+z2)/(re[0]+z1)),2.0));

    #/* calculate spline coefficients */
    spline (xs, ys, mn, yd1, yd2, y2out);
    x = zg/zgdif;
    y = [0.0]
    splint (xs, ys, y2out, mn, x, y);

    #/* temperature at altitude */
    tz[0] = 1.0 / y[0];
    if (xm!=0.0):
        #/* calaculate tropospheric / stratosphere density */
        glb = gsurf[0] / (pow((1.0 + z1/re[0]),2.0));
        gamm = xm * glb * zgdif / rgas;

        #/* Integrate temperature profile */
        yi = [0.0]
        splini(xs, ys, y2out, mn, x, yi);
        expl=gamm*yi[0];
        if (expl>50.0): # pragma: no cover
            expl=50.0;

        #/* Density at altitude */
        densm_tmp = densm_tmp * (t1 / tz[0]) * exp(-expl);
    
    if (xm==0.0):
        return tz[0];
    else:
        return densm_tmp;
1321970
def _dsinit(
       whichconst,
       cosim,  emsq,   argpo,   s1,     s2,
       s3,     s4,     s5,      sinim,  ss1,
       ss2,    ss3,    ss4,     ss5,    sz1,
       sz3,    sz11,   sz13,    sz21,   sz23,
       sz31,   sz33,   t,       tc,     gsto,
       mo,     mdot,   no,      nodeo,  nodedot,
       xpidot, z1,     z3,      z11,    z13,
       z21,    z23,    z31,     z33,    ecco,
       eccsq,  em,    argpm,  inclm, mm,
       nm,    nodem,
       irez,
       atime, d2201, d2211,  d3210, d3222,
       d4410, d4422, d5220,  d5232, d5421,
       d5433, dedt,  didt,   dmdt,
       dnodt, domdt, del1,   del2,  del3,
       xfact, xlamo, xli,    xni,
     ):

     q22    = 1.7891679e-6;
     q31    = 2.1460748e-6;
     q33    = 2.2123015e-7;
     root22 = 1.7891679e-6;
     root44 = 7.3636953e-9;
     root54 = 2.1765803e-9;
     rptim  = 4.37526908801129966e-3; # equates to 7.29211514668855e-5 rad/sec
     root32 = 3.7393792e-7;
     root52 = 1.1428639e-7;
     x2o3   = 2.0 / 3.0;
     znl    = 1.5835218e-4;
     zns    = 1.19459e-5;

     #  sgp4fix identify constants and allow alternate values
     xke = whichconst.xke

     #  -------------------- deep space initialization ------------
     irez = 0;
     if 0.0034906585 < nm < 0.0052359877:
         irez = 1;
     if 8.26e-3 <= nm <= 9.24e-3 and em >= 0.5:
         irez = 2;

     #  ------------------------ do solar terms -------------------
     ses  =  ss1 * zns * ss5;
     sis  =  ss2 * zns * (sz11 + sz13);
     sls  = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq);
     sghs =  ss4 * zns * (sz31 + sz33 - 6.0);
     shs  = -zns * ss2 * (sz21 + sz23);
     #  sgp4fix for 180 deg incl
     if inclm < 5.2359877e-2 or inclm > pi - 5.2359877e-2:
       shs = 0.0;
     if sinim != 0.0:
       shs = shs / sinim;
     sgs  = sghs - cosim * shs;

     #  ------------------------- do lunar terms ------------------
     dedt = ses + s1 * znl * s5;
     didt = sis + s2 * znl * (z11 + z13);
     dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq);
     sghl = s4 * znl * (z31 + z33 - 6.0);
     shll = -znl * s2 * (z21 + z23);
     #  sgp4fix for 180 deg incl
     if inclm < 5.2359877e-2 or inclm > pi - 5.2359877e-2:
         shll = 0.0;
     domdt = sgs + sghl;
     dnodt = shs;
     if sinim != 0.0:

         domdt = domdt - cosim / sinim * shll;
         dnodt = dnodt + shll / sinim;


     #  ----------- calculate deep space resonance effects --------
     dndt   = 0.0;
     theta  = (gsto + tc * rptim) % twopi
     em     = em + dedt * t;
     inclm  = inclm + didt * t;
     argpm  = argpm + domdt * t;
     nodem  = nodem + dnodt * t;
     mm     = mm + dmdt * t;
     

     #  -------------- initialize the resonance terms -------------
     if irez != 0:

         aonv = pow(nm / xke, x2o3);

         #  ---------- geopotential resonance for 12 hour orbits ------
         if irez == 2:

             cosisq = cosim * cosim;
             emo    = em;
             em     = ecco;
             emsqo  = emsq;
             emsq   = eccsq;
             eoc    = em * emsq;
             g201   = -0.306 - (em - 0.64) * 0.440;

             if em <= 0.65:

                 g211 =    3.616  -  13.2470 * em +  16.2900 * emsq;
                 g310 =  -19.302  + 117.3900 * em - 228.4190 * emsq +  156.5910 * eoc;
                 g322 =  -18.9068 + 109.7927 * em - 214.6334 * emsq +  146.5816 * eoc;
                 g410 =  -41.122  + 242.6940 * em - 471.0940 * emsq +  313.9530 * eoc;
                 g422 = -146.407  + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc;
                 g520 = -532.114  + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc;

             else:

                 g211 =   -72.099 +   331.819 * em -   508.738 * emsq +   266.724 * eoc;
                 g310 =  -346.844 +  1582.851 * em -  2415.925 * emsq +  1246.113 * eoc;
                 g322 =  -342.585 +  1554.908 * em -  2366.899 * emsq +  1215.972 * eoc;
                 g410 = -1052.797 +  4758.686 * em -  7193.992 * emsq +  3651.957 * eoc;
                 g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc;
                 if em > 0.715:
                     g520 =-5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc;
                 else:
                     g520 = 1464.74 -  4664.75 * em +  3763.64 * emsq;

             if em < 0.7:

                 g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21  * eoc;
                 g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc;
                 g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4  * eoc;

             else:

                 g533 =-37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc;
                 g521 =-51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc;
                 g532 =-40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc;

             sini2=  sinim * sinim;
             f220 =  0.75 * (1.0 + 2.0 * cosim+cosisq);
             f221 =  1.5 * sini2;
             f321 =  1.875 * sinim  *  (1.0 - 2.0 * cosim - 3.0 * cosisq);
             f322 = -1.875 * sinim  *  (1.0 + 2.0 * cosim - 3.0 * cosisq);
             f441 = 35.0 * sini2 * f220;
             f442 = 39.3750 * sini2 * sini2;
             f522 =  9.84375 * sinim * (sini2 * (1.0 - 2.0 * cosim- 5.0 * cosisq) +
                     0.33333333 * (-2.0 + 4.0 * cosim + 6.0 * cosisq) );
             f523 = sinim * (4.92187512 * sini2 * (-2.0 - 4.0 * cosim +
                    10.0 * cosisq) + 6.56250012 * (1.0+2.0 * cosim - 3.0 * cosisq));
             f542 = 29.53125 * sinim * (2.0 - 8.0 * cosim+cosisq *
                    (-12.0 + 8.0 * cosim + 10.0 * cosisq));
             f543 = 29.53125 * sinim * (-2.0 - 8.0 * cosim+cosisq *
                    (12.0 + 8.0 * cosim - 10.0 * cosisq));
             xno2  =  nm * nm;
             ainv2 =  aonv * aonv;
             temp1 =  3.0 * xno2 * ainv2;
             temp  =  temp1 * root22;
             d2201 =  temp * f220 * g201;
             d2211 =  temp * f221 * g211;
             temp1 =  temp1 * aonv;
             temp  =  temp1 * root32;
             d3210 =  temp * f321 * g310;
             d3222 =  temp * f322 * g322;
             temp1 =  temp1 * aonv;
             temp  =  2.0 * temp1 * root44;
             d4410 =  temp * f441 * g410;
             d4422 =  temp * f442 * g422;
             temp1 =  temp1 * aonv;
             temp  =  temp1 * root52;
             d5220 =  temp * f522 * g520;
             d5232 =  temp * f523 * g532;
             temp  =  2.0 * temp1 * root54;
             d5421 =  temp * f542 * g521;
             d5433 =  temp * f543 * g533;
             xlamo =  (mo + nodeo + nodeo-theta - theta) % twopi
             xfact =  mdot + dmdt + 2.0 * (nodedot + dnodt - rptim) - no;
             em    = emo;
             emsq  = emsqo;

         #  ---------------- synchronous resonance terms --------------
         if irez == 1:

             g200  = 1.0 + emsq * (-2.5 + 0.8125 * emsq);
             g310  = 1.0 + 2.0 * emsq;
             g300  = 1.0 + emsq * (-6.0 + 6.60937 * emsq);
             f220  = 0.75 * (1.0 + cosim) * (1.0 + cosim);
             f311  = 0.9375 * sinim * sinim * (1.0 + 3.0 * cosim) - 0.75 * (1.0 + cosim);
             f330  = 1.0 + cosim;
             f330  = 1.875 * f330 * f330 * f330;
             del1  = 3.0 * nm * nm * aonv * aonv;
             del2  = 2.0 * del1 * f220 * g200 * q22;
             del3  = 3.0 * del1 * f330 * g300 * q33 * aonv;
             del1  = del1 * f311 * g310 * q31 * aonv;
             xlamo = (mo + nodeo + argpo - theta) % twopi
             xfact = mdot + xpidot - rptim + dmdt + domdt + dnodt - no;

         #  ------------ for sgp4, initialize the integrator ----------
         xli   = xlamo;
         xni   = no;
         atime = 0.0;
         nm    = no + dndt;

     return (
       em,    argpm,  inclm, mm,
       nm,    nodem,
       irez, atime,
       d2201, d2211,  d3210, d3222,
       d4410, d4422, d5220,  d5232,
       d5421, d5433, dedt,  didt,
       dmdt,  dndt, dnodt, domdt,
       del1,   del2,  del3, xfact,
       xlamo, xli,    xni,
       )