Merge pull request #53 from mtakahiro/gsf-bpass

Gsf bpass

gsf/fitting.py

@@ -137,9 +137,9 @@ def __init__(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set:floa
                       'CAT_BB', 'CAT_BB_DUST', 'SNLIM',
                       'MORP', 'MORP_FILE', 
                       'SPEC_FILE', 'DIR_EXTR', 'MAGZP_SPEC', 'UNIT_SPEC', 'DIFF_CONV', 
-                      'CZ0', 'CZ1', 'CZ2', 'LINE', ],
+                      'CZ0', 'CZ1', 'CZ2', 'LINE', 'PA', ],
 
-            'Misc' : ['CONFIG', 'DIR_OUT', 'FILTER', 'SKIPFILT', 'FIR_FILTER']
+            'Misc' : ['CONFIG', 'DIR_OUT', 'FILTER', 'SKIPFILT', 'FIR_FILTER', 'DIR_FILT']
 
             }
 
@@ -600,7 +600,7 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set:
             self.fzmc = int(inputs['F_ZMC'])
         except:
             self.fzmc = 0
-            self.logger.warning('Cannot find ZMC. Set to %d.'%(self.fzmc))
+            self.logger.warning('Cannot find F_ZMC. Set to %d.'%(self.fzmc))
 
         # Metallicity
         self.has_ZFIX = False
@@ -667,6 +667,12 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set:
                 self.Zmax, self.Zmin = float(inputs['ZMAX']), float(inputs['ZMIN'])
                 con_z = np.where((Zbpass >= self.Zmin) & (Zbpass <= self.Zmax))
                 self.Zall = Zbpass[con_z]
+
+                if len(self.Zall) == 0:
+                    self.logger.warning('No metallicity can be found. Available Zs are:')
+                    self.logger.info(Zbpass)
+                    sys.exit()
+
                 self.delZ = 0.0001
                 self.Zmax,self.Zmin = np.max(self.Zall), np.min(self.Zall)
                 self.logger.info('Final list for log(Z_BPASS/Zsun) is : ')
@@ -785,15 +791,16 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set:
             elif self.dust_model == 4:
                 self.dust_model_name = 'KriekConroy'
             else:
-                self.logger.warning('Unknown index number for dust attenuation. Setting to Calzetti.')
+                self.logger.warning('Unknown index number for dust model. Setting to Calzetti.')
                 self.dust_model = 0
                 self.dust_model_name = 'Calz'
         except:
+            self.logger.warning('Index number for dust model (DUST_MODEL) cannot be found. Setting to Calzetti.')
             self.dust_model = 0
             self.dust_model_name = 'Calz'
 
         if self.verbose:
-            self.logger.info('Dust attenuation is set to %s\n'%self.dust_model_name)
+            self.logger.info('Dust model is set to %s\n'%self.dust_model_name)
 
         # If FIR data;
         try:
@@ -1145,7 +1152,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, #zliml=0.01, zlim
         x_cz = self.dict['x'][con_cz] # Observed range
         NR_cz = self.dict['NR'][con_cz]
         if len(NR_cz) == 0:
-            self.logger.error('No data point exists at SN>%.1f'%(snlim))
+            self.logger.error('No data point exists at SNR > `snlim` (%.1f)'%(snlim))
             self.logger.error('Decrease `snlim` or turn `f_exclude_negative` to False')
             return False
 
