minor

mtakahiro · Dec 15, 2022 · f11ad66 · f11ad66
1 parent a1bc2bf
commit f11ad66
Show file tree

Hide file tree

Showing 2 changed files with 39 additions and 14 deletions.
diff --git a/gsf/function.py b/gsf/function.py
@@ -489,16 +489,13 @@ def get_Fint(lmtmp, ftmp, lmin=1400, lmax=1500):
 
 
 def get_Fuv(lmtmp, ftmp, lmin=1400, lmax=1500):
-    '''
-    Purpose
-    -------
-    Get RF UV (or any wavelength) flux density.
+    '''Get RF UV (or any wavelength) flux density.
 
     Parameters
     ----------
-    lmtmp : 
+    lmtmp : float array
         Rest-frame wavelength, in AA.
-    ftmp : 
+    ftmp : float array
         Fnu
     
     Returns
@@ -548,12 +545,12 @@ def fnutonu(fnu, m0set=25.0, m0input=-48.6):
     return fnu_new
 
 
-def flamtonu(lam, flam, m0set=25.0):
+def flamtonu(lam, flam, m0set=25.0, m0=-48.6):
     '''
     Converts from Flam to Fnu, with mag zeropoint of m0set.
     
     '''
-    Ctmp = lam**2/c * 10**((48.6+m0set)/2.5) #/ delx_org
+    Ctmp = lam**2/c * 10**((m0set-m0)/2.5) #/ delx_org
     fnu  = flam * Ctmp
     return fnu
 

diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py
@@ -594,13 +594,45 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000,
                     ms[:] = mshdu['ms_'+str(zz)][:] # [:] is necessary.
                     Ls[:] = mshdu['Ls_'+str(zz)][:]
 
+                    # Distance;
+                    DL = MB.cosmo.luminosity_distance(zbest).value * MB.Mpc_cm # Luminositydistance in cm
+                    DL10 = MB.Mpc_cm / 1e6 * 10 # 10pc in cm
+
                     for ss in range(Na):
                         wave = lm0
+                        wavetmp = wave*(1.+zbest)
                         if delwave>0:
                             fint = interpolate.interp1d(lm0_orig, spechdu['fspec_'+str(zz)+'_'+str(ss)+'_'+str(pp)][::nthin], kind='nearest', fill_value="extrapolate")
-                            spec_mul[ss] = fint(lm0)
+                            spec_mul[ss,:] = fint(lm0)
                         else:
-                            spec_mul[ss] = spechdu['fspec_'+str(zz)+'_'+str(ss)+'_'+str(pp)][::nthin] # Lsun/A
+                            spec_mul[ss,:] = spechdu['fspec_'+str(zz)+'_'+str(ss)+'_'+str(pp)][::nthin] # Lsun/A
+
+                        """
+                        ###################
+                        # Check xi_ion
+                        ###################
+                        h = 6.626e-34 # J s
+                        # 1 J = 1e7 erg
+                        h *= 1e7 # erg s
+                        nu = c / wave # A/s / A = 1/s
+                        # con_lyc = (wave[:-1]<912.0) & (wave[:-1]>228.0)
+                        # nph = spec_mul[ss,:] / (h * nu) # Lsun/A / (erg s * 1/s)
+                        # nph *= MB.Lsun # 1/A
+                        # delwave_array = np.diff(wave)
+                        # nph_Lyc = np.nansum(nph[:-1][con_lyc]*delwave_array[con_lyc])
+                        con_lyc = (wave[:]<912.0) & (wave[:]>228.0)
+                        nph = spec_mul[ss,:] / (h * nu) # Lsun/s/A / (erg s * 1/s) = Lsun/s / A / erg
+                        nph *= MB.Lsun # 1/s/A
+                        nph_Lyc = np.nansum(nph[:][con_lyc]) # 1/s
+
+                        # UV Flux density;
+                        fnu = flamtonu(wave, spec_mul[ss,:], m0set=MB.m0set) * MB.Lsun # erg/A
+                         #/ (4. * np.pi * DL10**2) # fnu, in erg/s/cm2/Hz.
+                        Fuv = get_Fuv(wave, fnu, lmin=1250, lmax=1650) # erg/A
+                        # MUV = -2.5 * np.log10(Fuv) + MB.m0set 
+                        Luv = Fuv * (1650-1250) # erg
+                        xi_ion = nph_Lyc / Luv # 1/s / (erg)
+                        """
 
                         ###################
                         # IGM attenuation.
@@ -611,10 +643,6 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000,
                             MB.x_HI = x_HI
                             spec_mul[ss,:] = spec_av_tmp
 
-                        # Distance;
-                        DL = MB.cosmo.luminosity_distance(zbest).value * MB.Mpc_cm # Luminositydistance in cm
-                        wavetmp = wave*(1.+zbest)
-
                         # Flam to Fnu
                         spec_mul_nu[ss,:] = flamtonu(wave, spec_mul[ss,:], m0set=MB.m0set)
                         spec_mul_nu[ss,:] *= MB.Lsun / (4.*np.pi*DL**2/(1.+zbest))