fit PEP8 bugs

csst_ifs_gehong/map2d.py

@@ -2,11 +2,11 @@ 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):
+
+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):
     """
     Sersic2D: Caculate the surface brightness at given spatial position base on the Sersic model.
 
@@ -35,7 +35,6 @@ def Sersic2D(x, y, mag = 12.0, r_eff = 1.0, n = 2.0, ellip = 0.5,
 
     Returns
     -------
-    
         _description_
     """
 
@@ -48,7 +47,7 @@ def Sersic2D(x, y, mag = 12.0, r_eff = 1.0, n = 2.0, ellip = 0.5,
     cos_theta, sin_theta = np.cos(theta_radian), np.sin(theta_radian)
     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)
+    z = (abs(x_maj / a) ** 2 + abs(x_min / b) ** 2) ** (1 / 2)
     profile = np.exp(-bn * (z ** (1 / n) - 1))
     
     # Normalization
@@ -61,8 +60,9 @@ def Sersic2D(x, y, mag = 12.0, r_eff = 1.0, n = 2.0, ellip = 0.5,
     
     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):
+
+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: Caculate the velocity at given spatial position base on the rotating disk model. 
     The rotation curve is adpot as the tanh model. 
@@ -100,13 +100,14 @@ def VelMap2D(x, y, vmax = 200.0, rt = 1.0, ellip = 0.5,
     cos_theta, sin_theta = np.cos(theta_radian), np.sin(theta_radian)
     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)
+    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.0, r_eff = 1.0, gred = -1.0, ellip = 0.5, 
-              theta = 0.0, x_0 = 0.0, y_0 = 0.0):
+
+def GradMap2D(x, y, a0=10.0, r_eff=1.0, gred=-1.0, ellip=0.5, 
+              theta=0.0, x_0=0.0, y_0=0.0):
     """
     GradMap2D: Caculate the intensity at given spatial position base on the disk model. 
     The radial profile is adpot as the gradient model. 
@@ -138,7 +139,6 @@ def GradMap2D(x, y, a0 = 10.0, r_eff = 1.0, gred = -1.0, ellip = 0.5,
         _description_
     """
 
-
     # Convert the angle in degree to angle to radians
     theta_radian = theta / 180 * np.pi
 
@@ -147,11 +147,12 @@ def GradMap2D(x, y, a0 = 10.0, r_eff = 1.0, gred = -1.0, ellip = 0.5,
     cos_theta, sin_theta = np.cos(theta_radian), np.sin(theta_radian)
     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)
+    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):
@@ -166,11 +167,11 @@ class Map2d(object):
         self.xsamp = config.dpix
         self.ysamp = config.dpix
         startx = -(config.nx - 1) / 2.0 * self.xsamp
-        stopx  = (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)
+        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
@@ -178,8 +179,8 @@ class Map2d(object):
 
         self.nx = config.nx
         self.ny = config.ny
-        self.x  = x
-        self.y  = y
+        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]
@@ -207,7 +208,7 @@ class Map2d(object):
         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):
+    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
 
@@ -229,7 +230,7 @@ class Map2d(object):
         if (mag > 26) or (mag < 8):
             print("Notice: Your input integral magnitude of Sersic mode (mag) is > 26 mag or < 8 mag.")
         if (r_eff <= 0):
-            #raise Exception("Effective radius (r_eff) should be > 0 arcsec!")
+            # raise Exception("Effective radius (r_eff) should be > 0 arcsec!")
             print("Effective radius (r_eff) should be > 0 arcsec!")
         if (n < 0):
             raise Exception("Sersic index (n) should be > 0!")
@@ -240,17 +241,17 @@ class Map2d(object):
         if (theta > 180) or (theta < -180):
             print("Notice: Your input position angle (theta) is > 180 degree or < -180 degree.")
 
-        self.reff  = r_eff / self.xsamp
-        self.mag   = mag
-        self.n     = n
+        self.reff = r_eff / self.xsamp
+        self.mag = mag
+        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)
+        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):
+    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
 
@@ -268,25 +269,25 @@ class Map2d(object):
 
         # Check Input Parameters
         if (vmax <= 0):
