From Morrissey+ 2005, with the actual data coming from
        http://www.astro.caltech.edu/~capak/filters/. According to the latter,
        these are in QE units and thus need to be multiplied by the wavelength
        when integrating per-energy.
517881
def _load_data(self, band):
        
        # `band` should be 'nuv' or 'fuv'
        df = bandpass_data_frame('filter_galex_' + band + '.dat', 'wlen resp')
        df.resp *= df.wlen # QE -> EE response convention.
        return df
487715
def Δk(self):
        
        Δc, Δv = self.Δc, self.Δv
        cos = self._energy.trig('cos θ')
        return lambda e: (1 / 2) * ((Δc + Δv) + (Δc - Δv) * cos(e))
1078245
def energy_ES(q, v):
    """"""
    # Body 0 is the sun, Body 1 is the earth
    m0 = mass[0]
    m1 = mass[1]

    # Positions of sun and earth
    q0: np.ndarray = q[:, slices[0]]
    q1: np.ndarray = q[:, slices[1]]

    # Velocities of sun and earth
    v0: np.ndarray = v[:, slices[0]]
    v1: np.ndarray = v[:, slices[1]]

    # Kinetic energy is 1/2 mv^2
    T0: np.ndarray = 0.5 * m0 * np.sum(v0 * v0, axis=1)
    T1: np.ndarray = 0.5 * m1 * np.sum(v1 * v1, axis=1)
    T: np.ndarray = T0 + T1
    
    # Potential energy is -G m1 m2  / r
    dv_01 = q1 - q0
    r_01 = np.linalg.norm(dv_01, axis=1)
    U_01: np.ndarray = -G * m0 * m1 / r_01
    U: np.ndarray = U_01
    
    # Total energy H = T + U
    H = T + U
    
    return H, T, U
682372
def waveset(self):
        """"""
        w = get_waveset(self.model)
        if w is not None:
            utils.validate_wavelengths(w)
            w = w * self._internal_wave_unit
        return w
592288
def unitResponse(self, band):
        """"""
        #sumfilt(wave,0,band)
        # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ]
        total = band.trapezoidIntegration(band.wave, band.throughput)
        return 1.0/total
1349549
def smoothness(self):
        
        energies = self.energies
        try:
            sp = self.spline()
        except:
            print("Energy spline failed.")
            return None
        spline_energies = sp(range(len(energies)))
        diff = spline_energies - energies
        rms = np.sqrt(np.sum(np.square(diff)) / len(energies))
        return rms
