# Path: smpy/dust.py
# def Calzetti(ta_grid, wave, Av):
#     Att = np.ones(np.append(ta_grid.shape, len(wave)))
#     k = np.zeros_like(wave.value)
# 
#     w0 = wave <= 1200 * u.AA
#     w1 = wave < 6300 * u.AA
#     w2 = wave >= 6300 * u.AA
#     w_u = wave.to(u.um).value
# 
#     x1 = np.argmin(np.abs(wave - 1200 * u.AA))
#     x2 = np.argmin(np.abs(wave - 1250 * u.AA))
# 
#     k[w2] = 2.659 * (-1.857 + 1.040 / w_u[w2])
#     k[w1] = 2.659 * (-2.156 + (1.509 / w_u[w1]) - (0.198 / w_u[w1] ** 2) + (0.011 / w_u[w1] ** 3))
#     k[w0] = k[x1] + ((wave[w0] - 1200. * u.AA) * (k[x1] - k[x2]) / (wave[x1] - wave[x2]))
# 
#     k += 4.05
#     k[k < 0.] = 0.
# 
#     tv = Av * k / 4.05
#     Att *= np.power(10, -0.4 * tv)
#     return Att
# 
# Path: smpy/sfh.py
# def exponential(t, tau):
#     """ Exponential star formation history
# 
#     """
#     sfh = np.exp(-1 * t / tau) / np.abs(tau)
#     return sfh
# 
# Path: smpy/misc.py
# def tau_madau(wave, z):
#     """ Lyman limit system absorption following Madau 1995
#     
#     This version is in part based on my previous python implementation
#     and the simpler version shown here
#     https://github.com/spacetelescope/pysynphot/issues/77
#     
#     This version has been extended to include the additional
#     higher order lines calculated in EAZY. 
#     
#     Parameters:
#         wave (array): Array with astropy.units of length
#         z (float): redshift at which the source is being observed
#         
#     Returns:
#         LyL_absorption (array): Throughput along the line of sight
#             for wavelengths in wave from Lyman-alpha forest and 
#             Lyman-limit systems
#     
#     """
# 
#     lymanlim = 912. * u.AA
# 
#     lyn_w = np.array([1216, 1026., 972.8, 950., 938.1,
#                       931.0, 926.5, 923.4, 921.2,
#                       919.6, 918.4, 917.5, 916.7,
#                       916.1, 915.6, 915.2]) * u.AA
#     lyn_c = np.array([3.6e-3, 1.7e-3, 1.2e-3, 9.3e-4, 8.2e-4,
#                       7.5e-4, 7.1e-4, 6.8e-4, 6.6e-4,
#                       6.4e-4, 6.3e-4, 6.2e-4, 6.1e-4,
#                       6.0e-4, 6.0e-4, 6.0e-4])
#     tau = np.zeros(len(wave))
# 
#     for i in range(len(lyn_w)):
#         tau = np.where(wave <= lyn_w[i] * (1 + z), 
#                        tau + lyn_c[i] * (wave / lyn_w[i]) ** 3.46, 
#                        tau)
# 
#     wly = wave / lymanlim
#     wly3 = wly ** 3
#     pe_abs = (0.25 * wly3 * ((1 + z) ** 0.46 - wly ** 0.46) +
#               9.4 * wly ** 1.5 * ((1 + z) ** 0.18 - wly ** 0.18) -
#               0.7 * wly3 * (wly ** (-1.32) - (1 + z) ** (-1.32)) -
#               0.023 * ((1 + z) ** 1.68 - wly ** 1.68))
#     tau = np.where(wave <= lymanlim * (1 + z), tau + pe_abs, tau)
#     lyl_absorption = np.where(tau > 700., 0., np.exp(-tau))
#     return np.clip(lyl_absorption, 0., 1.)