-            #print("Notice: Your input maximum rotational velocity (vmax) is <= 0 km/s!")
+            # print("Notice: Your input maximum rotational velocity (vmax) is <= 0 km/s!")
             raise Exception("Maximum rotational velocity (vmax) should be >0 km/s!")
         if (rt <= 0):
-            #raise Exception("Turn-over radius (rt) should be > 0 arcsec!")
+            # raise Exception("Turn-over radius (rt) should be > 0 arcsec!")
             raise Exception("Turn-over radius (rt) should be > 0 arcsec!")
         if (ellip > 1) or (ellip < 0):
             raise Exception("Ellipcity (ellip) should be >= 0 and < 1!")
-            #print("Ellipcity (ellip) should be >= 0 and < 1!")
+            # print("Ellipcity (ellip) should be >= 0 and < 1!")
         if (theta > 180) or (theta < -180):
             print("Notice: Your input position angle (theta) is > 180 degree or < -180 degree.")
 
-        self.vmax  = vmax
-        self.rt    = rt / self.xsamp
+        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)
+        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):
+    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
 
@@ -306,20 +307,20 @@ class Map2d(object):
         # Check Input Parameters
         if (r_eff <= 0):
             raise Exception("Effective radius (r_eff) should be > 0 arcsec!")
-            #print("Effective radius (r_eff) should be > 0 arcsec!")
+            # print("Effective radius (r_eff) should be > 0 arcsec!")
         if (ellip > 1) or (ellip < 0):
-            #raise Exception("Ellipcity (ellip) should be >= 0 and < 1!")
+            # raise Exception("Ellipcity (ellip) should be >= 0 and < 1!")
             print("Notice: Ellipcity (ellip) should be >= 0 and < 1!")
         if (theta > 180) or (theta < -180):
             print("Notice: Your input position angle (theta) is > 180 degree or < -180 degree.")
 
-        self.a0    = a0
-        self.reff  = r_eff / self.xsamp
-        self.gred  = gred
+        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)
+        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):
         """
@@ -334,7 +335,8 @@ class Map2d(object):
             self.map = resize(image, (self.nx, self.ny))
         else:
             print("Input array should be a 2d-array!")
-            
+
+
 class StellarPopulationMap():
     """
     Class of 2D maps for the parameters of stellar population, such as 
@@ -358,55 +360,56 @@ class StellarPopulationMap():
     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):
+    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.nx = config.nx
+        self.ny = config.ny
+        self.dpix = config.dpix
         self.fov_x = config.fov_x
         self.fov_y = config.fov_y
         
-        if (sbright == None):
+        if (sbright is None):
             print('Input SurfaceBrightness Map is empty!')
         else:
             self.sbright = sbright.map
             self.mag = self.sbright - 2.5 * np.log10(self.dpix * self.dpix)
         
-        if (logage == None):
+        if (logage is None):
             print('Input Age Map is empty!')
         else:
-            self.logage  = logage.map
+            self.logage = logage.map
             self.age = 10 ** self.logage / 1e9
         
-        if (feh == None):
+        if (feh is None):
             print('Input Metallicity Map is empty!')
         else:
-            self.feh     = feh.map
+            self.feh = feh.map
         
-        if (vel == None):
+        if (vel is None):
             print('Input Velocity Map is empty!')
         else:
-            self.vel     = vel.map
+            self.vel = vel.map
         
-        if (vdisp == None):
+        if (vdisp is None):
             print('Input VelocityDispersion Map is empty!')
         else:
-            self.vdisp   = vdisp.map
+            self.vdisp = vdisp.map
             ind_overrange = (self.vdisp < 10)
             if len(self.vdisp[ind_overrange]) > 0:
                 print("Notice: Spaxel with <10km/s in the input vdisp map will be automatically adjusted to 10km/s.")
             self.vdisp[ind_overrange] = 10
         
-        if (ebv == None):
+        if (ebv is None):
             print('Input EBV Map is empty!')
         else:
-            self.ebv     = ebv.map
+            self.ebv = ebv.map
             ind_overrange = (self.ebv < 0)
             if len(self.ebv[ind_overrange]) > 0:
                 print("Notice: Spaxel with < 0 mag in the input ebv map will be automatically adjusted to 0 mag.")
             self.ebv[ind_overrange] = 0
         