113349
def optimise(self, data, initial_theta=None, full_output=False, **kwargs):
        

        data = self._format_data(data)

        if initial_theta is None:
            initial_theta = self.estimate(data)

        # Which parameters will be optimised, and which will be fixed?
        matched_channels, missing_channels, ignore_parameters \
            = self._match_channels_to_data(data)

        # For the Cannon model, the normalised fluxes are a multiplicative
        # operation, so we don't actually need to optimise the astrophysical
        # parameters, we only need to optimise the other ones. At each set of
        # non-astrophysical parameters we can determine the optimal
        # astrophysical parameters by a matrix inversion.
        parameters = [p for p in self.parameters if p not in ignore_parameters \
            and p not in self.grid_points.dtype.names]

        # What model wavelength ranges will be required?
        wavelengths_required = []
        for channel, spectrum in zip(matched_channels, data):
            if channel is None: continue
            z = initial_theta.get("z", initial_theta.get("z_{}".format(channel),
                0))
            wavelengths_required.append(
                [spectrum.disp[0] * (1 - z), spectrum.disp[-1] * (1 - z)])

        # (If necessary,..) train the Cannon model around the closest point.
        subset_bounds = self._initialise_approximator(closest_point=\
            [initial_theta[p] for p in self.grid_points.dtype.names], 
            wavelengths_required=wavelengths_required, **kwargs)
        
        # Get the optimisation keyword arguments.
        op_kwargs = self._configuration.get("optimise", {}).copy()
        op_kwargs.update(kwargs)

        # Get fixed keywords.
        fixed = op_kwargs.pop("fixed", {})
        if fixed is None: fixed = {}
        if fixed:
            # Remove non-parameters from the 'fixed' keywords.
            keys = set(fixed).intersection(self.parameters)
            ignored_keys = set(fixed).difference(parameters)
            if ignored_keys:
                logger.warn("Ignoring optimise.fixed keys because they are not"
                    " model parameters: {0}".format(", ".join(ignored_keys)))
            # If the 'fixed' value is provided, use that. Otherwise if it is
            # None then use the initial_theta value.
            fixed = dict(zip(keys, 
                [(fixed[k], initial_theta.get(k, None))[fixed[k] is None] \
                    for k in keys]))

            logger.info("Fixing keyword arguments (these will not be optimised)"
                ": {}".format(fixed))

        # Remove fixed parameters from the parameters to be optimised
        #parameters = list(set(parameters).difference(fixed))
        parameters = [p for p in parameters if p not in fixed]

        # Translate input bounds.
        nbs = (None, None) # No boundaries.
        input_bounds = op_kwargs.pop("bounds", {})
        # Because astrophysical parameters are not being 'optimised', we may not
        # actually need the subset_bounds.get call, but there may be fringe
        # cases that I cannot immediately think of.
        op_kwargs["bounds"] = [input_bounds.get(p, subset_bounds.get(p, nbs)) \
            for p in parameters]

        # Apply data masks now so we don't have to do it on the fly.
        masked_data, pixels_affected = self._apply_data_mask(data)

        # Prepare the convolution functions.
        self._create_convolution_functions(matched_channels, data, parameters)

        logger.info("Optimising parameters: {0}".format(", ".join(parameters)))
        logger.info("Optimisation keywords: {0}".format(op_kwargs))

        # Create the objective function. The objective function needs to
        # continuum-normalise the observed spectra, put it at rest, then solve
        # for the best astrophysical parameters, then *generate* that spectrum.
        debug = kwargs.get("debug", False)
        def nlp(t, return_labels=False):
            # Apply fixed keywords.
            theta = dict(zip(parameters, t))
            theta.update(fixed)

            observed_variances = []
            observed_intensities = []
            for channel, spectrum in zip(matched_channels, data):

                # For each channel:
                # 1) Put the observed spectrum at rest.
                z = theta.get("z", theta.get("z_{}".format(channel), 0))
                rest_observed_disp = spectrum.disp * (1 - z)

                # 2) Correct for continuum.
                j, coeffs = 0, []
                while theta.get("continuum_{0}_{1}".format(channel, j), None):
                    coeffs.append(theta["continuum_{0}_{1}".format(channel, j)])
                    j += 1

                # Remember: model.__call__ calculates continuum based on the
                # *observed* wavelength points, so here we do the same (e.g.,
                # not those that have potentially been corrected for redshift)
                if not coeffs: continuum = 1.0
                else:
                    continuum = np.polyval(coeffs[::-1], spectrum.disp) 
                
                # 3) Put the observed data onto the self.wavelengths scale.
                # [TODO] do the resampling correctly.
                rebinned_observed_intensities = np.interp(self.wavelengths,
                    rest_observed_disp, spectrum.flux / continuum,
                    left=np.nan, right=np.nan)

                rebinned_observed_variances = np.interp(self.wavelengths,
                    rest_observed_disp, spectrum.variance / continuum,
                    left=np.nan, right=np.nan)

                observed_variances.append(rebinned_observed_variances)
                observed_intensities.append(rebinned_observed_intensities)

            # [TODO] This may be the wrong thing to do.
            observed_variances \
                = np.nanmean(np.vstack(observed_variances), axis=0)
            observed_intensities \
                = np.nanmean(np.vstack(observed_intensities), axis=0)

            # Solve for the astrophysical parameters.
            try:
                labels = self._solve_labels(observed_intensities,
                    observed_variances)

            except:
                logger.exception("Could not determine labels:")
                if debug: raise
                return np.inf

            if return_labels:
                return labels

            # Put the optimal astrophysical parameters into our theta.
            theta.update(labels)

            # Make the log-probability call.
            # Note: The inference._ln_probability call takes a dictionary.
            return -inference._ln_probability(theta, self, data, debug,
                matched_channels=matched_channels)

        # Do the optimisation.
        p0 = np.array([initial_theta[p] for p in parameters])
        x_opt = optimise.minimise(nlp, p0, **op_kwargs)

        # Put the result into a usable form.
        result = dict(zip(parameters, x_opt))
        result.update(nlp(x_opt, return_labels=True))
        result.update(fixed)
        
        result = OrderedDict([(k, result[k]) \
            for k in self.parameters if k in result])

        if full_output:
            # Create model fluxes and calculate some metric.
            chi_sq, dof, model_fluxes = self._chi_sq(result, data)
            result = (result, chi_sq, dof, model_fluxes)

        # Remove any prepared convolution functions.
        self._destroy_convolution_functions()

        return result