@@ -2408,6 +2415,7 @@ def search_redshift(self, dict, xm_tmp, fm_tmp, zliml=0.01, zlimu=6.0, delzz=0.0
             fit_par_cz.add('C%d'%nn, value=1., min=0., max=1e5)
 
         def residual_z(pars,z):
+            ''''''
             vals  = pars.valuesdict()
             xm_s = xm_tmp * (1+z)
             fm_s = np.zeros(len(xm_tmp),float)
@@ -2416,7 +2424,6 @@ def residual_z(pars,z):
                 fm_s += fm_tmp[nn,:] * pars['C%d'%nn]
 
             fint = interpolate.interp1d(xm_s, fm_s, kind='nearest', fill_value="extrapolate")
-            #fm_int = np.interp(xobs, xm_s, fm_s)
             fm_int = fint(xobs)
 
             if fcon is None:
@@ -2427,7 +2434,6 @@ def residual_z(pars,z):
 
         # Start redshift search;
         for zz in range(len(zspace)):
-            # Best fit
             out_cz = minimize(residual_z, fit_par_cz, args=([zspace[zz]]), method=method)
             keys = fit_report(out_cz).split('\n')
 

gsf/function.py

@@ -739,24 +739,22 @@ def apply_dust(yy, xx, nr, Av, dust_model=0):
 	elif dust_model == 2: # LMC
 		yyd, xxd, nrd = dust_gen(xx, yy, Av, nr, Rv=4.05, gamma=-0.06, Eb=2.8)
 	elif dust_model == 3: # SMC
-		yyd, xxd, nrd = dust_smc(xx, yy, Av, nr, Rv=2.74, x0=4.703, gamma=1.212, lmlimu=3.115, f_Alam=False)
+		yyd, xxd, nrd = dust_smc(xx, yy, Av, nr, Rv=2.74, x0=4.703, gamma=1.212, f_Alam=False)
 	elif dust_model == 4: # Kriek&Conroy with gamma=-0.2
 		yyd, xxd, nrd = dust_kc(xx, yy, Av, nr, Rv=4.05, gamma=-0.2)
 	else:
 		yyd, xxd, nrd = dust_calz(xx, yy, Av, nr)
 	return yyd, xxd, nrd
 
 
-def dust_smc(lm, fl, Av, nr, Rv=2.74, x0=4.703, gamma=1.212, lmlimu=3.115, f_Alam=False):
+def dust_smc(lm, fl, Av, nr, Rv=2.74, x0=4.6, gamma=1.0, f_Alam=False):
     '''
-    For general purpose (Noll+09).
-    This function is much better than previous, but is hard to impliment for the current version.
-    A difference from dust_gen is Eb is defined as a function of gamma.
+    Dust law for SMC (Gordon+03).
 
     Parameters
     ----------
     lm : float array
-        wavelength, at RF.
+        wavelength in AA, at RF.
     fl : float array
         fnu
     Av : float
@@ -772,42 +770,42 @@ def dust_smc(lm, fl, Av, nr, Rv=2.74, x0=4.703, gamma=1.212, lmlimu=3.115, f_Ala
     # if any(np.diff(lm)<0):
     #     print('Something is wrong in lm: dust_smc of function.py')
 
-    lmm  = lm/10000. # in micron
-    con1 = (lmm<=0.63)
-    con2 = (lmm>0.63)  & (lmm<=lmlimu)
-    con3 = (lmm>lmlimu)
-
-    #nr0 = nr[con0]
-    nr1 = nr[con1]
-    nr2 = nr[con2]
-    nr3 = nr[con3]
-
-    #lmm0 = lmm[con0]
-    lmm1 = lmm[con1]
-    lmm2 = lmm[con2]
-    lmm3 = lmm[con3]
-
-    #fl0 = fl[con0]
-    fl1 = fl[con1]
-    fl2 = fl[con2]
-    fl3 = fl[con3]
-
-    nrd  = np.concatenate([nr1,nr2,nr3])
-    lmmc = np.concatenate([lmm1,lmm2,lmm3])
-    flc  = np.concatenate([fl1,fl2,fl3])
+    lmm  = lm/10000. # into micron
+    nrd = nr #np.concatenate([nr1,nr2,nr3])
+    lmmc = lmm #np.concatenate([lmm1,lmm2,lmm3])
+    flc = fl #np.concatenate([fl1,fl2,fl3])
 
     # Using average for SMC;
-    c1,c2,c3,c4 = -0.856, 1.038, 3.215, 0.107
+    c1,c2,c3,c4 = -4.959, 2.264, 0.389, 0.461
+    # SMC Wing Sample;
+    # c1,c2,c3,c4 = -0.856, 1.038, 3.215, 0.107
+
     x = 1./lmmc
-    Dx = x**2 / ((x**2-x0**2)**2 +x**2*gamma**2)
-    Fx = 0.5329 * (x - 5.9)**2 + 0.05644 * (x-5.9)**3
+    Dx = x**2 / ((x**2-x0**2)**2 + x**2*gamma**2)
+    Fx = 0.5392 * (x - 5.9)**2 + 0.05644 * (x-5.9)**3
     con_fx = (x<5.9)
     Fx[con_fx] = 0
 
-    EBlam_to_EB = c1+c2*x+c3*Dx+c4*Fx
+    EBlam_to_EB = c1 + c2*x + c3*Dx + c4*Fx
     Alam = Av / Rv * EBlam_to_EB
     fl_cor = flc[:] * 10**(-0.4*Alam[:])
 
+    if False:
+        import matplotlib.pyplot as plt
+        plt.close()
+        xs = np.arange(0.5, 10, 0.01)
+        x = xs
+        Dx = x**2 / ((x**2-x0**2)**2 + x**2*gamma**2)
+        Fx = 0.5392 * (x - 5.9)**2 + 0.05644 * (x-5.9)**3
+        con_fx = (x<5.9)
+        Fx[con_fx] = 0
+        EBlam_to_EB = c1 + c2*x + c3*Dx + c4*Fx
+        Av = 2.0
+        Alam = Av / Rv * EBlam_to_EB
+        plt.scatter(x, EBlam_to_EB, )
+        plt.show()
+        hoge
+
     if f_Alam:
         return fl_cor, lmmc*10000., nrd, Alam
     else:

gsf/maketmp_filt.py

@@ -1762,7 +1762,10 @@ def show_template_parameters(af):
     print('nimf     :',af['nimf'])
     print('age/Gyr  :',af['age'][:])
     print('logZ/Zsun:',af['Z'][:])
-    print('logU     :',np.arange(af['logUMIN'], af['logUMAX']+0.01, af['DELlogU']))
+    try:
+        print('logU     :',np.arange(af['logUMIN'], af['logUMAX']+0.01, af['DELlogU']))
+    except:
+        pass
     print('\n')
     return 
 

gsf/maketmp_z0.py

@@ -553,6 +553,17 @@ def make_tmp_z0_bpass(MB, lammin=100, lammax=160000, Zforce=None, Zsun=0.02):
                         'version_gsf': gsf.__version__
                     }
 
+                    tree.update({'age': MB.age})
+                    tree.update({'Z': MB.Zall})
+                    if fneb:
+                        tree.update({'logUMIN': MB.logUMIN})
+                        tree.update({'logUMAX': MB.logUMAX})
+                        tree.update({'DELlogU': MB.DELlogU})
+                    if MB.fagn:
+                        tree.update({'AGNTAUMIN': MB.AGNTAUMIN})
+                        tree.update({'AGNTAUMAX': MB.AGNTAUMAX})
+                        tree.update({'DELAGNTAU': MB.DELAGNTAU})
+
                     # ASDF
                     tree_spec.update({'wavelength': wave})
                     flagz = True
@@ -579,6 +590,7 @@ def make_tmp_z0_bpass(MB, lammin=100, lammax=160000, Zforce=None, Zsun=0.02):
                 col4 = fits.Column(name='tau_'+str(zz), format='E', unit='Gyr', array=tau_age)
                 tree_ML.update({'realtau_'+str(zz): ms})
 
+
     # Write;
     for aa in range(len(age)):
         tree.update({'realtau%d(Gyr)'%(aa): tau_age[aa]})