+
 class IonizedGasMap():
     """
     Class of 2D maps for the parameters of ionized gas, such as 
@@ -428,42 +431,42 @@ class IonizedGasMap():
     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):
+    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.nx = config.nx
+        self.ny = config.ny
+        self.dpix = config.dpix
+        self.fov_x = config.fov_x
+        self.fov_y = config.fov_y
         
-        if (halpha == None):
+        if (halpha is None):
             print('Input Halpha Map is empty!')
         else:
             self.halpha = halpha.map
 
-        if (zh == None):
+        if (zh is None):
             print('Input ZH Map is empty!')
         else:
-            self.zh     = zh.map
+            self.zh = zh.map
 
-        if (vel == None):
+        if (vel is None):
             print('Input Vel Map is empty!')
         else:
-            self.vel    = vel.map
+            self.vel = vel.map
 
-        if (vdisp == None):
+        if (vdisp is None):
             print('Input Vdisp Map is empty!')
             ind_overrange = (self.vdisp < 10)
             if len(self.vdisp[ind_overrange]) > 0:
                 print("Notice: Spaxel with <10km/s in the input vdisp map will be automatically adjusted to 10km/s.")
             self.vdisp[ind_overrange] = 10
         else:
-            self.vdisp  = vdisp.map
+            self.vdisp = vdisp.map
 
-        if (ebv == None):
+        if (ebv is None):
             print('Input EBV Map is empty!')
         else:
-            self.ebv    = ebv.map
+            self.ebv = ebv.map
             ind_overrange = (self.ebv < 0)
             if len(self.ebv[ind_overrange]) > 0:
                 print("Notice: Spaxel with < 0 mag in the input ebv map will be automatically adjusted to 0 mag.")

csst_ifs_gehong/spec1d.py

@@ -6,10 +6,12 @@ import numpy as np
 import astropy.units as u
 from astropy.io import fits
 from scipy.stats import norm
+from scipy import ndimage
 from scipy.interpolate import interp1d
 
 data_path = os.path.dirname(os.path.abspath(__file__))
 
+
 def readcol(filename, **kwargs):
     """
     readcol, taken from ppxf.
@@ -27,7 +29,7 @@ def readcol(filename, **kwargs):
     f = np.genfromtxt(filename, dtype=None, **kwargs)
 
     t = type(f[0])
-    if t == np.ndarray or t == np.void: # array or structured array
+    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.
@@ -38,10 +40,11 @@ def readcol(filename, **kwargs):
     # 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]
+        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. 
@@ -78,38 +81,39 @@ def log_rebin(lamRange, spec, oversample=False, velscale=None, flux=False):
     assert len(s) == 1, 'input spectrum must be a vector'
     n = s[0]
     if oversample:
-        m = int(n*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
+    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
+        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
+        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)
+    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])
+    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)
+    logLam = np.log(np.sqrt(newBorders[1:] * newBorders[:-1]) * dLam)
 
     return specNew, logLam, velscale
 
+
 def gaussian_filter1d(spec, sig):
     """
     One-dimensional Gaussian convolution
@@ -126,23 +130,24 @@ def gaussian_filter1d(spec, sig):
         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
+    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]
+        a[j, p:-p] = spec[j:n - m + j + 1]
 
-    gau = np.exp(-x2[:, None]/(2*sig**2))
+    gau = np.exp(-x2[:, None] / (2 * sig**2))
     gau /= np.sum(gau, 0)[None, :]  # Normalize kernel
 
-    conv_spectrum = np.sum(a*gau, 0)
+    conv_spectrum = np.sum(a * gau, 0)
 
     return conv_spectrum
 
-def calibrate(wave, flux, mag, filtername = 'SLOAN_SDSS.r'):
+
+def calibrate(wave, flux, mag, filtername='SLOAN_SDSS.r'):
     """
     Flux calibration of spectrum
 