import numpy as np
import copy
import os
import sys
import h5py
import types
import glob
import pkg_resources
import ast
from scipy.interpolate import griddata
from scipy import spatial
from astropy import units as u
from astropy import constants as c
from astropy import cosmology as cos
from astropy.utils.console import ProgressBar
from .dust import Calzetti
from .sfh import exponential
from .misc import tau_madau

        if hasattr(self.SSP, 'strm_arr') and hasattr(self.SSP, 'rmtm_arr'):
            self.strm_arr = np.array(self.SSP.strm_arr)
            self.rmtm_arr = np.array(self.SSP.rmtm_arr)

        self.iseds = np.array(self.SSP.iseds)

        # Find closest match for each tg value in ta - set tg to these values

        nebular_old = np.loadtxt(nebular_old_path, skiprows=1)
        self.neb_cont = nebular_old[:, 1]

        #self.neb_lines, self.neb_metal = nebular.inoue_lines(self.wave)
        #self.neb_cont = nebular.continuum(sek)

        self.neb_hlines = nebular_old[:, 2]
        self.neb_metal = nebular_old[:, 3:]
        self.neb_wave = nebular_old[:, 0]

        if None not in (age, sfh, dust, metal_ind):
            if f_esc == None:
                self.build(age, sfh, dust, metal_ind, sfh_law=sfh_law,
                           dust_model=dust_model,
                           neb_cont=neb_cont, neb_met=neb_met)

            else:
                self.build(age, sfh, dust, metal_ind,
                           f_esc, sfh_law, dust_model, neb_cont, neb_met)

    def build(self, age, sfh, dust, metal, fesc = 1.,

# Path: smpy/dust.py
# def Calzetti(ta_grid, wave, Av):
#     Att = np.ones(np.append(ta_grid.shape, len(wave)))
#     k = np.zeros_like(wave.value)
# 
#     w0 = wave <= 1200 * u.AA
#     w1 = wave < 6300 * u.AA
#     w2 = wave >= 6300 * u.AA
#     w_u = wave.to(u.um).value
# 
#     x1 = np.argmin(np.abs(wave - 1200 * u.AA))
#     x2 = np.argmin(np.abs(wave - 1250 * u.AA))
# 
#     k[w2] = 2.659 * (-1.857 + 1.040 / w_u[w2])
#     k[w1] = 2.659 * (-2.156 + (1.509 / w_u[w1]) - (0.198 / w_u[w1] ** 2) + (0.011 / w_u[w1] ** 3))
#     k[w0] = k[x1] + ((wave[w0] - 1200. * u.AA) * (k[x1] - k[x2]) / (wave[x1] - wave[x2]))
# 
#     k += 4.05
#     k[k < 0.] = 0.
# 
#     tv = Av * k / 4.05
#     Att *= np.power(10, -0.4 * tv)
#     return Att
# 
# Path: smpy/sfh.py
# def exponential(t, tau):
#     """ Exponential star formation history
# 
#     """
#     sfh = np.exp(-1 * t / tau) / np.abs(tau)
#     return sfh
# 
# Path: smpy/misc.py
# def tau_madau(wave, z):
#     """ Lyman limit system absorption following Madau 1995
#     
#     This version is in part based on my previous python implementation
#     and the simpler version shown here
#     https://github.com/spacetelescope/pysynphot/issues/77
#     
#     This version has been extended to include the additional
#     higher order lines calculated in EAZY. 
#     
#     Parameters:
#         wave (array): Array with astropy.units of length
#         z (float): redshift at which the source is being observed
#         
#     Returns:
#         LyL_absorption (array): Throughput along the line of sight
#             for wavelengths in wave from Lyman-alpha forest and 
#             Lyman-limit systems
#     
#     """
# 
#     lymanlim = 912. * u.AA
# 
#     lyn_w = np.array([1216, 1026., 972.8, 950., 938.1,
#                       931.0, 926.5, 923.4, 921.2,
#                       919.6, 918.4, 917.5, 916.7,
#                       916.1, 915.6, 915.2]) * u.AA
#     lyn_c = np.array([3.6e-3, 1.7e-3, 1.2e-3, 9.3e-4, 8.2e-4,
#                       7.5e-4, 7.1e-4, 6.8e-4, 6.6e-4,
#                       6.4e-4, 6.3e-4, 6.2e-4, 6.1e-4,
#                       6.0e-4, 6.0e-4, 6.0e-4])
#     tau = np.zeros(len(wave))
# 
#     for i in range(len(lyn_w)):
#         tau = np.where(wave <= lyn_w[i] * (1 + z), 
#                        tau + lyn_c[i] * (wave / lyn_w[i]) ** 3.46, 
#                        tau)
# 
#     wly = wave / lymanlim
#     wly3 = wly ** 3
#     pe_abs = (0.25 * wly3 * ((1 + z) ** 0.46 - wly ** 0.46) +
#               9.4 * wly ** 1.5 * ((1 + z) ** 0.18 - wly ** 0.18) -
#               0.7 * wly3 * (wly ** (-1.32) - (1 + z) ** (-1.32)) -
#               0.023 * ((1 + z) ** 1.68 - wly ** 1.68))
#     tau = np.where(wave <= lymanlim * (1 + z), tau + pe_abs, tau)
#     lyl_absorption = np.where(tau > 700., 0., np.exp(-tau))
#     return np.clip(lyl_absorption, 0., 1.)


# Path: smpy/smpy.py
import numpy as np
import copy
import os
import sys
import h5py
import types
import glob
import pkg_resources
import ast
from scipy.interpolate import griddata
from scipy import spatial
from astropy import units as u
from astropy import constants as c
from astropy import cosmology as cos
from astropy.utils.console import ProgressBar
from .dust import Calzetti
from .sfh import exponential
from .misc import tau_madau


        if hasattr(self.SSP, 'strm_arr') and hasattr(self.SSP, 'rmtm_arr'):
            self.strm_arr = np.array(self.SSP.strm_arr)
            self.rmtm_arr = np.array(self.SSP.rmtm_arr)

        self.iseds = np.array(self.SSP.iseds)

        # Find closest match for each tg value in ta - set tg to these values

        nebular_old = np.loadtxt(nebular_old_path, skiprows=1)
        self.neb_cont = nebular_old[:, 1]

        #self.neb_lines, self.neb_metal = nebular.inoue_lines(self.wave)
        #self.neb_cont = nebular.continuum(sek)

        self.neb_hlines = nebular_old[:, 2]
        self.neb_metal = nebular_old[:, 3:]
        self.neb_wave = nebular_old[:, 0]

        if None not in (age, sfh, dust, metal_ind):
            if f_esc == None:
                self.build(age, sfh, dust, metal_ind, sfh_law=sfh_law,
                           dust_model=dust_model,
                           neb_cont=neb_cont, neb_met=neb_met)

            else:
                self.build(age, sfh, dust, metal_ind,
                           f_esc, sfh_law, dust_model, neb_cont, neb_met)

    def build(self, age, sfh, dust, metal, fesc = 1.,

              sfh_law = exponential, dust_model = Calzetti,