# -*- coding: UTF-8 -*-
"""
Instrumental throughput models and calibration and synthetic photometry
routines
"""

from __future__ import absolute_import
from __future__ import unicode_literals
import numpy as np
import pysynphot as S

def synflux(spec, pb):
    """
    Compute the synthetic flux of spectrum ``spec`` through passband ``pb``

    Parameters
    ----------
    spec : :py:class:`pysynphot.ArraySpectrum`
        The spectrum. Must have ``dtype=[('wave', '<f8'), ('flux', '<f8')]``
    pb : :py:class:`pysynphot.ArrayBandpass` or :py:class:`pysynphot.obsbandpass.ObsModeBandpass`
        The passband data.
        Must have ``dtype=[('wave', '<f8'), ('throughput', '<f8')]``

    Returns
    -------
    flux : float
        The normalized flux of the spectrum through the passband

    Notes
    -----
        The passband is assumed to be dimensionless photon transmission
        efficiency.

        Routine is intended to be a mch faster implementation of
        :py:meth:`pysynphot.observation.Observation.effstim`, since it is called over and
        over by the samplers as a function of model parameters.

        Uses :py:func:`numpy.trapz` for interpolation.
    """
    flux = spec.sample(pb.wave)
    n = np.trapz(flux*pb.wave*pb.throughput, pb.wave)
    d = np.trapz(pb.wave*pb.throughput, pb.wave)
    out = n/d
    return out


def synphot(spec, pb, zp=0.):
    """
    Compute the synthetic magnitude of spectrum ``spec`` through passband ``pb``

    Parameters
    ----------
    spec : :py:class:`pysynphot.ArraySpectrum`
        The spectrum. Must have ``dtype=[('wave', '<f8'), ('flux', '<f8')]``
    pb : array-like
        The passband transmission.
        Must satisfy ``pb.shape == spec[ind].flux.shape``
    zp : float, optional
        The zeropoint to apply to the synthetic flux

    Returns
    -------
    mag : float
        The synthetic magnitude of the spectrum through the passband

    See Also
    --------
    :py:func:`source_synphot.passband.synflux`
    """
    flux = synflux(spec, pb)
    m = -2.5*np.log10(flux) + zp
    return m


def get_pb_zpt(pb, reference='AB', model_mag=None):
    """
    Determines a passband zeropoint for synthetic photometry, given a reference
    standard and its model magnitude in the passband

    Parameters
    ----------
    pb : :py:class:`pysynphot.ArrayBandpass` or :py:class:`pysynphot.obsbandpass.ObsModeBandpass`
        The passband data.
        Must have ``dtype=[('wave', '<f8'), ('throughput', '<f8')]``
    reference : str, optional
        The name of the reference spectrophotometric standard to use to determine the passband zeropoint.
        'AB' or 'Vega' (default: 'AB')
    model_mag : float, optional 
        The magnitude of the reference spectrophotometric standard in the passband.
        default = None

    Returns
    -------
    pbzp : float
        passband AB zeropoint 
        
    Raises
    ------
    RuntimeError
        If the passband and reference standard do not overlap
    ValueError
        If the reference model_mag is invalid
    """

    # setup the photometric system by defining the standard and corresponding magnitude system
    if reference.lower() not in ('vega', 'ab'):
        message = 'Do not understand mag system reference standard {}. Must be AB or Vega'.format(reference)
        raise RuntimeError(message)
    if reference.lower() == 'ab':
        refspec = S.FlatSpectrum(3631, fluxunits='jy')
        mag_type= 'abmag'
    else:
        refspec    = S.Vega
        mag_type= 'vegamag'

    refspec.convert('flam')

    if model_mag is None:
        ob = S.Observation(refspec, pb)
        model_mag = ob.effstim(mag_type)
    
    synphot_mag = synphot(refspec, pb, zp=0.)
    outzp = model_mag - synphot_mag

    return outzp

@@ -186,16 +186,33 @@
 iterals%0A
+from . import io%0A
 import n
@@ -243,16 +243,52 @@
 hot as S
+%0Afrom collections import OrderedDict
 %0A%0Adef sy
@@ -3820,17 +3820,16 @@
 hot_mag%0A
-%0A
     retu
@@ -3836,9 +3836,758 @@
 rn outzp
+%0A%0A%0Adef get_pbs(pbnames, model_mags, model='AB'):%0A%0A    pbs = OrderedDict()%0A    if np.isscalar(pbnames):%0A        pbnames = np.array(pbnames, ndmin=1)%0A    else:%0A        pbnames = np.array(pbnames).flatten()%0A    npb = len(pbnames)%0A%0A    if np.isscalar(model_mags):%0A        model_mags = np.repeat(model_mags, npb)%0A    else:%0A        model_mags = np.array(model_mags).flatten()%0A%0A    if len(model_mags) != npb:%0A        message = 'Model mags and pbnames must be 1-D arrays with the same shape'%0A        raise ValueError(message)%0A%0A    for i, pb, mag in enumerate(zip(pbnames, model_mags)):%0A        thispb, _ = io.get_passband(pb)%0A        thispbzp = passband.get_pb_zpt(sourcepb, model_mag=mag, reference=model)%0A        pbs%5Bpb%5D = thispb, thispbzp%0A%0A    return pbs
 %0A

add method to bullk load several passbands