@@ -164,11 +169,11 @@ def calibrate(wave, flux, mag, filtername = 'SLOAN_SDSS.r'):
     """
     # Loading response curve
     if filtername == '5100':
-        wave0 = np.linspace(3000,10000,7000)
+        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'
+        filter_file = data_path + '/data/filter/' + filtername + '.filter'
         wave0, response0 = readcol(filter_file)
     
     # Setting the response
@@ -183,7 +188,7 @@ def calibrate(wave, flux, mag, filtername = 'SLOAN_SDSS.r'):
     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
+    realflux = (mag * u.STmag).to(u.erg / u.s / u.cm**2 / u.AA).value
 
     # Normalization
     flux_ratio = realflux / preflux
@@ -194,7 +199,8 @@ def calibrate(wave, flux, mag, filtername = 'SLOAN_SDSS.r'):
 # ----------------
 # Reddening Module
 
-def Calzetti_Law(wave, Rv = 4.05):
+
+def Calzetti_Law(wave, Rv=4.05):
     """
     Dust Extinction Curve of Calzetti et al. (2000)
 
@@ -210,18 +216,19 @@ def Calzetti_Law(wave, Rv = 4.05):
     float
         Extinction value corresponding to the input wavelength
     """
-    wave_number = 1./(wave * 1e-4)
+    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
+    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
+    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):
+
+def reddening(wave, flux, ebv=0.0, law='calzetti', Rv=4.05):
     """
     Reddening an input spectra through a given reddening curve.
 
@@ -243,7 +250,7 @@ def reddening(wave, flux, ebv = 0.0, law = 'calzetti', Rv = 4.05):
     float array
         Flux of spectra after reddening
     """
-    curve = Calzetti_Law(wave, Rv = Rv)
+    curve = Calzetti_Law(wave, Rv=Rv)
     fluxNew = flux / (10. ** (0.4 * ebv * curve))
     return fluxNew
 
@@ -251,6 +258,7 @@ def reddening(wave, flux, ebv = 0.0, law = 'calzetti', Rv = 4.05):
 # Emission Lines
 ################
 
+
 def SingleEmissinoLine(wave, line_wave, FWHM_inst):
     """
     Profile of single emission line (Gaussian profile)
@@ -273,6 +281,7 @@ def SingleEmissinoLine(wave, line_wave, FWHM_inst):
     flux = norm.pdf(wave, line_wave, sigma)
     return flux
 
+
 class EmissionLineTemplate():
     """
     Template for the emission lines
@@ -290,7 +299,7 @@ class EmissionLineTemplate():
         by default 'hii'
     """
     
-    def __init__(self, config, lam_range = [500, 15000], dlam = 0.1, model = '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)
@@ -305,8 +314,8 @@ class EmissionLineTemplate():
             line_table = fits.open(flux_table_file)
 
             # Emission line list
-            line_list  = line_table[1].data
-            line_wave  = line_list['Wave']
+            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])
@@ -315,7 +324,7 @@ class EmissionLineTemplate():
             
             # Make parameter grid
             grid = line_table[2].data
-            self.logz_grid   = grid['logZ']
+            self.logz_grid = grid['logZ']
         
         # Narrow line region model of 
         if model == 'nlr':
@@ -325,8 +334,8 @@ class EmissionLineTemplate():
             line_table = fits.open(flux_table_file)
 
             # Emission line list
-            line_list  = line_table[1].data
-            line_wave  = line_list['Wave']
+            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])
@@ -335,7 +344,7 @@ class EmissionLineTemplate():
             
             # Make parameter grid
             grid = line_table[2].data
-            self.logz_grid   = grid['logZ']
+            self.logz_grid = grid['logZ']
         
         # Flux ratio
         flux_ratio = line_table[3].data
@@ -344,14 +353,15 @@ class EmissionLineTemplate():
         # 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)