609462
def _compute_fluxes(self):
        ''' '''
        self.emission = self._compute_emission()
        self.emission_sfc = self._compute_emission_sfc()
        fromspace = self._from_space()
        self.flux_down = self.trans.flux_down(fromspace, self.emission)
        self.flux_reflected_up = self.trans.flux_reflected_up(self.flux_down, self.albedo_sfc)
        # this ensure same dimensions as other fields
        self.flux_to_sfc = self.flux_down[..., -1, np.newaxis]
        self.flux_from_sfc = (self.emission_sfc +
                              self.flux_reflected_up[..., -1, np.newaxis])
        self.flux_up = self.trans.flux_up(self.flux_from_sfc,
                            self.emission + self.flux_reflected_up[...,0:-1])
        self.flux_net = self.flux_up - self.flux_down
        # absorbed radiation (flux convergence) in W / m**2 (per band)
        self.absorbed = np.diff(self.flux_net, axis=-1)
        self.absorbed_total = np.sum(self.absorbed, axis=-1)
        self.flux_to_space = self._compute_flux_top()
1276514
def sed_conversion(energy, model_unit, sed):
    

    model_pt = model_unit.physical_type

    ones = np.ones(energy.shape)

    if sed:
        # SED
        f_unit = u.Unit("erg/s")
        if model_pt == "power" or model_pt == "flux" or model_pt == "energy":
            sedf = ones
        elif "differential" in model_pt:
            sedf = energy ** 2
        else:
            raise u.UnitsError(
                "Model physical type ({0}) is not supported".format(model_pt),
                "Supported physical types are: power, flux, differential"
                " power, differential flux",
            )

        if "flux" in model_pt:
            f_unit /= u.cm ** 2
        elif "energy" in model_pt:
            # particle energy distributions
            f_unit = u.erg
    else:
        # Differential spectrum
        f_unit = u.Unit("1/(s TeV)")
        if "differential" in model_pt:
            sedf = ones
        elif model_pt == "power" or model_pt == "flux" or model_pt == "energy":
            # From SED to differential
            sedf = 1 / (energy ** 2)
        else:
            raise u.UnitsError(
                "Model physical type ({0}) is not supported".format(model_pt),
                "Supported physical types are: power, flux, differential"
                " power, differential flux",
            )

        if "flux" in model_pt:
            f_unit /= u.cm ** 2
        elif "energy" in model_pt:
            # particle energy distributions
            f_unit = u.Unit("1/TeV")

    log.debug(
        "Converted from {0} ({1}) into {2} ({3}) for sed={4}".format(
            model_unit, model_pt, f_unit, f_unit.physical_type, sed
        )
    )

    return f_unit, sedf
