0415 first upload

LICENSE

+Copyright (c) 2018 The Python Packaging Authority
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
\ No newline at end of file

PKG-INFO

+Metadata-Version: 2.1
+Name: csst-ifs-gehong
+Version: 1.0.0
+Home-page: https://csst-ifs-gehong.readthedocs.io/en/latest/
+Author: Shuai Feng
+Author-email: sfeng@hebtu.edu.cn
+License: MIT
+Keywords: CSST-IFS
+License-File: LICENSE

README.md

-# GEHONG
+# csst-ifs-gehong Package
 
-GEnerate tHe data Of iNtegral field spectrograph of Galaxy
+This is a Python package for modelling the data of intergral field spectrascopy mounted on the Chinese Space Station Telescopy (CSST-IFS). See more detail at [csst-ifs-gehong](https://csst-ifs-gehong.readthedocs.io/).
\ No newline at end of file

setup.cfg

+[metadata]
+description-file = README.md
+license_files = LICENSE
+
+[egg_info]
+tag_build = 
+tag_date = 0
+

setup.py

+from setuptools import setup, find_packages
+
+
+setup(
+    name='csst-ifs-gehong',
+    version='1.0.0',
+    license='MIT',
+    author="Shuai Feng",
+    author_email='sfeng@hebtu.edu.cn',
+    packages=find_packages('src'),
+    package_dir={'': 'src'},
+    url='https://csst-ifs-gehong.readthedocs.io/en/latest/',
+    keywords='CSST-IFS',
+    install_requires=[
+          'astropy',
+      ],
+
+)
\ No newline at end of file

src/csst_ifs_gehong.egg-info/PKG-INFO

+Metadata-Version: 2.1
+Name: csst-ifs-gehong
+Version: 1.0.0
+Home-page: https://csst-ifs-gehong.readthedocs.io/en/latest/
+Author: Shuai Feng
+Author-email: sfeng@hebtu.edu.cn
+License: MIT
+Keywords: CSST-IFS
+License-File: LICENSE

src/csst_ifs_gehong.egg-info/SOURCES.txt

+LICENSE
+README.md
+setup.cfg
+setup.py
+src/csst_ifs_gehong.egg-info/PKG-INFO
+src/csst_ifs_gehong.egg-info/SOURCES.txt
+src/csst_ifs_gehong.egg-info/dependency_links.txt
+src/csst_ifs_gehong.egg-info/requires.txt
+src/csst_ifs_gehong.egg-info/top_level.txt
+src/gehong/__init__.py
+src/gehong/config.py
+src/gehong/cube3d.py
+src/gehong/map2d.py
+src/gehong/spec1d.py
+test/test_cube3d.py
+test/test_map2d.py
+test/test_spec1d.py
\ No newline at end of file

src/csst_ifs_gehong.egg-info/dependency_links.txt

+

src/csst_ifs_gehong.egg-info/requires.txt

+astropy

src/csst_ifs_gehong.egg-info/top_level.txt

+gehong

src/gehong/config.py

+import numpy as np
+
+class config():
+    """
+    The configuration of spectral modeling. The default value is used for ETC calculation.
+
+    Parameters
+    ----------
+    wave_min : float, optional
+        Minimum value of wavelength coverage, by default 3500.0A
+    wave_max : float, optional
+        Minimum value of wavelength coverage, by default 10000.0A
+    dlam : float, optional
+        Wavelength width of each spaxel, by default 2.0A
+    inst_fwhm : float, optional
+        Spectral resolution, by default 0.1A
+    nx : int, optional
+        Number of spaxel in a spatial direction, by default 30
+    ny : int, optional
+        Number of spaxel in a spatial direction, by default 30
+    dpix : float, optional
+        Pixel size in the spatial direction, by default 0.2arcsec
+    """
+    def __init__(self, wave_min = 3500.0, wave_max = 10000.0, 
+                 dlam = 2.0, inst_fwhm = 0.1,
+                 nx = 30, ny = 30, dpix = 0.2):
+        self.dlam = dlam
+        self.wave = np.arange(wave_min, wave_max, dlam)
+        self.wave_min = wave_min
+        
+        self.inst_fwhm = inst_fwhm
+        self.nx = nx
+        self.ny = ny
+        self.dpix = dpix
+        self.fov_x = nx * dpix
+        self.fov_y = ny * dpix
\ No newline at end of file

src/gehong/cube3d.py

+import astropy.units as u
+import numpy as np
+from scipy.interpolate import interp1d
+from .spec1d import *
+import astropy.wcs
+
+class Cube3D():
+    """
+    Class of 3-dimentional spectral cube
+    """
+    def __init__(self, config, stellar_map = None, gas_map = None):
+        
+        self.config= config
+        
+        self.nx    = config.nx
+        self.ny    = config.ny
+        self.dpix  = config.dpix
+        self.fov_x = config.fov_x
+        self.fov_y = config.fov_y
+        
+        self.wave  = config.wave
+        self.nz    = len(self.wave)
+        self.wave0 = np.min(self.wave)
+        self.inst_fwhm = config.inst_fwhm
+        
+        self.flux  = np.zeros((self.nx,self.ny,self.nz))
+        
+        self.stellar_map = stellar_map
+        self.gas_map     = gas_map
+        
+    def make_cube(self, stellar_tem = None, hii_tem = None):
+        
+        for i in range(self.nx):
+            for j in range(self.ny):
+                if self.stellar_map is not None:
+                    ss = StellarContinuum(self.config, stellar_tem, mag = self.stellar_map.mag[i,j], 
+                                          age = self.stellar_map.age[i,j], feh = self.stellar_map.feh[i,j], 
+                                          vel = self.stellar_map.vel[i,j], vdisp = self.stellar_map.vdisp[i,j], 
+                                          ebv = self.stellar_map.ebv[i,j])
+                if self.gas_map is not None:
+                    gg = HII_Region(self.config, hii_tem, halpha = self.gas_map.halpha[i,j], 
+                                    logz = self.gas_map.zh[i,j], vel = self.gas_map.vel[i,j], 
+                                    vdisp = self.gas_map.vdisp[i,j], ebv = self.gas_map.ebv[i,j])
+                    self.flux[i,j,:] = ss.flux + gg.flux
+                else:
+                    self.flux[i,j,:] = ss.flux
+                    
+    def wcs_info(self):
+        wcs = fits.Header()
+        wcs_dict = {'CTYPE1': 'WAVE    ', 
+                    'CUNIT1': 'Angstrom', 
+                    'CDELT1': self.config.dlam, 
+                    'CRPIX1': 1, 
+                    'CRVAL1': np.min(self.wave), 
+                    'CTYPE2': 'RA---TAN', 
+                    'CUNIT2': 'deg', 
+                    'CDELT2': self.dpix / 3600., 
+                    'CRPIX2': np.round(self.ny / 2.), 
+                    'CRVAL2': 0.5, 
+                    'CTYPE3': 'DEC--TAN', 
+                    'CUNIT3': 'deg', 
+                    'CDELT3': self.dpix / 3600., 
+                    'CRPIX3': np.round(self.nx / 2.), 
+                    'CRVAL3': 1, 
+                    'BUNIT' : '10**(-17)*erg/s/cm**2/Angstrom'}
+        input_wcs = astropy.wcs.WCS(wcs_dict)
+        self.wcs_header = input_wcs.to_header()
+        
+    def insert_spec(self, spec, dx = 0, dy = 0):
+        x0 = np.int(np.round(self.config.nx / 2.))
+        y0 = np.int(np.round(self.config.ny / 2.))
+        self.flux[x0 + dx, y0 + dy, :] = self.flux[x0 + dx, y0 + dy, :] + spec.flux
+        
+    def savefits(self, filename, path = './'): 
+        
+        hdr = fits.Header()
+        
+        hdr['FILETYPE'] = 'SCICUBE'
+        hdr['CODE']     = 'CSST-IFS-GEHONG'
+        hdr['VERSION']  = '0.0.1'
+        
+        hdr['OBJECT']   = 'NGC1234'
+        hdr['RA']       = 0.0
+        hdr['DEC']      = 0.0
+        
+        hdu0 = fits.PrimaryHDU(header = hdr)
+        
+        self.wcs_info()
+        hdr = self.wcs_header
+        hdu1 = fits.ImageHDU(self.flux, header = hdr)
+        
+        # Output
+        hdulist = fits.HDUList([hdu0, hdu1])
+        hdulist.writeto(path + filename, overwrite = True)
\ No newline at end of file

src/gehong/map2d.py

+from __future__ import division
+
+import scipy.special as sp
+import numpy as np
+from astropy.io import fits
+from skimage.transform import resize
+
+def Sersic2D(x, y, mag = 12.0, r_eff = 1.0, n = 2.0, ellip = 0.5, 
+             theta = 0.0, x_0 = 0.0, y_0 = 0.0, pixelscale = 0.01):
+    """
+    Model of 2D - Sersic profile
+
+    Parameters
+    ----------
+    x : float array
+        _description_
+    y : _type_
+        _description_
+    mag : float, optional
+        Integral magnitude of sersic model, by default 12.0
+    r_eff : float, optional
+        Effective radius in pixel, by default 1.0
+    n : float, optional
+        Sersic index, by default 2.0
+    ellip : float, optional
+        Ellipticity, by default 0.5
+    theta : float, optional
+        Position angle in degree, by default 0.0
+    x_0 : float, optional
+        Offset of the center of Sersic model, by default 0.0
+    y_0 : float, optional
+        Offset of the center of Sersic model, by default 0.0
+    pixelscale : float, optional
+        Size of each pixel in arcsec^2, by default 0.01
+
+    Returns
+    -------
+    
+        _description_
+    """
+    # Produce Sersic profile
+    bn = sp.gammaincinv(2. * n, 0.5)
+    a, b = r_eff, (1 - ellip) * r_eff
+    cos_theta, sin_theta = np.cos(theta), np.sin(theta)
+    x_maj = (x - x_0) * cos_theta + (y - y_0) * sin_theta
+    x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta
+    z = (abs(x_maj / a) ** 2 + abs(x_min / b)  ** 2) ** (1 / 2)
+    profile = np.exp(-bn * (z ** (1 / n) - 1))
+    
+    # Normalization
+    integral = a * b * 2 * np.pi * n * np.exp(bn) / (bn ** (2 * n)) * sp.gamma(2 * n)
+    prof_norm = profile / integral * pixelscale
+    
+    # Calibration
+    total_flux = 10. ** ((22.5 - mag) * 0.4)
+    sb_mag = 22.5 - 2.5 * np.log10(prof_norm * total_flux / pixelscale)
+    
+    return sb_mag
+
+def VelMap2D(x, y, vmax = 200.0, rt = 1.0, ellip = 0.5, 
+             theta = 0.0, x_0 = 0.0, y_0 = 0.0):
+    """
+    VelMap2D _summary_
+
+    Parameters
+    ----------
+    x : _type_
+        _description_
+    y : _type_
+        _description_
+    vmax : int, optional
+        _description_, by default 200
+    rt : int, optional
+        _description_, by default 1
+    ellip : float, optional
+        _description_, by default 0.5
+    theta : int, optional
+        _description_, by default 0
+    x_0 : int, optional
+        _description_, by default 0
+    y_0 : int, optional
+        _description_, by default 0
+
+    Returns
+    -------
+    _type_
+        _description_
+    """
+    # Produce tanh profile
+    a, b = rt, (1 - ellip) * rt
+    cos_theta, sin_theta = np.cos(theta), np.sin(theta)
+    x_maj = (x - x_0) * cos_theta + (y - y_0) * sin_theta
+    x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta
+    z = (abs(x_maj / a) ** 2 + abs(x_min / b)  ** 2) ** (1 / 2)
+    profile = vmax * np.tanh(z) * ((x_maj / a) / z)
+    
+    return profile
+
+def GradMap2D(x, y, a0 = 10, r_eff = 1, gred = -1, ellip = 0.5, 
+              theta = 0, x_0 = 0, y_0 = 0):
+    """
+    GradMap2D _summary_
+
+    Parameters
+    ----------
+    x : _type_
+        _description_
+    y : _type_
+        _description_
+    a0 : int, optional
+        _description_, by default 10
+    r_eff : int, optional
+        _description_, by default 1
+    gred : int, optional
+        _description_, by default -1
+    ellip : float, optional
+        _description_, by default 0.5
+    theta : int, optional
+        _description_, by default 0
+    x_0 : int, optional
+        _description_, by default 0
+    y_0 : int, optional
+        _description_, by default 0
+
+    Returns
+    -------
+    _type_
+        _description_
+    """
+    # Produce gradiant profile
+    a, b = r_eff, (1 - ellip) * r_eff
+    cos_theta, sin_theta = np.cos(theta), np.sin(theta)
+    x_maj = (x - x_0) * cos_theta + (y - y_0) * sin_theta
+    x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta
+    z = (abs(x_maj / a) ** 2 + abs(x_min / b)  ** 2) ** (1 / 2)
+    profile = a0 + z * gred
+    
+    return profile
+
+class Map2d(object):
+
+    def __init__(self, config):
+        """
+        __init__ _summary_
+
+        Parameters
+        ----------
+        inst : _type_
+            _description_
+        """
+        self.xsamp = config.dpix
+        self.ysamp = config.dpix
+        startx = -(config.nx - 1) / 2.0 * self.xsamp
+        stopx  = (config.nx - 1) / 2.0 * self.xsamp
+        starty = -(config.ny - 1) / 2.0 * self.ysamp
+        stopy  = (config.ny - 1) / 2.0 * self.ysamp
+        xvals  = np.linspace(startx, stopx, num = config.nx)
+        yvals  = np.linspace(starty, stopy, num = config.ny)
+
+        ones = np.ones((config.ny, config.nx))
+        x = ones * xvals
+        y = np.flipud(ones * yvals.reshape(int(config.ny), 1))
+
+        self.nx = config.nx
+        self.ny = config.ny
+        self.x  = x
+        self.y  = y
+        self.row = xvals
+        # flip Y axis because we use Y increasing from bottom to top
+        self.col = yvals[::-1]
+
+    def shift_rotate(self, yoff, xoff, rot):
+        """
+        Return shifted/rotated (y, x) given offsets (yoff, xoff) and rotation, rot (degrees)
+
+        Parameters
+        ----------
+        yoff, xoff: float
+            yoff, xoff offsets in world coordinates
+        rot: float
+            rotation angle in degrees
+
+        Returns
+        -------
+        ysh_rot, xsh_rot: 2D numpy arrays
+            rotated and shifted copies of Grid.x and Grid.y
+        """
+        pa_radians = np.pi * rot / 180.0
+        xsh = self.x - xoff
+        ysh = self.y - yoff
+        xsh_rot = xsh * np.cos(pa_radians) + ysh * np.sin(pa_radians)
+        ysh_rot = -xsh * np.sin(pa_radians) + ysh * np.cos(pa_radians)
+        return ysh_rot, xsh_rot
+    
+    def sersic_map(self, mag = 12.0, r_eff = 2.0, n = 2.5, ellip = 0.5, theta = -50.0):
+        """
+        Generate 2D map of Sersic model
+
+        Parameters
+        ----------
+        mag : float, optional
+            Integral magnitude, by default 12.0
+        r_eff : float, optional
+            Effective radius in arcsec, by default 2.0
+        n : float, optional
+            Sersic index, by default 2.5
+        ellip : float, optional
+            Ellipcity, by default 0.5
+        theta : float, optional
+            Position angle in degree, by default -50.0
+        """
+        self.mag   = mag
+        self.reff  = r_eff / self.xsamp
+        self.n     = n
+        self.ellip = ellip
+        self.theta = theta
+        self.map = Sersic2D(self.x, self.y, mag = self.mag, 
+                            r_eff = self.reff, n = self.n, 
+                            ellip = self.ellip, theta = self.theta, 
+                            pixelscale = self.xsamp * self.ysamp)
+        
+    def tanh_map(self, vmax = 200.0, rt = 2.0, ellip = 0.5, theta = -50.0):
+        """
+        Generate 2D velocity map of rotating disk according to tanh rotation curve
+
+        Parameters
+        ----------
+        vmax : float, optional
+            Maximum rotational velocity, by default 200.0km/s
+        rt : float, optional
+            Turn-over radius of rotation curve, by default 2.0 arcsec
+        ellip : float, optional
+            Apparent ellipcity of rotating disk, by default 0.5
+        theta : float, optional
+            Position angle of rotating disk, by default -50.0
+        """
+        self.vmax  = vmax
+        self.rt    = rt / self.xsamp
+        self.ellip = ellip
+        self.theta = theta
+        self.map = VelMap2D(self.x, self.y, vmax = self.vmax, rt = self.rt, 
+                            ellip = self.ellip, theta = self.theta)
+    
+    def gred_map(self, a0 = 10, r_eff = 1, gred = -1, ellip = 0.5, theta = 0):
+        """
+        Generate 2D maps according to the radial gradient form
+
+        Parameters
+        ----------
+        a0 : float, optional
+            Amplitude at the central pixel, by default 10
+        r_eff : float, optional
+            Effective radius, by default 1
+        gred : float, optional
+            Gradient of radial profile, by default -1
+        ellip : float, optional
+            Ellipcity, by default 0.5
+        theta : int, optional
+            Position angle, by default 0
+        """
+        self.a0    = a0
+        self.reff  = r_eff / self.xsamp
+        self.gred  = gred
+        self.ellip = ellip
+        self.theta = theta
+        self.map = GradMap2D(self.x, self.y, a0 = self.a0, r_eff = self.reff, 
+                             gred = self.gred, ellip = self.ellip, theta = self.theta)
+        
+    def load_map(self, image):
+        """
+        Generate 2D map according to input image
+
+        Parameters
+        ----------
+        image : 2d numpy array
+            The 2d array to be loaded.
+        """
+        if np.ndim(image) == 2:
+            self.map = resize(image, (self.nx, self.ny))
+            
+class StellarPopulationMap():
+    """
+    Class of 2D maps for the parameters of stellar population, such as 
+    surface brightness, median age and metallicity of stellar population, 
+    velocity and velocity dispersion maps, and dust extinction.
+    
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    sbright : class, optional
+        Class of the map of surface brightness of stellar population, by default None
+    logage : class, optional
+        Class of the map of stellar age, by default None
+    feh : class, optional
+        Class of the map of stellar metellicity, by default None
+    vel : class, optional
+        Class of the map of stellar velocity, by default None
+    vdisp : class, optional
+        Class of the map of stellar velocity dispersion, by default None
+    ebv : class, optional
+        Class of the map of dust extinction, by default None
+    """
+    def __init__(self, config, sbright = None, logage = None, 
+                 feh = None, vel = None, vdisp = None, ebv = None):
+
+        self.nx    = config.nx
+        self.ny    = config.ny
+        self.dpix  = config.dpix
+        self.fov_x = config.fov_x
+        self.fov_y = config.fov_y
+        
+        self.sbright = sbright.map
+        self.logage  = logage.map
+        self.feh     = feh.map
+        self.vel     = vel.map
+        self.vdisp   = vdisp.map
+        self.ebv     = ebv.map
+        
+        self.mag = self.sbright - 2.5 * np.log10(self.dpix * self.dpix)
+        self.age = 10 ** self.logage / 1e9
+        
+        self.vdisp[self.vdisp < 10] = 10
+        self.ebv[self.ebv < 0]      = 0
+        
+class IonizedGasMap():
+    """
+    Class of 2D maps for the parameters of ionized gas, such as 
+    Halpha flux map, gas-phase metallicity map, 
+    velocity and velocity dispersion maps, and dust extinction.
+
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    halpha : class, optional
+        Class of the map of Halpha flux, by default None
+    zh : class, optional
+        Class of the map of gas-phase metallicity, by default None
+    vel : class, optional
+        Class of the map of gas velocity, by default None
+    vdisp : class, optional
+        Class of the map of gas velocity dispersion, by default None
+    ebv : class, optional
+        Class of the map of dust extinction, by default None
+    """
+    def __init__(self, config, halpha = None, zh = None, vel = None, vdisp = None, ebv = None):
+        
+        self.nx     = config.nx
+        self.ny     = config.ny
+        self.dpix   = config.dpix
+        self.fov_x  = config.fov_x
+        self.fov_y  = config.fov_y
+        
+        self.halpha = halpha.map
+        self.zh     = zh.map
+        self.vel    = vel.map
+        self.vdisp  = vdisp.map
+        self.ebv    = ebv.map
+        
+        self.vdisp[self.vdisp < 10] = 10
+        self.ebv[self.ebv < 0]      = 0
\ No newline at end of file

src/gehong/spec1d.py

+import os
+import glob
+import sys
+from os import path
+import numpy as np
+import astropy.units as u
+from astropy.io import fits
+from scipy.stats import norm
+from scipy.interpolate import interp1d
+
+data_path = os.getenv('GEHONG_DATA_PATH')
+
+def readcol(filename, **kwargs):
+    """
+    readcol, taken from ppxf.
+
+    Parameters
+    ----------
+    filename : string
+        The name of input ascii file
+
+    Returns
+    -------
+    float
+        The value of each columns
+    """
+    f = np.genfromtxt(filename, dtype=None, **kwargs)
+
+    t = type(f[0])
+    if t == np.ndarray or t == np.void: # array or structured array
+        f = map(np.array, zip(*f))
+
+    # In Python 3.x all strings (e.g. name='NGC1023') are Unicode strings by defauls.
+    # However genfromtxt() returns byte strings b'NGC1023' for non-numeric columns.
+    # To have the same behaviour in Python 3 as in Python 2, I convert the Numpy
+    # byte string 'S' type into Unicode strings, which behaves like normal strings.
+    # With this change I can read the string a='NGC1023' from a text file and the
+    # test a == 'NGC1023' will give True as expected.
+
+    if sys.version >= '3':
+        f = [v.astype(str) if v.dtype.char=='S' else v for v in f]
+
+    return f
+
+def log_rebin(lamRange, spec, oversample=False, velscale=None, flux=False):
+    """
+    Logarithmically rebin a spectrum, while rigorously conserving the flux. 
+    This function is taken from ppxf. 
+
+    Parameters
+    ----------
+    lamRange : array
+        Two elements vector containing the central wavelength
+        of the first and last pixels in the spectrum
+    spec : array
+        Input spectrum
+    oversample : bool, optional
+        Oversampling can be done, not to loose spectral resolution,
+        especally for extended wavelength ranges and to avoid aliasing, by default False
+    velscale : float, optional
+        velocity scale in km/s per pixels, by default None
+    flux : bool, optional
+        True to preserve total flux, by default False
+
+    Returns
+    -------
+    specNew : array
+        Output spectrum
+    logLam : array
+        Wavelength array in logarithm
+    velscale : array
+        velocity scale in km/s per pixels
+    """
+    lamRange = np.asarray(lamRange)
+    assert len(lamRange) == 2, 'lamRange must contain two elements'
+    assert lamRange[0] < lamRange[1], 'It must be lamRange[0] < lamRange[1]'
+    s = spec.shape
+    assert len(s) == 1, 'input spectrum must be a vector'
+    n = s[0]
+    if oversample:
+        m = int(n*oversample)
+    else:
+        m = int(n)
+
+    dLam = np.diff(lamRange)/(n - 1.)        # Assume constant dLam
+    lim = lamRange/dLam + [-0.5, 0.5]        # All in units of dLam
+    borders = np.linspace(*lim, num=n+1)     # Linearly
+    logLim = np.log(lim)
+
+    c = 299792.458                           # Speed of light in km/s
+    if velscale is None:                     # Velocity scale is set by user
+        velscale = np.diff(logLim)/m*c       # Only for output
+    else:
+        logScale = velscale/c
+        m = int(np.diff(logLim)/logScale)    # Number of output pixels
+        logLim[1] = logLim[0] + m*logScale
+
+    newBorders = np.exp(np.linspace(*logLim, num=m+1)) # Logarithmically
+    k = (newBorders - lim[0]).clip(0, n-1).astype(int)
+
+    specNew = np.add.reduceat(spec, k)[:-1]  # Do analytic integral
+    specNew *= np.diff(k) > 0    # fix for design flaw of reduceat()
+    specNew += np.diff((newBorders - borders[k])*spec[k])
+
+    if not flux:
+        specNew /= np.diff(newBorders)
+
+    # Output log(wavelength): log of geometric mean
+    logLam = np.log(np.sqrt(newBorders[1:]*newBorders[:-1])*dLam)
+
+    return specNew, logLam, velscale
+
+def gaussian_filter1d(spec, sig):
+    """
+    One-dimensional Gaussian convolution
+
+    Parameters
+    ----------
+    spec : float array
+        vector with the spectrum to convolve
+    sig : float
+        vector of sigma values (in pixels) for every pixel
+    Returns
+    -------
+    float array
+        Spectrum after convolution
+    """
+    sig = sig.clip(0.01)  # forces zero sigmas to have 0.01 pixels
+    p = int(np.ceil(np.max(3*sig)))
+    m = 2*p + 1  # kernel size
+    x2 = np.linspace(-p, p, m)**2
+
+    n = spec.size
+    a = np.zeros((m, n))
+    for j in range(m):   # Loop over the small size of the kernel
+        a[j, p:-p] = spec[j:n-m+j+1]
+
+    gau = np.exp(-x2[:, None]/(2*sig**2))
+    gau /= np.sum(gau, 0)[None, :]  # Normalize kernel
+
+    conv_spectrum = np.sum(a*gau, 0)
+
+    return conv_spectrum
+
+def calibrate(wave, flux, mag, filtername = 'SLOAN_SDSS.r'):
+    """
+    Flux calibration of spectrum
+
+    Parameters
+    ----------
+    wave : float array
+        Wavelength of input spectrum
+    flux : float array
+        Flux of input spectrum
+    mag : float
+        Magnitude used for flux calibration
+    filtername : str, optional
+        Filter band name, by default 'SLOAN_SDSS.r'
+
+    Returns
+    -------
+    float array
+        Spectrum after flux calibration
+    """
+    # Loading response curve
+    if filtername == '5100':
+        wave0 = np.linspace(3000,10000,7000)
+        response0 = np.zeros(7000)
+        response0[(wave0 > 5050) & (wave0 < 5150)] = 1.
+    else:
+        filter_file = data_path + '/data/filter/' + filtername+'.filter'
+        wave0, response0 = readcol(filter_file)
+    
+    # Setting the response
+    func = interp1d(wave0, response0)
+    response = np.copy(wave)
+    ind_extra = (wave > max(wave0)) | (wave < min(wave0))
+    response[ind_extra] = 0
+    ind_inside = (wave < max(wave0)) & (wave > min(wave0))
+    response[ind_inside] = func(wave[ind_inside])
+        
+    # Flux map of datacube for given filter band
+    preflux = np.sum(flux * response * np.mean(np.diff(wave))) / np.sum(response * np.mean(np.diff(wave)))
+    
+    # Real flux from magnitude for given filter
+    realflux = (mag * u.STmag).to(u.erg/u.s/u.cm**2/u.AA).value
+
+    # Normalization
+    flux_ratio = realflux / preflux
+    flux_calibrate = flux * flux_ratio * 1e17                        # Units: 10^-17 erg/s/A/cm^2
+    
+    return flux_calibrate
+
+# ----------------
+# Reddening Module
+
+def Calzetti_Law(wave, Rv = 4.05):
+    """
+    Dust Extinction Curve of Calzetti et al. (2000)
+
+    Parameters
+    ----------
+    wave : float, or float array
+        Wavelength
+    Rv : float, optional
+        Extinction coefficient, by default 4.05
+
+    Returns
+    -------
+    float
+        Extinction value corresponding to the input wavelength
+    """
+    wave_number = 1./(wave * 1e-4)
+    reddening_curve = np.zeros(len(wave))
+    
+    idx = (wave >= 1200) & (wave < 6300)
+    reddening_curve[idx] = 2.659 * ( -2.156 + 1.509 * wave_number[idx] - 0.198 * \
+                    (wave_number[idx] ** 2)) + 0.011 * (wave_number[idx] **3 ) + Rv
+                                    
+    idx = (wave >= 6300) & (wave <= 22000)
+    reddening_curve[idx] = 2.659 * ( -1.857 + 1.040 * wave_number[idx]) + Rv
+    return reddening_curve
+
+def reddening(wave, flux, ebv = 0.0, law = 'calzetti', Rv = 4.05):
+    """
+    Reddening an input spectra through a given reddening curve.
+
+    Parameters
+    ----------
+    wave : float array
+        Wavelength of input spectra
+    flux : float array
+        Flux of input spectra
+    ebv : float, optional
+        E(B-V) value, by default 0
+    law : str, optional
+        Extinction curve, by default 'calzetti'
+    Rv : float, optional
+        Extinction coefficient, by default 4.05
+
+    Returns
+    -------
+    float array
+        Flux of spectra after reddening
+    """
+    curve = Calzetti_Law(wave, Rv = Rv)
+    fluxNew = flux / (10. ** (0.4 * ebv * curve))
+    return fluxNew
+
+################
+# Emission Lines
+################
+
+def SingleEmissinoLine(wave, line_wave, FWHM_inst):
+    """
+    Profile of single emission line (Gaussian profile)
+
+    Parameters
+    ----------
+    wave : float array
+        Wavelength of spectrum
+    line_wave : float
+        Wavelength of emission line at the line center
+    FWHM_inst : float
+        Intrinsic broadening of emission line, units: A
+
+    Returns
+    -------
+    float array
+        Spectra of single emission line
+    """
+    sigma = FWHM_inst / 2.355 
+    flux = norm.pdf(wave, line_wave, sigma)
+    return flux
+
+class EmissionLineTemplate():
+    """
+    Template for the emission lines
+
+    Parameters
+    ----------
+    config : class
+        The class of configuration.
+    lam_range : list, optional
+        Wavelength range, by default [500, 15000]
+    dlam : float, optional
+        Wavelength width per pixel, by default 0.1A
+    model : str, optional
+        Emission line model, including 'hii' for HII region and 'nlr' for narrow line region of AGN, 
+        by default 'hii'
+    """
+    
+    def __init__(self, config, lam_range = [500, 15000], dlam = 0.1, model = 'hii'):
+        
+        self.lam_range = lam_range
+        self.wave = np.arange(lam_range[0], lam_range[1], 0.1)
+        self.FWHM_inst = config.inst_fwhm
+        self.model = model
+        
+        # HII region model of fsps-cloudy
+        if model == 'hii':
+            
+            # Loading emission line flux table
+            flux_table_file = data_path + '/data/fsps.nebular.fits'
+            line_table = fits.open(flux_table_file)
+
+            # Emission line list
+            line_list  = line_table[1].data
+            line_wave  = line_list['Wave']
+            line_names = line_list['Name']
+        
+            w = (line_wave > lam_range[0]) & (line_wave < lam_range[1])
+            self.line_names = line_names[w]
+            self.line_wave = line_wave[w]
+            
+            # Make parameter grid
+            grid = line_table[2].data
+            self.logz_grid   = grid['logZ']
+        
+        # Narrow line region model of 
+        if model == 'nlr':
+            
+            # Loading emission line flux table
+            flux_table_file = data_path + '/data/AGN.NLR.fits'
+            line_table = fits.open(flux_table_file)
+
+            # Emission line list
+            line_list  = line_table[1].data
+            line_wave  = line_list['Wave']
+            line_names = line_list['Name']
+        
+            w = (line_wave > lam_range[0]) & (line_wave < lam_range[1])
+            self.line_names = line_names[w]
+            self.line_wave = line_wave[w]
+            
+            # Make parameter grid
+            grid = line_table[2].data
+            self.logz_grid   = grid['logZ']
+        
+        # Flux ratio
+        flux_ratio = line_table[3].data
+        self.flux_ratio = flux_ratio
+        
+        # Make emission line
+        nline = len(line_wave)
+        for i in range(nline):
+            if i==0:
+                emission_line  = SingleEmissinoLine(self.wave, line_wave[i], self.FWHM_inst)
+                emission_lines = emission_line
+            else:
+                emission_line  = SingleEmissinoLine(self.wave, line_wave[i], self.FWHM_inst)
+                emission_lines = np.vstack((emission_lines, emission_line))
+        self.emission_lines = emission_lines.T
+
+class HII_Region():
+
+    """
+    Class for the spectra of HII region
+
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    temp : class
+        Class of emission line template
+    halpha : float, optional
+        Integral flux of Halpha emission line, by default 100 * 1e-17 erg/s/cm^2
+    logz : float, optional
+        Gas-phase metallicity, by default 0.0
+    vel : float, optional
+        Line of sight velocity, by default 100.0km/s
+    vdisp : float, optional
+        Velocity dispersion, by default 120.0km/s
+    ebv : float, optional
+        Dust extinction, by default 0.1
+
+    Raises
+    ------
+    ValueError
+        The value of logZ should be between -2 and 0.5.
+    """
+    
+    def __init__(self, config, temp, halpha = 100.0, logz = 0.0,
+                 vel = 100.0, vdisp = 120.0, ebv = 0.1):
+
+        if (logz > -2) & (logz < 0.5):
+            indz = np.argmin(np.abs(logz - temp.logz_grid))
+            flux_ratio = temp.flux_ratio[indz, :]
+        else:
+            raise ValueError('The value of logZ is not in the range of [-2, 0.5]!')
+        
+        # Make emission line spectra through adding emission lines                 
+        emlines = temp.emission_lines * flux_ratio
+        flux_combine = np.sum(emlines, axis = 1)
+        flux_calibrate = flux_combine * halpha      # Units: erg/s/A/cm^2
+        
+        # Dust attenuation
+        if np.isscalar(ebv):
+            flux_dust = reddening(temp.wave, flux_calibrate, ebv = ebv)
+            
+        # Broadening caused by Velocity Dispersion
+        velscale = 10
+        lam_range = [np.min(temp.wave), np.max(temp.wave)]
+        flux_logwave, logLam = log_rebin(lam_range, flux_dust, velscale=velscale)[:2]
+        
+        sigma_gas = vdisp / velscale                                # in pixel
+        sigma_LSF = temp.FWHM_inst / (np.exp(logLam)) * 3e5 / velscale          # in pixel
+        
+        if sigma_gas>0: 
+            sigma_dif = np.zeros(len(flux_logwave))
+            idx = (sigma_gas > sigma_LSF)
+            sigma_dif[idx] = np.sqrt(sigma_gas ** 2. - sigma_LSF[idx] ** 2.)
+            idx = (sigma_gas <= sigma_LSF)
+            sigma_dif[idx] = 0.1
+            flux_broad = gaussian_filter1d(flux_logwave, sigma_dif)
+        else:
+            flux_broad = flux_logwave
+            
+        # Redshift
+        redshift = vel / 3e5
+        wave_r = np.exp(logLam) * (1 + redshift)
+        flux_red = np.interp(config.wave, wave_r, flux_broad)
+            
+        self.wave = config.wave
+        self.flux = flux_red
+
+#############
+# AGN Spectra
+#############
+
+class AGN_NLR():
+
+    """
+    Class for narrow line region of AGN
+
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    temp : class
+        Class of emission line template
+    halpha : float, optional
+        Integral flux of Halpha emission line, by default 100 * 1e-17 erg/s/cm^2
+    logz : float, optional
+        Gas-phase metallicity, by default 0.0
+    vel : float, optional
+        Line of sight velocity, by default 100.0km/s
+    vdisp : float, optional
+        Velocity dispersion, by default 120.0km/s
+    ebv : float, optional
+        Dust extinction, by default 0.1
+
+    Raises
+    ------
+    ValueError
+        The value of logZ should be between -2 and 0.5.
+    """
+    
+    def __init__(self, config, temp, halpha = 100.0, logz = 0.0,
+                 vel = 100.0, vdisp = 120.0, ebv = 0.1):
+        
+        if (logz > -2) & (logz < 0.5):
+            indz = np.argmin(np.abs(logz - temp.logz_grid))
+            flux_ratio = temp.flux_ratio[indz, :]
+        else:
+            raise ValueError('The value of logZ is not in the range of [-2, 0.5]!')
+        
+        # Make emission line spectra through adding emission lines                 
+        emlines = temp.emission_lines * (flux_ratio / flux_ratio[6] )
+        flux_combine = np.sum(emlines, axis = 1)
+        flux_calibrate = flux_combine * halpha      # Units: 1e-17 erg/s/A/cm^2
+        
+        # Dust attenuation
+        if np.isscalar(ebv):
+            flux_dust = reddening(temp.wave, flux_calibrate, ebv = ebv)
+            
+        # Broadening caused by Velocity Dispersion
+        velscale = 10
+        lam_range = [np.min(temp.wave), np.max(temp.wave)]
+        flux_logwave, logLam = log_rebin(lam_range, flux_dust, velscale=velscale)[:2]
+        
+        sigma_gas = vdisp / velscale                                        # in pixel
+        sigma_LSF = temp.FWHM_inst / (np.exp(logLam)) * 3e5 / velscale      # in pixel
+        
+        if sigma_gas>0: 
+            sigma_dif = np.zeros(len(flux_logwave))
+            idx = (sigma_gas > sigma_LSF)
+            sigma_dif[idx] = np.sqrt(sigma_gas ** 2. - sigma_LSF[idx] ** 2.)
+            idx = (sigma_gas <= sigma_LSF)
+            sigma_dif[idx] = 0.1
+            flux_broad = gaussian_filter1d(flux_logwave, sigma_dif)
+            
+        # Redshift
+        redshift = vel / 3e5
+        wave_r = np.exp(logLam) * (1 + redshift)
+        flux_red = np.interp(config.wave, wave_r, flux_broad)
+        
+        self.wave = config.wave
+        self.flux = flux_red
+
+class AGN_BLR():
+
+    """
+    Class for the broad line region of AGN
+
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    hbeta_flux : float, optional
+        Integral flux of Hbeta broad line, by default 100 * 1e-17 erg/s/cm^2
+    hbeta_fwhm : float, optional
+        FWHM of Hbeta broad line, by default 2000.0km/s
+    vel : float, optional
+        Line of sight velocity, by default 100.0km/s
+    ebv : float, optional
+        Dust extinction, by default 0.1
+    lam_range : list, optional
+        Wavelength range, by default [500, 15000]
+    """
+
+    def __init__(self, config, hbeta_flux = 100.0, hbeta_fwhm = 2000.0, ebv = 0.1, 
+                 vel = 0., lam_range = [500, 15000]):
+
+        wave_rest = np.arange(lam_range[0], lam_range[1], 0.1)
+        
+        line_names = ['Hepsilon', 'Hdelta', 'Hgamma', 'Hbeta', 'Halpha']
+        line_waves = [3970.079, 4101.742, 4340.471, 4861.333, 6562.819]
+        line_ratio = [0.101, 0.208, 0.405, 1.000, 2.579]              # From Ilic et al. (2006)
+        
+        # Make emission lines
+        for i in range(len(line_names)):
+            if i==0:
+                emission_line  = SingleEmissinoLine(wave_rest, line_waves[i], 
+                                                    hbeta_fwhm / 3e5 * line_waves[i])
+                emission_lines = emission_line
+            else:
+                emission_line  = SingleEmissinoLine(wave_rest, line_waves[i], 
+                                                    hbeta_fwhm / 3e5 * line_waves[i])
+                emission_lines = np.vstack((emission_lines, emission_line))
+        emlines = emission_lines.T * line_ratio
+        flux_combine = np.sum(emlines, axis = 1)
+        
+        # Flux callibration
+        flux_calibrate = flux_combine * hbeta_flux      # Units: 1e-17 erg/s/A/cm^2
+        
+        # Dust attenuation
+        if np.isscalar(ebv):
+            flux_dust = reddening(wave_rest, flux_calibrate, ebv = ebv)
+        else:
+            flux_dust = flux_calibrate
+            
+        # Redshift
+        redshift = vel / 3e5
+        wave_r = wave_rest * (1 + redshift)
+        flux_red = np.interp(config.wave, wave_r, flux_dust)
+            
+        self.wave = config.wave
+        self.flux = flux_red
+
+class AGN_FeII():
+    """
+    Class for FeII emission lines of AGN
+
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    hbeta_broad : float, optional
+        Integral flux of Hbeta broad line, by default 100 * 1e-17 erg/s/cm^2
+    r4570 : float, optional
+        Flux ratio between Fe4570 flux and Hbeta broad line flux, by default 0.4
+    vel : float, optional
+        Line of sight velocity, by default 100.0km/s
+    ebv : float, optional
+        Dust extinction, by default 0.1
+    """
+    
+    def __init__(self, config, hbeta_broad = 100.0, r4570 = 0.4, ebv = 0.1, vel = 100.0):
+
+        filename = data_path + '/data/FeII.AGN.fits'
+        
+        # Loading FeII template
+        hdulist = fits.open(filename)
+        data = hdulist[1].data
+        wave_rest  = data['WAVE']
+        flux_model = data['FLUX']
+        
+        # Determine the flux of FeII
+        Fe4570_temp  = 100
+        Fe4570_model = hbeta_broad * r4570
+        Ratio_Fe4570 = Fe4570_model / Fe4570_temp
+        
+        # Flux calibration
+        flux_calibrate = flux_model * Ratio_Fe4570
+        
+        # Dust attenuation
+        if np.isscalar(ebv):
+            flux_dust = reddening(wave_rest, flux_calibrate, ebv = ebv)
+        else:
+            flux_dust = flux_calibrate
+             
+        # Redshift
+        redshift = vel / 3e5
+        wave_r   = wave_rest * (1 + redshift)
+        flux_red = np.interp(config.wave, wave_r, flux_dust)
+        
+        self.wave = config.wave
+        self.flux = flux_red
+        
+class AGN_Powerlaw():
+
+    """
+    The class of power-law spectrum of AGN
+
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    M5100 : float, optional
+        Magnitude of power law spectrum between 5050A and 5150A, by default 1000.0 * 1e-17 erg/s/cm^2
+    alpha : float, optional
+        Index of power law, by default -1.5
+    vel : float, optional
+        Line of sight velocity, by default 100.0km/s
+    Ebv : float, optional
+        Dust extinction, by default 0.1
+    """
+    
+    def __init__(self, config, m5100 = 1000.0, alpha = -1.5, vel = 100.0, ebv = 0.1):
+
+        wave_rest = np.linspace(1000,20000,10000)
+        flux      = wave_rest ** alpha
+        
+        # Flux calibration
+        flux_calibrate = calibrate(wave_rest, flux, m5100, filtername='5100')
+        
+        # Dust attenuation
+        if np.isscalar(ebv):
+            flux_dust = reddening(wave_rest, flux_calibrate, ebv = ebv)
+        else:
+            flux_dust = flux_calibrate
+            
+        # Redshift
+        redshift = vel / 3e5
+        wave_r   = wave_rest * (1 + redshift)
+        flux_red = np.interp(config.wave, wave_r, flux_dust)
+        
+        self.wave = config.wave
+        self.flux = flux_red
+        
+class AGN():
+
+    """
+    Class of singal spectra of AGN
+
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    nlr_template : class
+        Class of emission line template
+    bhmass : float, optional
+        Black hole mass used for calculating the luminosity of power law spectrum at 5100A, 
+        by default 1e6 solar mass
+    edd_ratio : float, optional
+        Eddinton ratio used for calculating the luminosity of power law spectrum at 5100A, by default 0.05
+    halpha_broad : float, optional
+        Integral flux of Halpha broad line, by default 100.0 * 1e-17 erg/s/cm^2
+    halpha_narrow : float, optional
+        Integral flux of Halpha narrow line, by default 100.0 * 1e-17 erg/s/cm^2
+    vdisp_broad : float, optional
+        Velocity dispersion of Halpha broad line, by default 5000.0km/s
+    vdisp_narrow : float, optional
+        Velocity dispersion of Halpha narrow line, by default 500.0km/s
+    vel : float, optional
+        Line of sight velocity, by default 1000.0km/s
+    logz : float, optional
+        Gas-phase metallicity of narrow line region, by default 0.0
+    ebv : float, optional
+        Dust extinction, by default 0.1
+    dist : float, optional
+        Luminosity distance of AGN, by default 20.0Mpc
+    """
+    def __init__(self, config, nlr_template, bhmass = 1e6, edd_ratio = 0.05, 
+                 halpha_broad = 100.0, halpha_narrow = 100.0, vdisp_broad = 5000.0, vdisp_narrow = 500.0, 
+                 vel = 1000.0, logz = 0.0, ebv = 0.1, dist = 20.0):
+        
+        NLR = AGN_NLR(config, nlr_template, halpha = halpha_narrow, logz = logz,
+                      vel = vel, vdisp = vdisp_narrow, ebv = ebv)
+        if halpha_broad > 0:
+            BLR = AGN_BLR(config, hbeta_flux = halpha_broad / 2.579, 
+                          hbeta_fwhm = vdisp_broad / 2.355, ebv = ebv, vel = vel)
+
+        m5100 = BHmass_to_M5100(bhmass, edd_ratio = edd_ratio, dist = dist)
+        PL  = AGN_Powerlaw(config, m5100 = m5100, ebv = ebv, vel = vel)
+        Fe  = AGN_FeII(config, hbeta_broad = halpha_broad / 2.579, ebv = ebv, vel = vel)
+        
+        self.wave = config.wave
+        self.flux = NLR.flux + PL.flux + Fe.flux
+        
+        if halpha_broad > 0:
+            self.flux = self.flux + BLR.flux
+    
+def BHmass_to_M5100(bhmass, edd_ratio = 0.05, dist = 21.0):
+    """
+    Caculate magnitude at 5100A according to the black hole mass
+
+    Parameters
+    ----------
+    bhmass : float
+        Black hole mass, unit: solar mass
+    edd_ratio : float, optional
+        Eddtington ratio, by default 0.05
+    dist : float, optional
+        Distance to the black hole, by default 21.0Mpc
+
+    Returns
+    -------
+    float
+        Magnitude at 5100A
+    """
+    
+    # Calculate bolometric luminosity
+    Ledd = 3e4 * bhmass
+    Lbol = Ledd * edd_ratio
+
+    # Convert bolometric luminosity to 5100A luminosity (Marconi et al. 2004)
+    L5100 = Lbol / 10.9
+    M5100 = 4.86 - 2.5 * np.log10(L5100)
+    m5100 = M5100 + 5. * np.log10(dist * 1e5)
+
+    return m5100
+
+#################
+# Stellar Spectra
+#################
+
+"""
+    Extract the age and metallicity from the name of a file of
+    the MILES library of Single Stellar Population models as
+    downloaded from http://miles.iac.es/ as of 2016
+
+    :param filename: string possibly including full path
+        (e.g. 'miles_library/Mun1.30Zm0.40T03.9811.fits')
+    :return: age (Gyr), [M/H]
+
+    """
+
+def age_metal(filename):
+    """
+    Extract the age and metallicity from the name of a file of
+    the MILES library of Single Stellar Population models.
+
+    Parameters
+    ----------
+    filename : string
+        Full path of template files
+
+    Returns
+    -------
+    age : float
+        Age of SSP (Gyr)
+    FeH : float
+        Metallicity of SSP
+
+    Raises
+    ------
+    ValueError
+        This is not a standard MILES filename
+    """
+    s = path.basename(filename)
+    age = float(s[s.find("T")+1:s.find("_iPp0.00_baseFe.fits")])
+    metal = s[s.find("Z")+1:s.find("T")]
+    if "m" in metal:
+        metal = -float(metal[1:])
+    elif "p" in metal:
+        metal = float(metal[1:])
+    else:
+        raise ValueError("This is not a standard MILES filename")
+
+    return age, metal
+
+class StellarContinuumTemplate(object):
+    """
+    Class of single stellar population template. 
+
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    velscale : array
+        velocity scale in km/s per pixels, by default 50.0km/s
+    pathname : string, optional
+        path with wildcards returning the list files to use, 
+        by default data_path+'/data/EMILES/Ech*_baseFe.fits'
+    normalize : bool, optional
+        Set to True to normalize each template to mean=1, by default False
+    """
+    def __init__(self, config, velscale = 50,
+                 pathname = data_path + '/data/EMILES/Ech*_baseFe.fits', 
+                 normalize = False):
+
+        FWHM_inst = config.inst_fwhm
+        
+        files = glob.glob(pathname)
+        assert len(files) > 0, "Files not found %s" % pathname
+
+        all = [age_metal(f) for f in files]
+        all_ages, all_metals = np.array(all).T
+        ages, metals = np.unique(all_ages), np.unique(all_metals)
+        n_ages, n_metal = len(ages), len(metals)
+
+        assert set(all) == set([(a, b) for a in ages for b in metals]), \
+            'Ages and Metals do not form a Cartesian grid'
+
+        # Extract the wavelength range and logarithmically rebin one spectrum
+        # to the same velocity scale of the SDSS galaxy spectrum, to determine
+        # the size needed for the array which will contain the template spectra.
+        hdu = fits.open(files[0])
+        ssp = hdu[0].data
+        h2 = hdu[0].header
+        lam_range_temp = h2['CRVAL1'] + np.array([0, h2['CDELT1']*(h2['NAXIS1']-1)])
+        sspNew, log_lam_temp = log_rebin(lam_range_temp, ssp, velscale=velscale)[:2]
+        #wave=((np.arange(hdr['NAXIS1'])+1.0)-hdr['CRPIX1'])*hdr['CDELT1']+hdr['CRVAL1']
+        
+        templates = np.empty((sspNew.size, n_ages, n_metal))
+        age_grid = np.empty((n_ages, n_metal))
+        metal_grid = np.empty((n_ages, n_metal))
+
+        # Convolve the whole Vazdekis library of spectral templates
+        # with the quadratic difference between the galaxy and the
+        # Vazdekis instrumental resolution. Logarithmically rebin
+        # and store each template as a column in the array TEMPLATES.
+
+        # Quadratic sigma difference in pixels Vazdekis --> galaxy
+        # The formula below is rigorously valid if the shapes of the
+        # instrumental spectral profiles are well approximated by Gaussians.
+        
+        # FWHM of Emiles templates
+        Emile_wave = np.exp(log_lam_temp)
+        Emile_FWHM = np.zeros(h2['NAXIS1'])
+        Emile_FWHM[np.where(Emile_wave < 3060)] = 3.
+        Emile_FWHM[np.where((Emile_wave >= 3060) & (Emile_wave < 3540))] = 3.
+        Emile_FWHM[np.where((Emile_wave >= 3540) & (Emile_wave < 8950))] = 2.5
+        Lwave = Emile_wave[np.where(Emile_wave >= 8950)]
+        Emile_FWHM[np.where(Emile_wave >= 8950)]=60*2.35/3.e5*Lwave  # sigma=60km/s at lambda > 8950
+        
+        LSF = Emile_FWHM
+        
+        FWHM_eff=Emile_FWHM.copy()   # combined FWHM from stellar library and instrument(input)
+        if np.isscalar(FWHM_inst):
+            FWHM_eff[Emile_FWHM < FWHM_inst] = FWHM_inst
+            LSF[Emile_FWHM < FWHM_inst] = FWHM_inst
+        else:
+            FWHM_eff[Emile_FWHM < FWHM_inst] = FWHM_inst[Emile_FWHM < FWHM_inst]
+            LSF[Emile_FWHM < FWHM_inst] = FWHM_inst[Emile_FWHM < FWHM_inst]
+        FWHM_dif  = np.sqrt(FWHM_eff**2 - Emile_FWHM**2)
+        sigma_dif = FWHM_dif/2.355/h2['CDELT1']   # Sigma difference in pixels
+
+        # Here we make sure the spectra are sorted in both [M/H] and Age
+        # along the two axes of the rectangular grid of templates.
+        for j, age in enumerate(ages):
+            for k, metal in enumerate(metals):
+                p = all.index((age, metal))
+                hdu = fits.open(files[p])
+                ssp = hdu[0].data
+                if np.isscalar(FWHM_dif):
+                    ssp = ndimage.gaussian_filter1d(ssp, sigma_dif)
+                else:
+                    ssp = gaussian_filter1d(ssp, sigma_dif)  # convolution with variable sigma
+                sspNew  = log_rebin(lam_range_temp, ssp, velscale=velscale)[0]
+                if normalize:
+                    sspNew /= np.mean(sspNew)
+                templates[:, j, k] = sspNew
+                age_grid[j, k]   = age
+                metal_grid[j, k] = metal
+
+        self.templates = templates/np.median(templates)  # Normalize by a scalar
+        self.log_lam_temp = log_lam_temp
+        self.age_grid = age_grid
+        self.metal_grid = metal_grid
+        self.n_ages = n_ages
+        self.n_metal = n_metal
+        self.LSF = log_rebin(lam_range_temp, LSF, velscale=velscale)[0]
+        self.velscale = velscale
+        
+    def fmass_ssp(self): 
+
+        isedpath = data_path + '/data/EMILES/model/'
+        massfile = isedpath + 'out_mass_CH_PADOVA00'
+        n_metal = self.n_metal
+        n_ages = self.n_ages
+        fage = self.age_grid[:,0]
+
+        fMs = np.zeros((n_ages,n_metal))
+        
+        Metal, Age, Ms = readcol(massfile, usecols=(2, 3, 6))
+        for i in range(n_metal):
+            for j in range(self.n_ages):
+                locmin = np.argmin(abs(Metal - self.metal_grid[j, i])) & np.argmin(abs(Age - self.age_grid[j, i]))
+                fMs[j,i] = Ms[locmin]
+
+        return fMs
+    
+class StellarContinuum():
+    """
+    The class of stellar continuum
+    
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    template : class
+        Class of single stellar population template
+    mag : float, optional
+        Magnitude in SDSS r-band, by default 15.0
+    age : float, optional
+        Median age of stellar continuum, by default 1.0Gyr
+    feh : float, optional
+        Metallicity of stellar continuum, by default 0.0
+    vel : float, optional
+        Line of sight velocity, by default 100.0km/s
+    vdisp : float, optional
+        Velocity dispersion, by default 120.0km/s
+    ebv : float, optional
+        Dust extinction, by default 0.1
+    """
+    def __init__(self, config, template, mag = 15.0, age = 1.0, feh = 0.0, 
+                 vel = 100.0, vdisp = 100.0, ebv = 0.1):
+
+        # -----------------
+        # Stellar Continuum
+        
+        SSP_temp = template.templates
+
+        # Select metal bins
+        metals = template.metal_grid[0,:]
+        minloc = np.argmin(abs(feh - metals))
+        tpls = SSP_temp[:, :, minloc]
+        fmass = template.fmass_ssp()[:, minloc]
+        
+        # Select age bins
+        Ages = template.age_grid[:,0]
+        minloc = np.argmin(abs(age-Ages))
+        Stellar = tpls[:, minloc]
+        
+        wave = np.exp(template.log_lam_temp)
+        
+        # Broadening caused by Velocity Dispersion
+        sigma_gal = vdisp / template.velscale                   # in pixel
+        sigma_LSF = template.LSF / template.velscale                 # in pixel
+        
+        if sigma_gal>0: 
+            sigma_dif = np.zeros(len(Stellar))
+            idx = (sigma_gal > sigma_LSF)
+            sigma_dif[idx] = np.sqrt(sigma_gal ** 2. - sigma_LSF[idx] ** 2.)
+            idx = (sigma_gal <= sigma_LSF)
+            sigma_dif[idx] = 0.1
+            flux0 = gaussian_filter1d(Stellar, sigma_dif)
+        
+        # Dust Reddening
+        if np.isscalar(ebv):
+            flux0 = reddening(wave, flux0, ebv = ebv)
+            
+        # Redshift
+        redshift = vel / 3e5
+        wave_r = wave * (1 + redshift)
+        
+        flux = np.interp(config.wave, wave_r, flux0)
+        
+        # Calibration
+        if np.isscalar(mag):
+            flux = calibrate(config.wave, flux, mag, filtername='SLOAN_SDSS.r')
+        
+        # Convert to input wavelength
+        self.wave = config.wave
+        self.flux = flux
+        
+#####################
+# Single Star Spectra
+#####################
+
+class SingleStarTemplate():
+    """
+    Class of single stellar template
+
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    velscale : float, option
+        velocity scale in km/s per pixels, by default 20.0km/s
+    """
+    def __init__(self, config, velscale = 20):
+
+        FWHM_inst = config.inst_fwhm
+        filename = data_path + '/data/Starlib.XSL.fits'
+        
+        hdulist = fits.open(filename)
+        lam  = hdulist[1].data['Wave']
+        flux = hdulist[2].data
+        par  = hdulist[3].data
+        
+        lam_range_temp = np.array([3500, 12000])
+        TemNew, log_lam_temp = log_rebin(lam_range_temp, flux[1, :], velscale = velscale)[:2]
+        
+        # FWHM of XLS templates
+        Temp_wave = np.exp(log_lam_temp)
+        Temp_FWHM = np.zeros(len(log_lam_temp))
+        Temp_FWHM[(Temp_wave < 5330)] = 13 * 2.35 / 3e5 * Temp_wave[(Temp_wave < 5330)]   # sigma = 13km/s at lambda <5330
+        Temp_FWHM[(Temp_wave >= 5330) & (Temp_wave < 9440)] = 11 * 2.35 / 3e5 * Temp_wave[(Temp_wave >= 5330) & (Temp_wave < 9440)]
+        # sigma = 13km/s at 5330 < lambda < 9440
+        Temp_FWHM[(Temp_wave >= 9440)] = 16 * 2.35 / 3e5 * Temp_wave[(Temp_wave >= 9440)] # sigma=16km/s at lambda > 9440
+        
+        LSF = Temp_FWHM
+        
+        FWHM_eff = Temp_FWHM.copy()   # combined FWHM from stellar library and instrument(input)
+        if np.isscalar(FWHM_inst):
+            FWHM_eff[Temp_FWHM < FWHM_inst] = FWHM_inst
+            LSF[Temp_FWHM < FWHM_inst]      = FWHM_inst
+        else:
+            FWHM_eff[Temp_FWHM < FWHM_inst] = FWHM_inst[Temp_FWHM < FWHM_inst]
+            LSF[Temp_FWHM < FWHM_inst]      = FWHM_inst[Temp_FWHM < FWHM_inst]
+        FWHM_dif  = np.sqrt(FWHM_eff ** 2 - Temp_FWHM ** 2)
+        sigma_dif = FWHM_dif / 2.355 / (lam[1] - lam[0])  # Sigma difference in pixels
+
+        temp = np.empty((TemNew.size, par.size))
+        for i in range(par.size):
+            temp0 = log_rebin(lam_range_temp, flux[i, :], velscale=velscale)[0]
+            if np.isscalar(FWHM_dif):
+                temp1 = ndimage.gaussian_filter1d(temp0, sigma_dif)
+            else:
+                temp1 = gaussian_filter1d(temp0, sigma_dif)             # convolution with variable sigma
+            tempNew = temp1 / np.mean(temp1)
+            temp[:, i] = tempNew
+            
+
+        self.templates    = temp
+        self.log_lam_temp = log_lam_temp
+        self.teff_grid = par['Teff']
+        self.feh_grid  = par['FeH']
+        self.logg_grid = par['logg']
+        self.LSF       = Temp_FWHM
+        self.velscale  = velscale
+        
+class SingleStar():
+
+    """
+    Class of single stelar spectrum
+
+    Parameters
+    ----------
+    config : class
+        Class of configuration
+    template : class
+        Class of single stellar population template
+    mag : float, optional
+        Magnitude in SDSS r-band, by default 15.0
+    Teff : float, optional
+        Effective tempreture, by default 10000.0K
+    FeH : float, optional
+        Metallicity of stellar, by default 0.0
+    vel : float, optional
+        Line of sight velocity, by default 100.0km/s
+    Ebv : float, optional
+        Dust extinction, by default 0.1
+    """
+
+    def __init__(self, config, template, mag = 15.0, teff = 10000.0, feh = 0.0, vel = 100.0, ebv = 0.0):
+    
+        StarTemp = template.templates
+        
+        # Select metal bins
+        idx_FeH = (np.abs(template.feh_grid - feh) < 0.5)
+        tpls = StarTemp[:, idx_FeH]
+        
+        # Select Teff bins
+        Teff_FeH = template.teff_grid[idx_FeH]
+        minloc   = np.argmin(abs(teff - Teff_FeH))
+        starspec = tpls[:, minloc]
+        
+        wave = np.exp(template.log_lam_temp)
+        
+        # Dust Reddening
+        if np.isscalar(ebv):
+            starspec = reddening(wave, starspec, ebv = ebv)
+            
+        # Redshift
+        redshift = vel / 3e5
+        wave_r = wave * (1 + redshift)
+        
+        flux = np.interp(config.wave, wave_r, starspec)
+        
+        # Calibration
+        if np.isscalar(mag):
+            flux = calibrate(config.wave, flux, mag, filtername='SLOAN_SDSS.r')
+        
+        # Convert to input wavelength
+        self.wave = config.wave
+        self.flux = flux
\ No newline at end of file