# Path: smpy/smpy.py
# class CSP:
# class Filter:
# class FileFilter(Filter):
# class TophatFilter(Filter):
# class LoadEAZYFilters:
# class FilterSet:
# class Observe:
# class ObserveToFile:
#     def __init__(self, ssp, age=None, sfh=None, dust=None, metal_ind=None,
#                  f_esc=None, sfh_law=exponential, dust_model=Calzetti,
#                  neb_cont=True, neb_met=True):
#     def build(self, age, sfh, dust, metal, fesc = 1.,
#               sfh_law = exponential, dust_model = Calzetti,
#               neb_dust_weight = 1.,
#               neb_cont = True, neb_met = True,
#               timesteps = 500, verbose = False):
#     def calc_lyman(wave, seds):
#     def __add__(self, other):
#     def __iadd__(self, other):
#     def __mul__(self, other):
#     def __imul__(self, other):
#     def __div__(self, other):
#     def __idiv__(self, other):
#     def __rmul__(self, other):
#     def addEmissionLine(self, wavelength, EqW):
#     def save_to_hdf(self, hdfpath):
#     def load_from_hdf(self, hdfpath):
#     def __init__(self):
#     def __init__(self, filepath, maxbins=500):
#     def __init__(self, centre, width, steps=200):
#     def __init__(self, path):
#     def search(self, searchterm, error=0.1):
#     def __init__(self, path=None):
#     def addFileFilter(self, path):
#     def addTophatFilter(self, centre, width, steps=200):
#     def addEAZYFilter(self, EAZYset, N):
#     def __init__(self, SED, Filters, redshift, force_age=True,
#                  madau=True, units=u.uJy, verbose=False):
#     def calcflux(self, SED, filt, z, dl, units):
#     def __getitem__(self, items):
#     def __init__(self):
#     def build(self, SED, Filters, redshift, savepath,
#               force_age=True, madau=True, units=u.uJy,
#               verbose=False, clobber=True):
#     def load(self, loadpath):
#     def calcflux(self, SED, wf, tp, z, dl, units):
#     def __getitem__(self, items):
#             X = tri_grid.transform[ss, :2]
#             Y = points - tri_grid.transform[ss, 2]
#         S = np.array(seds.shape)
#         S[-1] = n+1
# 
# Path: smpy/ssp.py
# class Ised:
# class SSP:
# class BC(SSP):
# class BPASS(SSP):
# class BPASS2(SSP):
#     def __init__(self, path):
#     def read_ised(self, filename):
#     def __init__(self, path=None):
#     def __init__(self, path='../ssp/bc03/salpeter/lr/', verbose=False):
#     def __init__(self, path='../ssp/bpass_v1.1/SEDS/', nsteps = 3000):
#     def __init__(self, path='../ssp/bpass_v1.1/SEDS/', nsteps = 3000):
# def rebin_sed(wave_in, sed_in, wave_out):
#             G = array.array('f')
#             H = array.array('f')
#                 F = array.array('i')
#                 G = array.array('f')
#                 H = array.array('f')
# 
# Path: smpy/sfh.py
# def dblpower(t, sl):
#     """ Power-law star formation history
#     """
#     alpha, beta, tau = sl
#     sfh = (np.power(np.abs(t)/(tau*np.max(t)), alpha) + np.power(np.abs(t)/(tau*np.max(t)), beta))**-1
#     return sfh


import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
from smpy import smpy as S
from smpy import ssp as B
from smpy.sfh import dblpower

BC = B.BC('data/ssp/bc03/chab/lr/')
models = S.CSP(BC)

tau = 1*u.Gyr

age = 3*u.Gyr
sfh = [(1., -1., 0.3), (1., -1., 0.5)]
dust = 0.
metal = 1.

#BC = B.BC('')