1270774
def perform_es_corr(self, lattice, q, step=1e-4):
        
        logger.info("Running Freysoldt 2011 PC calculation (should be " "equivalent to sxdefectalign)")
        logger.debug("defect lattice constants are (in angstroms)" + str(lattice.abc))

        [a1, a2, a3] = ang_to_bohr * np.array(lattice.get_cartesian_coords(1))
        logging.debug("In atomic units, lat consts are (in bohr):" + str([a1, a2, a3]))
        vol = np.dot(a1, np.cross(a2, a3))  # vol in bohr^3

        def e_iso(encut):
            gcut = eV_to_k(encut)  # gcut is in units of 1/A
            return scipy.integrate.quad(lambda g: self.q_model.rho_rec(g * g)**2, step, gcut)[0] * (q**2) / np.pi

        def e_per(encut):
            eper = 0
            for g2 in generate_reciprocal_vectors_squared(a1, a2, a3, encut):
                eper += (self.q_model.rho_rec(g2)**2) / g2
            eper *= (q**2) * 2 * round(np.pi, 6) / vol
            eper += (q**2) * 4 * round(np.pi, 6) \
                * self.q_model.rho_rec_limit0 / vol
            return eper

        eiso = converge(e_iso, 5, self.madetol, self.energy_cutoff)
        logger.debug("Eisolated : %f", round(eiso, 5))

        eper = converge(e_per, 5, self.madetol, self.energy_cutoff)

        logger.info("Eperiodic : %f hartree", round(eper, 5))
        logger.info("difference (periodic-iso) is %f hartree", round(eper - eiso, 6))
        logger.info("difference in (eV) is %f", round((eper - eiso) * hart_to_ev, 4))

        es_corr = round((eiso - eper) / self.dielectric * hart_to_ev, 6)
        logger.info("Defect Correction without alignment %f (eV): ", es_corr)
        return es_corr
199002
def unitResponse(self,band):
        """"""
        #sumfilt(wave,1,band)
        # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ]
        wave=band.wave
        total = band.trapezoidIntegration(wave,band.throughput*wave)
        modtot = total/(H*C)
        return 2.5*math.log10(modtot) + STZERO
1128232
def Kerr_factor(final_mass, distance):
    

    # Convert solar masses to meters
    mass = final_mass * lal.MSUN_SI * lal.G_SI / lal.C_SI**2
    # Convert Mpc to meters
    dist = distance * 1e6 * lal.PC_SI

    return mass / dist
1302578
def ToInternal(self):
        
        self.validate_units()

        savewunits = self.waveunits
        savefunits = self.fluxunits

        if hasattr(self, 'primary_area'):
            area = self.primary_area
        else:
            area = None

        angwave = self.waveunits.Convert(self.GetWaveSet(), 'angstrom')
        phoflux = self.fluxunits.Convert(angwave, self._fluxtable, 'photlam',
                                         area=area)

        self._wavetable = angwave.copy()
        self._fluxtable = phoflux.copy()

        self.waveunits = savewunits
        self.fluxunits = savefunits