+            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_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():
 
     """
@@ -380,29 +390,28 @@ class HII_Region():
         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):
+    def __init__(self, config, temp, halpha=100.0, logz=0.0, vel=100.0, vdisp=120.0, ebv=0.1):
         
         # Check Input Parameters
-        #if (halpha < 0):
+        # if (halpha < 0):
         #    print("Notice: Your input Halpha flux (halpha) is < 0 erg/s/A/cm^2. ")
         if (logz >= -2) & (logz <= 0.5):
             indz = np.argmin(np.abs(logz - temp.logz_grid))
             flux_ratio = temp.flux_ratio[indz, :]
         else:
             raise ValueError('Input gas-phase metallicity (logz) should be [-2, 0.5]!')
-        #if (ebv < 0):
+        # if (ebv < 0):
         #    ebv = 0
         #    print('Notice: Your input dust extinction (ebv) is < 0 mag, which will be automatically adjusted to 0 mag. ')
         
         # Make emission line spectra through adding emission lines                 
         emlines = temp.emission_lines * flux_ratio
-        flux_combine = np.sum(emlines, axis = 1)
+        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)
+            flux_dust = reddening(temp.wave, flux_calibrate, ebv=ebv)
             
         # Broadening caused by Velocity Dispersion
         velscale = 10
@@ -412,7 +421,7 @@ class HII_Region():
         sigma_gas = vdisp / velscale                                # in pixel
         sigma_LSF = temp.FWHM_inst / (np.exp(logLam)) * 3e5 / velscale          # in pixel
         
-        if sigma_gas>0: 
+        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.)
@@ -435,6 +444,7 @@ class HII_Region():
 # AGN Spectra
 #############
 
+
 class AGN_NLR():
 
     """
@@ -463,8 +473,7 @@ class AGN_NLR():
         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):
+    def __init__(self, config, temp, halpha=100.0, logz=0.0, vel=100.0, vdisp=120.0, ebv=0.1):
         
         if (logz >= -2.3) & (logz <= 0.54):
             indz = np.argmin(np.abs(logz - temp.logz_grid))
@@ -473,13 +482,13 @@ class AGN_NLR():
             raise ValueError('The value of logZ is not in the range of [-2.3, 0.54]')
         
         # Make emission line spectra through adding emission lines                 
-        emlines = temp.emission_lines * (flux_ratio / flux_ratio[6] )
-        flux_combine = np.sum(emlines, axis = 1)
+        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)
+            flux_dust = reddening(temp.wave, flux_calibrate, ebv=ebv)
             
         # Broadening caused by Velocity Dispersion
         velscale = 10
@@ -489,7 +498,7 @@ class AGN_NLR():
         sigma_gas = vdisp / velscale                                        # in pixel
         sigma_LSF = temp.FWHM_inst / (np.exp(logLam)) * 3e5 / velscale      # in pixel
         
-        if sigma_gas>0: 
+        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.)
@@ -505,6 +514,7 @@ class AGN_NLR():
         self.wave = config.wave
         self.flux = flux_red
 
+
 class AGN_BLR():
 
     """
@@ -526,8 +536,8 @@ class AGN_BLR():
         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]):
+    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)
         
@@ -537,23 +547,23 @@ class AGN_BLR():
         
         # 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])
+            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_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_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)
+            flux_dust = reddening(wave_rest, flux_calibrate, ebv=ebv)
         else:
             flux_dust = flux_calibrate
             
@@ -565,6 +575,7 @@ class AGN_BLR():
         self.wave = config.wave
         self.flux = flux_red
 
+
 class AGN_FeII():
     """
     Class for FeII emission lines of AGN
@@ -583,18 +594,18 @@ class AGN_FeII():
         Dust extinction, by default 0.1
     """
     
-    def __init__(self, config, hbeta_broad = 100.0, r4570 = 0.4, ebv = 0.1, vel = 100.0):
+    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']
+        wave_rest = data['WAVE']
         flux_model = data['FLUX']
         
         # Determine the flux of FeII
-        Fe4570_temp  = 100
+        Fe4570_temp = 100
         Fe4570_model = hbeta_broad * r4570
         Ratio_Fe4570 = Fe4570_model / Fe4570_temp
         
@@ -603,18 +614,19 @@ class AGN_FeII():
         
         # Dust attenuation
         if np.isscalar(ebv):
-            flux_dust = reddening(wave_rest, flux_calibrate, ebv = 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)
+        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():
 
     """
@@ -634,28 +646,29 @@ class AGN_Powerlaw():
         Dust extinction, by default 0.1
     """
     
-    def __init__(self, config, m5100 = 1000.0, alpha = -1.5, vel = 100.0, ebv = 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
+        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)
+            flux_dust = reddening(wave_rest, flux_calibrate, ebv=ebv)
         else:
             flux_dust = flux_calibrate
             
         # Redshift
         redshift = vel / 3e5
-        wave_r   = wave_rest * (1 + redshift)
+        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():
 
     """