def make_dbl_sfhs(nalpha, nbeta, ntau,
                  lalpha_min=-1., lalpha_max=3.,
                  lbeta_min=-1., lbeta_max=3.,
                  tau_min=0.1, tau_max=1.0):

    alphas, betas, taus = np.meshgrid(np.logspace(lalpha_min, lalpha_max, nalpha),

# Path: smpy/smpy.py
# class CSP:
# class Filter:
# class FileFilter(Filter):
# class TophatFilter(Filter):
# class LoadEAZYFilters:
# class FilterSet:
# class Observe:
# class ObserveToFile:
#     def __init__(self, ssp, age=None, sfh=None, dust=None, metal_ind=None,
#                  f_esc=None, sfh_law=exponential, dust_model=Calzetti,
#                  neb_cont=True, neb_met=True):
#     def build(self, age, sfh, dust, metal, fesc = 1.,
#               sfh_law = exponential, dust_model = Calzetti,
#               neb_dust_weight = 1.,
#               neb_cont = True, neb_met = True,
#               timesteps = 500, verbose = False):
#     def calc_lyman(wave, seds):
#     def __add__(self, other):
#     def __iadd__(self, other):
#     def __mul__(self, other):
#     def __imul__(self, other):
#     def __div__(self, other):
#     def __idiv__(self, other):
#     def __rmul__(self, other):
#     def addEmissionLine(self, wavelength, EqW):
#     def save_to_hdf(self, hdfpath):
#     def load_from_hdf(self, hdfpath):
#     def __init__(self):
#     def __init__(self, filepath, maxbins=500):
#     def __init__(self, centre, width, steps=200):
#     def __init__(self, path):
#     def search(self, searchterm, error=0.1):
#     def __init__(self, path=None):
#     def addFileFilter(self, path):
#     def addTophatFilter(self, centre, width, steps=200):
#     def addEAZYFilter(self, EAZYset, N):
#     def __init__(self, SED, Filters, redshift, force_age=True,
#                  madau=True, units=u.uJy, verbose=False):
#     def calcflux(self, SED, filt, z, dl, units):
#     def __getitem__(self, items):
#     def __init__(self):
#     def build(self, SED, Filters, redshift, savepath,
#               force_age=True, madau=True, units=u.uJy,
#               verbose=False, clobber=True):
#     def load(self, loadpath):
#     def calcflux(self, SED, wf, tp, z, dl, units):
#     def __getitem__(self, items):
#             X = tri_grid.transform[ss, :2]
#             Y = points - tri_grid.transform[ss, 2]
#         S = np.array(seds.shape)
#         S[-1] = n+1
# 
# Path: smpy/ssp.py
# class Ised:
# class SSP:
# class BC(SSP):
# class BPASS(SSP):
# class BPASS2(SSP):
#     def __init__(self, path):
#     def read_ised(self, filename):
#     def __init__(self, path=None):
#     def __init__(self, path='../ssp/bc03/salpeter/lr/', verbose=False):
#     def __init__(self, path='../ssp/bpass_v1.1/SEDS/', nsteps = 3000):
#     def __init__(self, path='../ssp/bpass_v1.1/SEDS/', nsteps = 3000):
# def rebin_sed(wave_in, sed_in, wave_out):
#             G = array.array('f')
#             H = array.array('f')
#                 F = array.array('i')
#                 G = array.array('f')
#                 H = array.array('f')
# 
# Path: smpy/sfh.py
# def dblpower(t, sl):
#     """ Power-law star formation history
#     """
#     alpha, beta, tau = sl
#     sfh = (np.power(np.abs(t)/(tau*np.max(t)), alpha) + np.power(np.abs(t)/(tau*np.max(t)), beta))**-1
#     return sfh


# Path: scripts/dblpower_test.py
import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
from smpy import smpy as S
from smpy import ssp as B
from smpy.sfh import dblpower


BC = B.BC('data/ssp/bc03/chab/lr/')
models = S.CSP(BC)

tau = 1*u.Gyr

age = 3*u.Gyr
sfh = [(1., -1., 0.3), (1., -1., 0.5)]
dust = 0.
metal = 1.

#BC = B.BC('')

def make_dbl_sfhs(nalpha, nbeta, ntau,
                  lalpha_min=-1., lalpha_max=3.,
                  lbeta_min=-1., lbeta_max=3.,
                  tau_min=0.1, tau_max=1.0):

    alphas, betas, taus = np.meshgrid(np.logspace(lalpha_min, lalpha_max, nalpha),

                                      -1*np.logspace(lbeta_min, lbeta_max, nbeta),