255315
def _get_binned_arrays(self, wavelengths, flux_unit, area=None,
                           vegaspec=None):
        """"""
        x = self._validate_binned_wavelengths(wavelengths)
        y = self.sample_binned(wavelengths=x, flux_unit=flux_unit, area=area,
                               vegaspec=vegaspec)

        if isinstance(wavelengths, u.Quantity):
            w = x.to(wavelengths.unit, u.spectral())
        else:
            w = x

        return w, y
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
585271
def optimise(self, data, initial_theta=None, full_output=False, **kwargs):
        

        data = self._format_data(data)

        if initial_theta is None:
            initial_theta = self.estimate(data)

        # Which parameters will be optimised, and which will be fixed?
        matched_channels, missing_channels, ignore_parameters \
            = self._match_channels_to_data(data)
        #parameters = set(self.parameters).difference(ignore_parameters)
        parameters = [p for p in self.parameters if p not in ignore_parameters]

        # What model wavelength ranges will be required?
        wavelengths_required = []
        for channel, spectrum in zip(matched_channels, data):
            if channel is None: continue
            z = initial_theta.get("z",
                initial_theta.get("z_{}".format(channel), 0))
            wavelengths_required.append(
                [spectrum.disp[0] * (1 - z), spectrum.disp[-1] * (1 - z)])

        # Create the spectrum approximator/interpolator.
        closest_point = [initial_theta[p] for p in self.grid_points.dtype.names]
        subset_bounds = self._initialise_approximator(
            closest_point=closest_point, 
            wavelengths_required=wavelengths_required, **kwargs)
        
        # Get the optimisation keyword arguments.
        op_kwargs = self._configuration.get("optimise", {}).copy()
        op_kwargs.update(kwargs)

        # Get fixed keywords.
        fixed = op_kwargs.pop("fixed", {})
        if fixed:
            # Remove non-parameters from the 'fixed' keywords.
            keys = set(fixed).intersection(parameters)
            # If the 'fixed' value is provided, use that. Otherwise if it is
            # None then use the initial_theta value.
            fixed = dict(zip(keys, 
                [(fixed[k], initial_theta.get(k, None))[fixed[k] is None] \
                    for k in keys]))

            logger.info("Fixing keyword arguments (these will not be optimised)"
                ": {}".format(fixed))

        # Remove fixed parameters from the parameters to be optimised
        #parameters = list(set(parameters).difference(fixed))
        parameters = [p for p in parameters if p not in fixed]

        # Translate input bounds.
        nbs = (None, None) # No boundaries.
        input_bounds = op_kwargs.pop("bounds", {})
        op_kwargs["bounds"] = [input_bounds.get(p, subset_bounds.get(p, nbs)) \
            for p in parameters]
        
        # Apply data masks now so we don't have to do it on the fly.
        masked_data, pixels_affected = self._apply_data_mask(data)

        # Prepare the convolution functions.
        self._create_convolution_functions(matched_channels, data, parameters,
            fixed_parameters=fixed)

        logger.info("Optimising parameters: {0}".format(", ".join(parameters)))
        logger.info("Optimisation keywords: {0}".format(op_kwargs))

        # Create the objective function.
        debug = kwargs.get("debug", False)
        def nlp(theta):
            # Apply fixed keywords
            t_ = theta.copy()
            p_ = [] + parameters
            for parameter, value in fixed.iteritems():
                p_.append(parameter)
                t_ = np.append(theta, value)

            return -inference.ln_probability(t_, p_, self, data, debug,
                matched_channels=matched_channels)

        # Do the optimisation.
        p0 = np.array([initial_theta[p] for p in parameters])

        x_opt = op.minimise(nlp, p0, **op_kwargs)

        # Put the result into a usable form.
        x_opt_theta = OrderedDict(zip(parameters, x_opt))
        x_opt_theta.update(fixed)

        if full_output:
            # Create model fluxes and calculate some metric.
            chi_sq, dof, model_fluxes = self._chi_sq(x_opt_theta, data)


            # Remove any prepared convolution functions.
            self._destroy_convolution_functions()
    
            return (x_opt_theta, chi_sq, dof, model_fluxes)

        # Remove any prepared convolution functions.
        self._destroy_convolution_functions()

        return x_opt_theta
1230659
def capacitance(self, filter_order=None, window_size=None, tol=0.05):
        '''
        
        '''
        C = np.ma.masked_invalid(1.0 / (2.0 * math.pi * self.frequency *
                   self.Z_device(filter_order=filter_order,
                                 window_size=window_size, tol=tol)))
        C.fill_value = np.nan
        C.data[C.mask] = C.fill_value
        return C