@@ -689,12 +702,12 @@ class AGN():
     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):
+    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):
         
         # Check Input Parameters
-        if (bhmass <0):
+        if (bhmass < 0):
             raise Exception("Input black hole mass (bhmass) should be >= 0 solarmass.")
         if (edd_ratio < 0):
             raise Exception("Input Eddington ratio (edd_ratio) should be >= 0.")
@@ -716,22 +729,23 @@ class AGN():
         self.flux = self.wave * 0
         
         if (vdisp_narrow > 0) & (halpha_narrow > 0): 
-            NLR = AGN_NLR(config, nlr_template, halpha = halpha_narrow, logz = logz,
-                        vel = vel, vdisp = vdisp_narrow, ebv = ebv)
+            NLR = AGN_NLR(config, nlr_template, halpha=halpha_narrow, logz=logz,
+                          vel=vel, vdisp=vdisp_narrow, ebv=ebv)
             self.flux = self.flux + NLR.flux
 
         if (halpha_broad > 0) & (vdisp_broad > 0):
-            BLR = AGN_BLR(config, hbeta_flux = halpha_broad / 2.579, 
-                          hbeta_fwhm = vdisp_broad / 2.355, ebv = ebv, vel = vel)
+            BLR = AGN_BLR(config, hbeta_flux=halpha_broad / 2.579, 
+                          hbeta_fwhm=vdisp_broad / 2.355, ebv=ebv, vel=vel)
             self.flux = self.flux + BLR.flux
 
         if (bhmass > 0) & (edd_ratio > 0):
-            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)
+            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.flux = self.flux + PL.flux + Fe.flux
     
-def BHmass_to_M5100(bhmass, edd_ratio = 0.05, dist = 21.0):
+
+def BHmass_to_M5100(bhmass, edd_ratio=0.05, dist=21.0):
     """
     Caculate magnitude at 5100A according to the black hole mass
 
@@ -765,16 +779,6 @@ def BHmass_to_M5100(bhmass, edd_ratio = 0.05, dist = 21.0):
 # 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):
     """
@@ -799,8 +803,8 @@ def age_metal(filename):
         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")]
+    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:
@@ -810,6 +814,7 @@ def age_metal(filename):
 
     return age, metal
 
+
 class StellarContinuumTemplate(object):
     """
     Class of single stellar population template. 
@@ -826,9 +831,9 @@ class StellarContinuumTemplate(object):
     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):
+    def __init__(self, config, velscale=50,
+                 pathname=data_path + '/data/EMILES/Ech*_baseFe.fits', 
+                 normalize=False):
 
         FWHM_inst = config.inst_fwhm
         
@@ -849,9 +854,9 @@ class StellarContinuumTemplate(object):
         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)])
+        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']
+        # 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))
@@ -873,19 +878,19 @@ class StellarContinuumTemplate(object):
         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
+        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)
+        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
+        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.
@@ -898,14 +903,14 @@ class StellarContinuumTemplate(object):
                     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]
+                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
+                age_grid[j, k] = age
                 metal_grid[j, k] = metal
 
-        self.templates = templates/np.median(templates)  # Normalize by a scalar
+        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
@@ -920,18 +925,19 @@ class StellarContinuumTemplate(object):
         massfile = isedpath + 'out_mass_CH_PADOVA00'
         n_metal = self.n_metal
         n_ages = self.n_ages
-        fage = self.age_grid[:,0]
+        # fage = self.age_grid[:,0]
 
-        fMs = np.zeros((n_ages,n_metal))
+        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]
+                fMs[j, i] = Ms[locmin]
 
         return fMs
     
+
 class StellarContinuum():
     """
     The class of stellar continuum
@@ -955,8 +961,8 @@ class StellarContinuum():
     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):
+    def __init__(self, config, template, mag=15.0, age=1.0, feh=0.0, 
+                 vel=100.0, vdisp=100.0, ebv=0.1):
         
         # Check Input Parameters
         if (mag > 26) or (mag < 8):
