add justification of input pars range

csst_ifs_gehong/map2d.py

@@ -8,14 +8,14 @@ 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
+    Sersic2D: Caculate the surface brightness at given spatial position base on the Sersic model.
 
     Parameters
     ----------
     x : float array
-        _description_
-    y : _type_
-        _description_
+        Coordinate of spatial position in the x-axis
+    y : float array
+        Coordinate of spatial position in the y-axis
     mag : float, optional
         Integral magnitude of sersic model, by default 12.0
     r_eff : float, optional
@@ -27,9 +27,9 @@ def Sersic2D(x, y, mag = 12.0, r_eff = 1.0, n = 2.0, ellip = 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
+        Offset of the center of Sersic model in x-axis, by default 0.0
     y_0 : float, optional
-        Offset of the center of Sersic model, by default 0.0
+        Offset of the center of Sersic model in y-axis, by default 0.0
     pixelscale : float, optional
         Size of each pixel in arcsec^2, by default 0.01
 
@@ -38,10 +38,14 @@ def Sersic2D(x, y, mag = 12.0, r_eff = 1.0, n = 2.0, ellip = 0.5,
     
         _description_
     """
+
+    # Convert the angle in degree to angle to radians
+    theta_radian = theta / 180 * np.pi
+
     # 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)
+    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)
@@ -60,35 +64,40 @@ def Sersic2D(x, y, mag = 12.0, r_eff = 1.0, n = 2.0, ellip = 0.5,
 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_
+    VelMap2D: Caculate the velocity at given spatial position base on the rotating disk model. 
+    The rotation curve is adpot as the tanh model. 
 
     Parameters
     ----------
-    x : _type_
-        _description_
-    y : _type_
-        _description_
-    vmax : int, optional
-        _description_, by default 200
-    rt : int, optional
-        _description_, by default 1
+    x : float array
+        Coordinate of spatial position in the x-axis
+    y : float array
+        Coordinate of spatial position in the y-axis
+    vmax : float, optional
+        Maximum rotational velocity of rotation curve, by default 200.0
+    rt : float, optional
+        Turn-over radius of rotation curve, by default 1.0
     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
+        Ellipicity of disk model, by default 0.5
+    theta : float, optional
+        Position angle of disk model, by default 0.0
+    x_0 : float, optional
+        Offset of the center of disk model in the x-axis, by default 0.0
+    y_0 : float, optional
+        Offset of the center of disk model in the y-axis, by default 0.0
 
     Returns
     -------
     _type_
         _description_
     """
+
+    # Convert the angle in degree to angle to radians
+    theta_radian = theta / 180 * np.pi
+
     # Produce tanh profile
     a, b = rt, (1 - ellip) * rt
-    cos_theta, sin_theta = np.cos(theta), np.sin(theta)
+    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)
@@ -96,40 +105,46 @@ def VelMap2D(x, y, vmax = 200.0, rt = 1.0, ellip = 0.5,
     
     return profile
 
-def GradMap2D(x, y, a0 = 10, r_eff = 1, gred = -1, ellip = 0.5, 
-              theta = 0, x_0 = 0, y_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 _summary_
+    GradMap2D: Caculate the intensity at given spatial position base on the disk model. 
+    The radial profile is adpot as the gradient model. 
 
     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
+    x : float array
+        Coordinate of spatial position in the x-axis
+    y : float array
+        Coordinate of spatial position in the y-axis
+    a0 : float, optional
+        Intensity at the center, by default 10.0
+    r_eff : float, optional
+        Effective radius, by default 1.0
+    gred : float, optional
+        Gradient, by default -1.0
     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
+        Ellipicity of disk model, by default 0.5
+    theta : float, optional
+        Position angle of disk model, by default 0.0
+    x_0 : float, optional
+        Offset of the center of disk model in the x-axis, by default 0.0
+    y_0 : float, optional
+        Offset of the center of disk model in the y-axis, by default 0.0
 
     Returns
     -------
     _type_
         _description_
     """
+
+
+    # Convert the angle in degree to angle to radians
+    theta_radian = theta / 180 * np.pi
+
     # Produce gradiant profile
     a, b = r_eff, (1 - ellip) * r_eff
-    cos_theta, sin_theta = np.cos(theta), np.sin(theta)
+    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)
@@ -200,7 +215,7 @@ class Map2d(object):
         ----------
         mag : float, optional
             Integral magnitude, by default 12.0
-        r_eff : float, optional
+        reff : float, optional
             Effective radius in arcsec, by default 2.0
         n : float, optional
             Sersic index, by default 2.5
@@ -209,8 +224,21 @@ class Map2d(object):
         theta : float, optional
             Position angle in degree, by default -50.0
         """
-        self.mag   = mag
+
+        # Check Input Parameters
+        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!")
+        if (n <= 0):
+            raise Exception("Sersic index (n) should be > 0!")
+        if (ellip >= 1) or (ellip < 0):
+            raise Exception("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.reff  = r_eff / self.xsamp
+        self.mag   = mag
         self.n     = n
         self.ellip = ellip
         self.theta = theta
@@ -234,6 +262,17 @@ class Map2d(object):
         theta : float, optional
             Position angle of rotating disk, by default -50.0
         """
+
+        # Check Input Parameters
+        if (vmax <= 0):
+            print("Notice: Your input maximum rotational velocity (vmax) is <= 0 km/s!")
+        if (rt <= 0):
+            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!")
+        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.ellip = ellip
@@ -248,7 +287,7 @@ class Map2d(object):
         Parameters
         ----------
         a0 : float, optional
-            Amplitude at the central pixel, by default 10
+            Amplitude at the center, by default 10
         r_eff : float, optional
             Effective radius, by default 1
         gred : float, optional
@@ -258,6 +297,14 @@ class Map2d(object):
         theta : int, optional
             Position angle, by default 0
         """
+        # Check Input Parameters
+        if (r_eff <= 0):
+            raise Exception("Effective radius (r_eff) should be > 0 arcsec!")
+        if (ellip >= 1) or (ellip < 0):
+            raise Exception("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
@@ -320,8 +367,16 @@ class StellarPopulationMap():
         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
+        # Check Input Maps
+        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
+
+        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():
     """
@@ -358,5 +413,16 @@ class IonizedGasMap():
         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
+        #self.vdisp[self.vdisp < 10] = 10
+        #self.ebv[self.ebv < 0]      = 0
+
+        # Check Input Maps
+        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
+
+        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
\ No newline at end of file

csst_ifs_gehong/spec1d.py

@@ -382,12 +382,18 @@ class HII_Region():
     
     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):
+        #    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('The value of logZ is not in the range of [-2, 0.5]!')
+            raise ValueError('Input gas-phase metallicity (logz) should be [-2, 0.5]!')
+        #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
@@ -414,6 +420,7 @@ class HII_Region():
             sigma_dif[idx] = 0.1
             flux_broad = gaussian_filter1d(flux_logwave, sigma_dif)
         else:
+            print('Notice: Your input velocity dispersion (vdisp) is <= 0 km/s. The nagtive value will be automatically adjusted to 0 km/s. ')
             flux_broad = flux_logwave
             
         # Redshift
@@ -459,11 +466,11 @@ class AGN_NLR():
     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):
+        if (logz >= -2.3) & (logz <= 0.54):
             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]!')
+            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] )
@@ -686,21 +693,43 @@ class AGN():
                  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)
+        # Check Input Parameters
+        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.")
+        if (halpha_broad < 0):
+            print("Notice: Your input braod Halpha emission line flux (halpha_broad) is <0 erg/s/A/cm^2, which will be automatically adjusted to 0.")
+        if (vdisp_broad < 0):
+            raise Exception('Input velocity dispersion of broad lines (vdisp_broad) should be >= 0 km/s')
+        if (vdisp_narrow < 0):
+            raise Exception('Input velocity dispersion of narrow lines (vdisp_narrow) should be >= 0 km/s')
+        if (logz > 0.54):
+            logz = 0.54
+            print("Notice: Your input metallicity of narrow lines (logz) is > 0.54, which will be automatically adjusted to 0.54")
+        if (logz < -2.3):
+            raise Exception('Input metallicity of narrow lines (logz) should be < -2.3')
+        if (dist <= 0):
+            raise Exception('Input luminosity distance (dist) should be > 0 Mpc')
         
         self.wave = config.wave
-        self.flux = NLR.flux + PL.flux + Fe.flux
+        self.flux = self.wave * 0
         
-        if halpha_broad > 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)
+            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)
             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)
+            self.flux = self.flux + PL.flux + Fe.flux
     
 def BHmass_to_M5100(bhmass, edd_ratio = 0.05, dist = 21.0):
     """
@@ -928,10 +957,23 @@ class StellarContinuum():
     """
     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):
+            print("Notice: Your input magnitude (mag) is > 26 mag or < 8 mag.")
+        if (age < 0.06) or (age > 17.8):
+            print("Notice: Your input median age (age) is beyond the range of stellar template [0.06, 17.8], which will be automatically adjusted to the upper/lower limit.")
+        if (feh > 0.22):
+            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):
+        #    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
@@ -958,6 +1000,9 @@ class StellarContinuum():
             idx = (sigma_gal <= sigma_LSF)
             sigma_dif[idx] = 0.1
             flux0 = gaussian_filter1d(Stellar, sigma_dif)
+        else:
+            if sigma_gal < 0:
+                print('Notice: Your input velocity dispersion (vdisp) is < 0 km/s, which will be automatically adjusted to 0 km/s. ')
         
         # Dust Reddening
         if np.isscalar(ebv):
@@ -1068,6 +1113,17 @@ class SingleStar():
     """
 
     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):
+            print("Notice: Your input magnitude (mag) is > 26 mag or < 8 mag.")
+        if (teff < 2820) or (teff > 20129):
+            print("Notice: Your input effective tempreture (teff) is beyond the range of stellar template [2820, 20129], which will be automatically adjusted to the upper/lower limit.")
+        if (feh > 0.81):
+            print("Notice: Your input metallicity (feh) is beyond the range of stellar template [-2.69, 0.81], which will be automatically adjusted to the upper limit.")
+        if (feh < -2.69):
+            raise Exception('Your input metallicity (feh) is beyond the range of stellar template [-2.69, 0.81]')
+        
     
         StarTemp = template.templates