@@ -967,24 +973,23 @@ class StellarContinuum():
             print("Notice: Your input median metallicity (feh) is beyond the range of stellar template [-2.32, 0.22], which will be automatically adjusted to the upper limit.")
         if (feh < -2.32):
             raise Exception('Your input median metallicity (feh) is beyond the range of stellar template [-2.32, 0.22]')
-        #if (ebv < 0):
+        # if (ebv < 0):
         #    ebv = 0
         #    print('Notice: Your input dust extinction (ebv) is < 0 mag, which will be automatically adjusted to 0 mag. ')
-        
 
         # -----------------
         # Stellar Continuum
         SSP_temp = template.templates
 
         # Select metal bins
-        metals = template.metal_grid[0,:]
+        metals = template.metal_grid[0, :]
         minloc = np.argmin(abs(feh - metals))
         tpls = SSP_temp[:, :, minloc]
-        fmass = template.fmass_ssp()[:, minloc]
+        # fmass = template.fmass_ssp()[:, minloc]
         
         # Select age bins
-        Ages = template.age_grid[:,0]
-        minloc = np.argmin(abs(age-Ages))
+        Ages = template.age_grid[:, 0]
+        minloc = np.argmin(abs(age - Ages))
         Stellar = tpls[:, minloc]
         
         wave = np.exp(template.log_lam_temp)
@@ -993,7 +998,7 @@ class StellarContinuum():
         sigma_gal = vdisp / template.velscale                   # in pixel
         sigma_LSF = template.LSF / template.velscale                 # in pixel
         
-        if sigma_gal>0: 
+        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.)
@@ -1007,7 +1012,7 @@ class StellarContinuum():
 
         # Dust Reddening
         if np.isscalar(ebv):
-            flux0 = reddening(wave, flux0, ebv = ebv)
+            flux0 = reddening(wave, flux0, ebv=ebv)
             
         # Redshift
         redshift = vel / 3e5
@@ -1027,6 +1032,7 @@ class StellarContinuum():
 # Single Star Spectra
 #####################
 
+
 class SingleStarTemplate():
     """
     Class of single stellar template
@@ -1038,18 +1044,18 @@ class SingleStarTemplate():
     velscale : float, option
         velocity scale in km/s per pixels, by default 20.0km/s
     """
-    def __init__(self, config, velscale = 20):
+    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']
+        lam = hdulist[1].data['Wave']
         flux = hdulist[2].data
-        par  = hdulist[3].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]
+        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)
@@ -1057,18 +1063,18 @@ class SingleStarTemplate():
         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
+        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
+            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)
+            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))
@@ -1080,16 +1086,16 @@ class SingleStarTemplate():
                 temp1 = gaussian_filter1d(temp0, sigma_dif)             # convolution with variable sigma
             tempNew = temp1 / np.mean(temp1)
             temp[:, i] = tempNew
-            
 
-        self.templates    = temp
+        self.templates = temp
         self.log_lam_temp = log_lam_temp
         self.teff_grid = par['Teff']
-        self.feh_grid  = par['FeH']
+        self.feh_grid = par['FeH']
         self.logg_grid = par['logg']
-        self.LSF       = Temp_FWHM
-        self.velscale  = velscale
+        self.LSF = Temp_FWHM
+        self.velscale = velscale
         
+
 class SingleStar():
 
     """
@@ -1113,7 +1119,7 @@ class SingleStar():
         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):
+    def __init__(self, config, template, mag=15.0, teff=10000.0, feh=0.0, vel=100.0, ebv=0.0):
 
         # Check Input Parameters
         if (mag > 26) or (mag < 8):
@@ -1130,7 +1136,6 @@ class SingleStar():
         if (feh < -2.69):
             raise Exception('Your input metallicity (feh) is beyond the range of stellar template [-2.69, 0.81]')
         
-    
         StarTemp = template.templates
         
         # Select metal bins
@@ -1139,14 +1144,14 @@ class SingleStar():
         
         # Select Teff bins
         Teff_FeH = template.teff_grid[idx_FeH]
-        minloc   = np.argmin(abs(teff - Teff_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)
+            starspec = reddening(wave, starspec, ebv=ebv)
             
         # Redshift
         redshift = vel / 3e5