From 9b7bacdab7ac92524c3aaa1798d4418abb145e04 Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Sun, 19 Jun 2022 16:47:06 -0700 Subject: [PATCH 01/14] grism plot --- gsf/maketmp_filt.py | 37 +++++++++++++++++-------------------- gsf/plot_sed.py | 14 +++++++++----- gsf/zfit.py | 2 ++ 3 files changed, 28 insertions(+), 25 deletions(-) diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py index 395e7a1..d7462ea 100644 --- a/gsf/maketmp_filt.py +++ b/gsf/maketmp_filt.py @@ -187,27 +187,24 @@ def get_LSF(inputs, DIR_EXTR, ID, lm, c=3e18): ''' Amp = 0 f_morp = False - try: - if inputs['MORP'] == 'moffat' or inputs['MORP'] == 'gauss': - f_morp = True - try: - mor_file = inputs['MORP_FILE'].replace('$ID','%s'%(ID)) - fm = ascii.read(DIR_EXTR + mor_file) - Amp = fm['A'] - gamma = fm['gamma'] - if inputs['MORP'] == 'moffat': - alp = fm['alp'] - else: - alp = 0 - except Exception: - print('Error in reading morphology params.') - print('No morphology convolution.') - pass - else: - print('MORP Keywords does not match.') + if inputs['MORP'] == 'moffat' or inputs['MORP'] == 'gauss': + f_morp = True + try: + mor_file = inputs['MORP_FILE'].replace('$ID','%s'%(ID)) + fm = ascii.read(DIR_EXTR + mor_file) + Amp = fm['A'] + gamma = fm['gamma'] + if inputs['MORP'] == 'moffat': + alp = fm['alp'] + else: + alp = 0 + except Exception: + print('Error in reading morphology params.') print('No morphology convolution.') - except: - pass + pass + else: + print('MORP Keywords does not match.') + print('No morphology convolution.') ############################ # Template convolution; diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py index ae50467..9a8df1c 100644 --- a/gsf/plot_sed.py +++ b/gsf/plot_sed.py @@ -52,11 +52,11 @@ def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=Fal from astropy.io import ascii import time - #if False:#f_silence: - # import matplotlib - # matplotlib.use("Agg") - #else: - # matplotlib.use("MacOSX") + if False:#f_silence: + import matplotlib + matplotlib.use("Agg") + else: + matplotlib.use("MacOSX") def gaus(x,a,x0,sigma): return a*exp(-(x-x0)**2/(2*sigma**2)) @@ -263,6 +263,10 @@ def gaus(x,a,x0,sigma): ax2t = ax1.inset_axes((1-xsize-0.01,1-ysize-0.01,xsize,ysize)) if MB.f_dust: ax3t = ax1.inset_axes((0.7,.35,.28,.25)) + + f_plot_resid = False + print('Grism data. f_plot_resid is turned off.') + else: if f_plot_resid: fig_mosaic = """ diff --git a/gsf/zfit.py b/gsf/zfit.py index 21f5e73..c2e60c1 100644 --- a/gsf/zfit.py +++ b/gsf/zfit.py @@ -93,6 +93,8 @@ def residual_z(pars): ############################### def lnprob_cz(pars): + ''' + ''' resid = residual_z(pars) # i.e. (data - model) * wht z = pars['z'] s_z = 1 #pars['f_cz'] From 45aa2aeca2cd088494b191a4fbca8083a547bb54 Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Sun, 19 Jun 2022 16:57:58 -0700 Subject: [PATCH 02/14] grism cycle --- gsf/fitting.py | 3 +-- gsf/maketmp_filt.py | 47 ++++++++++++++++++++++----------------------- 2 files changed, 24 insertions(+), 26 deletions(-) diff --git a/gsf/fitting.py b/gsf/fitting.py index ceb91ee..47f8af7 100644 --- a/gsf/fitting.py +++ b/gsf/fitting.py @@ -1932,14 +1932,13 @@ def __init__(self, flatchain, var_names, params_value, res): print('\n\n') flag_gen = raw_input('Do you want to make templates with recommended redshift, Cz0, and Cz1 , %.5f %.5f %.5f? ([y]/n) '%(self.zrecom, self.Czrec0, self.Czrec1)) if flag_gen == 'y' or flag_gen == '': - self.zprev = self.zgal # Input redshift for previous run + self.zprev = self.zgal # Input redshift for previous run self.zgal = self.zrecom # Recommended redshift from previous run self.Cz0 = self.Czrec0 self.Cz1 = self.Czrec1 return True else: print('\n\n') - print('There is nothing to do.') print('Terminating process.') print('\n\n') return False diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py index d7462ea..05fa46a 100644 --- a/gsf/maketmp_filt.py +++ b/gsf/maketmp_filt.py @@ -408,7 +408,7 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, lm = np.zeros(np.sum(ninp0[:]),dtype='float') fobs = np.zeros(np.sum(ninp0[:]),dtype='float') eobs = np.zeros(np.sum(ninp0[:]),dtype='float') - fgrs = np.zeros(np.sum(ninp0[:]),dtype='int') # FLAG for G102/G141. + fgrs = np.zeros(np.sum(ninp0[:]),dtype='int') # FLAG for each grism. for ff, spec_file in enumerate(spec_files): try: fd0 = np.loadtxt(DIR_EXTR + spec_file, comments='#') @@ -835,17 +835,13 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, ########################################## fw = open(DIR_TMP + 'spec_obs_' + ID + '.cat', 'w') fw.write('# BB data (>%d) in this file are not used in fitting.\n'%(ncolbb)) + for ii in range(len(lm)): - if fgrs[ii]==0: # G102 - if lm[ii]/(1.+zbest) > lamliml and lm[ii]/(1.+zbest) < lamlimu: - fw.write('%d %.5f %.5e %.5e\n'%(ii, lm[ii], fobs[ii], eobs[ii])) - else: - fw.write('%d %.5f 0 1000\n'%(ii, lm[ii])) - elif fgrs[ii]==1: # G141 - if lm[ii]/(1.+zbest) > lamliml and lm[ii]/(1.+zbest) < lamlimu: - fw.write('%d %.5f %.5e %.5e\n'%(ii+1000, lm[ii], fobs[ii], eobs[ii])) - else: - fw.write('%d %.5f 0 1000\n'%(ii+1000, lm[ii])) + g_offset = 1000 * fgrs[ii] + if lm[ii]/(1.+zbest) > lamliml and lm[ii]/(1.+zbest) < lamlimu: + fw.write('%d %.5f %.5e %.5e\n'%(ii+g_offset, lm[ii], fobs[ii], eobs[ii])) + else: + fw.write('%d %.5f 0 1000\n'%(ii+g_offset, lm[ii])) for ii in range(len(ltmpbb[0,:])): if SFILT[ii] in SKIPFILT:# data point to be skiped; @@ -1012,8 +1008,8 @@ def maketemp_tau(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tau_ ninp0 = np.zeros(len(spec_files), dtype='int') for ff, spec_file in enumerate(spec_files): try: - fd0 = np.loadtxt(DIR_EXTR + spec_file, comments='#') - lm0tmp= fd0[:,0] + fd0 = np.loadtxt(DIR_EXTR + spec_file, comments='#') + lm0tmp = fd0[:,0] fobs0 = fd0[:,1] eobs0 = fd0[:,2] ninp0[ff] = len(lm0tmp)#[con_tmp]) @@ -1021,14 +1017,14 @@ def maketemp_tau(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tau_ print('File, %s/%s, cannot be open.'%(DIR_EXTR,spec_file)) pass # Constructing arrays. - lm = np.zeros(np.sum(ninp0[:]),dtype='float') + lm = np.zeros(np.sum(ninp0[:]),dtype='float') fobs = np.zeros(np.sum(ninp0[:]),dtype='float') eobs = np.zeros(np.sum(ninp0[:]),dtype='float') - fgrs = np.zeros(np.sum(ninp0[:]),dtype='int') # FLAG for G102/G141. + fgrs = np.zeros(np.sum(ninp0[:]),dtype='int') # FLAG for each grism. for ff, spec_file in enumerate(spec_files): try: - fd0 = np.loadtxt(DIR_EXTR + spec_file, comments='#') - lm0tmp= fd0[:,0] + fd0 = np.loadtxt(DIR_EXTR + spec_file, comments='#') + lm0tmp = fd0[:,0] fobs0 = fd0[:,1] eobs0 = fd0[:,2] for ii1 in range(ninp0[ff]): @@ -1037,7 +1033,7 @@ def maketemp_tau(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tau_ else: ii = ii1 + np.sum(ninp0[:ff]) fgrs[ii] = ff - lm[ii] = lm0tmp[ii1] + lm[ii] = lm0tmp[ii1] fobs[ii] = fobs0[ii1] eobs[ii] = eobs0[ii1] f_spec = True @@ -1474,17 +1470,20 @@ def maketemp_tau(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tau_ ########################################## fw = open(DIR_TMP + 'spec_obs_' + ID + '.cat', 'w') fw.write('# BB data (>%d) in this file are not used in fitting.\n'%(ncolbb)) + for ii in range(len(lm)): - if fgrs[ii]==0: # G102 - if lm[ii]/(1.+zbest) > lamliml and lm[ii]/(1.+zbest) < lamlimu: - fw.write('%d %.5f %.5e %.5e\n'%(ii, lm[ii], fobs[ii], eobs[ii])) - else: - fw.write('%d %.5f 0 1000\n'%(ii, lm[ii])) - elif fgrs[ii]==1: # G141 + g_offset = 1000 * fgrs[ii] + if lm[ii]/(1.+zbest) > lamliml and lm[ii]/(1.+zbest) < lamlimu: + fw.write('%d %.5f %.5e %.5e\n'%(ii+g_offset, lm[ii], fobs[ii], eobs[ii])) + else: + fw.write('%d %.5f 0 1000\n'%(ii+g_offset, lm[ii])) + ''' + elif fgrs[ii]==1: # grism 2 if lm[ii]/(1.+zbest) > lamliml and lm[ii]/(1.+zbest) < lamlimu: fw.write('%d %.5f %.5e %.5e\n'%(ii+1000, lm[ii], fobs[ii], eobs[ii])) else: fw.write('%d %.5f 0 1000\n'%(ii+1000, lm[ii])) + ''' for ii in range(len(ltmpbb[0,:])): if SFILT[ii] in SKIPFILT:# data point to be skiped; From 7607ff4cb8a29d0a171adb2e01756e3f1d111cd7 Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Mon, 20 Jun 2022 17:10:23 -0700 Subject: [PATCH 03/14] third grism --- docs/parameters.rst | 2 +- gsf/fitting.py | 80 ++++++++++++++------ gsf/plot_sed.py | 179 +++++++++++++++++++++++++------------------- gsf/writing.py | 10 ++- gsf/zfit.py | 28 ++++--- 5 files changed, 183 insertions(+), 116 deletions(-) diff --git a/docs/parameters.rst b/docs/parameters.rst index 6e13098..10b60ee 100644 --- a/docs/parameters.rst +++ b/docs/parameters.rst @@ -81,7 +81,7 @@ An example config file can be generated by executing the following command (TBD) - Comma separated list for spectral files. Path is relative to DIR_EXTR. * - MORP - str - - Profile shape, if grism spectra are provided. Will be used to convolve templates. "moffat", "gauss", or none. + - "moffat", "gauss", or none. Profile shape, if grism spectra are provided. Will be used to convolve templates. * - MORP_FILE - str - Ascii file for morphology parameters. diff --git a/gsf/fitting.py b/gsf/fitting.py index 47f8af7..f10d447 100644 --- a/gsf/fitting.py +++ b/gsf/fitting.py @@ -41,6 +41,8 @@ LW = [2800, 3347, 3727, 3799, 3836, 3869, 4102, 4341, 4861, 4960, 5008, 5175, 6563, 6717, 6731] fLW = np.zeros(len(LW), dtype='int') +NRbb_lim = 10000 # BB data is associated with ids greater than this number. + class Mainbody(): ''' The list of (possible) `Mainbody` attributes is given below: @@ -178,6 +180,7 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: # Scaling for grism; self.Cz0 = float(inputs['CZ0']) self.Cz1 = float(inputs['CZ1']) + self.Cz2 = float(inputs['CZ2']) try: self.DIR_EXTR = inputs['DIR_EXTR'] @@ -604,7 +607,7 @@ def get_lines(self, LW0): return LW, fLW - def read_data(self, Cz0:float, Cz1:float, zgal:float, add_fir:bool=False, idman=None): + def read_data(self, Cz0:float, Cz1:float, Cz2:float, zgal:float, add_fir:bool=False, idman=None): ''' Parameters ---------- @@ -621,7 +624,7 @@ def read_data(self, Cz0:float, Cz1:float, zgal:float, add_fir:bool=False, idman= ----- This modeule can be used for any SFHs. ''' - print('READ data with Cz0=%.2f, Cz0=%.2f, zgal=%.2f'%(Cz0, Cz1, zgal)) + print('READ data with Cz0=%.2f, Cz1=%.2f, Cz2=%.2f, zgal=%.2f'%(Cz0, Cz1, Cz2, zgal)) ############## # Spectrum @@ -636,10 +639,14 @@ def read_data(self, Cz0:float, Cz1:float, zgal:float, add_fir:bool=False, idman= xx0 = x[con0] fy0 = fy00[con0] * Cz0 ey0 = ey00[con0] * Cz0 - con1 = (NR>=1000) & (NR<10000) + con1 = (NR>=1000) & (NR<2000) xx1 = x[con1] fy1 = fy00[con1] * Cz1 ey1 = ey00[con1] * Cz1 + con2 = (NR>=2000) & (NRsnlim) + con_cz = (self.dict['NR']snlim) if len(self.dict['fy'][con_cz])==0: if f_bb_zfit: con_cz = (sn>snlim) @@ -957,10 +969,12 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, zlimu z_cz = np.percentile(res_cz.flatchain['z'], [16,50,84]) scl_cz0 = np.percentile(res_cz.flatchain['Cz0'], [16,50,84]) scl_cz1 = np.percentile(res_cz.flatchain['Cz1'], [16,50,84]) + scl_cz2 = np.percentile(res_cz.flatchain['Cz2'], [16,50,84]) zrecom = z_cz[1] Czrec0 = scl_cz0[1] Czrec1 = scl_cz1[1] + Czrec2 = scl_cz2[1] # Switch to peak redshift: # find minimum and maximum of xticks, so we know @@ -969,7 +983,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, zlimu xmin, xmax = self.zgal-0.2, self.zgal+0.2 lnspc = np.linspace(xmin, xmax, len(ser)) print('\n\n') - print('Recommended redshift, Cz0 and Cz1, %.5f %.5f %.5f, with chi2/nu=%.3f'%(zrecom, Czrec0, Czrec1, fitc_cz[1])) + print('Recommended redshift, Cz0, Cz1, and Cz2, %.5f %.5f %.5f %.5f, with chi2/nu=%.3f'%(zrecom, Czrec0, Czrec1, Czrec2, fitc_cz[1])) print('\n\n') fit_label = 'Proposed model' @@ -988,8 +1002,10 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, zlimu zrecom = z_cz[1] scl_cz0 = [1.,1.,1.] scl_cz1 = [1.,1.,1.] + scl_cz2 = [1.,1.,1.] Czrec0 = scl_cz0[1] Czrec1 = scl_cz1[1] + Czrec2 = scl_cz2[1] res_cz = None # If this label is being used, it means that the fit is failed. @@ -1044,8 +1060,6 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, zlimu if len(fm_tmp) == len(self.dict['xbb']): # BB only; data_obsbb[:,2] = fm_tmp data_obsbb_sort = np.sort(data_obsbb, axis=0) - # y-Scale seems to be wrong? - #plt.plot(data_obsbb_sort[:,0], data_obsbb_sort[:,1], marker='.', color='r', ms=10, linestyle='', linewidth=0, zorder=4, label='Obs.(BB)') if len(fm_tmp) == len(self.dict['xbb']): # BB only; plt.scatter(data_obsbb_sort[:,0], data_obsbb_sort[:,2], color='none', marker='d', s=50, edgecolor='gray', zorder=4, label='Current model ($z=%.5f$)'%(self.zgal)) @@ -1076,6 +1090,8 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, zlimu eC0sigma = ((scl_cz0[2]-scl_cz0[0])/2.)/self.Cz0 C1sigma = np.abs(Czrec1-self.Cz1)/self.Cz1 eC1sigma = ((scl_cz1[2]-scl_cz1[0])/2.)/self.Cz1 + C2sigma = np.abs(Czrec2-self.Cz2)/self.Cz2 + eC2sigma = ((scl_cz2[2]-scl_cz2[0])/2.)/self.Cz2 print('\n##############################################################') print('Input redshift is %.3f per cent agreement.'%((1.-zsigma)*100)) @@ -1084,10 +1100,13 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, zlimu print('Error is %.3f per cent.'%(eC0sigma*100)) print('Input Cz1 is %.3f per cent agreement.'%((1.-C1sigma)*100)) print('Error is %.3f per cent.'%(eC1sigma*100)) + print('Input Cz2 is %.3f per cent agreement.'%((1.-C2sigma)*100)) + print('Error is %.3f per cent.'%(eC2sigma*100)) print('##############################################################\n') plt.show() - flag_z = raw_input('Do you want to continue with the input redshift, Cz0 and Cz1, %.5f %.5f %.5f? ([y]/n/m) '%(self.zgal, self.Cz0, self.Cz1)) + flag_z = raw_input('Do you want to continue with the input redshift, Cz0, Cz1 and Cz2, %.5f %.5f %.5f %.5f? ([y]/n/m) '%\ + (self.zgal, self.Cz0, self.Cz1, self.Cz2)) else: flag_z = 'y' @@ -1095,9 +1114,11 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, zlimu self.zrecom = zrecom self.Czrec0 = Czrec0 * self.Cz0 self.Czrec1 = Czrec1 * self.Cz1 + self.Czrec2 = Czrec2 * self.Cz2 self.z_cz = z_cz self.scl_cz0 = scl_cz0 self.scl_cz1 = scl_cz1 + self.scl_cz2 = scl_cz2 self.res_cz = res_cz return flag_z @@ -1122,7 +1143,9 @@ def get_zdist(self, f_interact=False): yy = np.arange(0,np.max(n),1) xx = yy * 0 + self.z_cz[1] ax1.plot(xx,yy,linestyle='-',linewidth=1,color='orangered',\ - label='$z=%.5f_{-%.5f}^{+%.5f}$\n$C_z0=%.3f$\n$C_z1=%.3f$'%(self.z_cz[1],self.z_cz[1]-self.z_cz[0],self.z_cz[2]-self.z_cz[1], self.Cz0, self.Cz1)) + label='$z=%.5f_{-%.5f}^{+%.5f}$\n$C_z0=%.3f$\n$C_z1=%.3f$\n$C_z2=%.3f$'%\ + (self.z_cz[1],self.z_cz[1]-self.z_cz[0],self.z_cz[2]-self.z_cz[1], self.Cz0, self.Cz1, self.Cz2)) + xx = yy * 0 + self.z_cz[0] ax1.plot(xx,yy,linestyle='--',linewidth=1,color='orangered') xx = yy * 0 + self.z_cz[2] @@ -1183,7 +1206,7 @@ def add_param(self, fit_params, sigz=1.0, zmin=None, zmax=None): fit_params.add('MDUST', value=9, min=0, max=15) self.ndim += 1 - self.dict = self.read_data(self.Cz0, self.Cz1, self.zgal, add_fir=self.f_dust) + self.dict = self.read_data(self.Cz0, self.Cz1, self.Cz2, self.zgal, add_fir=self.f_dust) f_add = True # Nebular; ver1.6 @@ -1410,7 +1433,7 @@ def prepare_class(self, add_fir=None): # # Observed Data; # - self.dict = self.read_data(self.Cz0, self.Cz1, self.zgal, add_fir=add_fir) + self.dict = self.read_data(self.Cz0, self.Cz1, self.Cz2, self.zgal, add_fir=add_fir) # Set parameters; self.set_param() @@ -1919,10 +1942,17 @@ def __init__(self, flatchain, var_names, params_value, res): else: Czrec1 = self.Cz1 + Czrec2 = raw_input('What is your manual input for Cz2? [%.3f] '%(self.Cz2)) + if Czrec2 != '': + Czrec2 = float(Czrec2) + else: + Czrec2 = self.Cz2 + self.zprev = self.zgal # Input redshift for previous run self.zgal = zrecom # Recommended redshift from previous run self.Cz0 = Czrec0 self.Cz1 = Czrec1 + self.Cz2 = Czrec2 print('\n\n') print('Generating model templates with input redshift and Scale.') print('\n\n') @@ -1930,12 +1960,14 @@ def __init__(self, flatchain, var_names, params_value, res): else: print('\n\n') - flag_gen = raw_input('Do you want to make templates with recommended redshift, Cz0, and Cz1 , %.5f %.5f %.5f? ([y]/n) '%(self.zrecom, self.Czrec0, self.Czrec1)) + flag_gen = raw_input('Do you want to make templates with recommended redshift, Cz0, Cz1, and Cz2 , %.5f %.5f %.5f %.5f? ([y]/n) '%\ + (self.zrecom, self.Czrec0, self.Czrec1, self.Czrec2)) if flag_gen == 'y' or flag_gen == '': self.zprev = self.zgal # Input redshift for previous run self.zgal = self.zrecom # Recommended redshift from previous run self.Cz0 = self.Czrec0 self.Cz1 = self.Czrec1 + self.Cz2 = self.Czrec2 return True else: print('\n\n') diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py index 9a8df1c..109d42c 100644 --- a/gsf/plot_sed.py +++ b/gsf/plot_sed.py @@ -19,7 +19,7 @@ def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, mmax=300, dust_model=0, DIR_TMP='./templates/', f_label=False, f_bbbox=False, verbose=False, f_silence=True, - f_fill=False, f_fancyplot=False, f_Alog=True, dpi=300, f_plot_filter=True, f_plot_resid=False): + f_fill=False, f_fancyplot=False, f_Alog=True, dpi=300, f_plot_filter=True, f_plot_resid=False, NRbb_lim=10000): ''' Parameters ---------- @@ -182,6 +182,7 @@ def gaus(x,a,x0,sigma): Cz0 = hdul[0].header['Cz0'] Cz1 = hdul[0].header['Cz1'] + Cz2 = hdul[0].header['Cz2'] zbes = zp50 zscl = (1.+zbes) @@ -189,9 +190,9 @@ def gaus(x,a,x0,sigma): # Data taken from ############################### if MB.f_dust: - MB.dict = MB.read_data(Cz0, Cz1, zbes, add_fir=True) + MB.dict = MB.read_data(Cz0, Cz1, Cz2, zbes, add_fir=True) else: - MB.dict = MB.read_data(Cz0, Cz1, zbes) + MB.dict = MB.read_data(Cz0, Cz1, Cz2, zbes) NR = MB.dict['NR'] x = MB.dict['x'] @@ -202,11 +203,16 @@ def gaus(x,a,x0,sigma): xg0 = x[con0] fg0 = fy[con0] eg0 = ey[con0] - con1 = (NR>=1000) & (NR<10000) + con1 = (NR>=1000) & (NR<2000) xg1 = x[con1] fg1 = fy[con1] eg1 = ey[con1] - if len(xg0)>0 or len(xg1)>0: + con2 = (NR>=1000) & (NR0 or len(xg1)>0 or len(xg2)>0: f_grsm = True else: f_grsm = False @@ -576,14 +582,13 @@ def gaus(x,a,x0,sigma): ########################## if f_grsm: conspec = (NR<10000) #& (fy/ey>1) - #ax2t.fill_between(xg1, (fg1-eg1) * c/np.square(xg1)/d, (fg1+eg1) * c/np.square(xg1)/d, lw=0, color='#DF4E00', zorder=10, alpha=0.7, label='') - #ax2t.fill_between(xg0, (fg0-eg0) * c/np.square(xg0)/d, (fg0+eg0) * c/np.square(xg0)/d, lw=0, color='royalblue', zorder=10, alpha=0.2, label='') - ax2t.errorbar(xg1, fg1 * c/np.square(xg1)/d, yerr=eg1 * c/np.square(xg1)/d, lw=0.5, color='#DF4E00', zorder=10, alpha=1., label='', capsize=0) + ax2t.errorbar(xg2, fg2 * c/np.square(xg2)/d, yerr=eg2 * c/np.square(xg2)/d, lw=0.5, color='#DF4E00', zorder=10, alpha=1., label='', capsize=0) + ax2t.errorbar(xg1, fg1 * c/np.square(xg1)/d, yerr=eg1 * c/np.square(xg1)/d, lw=0.5, color='g', zorder=10, alpha=1., label='', capsize=0) ax2t.errorbar(xg0, fg0 * c/np.square(xg0)/d, yerr=eg0 * c/np.square(xg0)/d, lw=0.5, linestyle='', color='royalblue', zorder=10, alpha=1., label='', capsize=0) - xgrism = np.concatenate([xg0,xg1]) - fgrism = np.concatenate([fg0,fg1]) - egrism = np.concatenate([eg0,eg1]) + xgrism = np.concatenate([xg0,xg1,xg2]) + fgrism = np.concatenate([fg0,fg1,fg2]) + egrism = np.concatenate([eg0,eg1,eg2]) con4000b = (xgrism/zscl>3400) & (xgrism/zscl<3800) & (fgrism>0) & (egrism>0) con4000r = (xgrism/zscl>4200) & (xgrism/zscl<5000) & (fgrism>0) & (egrism>0) print('Median SN at 3400-3800 is;', np.median((fgrism/egrism)[con4000b])) @@ -1135,6 +1140,9 @@ def func_tmp(xint,eobs,fmodel): tree_spec.update({'fg1_obs': fg1 * c/np.square(xg1)/d}) tree_spec.update({'eg1_obs': eg1 * c/np.square(xg1)/d}) tree_spec.update({'wg1_obs': xg1}) + tree_spec.update({'fg2_obs': fg2 * c/np.square(xg2)/d}) + tree_spec.update({'eg2_obs': eg2 * c/np.square(xg2)/d}) + tree_spec.update({'wg2_obs': xg2}) af = asdf.AsdfFile(tree_spec) af.write_to(MB.DIR_OUT + 'gsf_spec_%s.asdf'%(ID), all_array_compression='zlib') @@ -1280,7 +1288,7 @@ def func_tmp(xint,eobs,fmodel): def plot_sed_tau(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, mmax=300, dust_model=0, DIR_TMP='./templates/', f_label=False, f_bbbox=False, verbose=False, f_silence=True, - f_fill=False, f_fancyplot=False, f_Alog=True, dpi=300, f_plot_filter=True, f_plot_resid=False): + f_fill=False, f_fancyplot=False, f_Alog=True, dpi=300, f_plot_filter=True, f_plot_resid=False, NRbb_lim=10000): ''' Parameters ---------- @@ -1481,11 +1489,15 @@ def gaus(x,a,x0,sigma): xg0 = x[con0] fg0 = fy[con0] #* Cz0 eg0 = ey[con0] #* Cz0 - con1 = (NR>=1000) & (NR<10000) #& (fy/ey>SNlim) + con1 = (NR>=1000) & (NR<2000) #& (fy/ey>SNlim) xg1 = x[con1] fg1 = fy[con1] #* Cz1 eg1 = ey[con1] #* Cz1 - if len(xg0)>0 or len(xg1)>0: + con2 = (NR>=2000) & (NRSNlim) + xg2 = x[con2] + fg2 = fy[con2] #* Cz1 + eg2 = ey[con2] #* Cz1 + if len(xg0)>0 or len(xg1)>0 or len(xg2)>0: f_grsm = True else: f_grsm = False @@ -1822,14 +1834,13 @@ def gaus(x,a,x0,sigma): ########################## if f_grsm: conspec = (NR<10000) #& (fy/ey>1) - #ax2t.fill_between(xg1, (fg1-eg1) * c/np.square(xg1)/d, (fg1+eg1) * c/np.square(xg1)/d, lw=0, color='#DF4E00', zorder=10, alpha=0.7, label='') - #ax2t.fill_between(xg0, (fg0-eg0) * c/np.square(xg0)/d, (fg0+eg0) * c/np.square(xg0)/d, lw=0, color='royalblue', zorder=10, alpha=0.2, label='') - ax2t.errorbar(xg1, fg1 * c/np.square(xg1)/d, yerr=eg1 * c/np.square(xg1)/d, lw=0.5, color='#DF4E00', zorder=10, alpha=1., label='', capsize=0) + ax2t.errorbar(xg2, fg2 * c/np.square(xg2)/d, yerr=eg2 * c/np.square(xg2)/d, lw=0.5, color='#DF4E00', zorder=10, alpha=1., label='', capsize=0) + ax2t.errorbar(xg1, fg1 * c/np.square(xg1)/d, yerr=eg1 * c/np.square(xg1)/d, lw=0.5, color='g', zorder=10, alpha=1., label='', capsize=0) ax2t.errorbar(xg0, fg0 * c/np.square(xg0)/d, yerr=eg0 * c/np.square(xg0)/d, lw=0.5, linestyle='', color='royalblue', zorder=10, alpha=1., label='', capsize=0) - xgrism = np.concatenate([xg0,xg1]) - fgrism = np.concatenate([fg0,fg1]) - egrism = np.concatenate([eg0,eg1]) + xgrism = np.concatenate([xg0,xg1,xg2]) + fgrism = np.concatenate([fg0,fg1,fg2]) + egrism = np.concatenate([eg0,eg1,eg1]) con4000b = (xgrism/zscl>3400) & (xgrism/zscl<3800) & (fgrism>0) & (egrism>0) con4000r = (xgrism/zscl>4200) & (xgrism/zscl<5000) & (fgrism>0) & (egrism>0) print('Median SN at 3400-3800 is;', np.median((fgrism/egrism)[con4000b])) @@ -2361,6 +2372,9 @@ def func_tmp(xint,eobs,fmodel): tree_spec.update({'fg1_obs': fg1 * c/np.square(xg1)/d}) tree_spec.update({'eg1_obs': eg1 * c/np.square(xg1)/d}) tree_spec.update({'wg1_obs': xg1}) + tree_spec.update({'fg2_obs': fg2 * c/np.square(xg2)/d}) + tree_spec.update({'eg2_obs': eg2 * c/np.square(xg2)/d}) + tree_spec.update({'wg2_obs': xg2}) af = asdf.AsdfFile(tree_spec) af.write_to(MB.DIR_OUT + 'gsf_spec_%s.asdf'%(ID), all_array_compression='zlib') @@ -2512,7 +2526,8 @@ def plot_filter(MB, ax, ymax, scl=0.3, cmap='gist_rainbow', alp=0.4, ind_remove= return ax -def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax=300, TMIN=0.0001, tau_lim=0.01, f_plot_filter=True, scale=1e-19): +def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax=300, TMIN=0.0001, tau_lim=0.01, f_plot_filter=True, + scale=1e-19, NRbb_lim=10000): ''' Purpose ------- @@ -2588,14 +2603,15 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax=30 Z84[aa] = hdul[1].data['Z'+str(aa)][2] #NZbest[aa]= bfnc.Z2NZ(Z50[aa]) - ZZ50 = np.sum(Z50*A50)/np.sum(A50) # Light weighted Z. - chi = hdul[1].data['chi'][0] - chin = hdul[1].data['chi'][1] - fitc = chin - Cz0 = hdul[0].header['Cz0'] - Cz1 = hdul[0].header['Cz1'] - zbes = hdul[0].header['z'] - zscl = (1.+zbes) + ZZ50 = np.sum(Z50*A50)/np.sum(A50) # Light weighted Z. + chi = hdul[1].data['chi'][0] + chin = hdul[1].data['chi'][1] + fitc = chin + Cz0 = hdul[0].header['Cz0'] + Cz1 = hdul[0].header['Cz1'] + Cz2 = hdul[0].header['Cz2'] + zbes = hdul[0].header['z'] + zscl = (1.+zbes) # plot Configuration if MB.fzmc == 1: @@ -2628,46 +2644,54 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax=30 ax2 = fig.add_axes([0.05,0.07,0.37,0.23]) if MB.f_dust: - MB.dict = MB.read_data(MB.Cz0, MB.Cz1, MB.zgal, add_fir=True) + MB.dict = MB.read_data(MB.Cz0, MB.Cz1, MB.Cz2, MB.zgal, add_fir=True) else: - MB.dict = MB.read_data(MB.Cz0, MB.Cz1, MB.zgal) + MB.dict = MB.read_data(MB.Cz0, MB.Cz1, MB.Cz2, MB.zgal) # Get data points; - NRbb = MB.dict['NRbb'] - xbb = MB.dict['xbb'] - fybb = MB.dict['fybb'] - eybb = MB.dict['eybb'] - exbb = MB.dict['exbb'] - snbb = fybb/eybb + NRbb = MB.dict['NRbb'] + xbb = MB.dict['xbb'] + fybb = MB.dict['fybb'] + eybb = MB.dict['eybb'] + exbb = MB.dict['exbb'] + snbb = fybb/eybb # Get spec data points; dat = np.loadtxt(DIR_TMP + 'spec_obs_' + ID + '.cat', comments='#') - NR = dat[:, 0] - x = dat[:, 1] - fy00 = dat[:, 2] - ey00 = dat[:, 3] + NR = dat[:,0] + x = dat[:,1] + fy00 = dat[:,2] + ey00 = dat[:,3] con0 = (NR<1000) #& (fy/ey>SNlim) xg0 = x[con0] fg0 = fy00[con0] * Cz0 eg0 = ey00[con0] * Cz0 - con1 = (NR>=1000) & (NR<10000) #& (fy/ey>SNlim) + con1 = (NR>=1000) & (NR<2000) #& (fy/ey>SNlim) xg1 = x[con1] fg1 = fy00[con1] * Cz1 eg1 = ey00[con1] * Cz1 - con2 = (NR>=10000)#& (fy/ey>SNlim) + con2 = (NR>=2000) & (NRSNlim) xg2 = x[con2] - fg2 = fy00[con2] - eg2 = ey00[con2] + fg2 = fy00[con2] * Cz2 + eg2 = ey00[con2] * Cz2 + + con_bb = (NR>=NRbb_lim) + xg_bb = x[con_bb] + fg_bb = fy00[con_bb] + eg_bb = ey00[con_bb] + fy01 = np.append(fg0,fg1) ey01 = np.append(eg0,eg1) + fy02 = np.append(fy01,fg2) + ey02 = np.append(ey01,eg2) - fy = np.append(fy01,fg2) - ey = np.append(ey01,eg2) + fy = np.append(fy01,fg_bb) + ey = np.append(ey01,eg_bb) wht = 1./np.square(ey) # BB photometry - conspec = (NR<10000) + conspec = (NR sigma) ax0.errorbar(xbb[conbb], fybb[conbb] * c / np.square(xbb[conbb]) / d, yerr=eybb[conbb]*c/np.square(xbb[conbb])/d, color='k', linestyle='', linewidth=0.5, zorder=4) @@ -2690,7 +2714,7 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax=30 burnin = data['burnin'] nmc = data['niter'] nwalk = data['nwalkers'] - Nburn = burnin #*20 + Nburn = burnin samples = data['chain'][:] except: if verbose: print(' = > NO keys of ndim and burnin found in cpkl, use input keyword values') @@ -3035,7 +3059,8 @@ def density_estimation(m1, m2): # For the last one ax0.plot(xg0, fg0 * c / np.square(xg0) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='royalblue', label='') - ax0.plot(xg1, fg1 * c / np.square(xg1) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='#DF4E00', label='') + ax0.plot(xg1, fg1 * c / np.square(xg1) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='g', label='') + ax0.plot(xg2, fg2 * c / np.square(xg2) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='#DF4E00', label='') ax0.set_xlim(2200, 88000) ax0.set_xscale('log') @@ -3093,7 +3118,7 @@ def density_estimation(m1, m2): def plot_corner_physparam_cumulative_frame(ID, Zall=np.arange(-1.2,0.4249,0.05), age=[0.01, 0.1, 0.3, 0.7, 1.0, 3.0], \ - tau0=[0.1,0.2,0.3], fig=None, dust_model=0, out_ind=0, snlimbb=1.0, DIR_OUT='./'): + tau0=[0.1,0.2,0.3], fig=None, dust_model=0, out_ind=0, snlimbb=1.0, DIR_OUT='./', NRbb_lim=10000): ''' Creat "cumulative" png for gif image. @@ -3159,6 +3184,7 @@ def plot_corner_physparam_cumulative_frame(ID, Zall=np.arange(-1.2,0.4249,0.05), fitc = chin Cz0 = hdul[0].header['Cz0'] Cz1 = hdul[0].header['Cz1'] + Cz2 = hdul[0].header['Cz2'] zbes = hdul[0].header['z'] zscl = (1.+zbes) @@ -3198,6 +3224,7 @@ def plot_corner_physparam_cumulative_frame(ID, Zall=np.arange(-1.2,0.4249,0.05), # For spec plot ax0 = fig.add_axes([0.62,0.61,0.37,0.33]) + ############################### # Data taken from ############################### @@ -3212,29 +3239,38 @@ def plot_corner_physparam_cumulative_frame(ID, Zall=np.arange(-1.2,0.4249,0.05), xg0 = x[con0] fg0 = fy00[con0] * Cz0 eg0 = ey00[con0] * Cz0 - con1 = (NR>=1000) & (NR<10000) #& (fy/ey>SNlim) + con1 = (NR>=1000) & (NR<2000) #& (fy/ey>SNlim) xg1 = x[con1] fg1 = fy00[con1] * Cz1 eg1 = ey00[con1] * Cz1 - con2 = (NR>=10000)#& (fy/ey>SNlim) + con2 = (NR>=2000) & (NRSNlim) xg2 = x[con2] - fg2 = fy00[con2] - eg2 = ey00[con2] + fg2 = fy00[con2] * Cz2 + eg2 = ey00[con2] * Cz2 + + con_bb = (NR>=NRbb_lim)#& (fy/ey>SNlim) + xg_bb = x[con_bb] + fg_bb = fy00[con_bb] + eg_bb = ey00[con_bb] + fy01 = np.append(fg0,fg1) - fy = np.append(fy01,fg2) ey01 = np.append(eg0,eg1) - ey = np.append(ey01,eg2) + fy02 = np.append(fy01,fg2) + ey02 = np.append(ey01,eg2) + + fy = np.append(fy02,fg_bb) + ey = np.append(ey02,eg_bb) wht=1./np.square(ey) dat = np.loadtxt(DIR_TMP + 'bb_obs_' + ID + '.cat', comments='#') - NRbb = dat[:, 0] - xbb = dat[:, 1] - fybb = dat[:, 2] - eybb = dat[:, 3] - exbb = dat[:, 4] + NRbb = dat[:,0] + xbb = dat[:,1] + fybb = dat[:,2] + eybb = dat[:,3] + exbb = dat[:,4] snbb = fybb/eybb - conspec = (NR<10000) #& (fy/ey>1) + conspec = (NR1) #ax0.plot(xg0, fg0 * c / np.square(xg0) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='royalblue', label='') #ax0.plot(xg1, fg1 * c / np.square(xg1) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='#DF4E00', label='') conbb = (fybb/eybb>snlimbb) @@ -3252,7 +3288,6 @@ def plot_corner_physparam_cumulative_frame(ID, Zall=np.arange(-1.2,0.4249,0.05), ax0.set_xscale('log') ax0.set_ylim(-0.05, ymax) - DIR_TMP = './templates/' #################### # MCMC corner plot. @@ -3323,7 +3358,7 @@ def plot_corner_physparam_cumulative_frame(ID, Zall=np.arange(-1.2,0.4249,0.05), ax0.plot(x0, y0 * c/ np.square(x0) / d, '--', lw=0.1, color=col[ii], zorder=-1, label='', alpha=0.1) else: y0_r, x0_tmp = fnc.tmp03(AA_tmp, Avtmp[kk], ii, ZZ_tmp, zbes, lib_all, tau0=tau0) - y0p, x0p = fnc.tmp03(AA_tmp, Avtmp[kk], ii, ZZ_tmp, zbes, lib, tau0=tau0) + y0p, x0p = fnc.tmp03(AA_tmp, Avtmp[kk], ii, ZZ_tmp, zbes, lib, tau0=tau0) ysump += y0p #* 1e18 ysum += y0_r #* 1e18 if AA_tmp/Asum > flim: @@ -3344,9 +3379,7 @@ def plot_corner_physparam_cumulative_frame(ID, Zall=np.arange(-1.2,0.4249,0.05), Ttmp[kk] = np.log10(Ttmp[kk]) - NPAR = [lmtmp[:kk+1], Ttmp[:kk+1], Avtmp[:kk+1], Ztmp[:kk+1]] - #NPARmin = [np.log10(M16)-0.1, -0.4, 0, -0.6] - #NPARmax = [np.log10(M84)+0.1, 0.5, 2., 0.5] + NPAR = [lmtmp[:kk+1], Ttmp[:kk+1], Avtmp[:kk+1], Ztmp[:kk+1]] NPARmin = [np.log10(M16)-0.1, -0.4, Av16-0.1, -0.5] NPARmax = [np.log10(M84)+0.1, 0.5, Av84+0.1, 0.5] @@ -3361,14 +3394,7 @@ def plot_corner_physparam_cumulative_frame(ID, Zall=np.arange(-1.2,0.4249,0.05), bins1 = np.arange(x1min, x1max + binwidth1, binwidth1) ax.hist(NPAR[i], bins=bins1, orientation='vertical', color='b', histtype='stepfilled', alpha=0.6) ax.set_xlim(x1min, x1max) - #print(x, x1min, x1max) - #ax2.scatter(np.log10(Ttmp), np.log10(Avtmp), c='r', s=1, marker='.', alpha=0.1) - #ax3.scatter(np.log10(Ztmp), np.log10(Avtmp), c='r', s=1, marker='.', alpha=0.1) - #ax.set_xlabel('$\log T_*$/Gyr', fontsize=12) - #ax.set_ylabel('$\log Z_*/Z_\odot$', fontsize=12) ax.set_yticklabels([]) - #ax.set_xticklabels([]) - #ax.set_title('%s'%(Par[i]), fontsize=12) if i == K-1: ax.set_xlabel('%s'%(Par[i]), fontsize=12) if i < K-1: @@ -3377,14 +3403,11 @@ def plot_corner_physparam_cumulative_frame(ID, Zall=np.arange(-1.2,0.4249,0.05), # Scatter and contour for i, x in enumerate(Par): for j, y in enumerate(Par): - #print(i,j,Par[j], Par[i]) if i > j: ax = axes[i, j] ax.scatter(NPAR[j], NPAR[i], c='b', s=1, marker='o', alpha=0.01) ax.set_xlabel('%s'%(Par[j]), fontsize=12) - #x1min, x1max = np.min(NPAR[j]), np.max(NPAR[j]) - #y1min, y1max = np.min(NPAR[i]), np.max(NPAR[i]) x1min, x1max = NPARmin[j], NPARmax[j] y1min, y1max = NPARmin[i], NPARmax[i] ax.set_xlim(x1min, x1max) diff --git a/gsf/writing.py b/gsf/writing.py index 6659f3a..2eac007 100644 --- a/gsf/writing.py +++ b/gsf/writing.py @@ -1,6 +1,7 @@ import numpy as np from astropy.io import fits import asdf +import gsf from .function import filconv, calc_Dn4 @@ -20,15 +21,18 @@ def get_param(self, res, fitc, tcalc=1., burnin=-1): zrecom = self.zgal Czrec0 = self.Cz0 Czrec1 = self.Cz1 + Czrec2 = self.Cz2 try: z_cz = self.z_cz scl_cz0 = self.scl_cz0 scl_cz1 = self.scl_cz1 + scl_cz2 = self.scl_cz2 except: # When redshiftfit is skipped. z_cz = np.asarray([self.zgal,self.zgal,self.zgal]) scl_cz0 = np.asarray([self.Cz0,self.Cz0,self.Cz0]) scl_cz1 = np.asarray([self.Cz1,self.Cz1,self.Cz1]) + scl_cz2 = np.asarray([self.Cz2,self.Cz2,self.Cz2]) tau0 = self.tau0 ID0 = self.ID @@ -163,6 +167,7 @@ def get_param(self, res, fitc, tcalc=1., burnin=-1): prihdr['ID'] = ID0 prihdr['Cz0'] = Czrec0 prihdr['Cz1'] = Czrec1 + prihdr['Cz2'] = Czrec2 prihdr['z'] = zrecom prihdr['zmc'] = zmc[1] prihdr['SN'] = SN @@ -174,7 +179,6 @@ def get_param(self, res, fitc, tcalc=1., burnin=-1): prihdr['bic'] = res.bic prihdr['nmc'] = nmc prihdr['nwalker'] = nwalker - import gsf prihdr['version'] = gsf.__version__ prihdu = fits.PrimaryHDU(header=prihdr) @@ -265,6 +269,10 @@ def get_param(self, res, fitc, tcalc=1., burnin=-1): col50 = fits.Column(name='Cscale1', format='E', unit='', array=scl_cz1[:]) col01.append(col50) + # C2 scale + col50 = fits.Column(name='Cscale2', format='E', unit='', array=scl_cz2[:]) + col01.append(col50) + col50 = fits.Column(name='logf', format='E', unit='', array=logf) col01.append(col50) diff --git a/gsf/zfit.py b/gsf/zfit.py index c2e60c1..7086d32 100644 --- a/gsf/zfit.py +++ b/gsf/zfit.py @@ -1,6 +1,6 @@ def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, zprior, prior, NR, zliml, zlimu, \ - nmc_cz=100, nwalk_cz=10, nthin=5, f_line_check=False, f_vary=True): + nmc_cz=100, nwalk_cz=10, nthin=5, f_line_check=False, f_vary=True, NRbb_lim=10000): ''' Purpose ------- @@ -49,6 +49,7 @@ def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, zprior, prior, NR, z fit_par_cz.add('z', value=zbest, min=zliml, max=zlimu, vary=f_vary) fit_par_cz.add('Cz0', value=1, min=0.5, max=1.5) fit_par_cz.add('Cz1', value=1, min=0.5, max=1.5) + fit_par_cz.add('Cz2', value=1, min=0.5, max=1.5) ############################## def residual_z(pars): @@ -56,6 +57,7 @@ def residual_z(pars): z = vals['z'] Cz0s = vals['Cz0'] Cz1s = vals['Cz1'] + Cz2s = vals['Cz2'] xm_s = xm_tmp * (1+z) fint = interpolate.interp1d(xm_s, fm_tmp, kind='nearest', fill_value="extrapolate") @@ -64,17 +66,23 @@ def residual_z(pars): con0 = (NR<1000) fy0 = fobs[con0] * Cz0s ey0 = eobs[con0] * Cz0s - con1 = (NR>=1000) & (NR<10000) + con1 = (NR>=1000) & (NR<2000) fy1 = fobs[con1] * Cz1s ey1 = eobs[con1] * Cz1s - con2 = (NR>=10000) # BB - fy2 = fobs[con2] - ey2 = eobs[con2] + con2 = (NR>=2000) & (NR=NRbb_lim) # BB + fy_bb = fobs[con_bb] + ey_bb = eobs[con_bb] fy01 = np.append(fy0,fy1) - fcon = np.append(fy01,fy2) ey01 = np.append(ey0,ey1) - eycon = np.append(ey01,ey2) + fy02 = np.append(fy01,fy2) + ey02 = np.append(ey01,ey2) + + fcon = np.append(fy02,fy_bb) + eycon = np.append(ey02,ey_bb) wht = 1./np.square(eycon) if f_line_check: @@ -110,12 +118,9 @@ def lnprob_cz(pars): respr = np.log(prior[nzz]) resid += np.log(2 * np.pi * s_z**2) return -0.5 * np.sum(resid) + respr - ################################# - # - # Best fit - # + # Get Best fit out_cz = minimize(residual_z, fit_par_cz, method='nelder') keys = fit_report(out_cz).split('\n') for key in keys: @@ -127,7 +132,6 @@ def lnprob_cz(pars): rcsq = float(skey[7]) fitc_cz = [csq, rcsq] # Chi2, Reduced-chi2 - mini_cz = Minimizer(lnprob_cz, out_cz.params) res_cz = mini_cz.emcee(burn=int(nmc_cz/2), steps=nmc_cz, thin=nthin, nwalkers=nwalk_cz, params=out_cz.params, is_weighted=True) From ed218b77e9b5d573579ed2d8d5f22142dbee901b Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Wed, 22 Jun 2022 09:51:27 -0700 Subject: [PATCH 04/14] minor to grism --- gsf/fitting.py | 14 ++++++++++++-- gsf/maketmp_filt.py | 18 +++++++++--------- gsf/plot_sed.py | 45 +++++++++++++++++++++++---------------------- 3 files changed, 44 insertions(+), 33 deletions(-) diff --git a/gsf/fitting.py b/gsf/fitting.py index f10d447..cfec5dd 100644 --- a/gsf/fitting.py +++ b/gsf/fitting.py @@ -1124,7 +1124,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, return flag_z - def get_zdist(self, f_interact=False): + def get_zdist(self, f_interact=False, f_ascii=True): ''' Saves a plot of z-distribution. @@ -1143,9 +1143,17 @@ def get_zdist(self, f_interact=False): yy = np.arange(0,np.max(n),1) xx = yy * 0 + self.z_cz[1] ax1.plot(xx,yy,linestyle='-',linewidth=1,color='orangered',\ - label='$z=%.5f_{-%.5f}^{+%.5f}$\n$C_z0=%.3f$\n$C_z1=%.3f$\n$C_z2=%.3f$'%\ + label='$z=%.5f_{-%.5f}^{+%.5f}$\n$C_{z0}=%.3f$\n$C_{z1}=%.3f$\n$C_{z2}=%.3f$'%\ (self.z_cz[1],self.z_cz[1]-self.z_cz[0],self.z_cz[2]-self.z_cz[1], self.Cz0, self.Cz1, self.Cz2)) + if f_ascii: + file_ascii_out = self.DIR_OUT + 'zprob_' + self.ID + '.txt' + fw_ascii = open(file_ascii_out,'w') + fw_ascii.write('# z pz\n') + for ii in range(len(xx)): + fw_ascii.write('%.3f %.3f\n'%(xx[ii],yy[ii])) + fw_ascii.close() + xx = yy * 0 + self.z_cz[0] ax1.plot(xx,yy,linestyle='--',linewidth=1,color='orangered') xx = yy * 0 + self.z_cz[2] @@ -1159,6 +1167,8 @@ def get_zdist(self, f_interact=False): # Label: ax1.set_xlabel('Redshift') ax1.set_ylabel('$dn/dz$') + zp_min,zp_max = self.z_cz[0] - (self.z_cz[1]-self.z_cz[0])*3, self.z_cz[2] + (self.z_cz[2]-self.z_cz[1])*3 + ax1.set_xlim(zp_min,zp_max) ax1.legend(loc=0) # Save: diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py index 05fa46a..b5d6931 100644 --- a/gsf/maketmp_filt.py +++ b/gsf/maketmp_filt.py @@ -211,21 +211,22 @@ def get_LSF(inputs, DIR_EXTR, ID, lm, c=3e18): ############################ try: sig_temp = float(inputs['SIG_TEMP']) + print('Template is set to %.1f km/s.'%(sig_temp)) except: sig_temp = 50. print('Template resolution is unknown.') print('Set to %.1f km/s.'%(sig_temp)) - dellam = lm[1] - lm[0] # AA/pix - R_temp = c/(sig_temp*1e3*1e10) - sig_temp_pix = np.median(lm) / R_temp / dellam # delta v in pixel; - # + iixlam = np.argmin(np.abs(lm-4000)) # Around 4000 AA + dellam = lm[iixlam+1] - lm[iixlam] # AA/pix + R_temp = c / (sig_temp*1e3*1e10) + sig_temp_pix = np.median(lm) / R_temp / dellam # delta v in pixel; sig_inst = 0 #65 #km/s for Manga # If grism; if f_morp: - print('Templates convolution (intrinsic morphology).') - if gamma>sig_temp_pix: + print('\nStarting templates convolution (intrinsic morphology).') + if gamma>sig_temp_pix:# and False: sig_conv = np.sqrt(gamma**2-sig_temp_pix**2) else: sig_conv = 0 @@ -250,7 +251,6 @@ def get_LSF(inputs, DIR_EXTR, ID, lm, c=3e18): try: vdisp = float(inputs['VDISP']) dellam = lm[1] - lm[0] # AA/pix - #R_disp = c/(vdisp*1e3*1e10) R_disp = c/(np.sqrt(vdisp**2-sig_inst**2)*1e3*1e10) vdisp_pix = np.median(lm) / R_disp / dellam # delta v in pixel; print('Templates are convolved at %.2f km/s.'%(vdisp)) @@ -289,7 +289,7 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, flag for adding nebular emissionself. tmp_norm : float Normalization of the stored templated. i.e. each template is in units of tmp_norm [Lsun]. - ''' + ''' inputs = MB.inputs ID = MB.ID @@ -592,7 +592,6 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, # Flam to Fnu spec_mul_nu[ss,:] = flamtonu(wave, spec_mul[ss,:], m0set=MB.m0set) - spec_mul_nu[ss,:] *= Lsun/(4.*np.pi*DL**2/(1.+zbest)) # (1.+zbest) takes acount of the change in delta lam by redshifting. # Note that this is valid only when F_nu. @@ -618,6 +617,7 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, try: spec_mul_nu_conv[ss,:] = convolve(spec_mul_nu[ss], LSF, boundary='extend') except: + print('Error. No convolution is happening...') spec_mul_nu_conv[ss,:] = spec_mul_nu[ss] if zz==0 and ss==0: print('Kernel is too small. No convolution.') diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py index 109d42c..5e7e9cd 100644 --- a/gsf/plot_sed.py +++ b/gsf/plot_sed.py @@ -737,7 +737,6 @@ def gaus(x,a,x0,sigma): # Update Progress Bar printProgressBar(kk, mmax, prefix = 'Progress:', suffix = 'Complete', length = 40) - # # Plot Median SED; # @@ -756,12 +755,15 @@ def gaus(x,a,x0,sigma): if f_grsm: from astropy.convolution import convolve from .maketmp_filt import get_LSF - LSF, lmtmp = get_LSF(MB.inputs, MB.DIR_EXTR, ID, x1_tot[::nstep_plot], c=3e18) - spec_grsm16 = convolve(ytmp16[::nstep_plot], LSF, boundary='extend') - spec_grsm50 = convolve(ytmp50[::nstep_plot], LSF, boundary='extend') - spec_grsm84 = convolve(ytmp84[::nstep_plot], LSF, boundary='extend') - ax2t.plot(x1_tot[::nstep_plot], spec_grsm50, '-', lw=0.5, color='gray', zorder=3., alpha=1.0) - + LSF, _ = get_LSF(MB.inputs, MB.DIR_EXTR, ID, x1_tot[:], c=3e18) + spec_grsm16 = convolve(ytmp16[:], LSF, boundary='extend') + spec_grsm50 = convolve(ytmp50[:], LSF, boundary='extend') + spec_grsm84 = convolve(ytmp84[:], LSF, boundary='extend') + if True: + ax2t.plot(x1_tot[:], ytmp50, '-', lw=0.5, color='gray', zorder=3., alpha=1.0) + else: + ax2t.plot(x1_tot[:], spec_grsm50, '-', lw=0.5, color='gray', zorder=3., alpha=1.0) + # Attach the data point in MB; MB.sed_wave_obs = xbb MB.sed_flux_obs = fybb * c / np.square(xbb) / d @@ -779,9 +781,7 @@ def gaus(x,a,x0,sigma): ysumtmp2 = ytmp[:, ::nstep_plot] * 0 ysumtmp2_prior = ytmp[0, ::nstep_plot] * 0 for ss in range(len(age)): - ii = int(len(nage) - ss - 1) # from old to young templates. - #ysumtmp += np.percentile(ytmp_each[:, ::nstep_plot, ii], 50, axis=0) - #ax1.plot(x1_tot[::nstep_plot], ysumtmp, linestyle='--', lw=.5, color=col[ii], alpha=0.5) + ii = int(len(nage) - ss - 1) # !! Take median after summation; ysumtmp2[:,:len(xm_tmp)] += ytmp_each[:, ::nstep_plot, ii] if f_fill: @@ -989,11 +989,11 @@ def func_tmp(xint,eobs,fmodel): fybb = np.append(fybb,fybbd) eybb = np.append(eybb,eybbd) - col5 = fits.Column(name='wave_obs', format='E', unit='AA', array=xbb) + col5 = fits.Column(name='wave_obs', format='E', unit='AA', array=xbb) col00.append(col5) - col6 = fits.Column(name='f_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=fybb[:] * c / np.square(xbb[:]) / d) + col6 = fits.Column(name='f_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=fybb[:] * c / np.square(xbb[:]) / d) col00.append(col6) - col7 = fits.Column(name='e_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=eybb[:] * c / np.square(xbb[:]) / d) + col7 = fits.Column(name='e_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=eybb[:] * c / np.square(xbb[:]) / d) col00.append(col7) hdr = fits.Header() @@ -1161,13 +1161,17 @@ def func_tmp(xint,eobs,fmodel): %(ID, zbes, float(fd['Mstel_50']), float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50']), fin_chi2) ylabel = ymax*0.25 - #print(ax1.get_ylim()) - ax1.text(0.77, 0.7, label,\ - fontsize=9, bbox=dict(facecolor='w', alpha=0.7), zorder=10, - ha='left', va='center', transform=ax1.transAxes) - - ####################################### + if f_grsm: + ax1.text(0.02, 0.7, label,\ + fontsize=9, bbox=dict(facecolor='w', alpha=0.7), zorder=10, + ha='left', va='center', transform=ax1.transAxes) + else: + ax1.text(0.77, 0.7, label,\ + fontsize=9, bbox=dict(facecolor='w', alpha=0.7), zorder=10, + ha='left', va='center', transform=ax1.transAxes) + ax1.xaxis.labelpad = -3 + if f_grsm: if np.max(xg0)<23000: # E.g. WFC3, NIRISS grisms conlim = (x0>10000) & (x0<25000) @@ -1219,9 +1223,6 @@ def func_tmp(xint,eobs,fmodel): ax3t.set_xscale('log') ax3t.set_xticks([100000, 1000000, 10000000]) ax3t.set_xticklabels(['10', '100', '1000']) - #ax3t.set_xlim(7e6, 3e7) - #ax3t.set_ylim(1e-4, 0.1) - #ax3t.set_yscale('log') ############### # Line name From 47bbfc61eb43135cb0848e2ace4324a9724e819c Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Wed, 22 Jun 2022 11:40:49 -0700 Subject: [PATCH 05/14] bug in tau model fixed --- gsf/fitting.py | 1 + gsf/function_class.py | 1 - gsf/maketmp_filt.py | 126 +++++------------------------------------- gsf/plot_sed.py | 32 +++++------ gsf/writing.py | 17 ++++-- 5 files changed, 44 insertions(+), 133 deletions(-) diff --git a/gsf/fitting.py b/gsf/fitting.py index cfec5dd..6025624 100644 --- a/gsf/fitting.py +++ b/gsf/fitting.py @@ -1390,6 +1390,7 @@ def prepare_class(self, add_fir=None): # Load Spectral library; self.lib = self.fnc.open_spec_fits(fall=0) self.lib_all = self.fnc.open_spec_fits(fall=1) + if self.f_dust: self.lib_dust = self.fnc.open_spec_dust_fits(fall=0) self.lib_dust_all = self.fnc.open_spec_dust_fits(fall=1) diff --git a/gsf/function_class.py b/gsf/function_class.py index baa831f..7cc7b91 100644 --- a/gsf/function_class.py +++ b/gsf/function_class.py @@ -776,7 +776,6 @@ def open_spec_fits(self, fall:int=0, orig:bool=False): hdu0 = self.MB.af['spec_full'] DIR_TMP = self.DIR_TMP - NZ = len(ZZ) NT = self.MB.ntau NA = self.MB.nage diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py index b5d6931..6d8e8a8 100644 --- a/gsf/maketmp_filt.py +++ b/gsf/maketmp_filt.py @@ -644,11 +644,10 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tree_spec_full.update({'wavelength':wavetmp}) tree_spec_full.update({'colnum':nd}) - spec_ap = np.append(ftmp_nu_int[ss,:], ftmpbb[ss,:]) - # ASDF - tree_spec.update({'fspec_'+str(zz)+'_'+str(ss)+'_'+str(pp): spec_ap}) - # ASDF + # ??? + #spec_ap = np.append(ftmp_nu_int[ss,:], ftmpbb[ss,:]) + #tree_spec.update({'fspec_'+str(zz)+'_'+str(ss)+'_'+str(pp): spec_ap}) tree_spec_full.update({'fspec_orig_'+str(zz)+'_'+str(ss)+'_'+str(pp): spec_mul_nu[ss,:]}) tree_spec_full.update({'fspec_'+str(zz)+'_'+str(ss)+'_'+str(pp): spec_mul_nu_conv[ss,:]}) @@ -1101,7 +1100,7 @@ def maketemp_tau(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tau_ ebb = np.zeros(len(SFILT), dtype='float') for ii in range(len(SFILT)): fbb[ii] = 0 - ebb[ii] = -99 #1000 + ebb[ii] = -99 # Dust flux; if MB.f_dust: @@ -1126,93 +1125,10 @@ def maketemp_tau(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tau_ Amp = 0 f_morp = False if f_spec: - try: - if inputs['MORP'] == 'moffat' or inputs['MORP'] == 'gauss': - f_morp = True - try: - mor_file = inputs['MORP_FILE'].replace('$ID','%s'%(ID)) - #fm = np.loadtxt(DIR_EXTR + mor_file, comments='#') - fm = ascii.read(DIR_EXTR + mor_file) - Amp = fm['A'] - gamma = fm['gamma'] - if inputs['MORP'] == 'moffat': - alp = fm['alp'] - else: - alp = 0 - except Exception: - print('Error in reading morphology params.') - print('No morphology convolution.') - pass - else: - print('MORP Keywords does not match.') - print('No morphology convolution.') - except: - pass - - ############################ - # Template convolution; - ############################ - try: - sig_temp = float(inputs['SIG_TEMP']) - except: - sig_temp = 50. - print('Template resolution is unknown.') - print('Set to %.1f km/s.'%(sig_temp)) - dellam = lm[1] - lm[0] # AA/pix - R_temp = c/(sig_temp*1e3*1e10) - sig_temp_pix = np.median(lm) / R_temp / dellam # delta v in pixel; - - # - sig_inst = 0 #65 #km/s for Manga - - # If grism; - if f_morp: - print('Templates convolution (intrinsic morphology).') - if gamma>sig_temp_pix: - sig_conv = np.sqrt(gamma**2-sig_temp_pix**2) - else: - sig_conv = 0 - print('Template resolution is broader than Morphology.') - print('No convolution is applied to templates.') - - xMof = np.arange(-5, 5.1, .1) # dimension must be even. - if inputs['MORP'] == 'moffat' and Amp>0 and alp>0: - LSF = moffat(xMof, Amp, 0, np.sqrt(gamma**2-sig_temp_pix**2), alp) - print('Template convolution with Moffat.') - elif inputs['MORP'] == 'gauss': - sigma = gamma - LSF = gauss(xMof, Amp, np.sqrt(sigma**2-sig_temp_pix**2)) - print('Template convolution with Gaussian.') - print('params is sigma;',sigma) - else: - print('Something is wrong with the convolution file. Exiting.') - return False - - else: # For slit spectroscopy. To be updated... - print('Templates convolution (intrinsic velocity).') - try: - vdisp = float(inputs['VDISP']) - dellam = lm[1] - lm[0] # AA/pix - #R_disp = c/(vdisp*1e3*1e10) - R_disp = c/(np.sqrt(vdisp**2-sig_inst**2)*1e3*1e10) - vdisp_pix = np.median(lm) / R_disp / dellam # delta v in pixel; - print('Templates are convolved at %.2f km/s.'%(vdisp)) - if vdisp_pix-sig_temp_pix>0: - sig_conv = np.sqrt(vdisp_pix**2-sig_temp_pix**2) - else: - sig_conv = 0 - except: - vdisp = 0. - print('Templates are not convolved.') - sig_conv = 0 #np.sqrt(sig_temp_pix**2) - pass - xMof = np.arange(-5, 5.1, .1) # dimension must be even. - Amp = 1. - LSF = gauss(xMof, Amp, sig_conv) + LSF, lm = get_LSF(inputs, DIR_EXTR, ID, lm) else: lm = [] - #################################### # Start generating templates #################################### @@ -1276,30 +1192,25 @@ def maketemp_tau(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tau_ spec_av_tmp = spec_mul[ss,:] spec_mul_nu[ss,:] = flamtonu(wave, spec_av_tmp, m0set=MB.m0set) - # (1.+zbest) takes acount of the change in delta lam by redshifting. - # Note that this is valid only when F_nu. - # When Flambda, /(1.+zbest) will be *(1.+zbest). + + spec_mul_nu[ss,:] *= Lsun/(4.*np.pi*DL**2/(1.+zbest)) + spec_mul_nu[ss,:] *= (1./Ls[ss])*tmp_norm # in unit of erg/s/Hz/cm2/ms[ss]. + ms[ss] *= (1./Ls[ss])*tmp_norm # M/L; 1 unit template has this mass in [Msolar]. if len(lm)>0: try: - spec_mul_nu_conv[ss,:] = convolve(spec_mul_nu[ss], LSF, boundary='extend') + spec_mul_nu_conv[ss,:] = convolve(spec_mul_nu[ss,:], LSF, boundary='extend') except: - spec_mul_nu_conv[ss,:] = spec_mul_nu[ss] + spec_mul_nu_conv[ss,:] = spec_mul_nu[ss,:] if zz==0 and ss==0: print('Kernel is too small. No convolution.') else: - spec_mul_nu_conv[ss,:] = spec_mul_nu[ss] - - spec_mul_nu_conv[ss,:] *= Lsun/(4.*np.pi*DL**2/(1.+zbest)) - spec_mul_nu_conv[ss,:] *= (1./Ls[ss])*tmp_norm # in unit of erg/s/Hz/cm2/ms[ss]. - ms[ss] *= (1./Ls[ss])*tmp_norm # M/L; 1 unit template has this mass in [Msolar]. + spec_mul_nu_conv[ss,:] = spec_mul_nu[ss,:] if f_spec: ftmp_nu_int[ss,:] = data_int(lm, wavetmp, spec_mul_nu_conv[ss,:]) # Register filter response; - #if ss == 0 and tt == 0 and zz == 0: - # filconv(SFILT, wavetmp, spec_mul_nu_conv[ss,:], DIR_FILT, fw=True, MB=MB, f_regist=True) ltmpbb[ss,:], ftmpbb[ss,:] = filconv(SFILT, wavetmp, spec_mul_nu_conv[ss,:], DIR_FILT, MB=MB, f_regist=False) ########################################## @@ -1323,11 +1234,12 @@ def maketemp_tau(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tau_ tree_spec_full.update({'colnum':nd}) # ASDF + # ??? spec_ap = np.append(ftmp_nu_int[ss,:], ftmpbb[ss,:]) tree_spec.update({'fspec_'+str(zz)+'_'+str(tt)+'_'+str(ss): spec_ap}) + tree_spec_full.update({'fspec_orig_'+str(zz)+'_'+str(tt)+'_'+str(ss): spec_mul_nu[ss,:]}) tree_spec_full.update({'fspec_'+str(zz)+'_'+str(tt)+'_'+str(ss): spec_mul_nu_conv[ss,:]}) - # Nebular library; if fneb == 1 and MB.f_bpass==0 and ss==0 and tt==0: if zz==0: @@ -1374,15 +1286,12 @@ def maketemp_tau(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tau_ spec_neb_ap = np.append(ftmp_neb_nu_int[zz,uu,:], ftmpbb_neb[zz,uu,:]) tree_spec.update({'fspec_nebular_Z%d_logU%d'%(zz,uu): spec_neb_ap}) - ######################### # Summarize the ML ######################### # ASDF tree_ML.update({'ML_'+str(zz)+'_'+str(tt): ms}) - - ######################### # Summarize the templates ######################### @@ -1477,13 +1386,6 @@ def maketemp_tau(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, tau_ fw.write('%d %.5f %.5e %.5e\n'%(ii+g_offset, lm[ii], fobs[ii], eobs[ii])) else: fw.write('%d %.5f 0 1000\n'%(ii+g_offset, lm[ii])) - ''' - elif fgrs[ii]==1: # grism 2 - if lm[ii]/(1.+zbest) > lamliml and lm[ii]/(1.+zbest) < lamlimu: - fw.write('%d %.5f %.5e %.5e\n'%(ii+1000, lm[ii], fobs[ii], eobs[ii])) - else: - fw.write('%d %.5f 0 1000\n'%(ii+1000, lm[ii])) - ''' for ii in range(len(ltmpbb[0,:])): if SFILT[ii] in SKIPFILT:# data point to be skiped; diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py index 5e7e9cd..737be6f 100644 --- a/gsf/plot_sed.py +++ b/gsf/plot_sed.py @@ -1430,12 +1430,10 @@ def gaus(x,a,x0,sigma): Z50 = np.zeros(len(age), dtype='float') Z16 = np.zeros(len(age), dtype='float') Z84 = np.zeros(len(age), dtype='float') - #NZbest = np.zeros(len(age), dtype='int') for aa in range(len(age)): Z16[aa] = hdul[1].data['Z'+str(aa)][0] Z50[aa] = hdul[1].data['Z'+str(aa)][1] Z84[aa] = hdul[1].data['Z'+str(aa)][2] - #NZbest[aa]= bfnc.Z2NZ(Z50[aa]) vals['Z'+str(aa)] = Z50[aa] # Light weighted Z. @@ -1470,6 +1468,7 @@ def gaus(x,a,x0,sigma): Cz0 = hdul[0].header['Cz0'] Cz1 = hdul[0].header['Cz1'] + Cz2 = hdul[0].header['Cz2'] zbes = zp50 zscl = (1.+zbes) @@ -1477,9 +1476,9 @@ def gaus(x,a,x0,sigma): # Data taken from ############################### if MB.f_dust: - MB.dict = MB.read_data(Cz0, Cz1, zbes, add_fir=True) + MB.dict = MB.read_data(Cz0, Cz1, Cz2, zbes, add_fir=True) else: - MB.dict = MB.read_data(Cz0, Cz1, zbes) + MB.dict = MB.read_data(Cz0, Cz1, Cz2, zbes) NR = MB.dict['NR'] x = MB.dict['x'] @@ -1488,16 +1487,16 @@ def gaus(x,a,x0,sigma): con0 = (NR<1000) xg0 = x[con0] - fg0 = fy[con0] #* Cz0 - eg0 = ey[con0] #* Cz0 + fg0 = fy[con0] + eg0 = ey[con0] con1 = (NR>=1000) & (NR<2000) #& (fy/ey>SNlim) xg1 = x[con1] - fg1 = fy[con1] #* Cz1 - eg1 = ey[con1] #* Cz1 + fg1 = fy[con1] + eg1 = ey[con1] con2 = (NR>=2000) & (NRSNlim) xg2 = x[con2] - fg2 = fy[con2] #* Cz1 - eg2 = ey[con2] #* Cz1 + fg2 = fy[con2] + eg2 = ey[con2] if len(xg0)>0 or len(xg1)>0 or len(xg2)>0: f_grsm = True else: @@ -1509,11 +1508,11 @@ def gaus(x,a,x0,sigma): wht[con_wht] = 1./np.square(ey[con_wht]) # BB data points; - NRbb = MB.dict['NRbb'] #dat[:, 0] - xbb = MB.dict['xbb'] #dat[:, 1] - fybb = MB.dict['fybb'] #dat[:, 2] - eybb = MB.dict['eybb'] #dat[:, 3] - exbb = MB.dict['exbb'] #dat[:, 4] + NRbb = MB.dict['NRbb'] + xbb = MB.dict['xbb'] + fybb = MB.dict['fybb'] + eybb = MB.dict['eybb'] + exbb = MB.dict['exbb'] snbb = fybb/eybb ###################### @@ -1552,6 +1551,7 @@ def gaus(x,a,x0,sigma): ax2t = ax1.inset_axes((1-xsize-0.01,1-ysize-0.01,xsize,ysize)) if f_dust: ax3t = ax1.inset_axes((0.7,.35,.28,.25)) + f_plot_resid = False else: if f_plot_resid: fig_mosaic = """ @@ -1841,7 +1841,7 @@ def gaus(x,a,x0,sigma): xgrism = np.concatenate([xg0,xg1,xg2]) fgrism = np.concatenate([fg0,fg1,fg2]) - egrism = np.concatenate([eg0,eg1,eg1]) + egrism = np.concatenate([eg0,eg1,eg2]) con4000b = (xgrism/zscl>3400) & (xgrism/zscl<3800) & (fgrism>0) & (egrism>0) con4000r = (xgrism/zscl>4200) & (xgrism/zscl<5000) & (fgrism>0) & (egrism>0) print('Median SN at 3400-3800 is;', np.median((fgrism/egrism)[con4000b])) diff --git a/gsf/writing.py b/gsf/writing.py index 2eac007..5c78752 100644 --- a/gsf/writing.py +++ b/gsf/writing.py @@ -115,10 +115,19 @@ def get_param(self, res, fitc, tcalc=1., burnin=-1): except: pass - AGEmc[aa,:] = np.percentile(res.flatchain['AGE'+str(aa)][burnin:], [16,50,84]) - TAUmc[aa,:] = np.percentile(res.flatchain['TAU'+str(aa)][burnin:], [16,50,84]) + # Do debug.. + try: + AGEmc[aa,:] = np.percentile(res.flatchain['AGE'+str(aa)][burnin:], [16,50,84]) + except: + AGEFIX = self.AGEFIX + AGEmc[aa,:] = [AGEFIX, AGEFIX, AGEFIX] + + try: + TAUmc[aa,:] = np.percentile(res.flatchain['TAU'+str(aa)][burnin:], [16,50,84]) + except: + TAUFIX = self.TAUFIX + TAUmc[aa,:] = [TAUFIX, TAUFIX, TAUFIX] - # msmc = np.percentile(msmc0, [16,50,84]) try: Avb = res.params['Av'].value @@ -153,7 +162,7 @@ def get_param(self, res, fitc, tcalc=1., burnin=-1): esp = fds[:,3] consp = (nrs<10000) & (lams/(1.+zrecom)>3600) & (lams/(1.+zrecom)<4200) & (esp>0) - NSN = len(fsp[consp]) + NSN = len(fsp[consp]) if NSN>0: SN = np.median((fsp/esp)[consp]) else: From 1fa2d3283628090becddf10f7181a89941a11501 Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Thu, 23 Jun 2022 17:15:07 -0700 Subject: [PATCH 06/14] bug in plot_tau fixed --- gsf/fitting.py | 22 +++++++++++++++++----- gsf/maketmp_filt.py | 14 ++++++++------ gsf/plot_sed.py | 27 +++++++++++++++++++++------ 3 files changed, 46 insertions(+), 17 deletions(-) diff --git a/gsf/fitting.py b/gsf/fitting.py index 6025624..2ddb6ed 100644 --- a/gsf/fitting.py +++ b/gsf/fitting.py @@ -177,15 +177,18 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: # # If Mdyn is included. # self.af = asdf.open(self.DIR_TMP + 'spec_all_' + self.ID + '_PA' + self.PA + '.asdf') - # Scaling for grism; - self.Cz0 = float(inputs['CZ0']) - self.Cz1 = float(inputs['CZ1']) - self.Cz2 = float(inputs['CZ2']) try: self.DIR_EXTR = inputs['DIR_EXTR'] + # Scaling for grism; + self.Cz0 = float(inputs['CZ0']) + self.Cz1 = float(inputs['CZ1']) + self.Cz2 = float(inputs['CZ2']) except: self.DIR_EXTR = False + self.Cz0 = 1 + self.Cz1 = 1 + self.Cz2 = 1 # BPASS Binary template try: @@ -318,12 +321,21 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: self.agemin = float(inputs['AGEMIN']) self.delage = float(inputs['DELAGE']) + if self.agemax-self.agemin0) & (fy/ey<=SNlim) & (f_ex == 0) From c34cbf0e5942f10f18600a73bd7710a6a09a49a7 Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Sat, 2 Jul 2022 16:52:54 -0700 Subject: [PATCH 07/14] Update function_igm.py --- gsf/function_igm.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/gsf/function_igm.py b/gsf/function_igm.py index 829215a..3195d40 100644 --- a/gsf/function_igm.py +++ b/gsf/function_igm.py @@ -77,17 +77,20 @@ def get_sig_lya(lam_o, z_s, T=1e4, c=3e18): def get_nH(z): ''' + Purpose + ------- + Get HI density by using Cen & Haiman 2000. + Returns ------- - HI density in IGM. + HI density in IGM, in cm^-3 ''' try: nH = np.zeros(len(z),dtype='float') except: nH = 0 - # From Cen & Haiman 2000 - nH = 8.5e-5 * ((1.+z)/8)**3 # in cm^-3 + nH = 8.5e-5 * ((1.+z)/8)**3 return nH From 2282299610f897bf1f127a61dc375045e4ba2e9a Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Sat, 2 Jul 2022 16:53:39 -0700 Subject: [PATCH 08/14] Update function_igm.py --- gsf/function_igm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gsf/function_igm.py b/gsf/function_igm.py index 3195d40..1dbf47f 100644 --- a/gsf/function_igm.py +++ b/gsf/function_igm.py @@ -98,7 +98,7 @@ def get_column(zin, cosmo, Mpc_cm=3.08568025e+24, z_r=6.0, delz=0.1): ''' Returns ------- - HI column density of IGM at zin, in cm^-3. + HI column density of IGM at zin, in cm^-2. ''' z = np.arange(z_r, zin, delz) try: From 7ace9c04937bd9f14928dab8aaf0a7db22365e45 Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Sat, 2 Jul 2022 21:01:45 -0700 Subject: [PATCH 09/14] uvbeta --- gsf/function.py | 25 +++++++++++++++++++++++++ gsf/plot_sed.py | 36 +++++++++++++++++++++++++++++------- gsf/plot_sfh.py | 3 +++ 3 files changed, 57 insertions(+), 7 deletions(-) diff --git a/gsf/function.py b/gsf/function.py index d6f16c2..0d9edc3 100644 --- a/gsf/function.py +++ b/gsf/function.py @@ -17,6 +17,31 @@ fLW = np.zeros(len(LW0), dtype='int') # flag. +def get_uvbeta(lm, flam, zbes, lam_blue=1650, lam_red=2300): + ''' + Purpose + ------- + get UV beta_lambda slope. + + Parameters + ---------- + lm : float array + in lambda + flam : float array + in flambda + ''' + con_uv = (lm/(1.+zbes)>lam_blue) & (lm/(1.+zbes) Date: Wed, 20 Jul 2022 12:38:17 -0700 Subject: [PATCH 10/14] minor --- gsf/function.py | 30 ++++++++++++++++++++++++++++++ gsf/gsf.py | 4 ++-- gsf/plot_sed.py | 18 +++++++++++++----- 3 files changed, 45 insertions(+), 7 deletions(-) diff --git a/gsf/function.py b/gsf/function.py index 0d9edc3..6475a95 100644 --- a/gsf/function.py +++ b/gsf/function.py @@ -5,6 +5,7 @@ import pickle as cPickle import os import scipy.interpolate as interpolate +import logging c = 3.e18 # A/s #d = 10**(73.6/2.5) # From [ergs/s/cm2/A] to [ergs/s/cm2/Hz] @@ -17,6 +18,19 @@ fLW = np.zeros(len(LW0), dtype='int') # flag. +def str2bool(v): + ''' + ''' + if isinstance(v, bool): + return v + if v.lower() in ('yes', 'true', 't', 'y', '1'): + return True + elif v.lower() in ('no', 'false', 'f', 'n', '0'): + return False + else: + raise argparse.ArgumentTypeError('Boolean value expected.') + + def get_uvbeta(lm, flam, zbes, lam_blue=1650, lam_red=2300): ''' Purpose @@ -548,6 +562,22 @@ def fnutolam(lam, fnu, m0set=25.0): return flam +def delta(x,A): + ''' + ''' + yy = np.zeros(len(x),float) + iix = np.argmin(np.abs(x)) + if len(x)%2 == 0: + logger= logging.getLogger( __name__ ) + print(logger) + print('Input array has an even number.') + yy[iix] = A/2. + yy[iix+1] = A/2. + else: + yy[iix] = A + return yy + + def gauss(x,A,sig): return A * np.exp(-0.5*x**2/sig**2) diff --git a/gsf/gsf.py b/gsf/gsf.py index e4c8c6f..92883fd 100644 --- a/gsf/gsf.py +++ b/gsf/gsf.py @@ -109,7 +109,7 @@ def run_gsf_template(inputs, fplt=0, tau_lim=0.001, idman=None, nthin=1, delwave def run_gsf_all(parfile, fplt, cornerplot=True, f_Alog=True, idman=None, zman=None, f_label=True, f_symbol=True, f_SFMS=True, f_fill=True, save_sed=True, figpdf=False, mmax=300, skip_sfh=False, f_fancyplot=False, skip_zhist=False, tau_lim=0.001, tset_SFR_SED=0.1, f_shuffle=False, amp_shuffle=1e-2, Zini=None, - nthin=1, delwave=1): + nthin=1, delwave=1, f_plot_resid=False): ''' Purpose ------- @@ -234,7 +234,7 @@ def run_gsf_all(parfile, fplt, cornerplot=True, f_Alog=True, idman=None, zman=No plot_sed(MB, fil_path=MB.DIR_FILT, figpdf=figpdf, save_sed=save_sed, mmax=mmax, dust_model=MB.dust_model, DIR_TMP=MB.DIR_TMP, f_label=f_label, f_fill=f_fill, - f_fancyplot=f_fancyplot, f_plot_resid=True) + f_fancyplot=f_fancyplot, f_plot_resid=f_plot_resid) ''' if fplt == 4: diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py index 5ee579d..bcae821 100644 --- a/gsf/plot_sed.py +++ b/gsf/plot_sed.py @@ -757,11 +757,19 @@ def gaus(x,a,x0,sigma): if f_grsm: from astropy.convolution import convolve from .maketmp_filt import get_LSF - LSF, _ = get_LSF(MB.inputs, MB.DIR_EXTR, ID, x1_tot[:], c=3e18) + LSF, _ = get_LSF(MB.inputs, MB.DIR_EXTR, ID, x1_tot[:]/(1.+zbes), c=3e18) spec_grsm16 = convolve(ytmp16[:], LSF, boundary='extend') spec_grsm50 = convolve(ytmp50[:], LSF, boundary='extend') spec_grsm84 = convolve(ytmp84[:], LSF, boundary='extend') - if True: + ''' + plt.close() + plt.plot(x1_tot, ytmp50, color='r') + plt.plot(x1_tot, spec_grsm50, color='gray') + print(spec_grsm50) + plt.xlim(8000,20000) + plt.show() + ''' + if False: ax2t.plot(x1_tot[:], ytmp50, '-', lw=0.5, color='gray', zorder=3., alpha=1.0) else: ax2t.plot(x1_tot[:], spec_grsm50, '-', lw=0.5, color='gray', zorder=3., alpha=1.0) @@ -1174,11 +1182,11 @@ def func_tmp(xint,eobs,fmodel): if f_grsm: ax1.text(0.02, 0.7, label,\ - fontsize=9, bbox=dict(facecolor='w', alpha=0.7), zorder=10, + fontsize=6, bbox=dict(facecolor='w', alpha=0.7), zorder=10, ha='left', va='center', transform=ax1.transAxes) else: - ax1.text(0.77, 0.7, label,\ - fontsize=9, bbox=dict(facecolor='w', alpha=0.7), zorder=10, + ax1.text(0.02, 0.7, label,\ + fontsize=6, bbox=dict(facecolor='w', alpha=0.7), zorder=10, ha='left', va='center', transform=ax1.transAxes) ax1.xaxis.labelpad = -3 From 53418c1fdd757afeafc6fb7f422c2b29b668de16 Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Fri, 22 Jul 2022 15:36:18 -0700 Subject: [PATCH 11/14] minor --- gsf/fitting.py | 21 +++++++++++++-------- gsf/plot_sed.py | 22 ++++++++++++++++++++-- 2 files changed, 33 insertions(+), 10 deletions(-) diff --git a/gsf/fitting.py b/gsf/fitting.py index 2ddb6ed..2c7ec7c 100644 --- a/gsf/fitting.py +++ b/gsf/fitting.py @@ -115,8 +115,12 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: print('\nFitting : %s\n'%self.ID) # Read catalog; - self.CAT_BB = inputs['CAT_BB'] - self.fd_cat = ascii.read(self.CAT_BB) + try: + self.CAT_BB = inputs['CAT_BB'] + self.fd_cat = ascii.read(self.CAT_BB) + except: + self.CAT_BB = None + self.fd_cat = None if zman != None: self.zgal = zman @@ -176,8 +180,6 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: #if self.f_Mdyn: # # If Mdyn is included. # self.af = asdf.open(self.DIR_TMP + 'spec_all_' + self.ID + '_PA' + self.PA + '.asdf') - - try: self.DIR_EXTR = inputs['DIR_EXTR'] # Scaling for grism; @@ -859,7 +861,8 @@ def residual_z(pars,z): def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, - zlimu=None, snlim=0, priors=None, f_bb_zfit=True, f_line_check=False, f_norm=True): + zlimu=None, snlim=0, priors=None, f_bb_zfit=True, f_line_check=False, + f_norm=True, f_lambda=False): ''' Find the best-fit redshift, before going into a big fit, through an interactive inspection. This module is effective only when spec data is provided. @@ -894,9 +897,11 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, self.nmc_cz = int(self.inputs['NMCZ']) # For z prior. - zliml = self.zgal - 0.5 - if zlimu == None: - zlimu = self.zgal + 0.5 + #zliml = self.zgal - 0.5 + #if zlimu == None: + # zlimu = self.zgal + 0.5 + zliml = self.zmcmin + zlimu = self.zmcmax # Observed data; sn = self.dict['fy'] / self.dict['ey'] diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py index bcae821..13287ff 100644 --- a/gsf/plot_sed.py +++ b/gsf/plot_sed.py @@ -618,6 +618,7 @@ def gaus(x,a,x0,sigma): # Saved template; ytmp = np.zeros((mmax,len(ysum)), dtype='float') ytmp_each = np.zeros((mmax,len(ysum),len(age)), dtype='float') + ytmp_nl = np.zeros((mmax,len(ysum)), dtype='float') # no line ytmpmax = np.zeros(len(ysum), dtype='float') ytmpmin = np.zeros(len(ysum), dtype='float') @@ -662,7 +663,8 @@ def gaus(x,a,x0,sigma): if ss == MB.aamin[0]: mod0_tmp, xm_tmp = fnc.tmp03(AA_tmp, Av_tmp, ss, ZZ_tmp, zmc, lib_all) - fm_tmp = mod0_tmp + fm_tmp = mod0_tmp.copy() + fm_tmp_nl = mod0_tmp.copy() if MB.fneb: Aneb_tmp = 10**samples['Aneb'][nr] if not MB.logUFIX == None: @@ -671,9 +673,13 @@ def gaus(x,a,x0,sigma): logU_tmp = samples['logU'][nr] mod0_tmp, xm_tmp = fnc.tmp03_neb(Aneb_tmp, Av_tmp, logU_tmp, ss, ZZ_tmp, zmc, lib_neb_all) fm_tmp += mod0_tmp + # Make no emission line template; + mod0_tmp_nl, xm_tmp_nl = fnc.tmp03_neb(0, Av_tmp, logU_tmp, ss, ZZ_tmp, zmc, lib_neb_all) + fm_tmp_nl += mod0_tmp_nl else: mod0_tmp, xx_tmp = fnc.tmp03(AA_tmp, Av_tmp, ss, ZZ_tmp, zmc, lib_all) fm_tmp += mod0_tmp + fm_tmp_nl += mod0_tmp # Each; ytmp_each[kk,:,ss] = mod0_tmp[:] * c / np.square(xm_tmp[:]) / d @@ -709,12 +715,15 @@ def gaus(x,a,x0,sigma): ytmp_dust[kk,:] = model_dust * c/np.square(x1_dust)/d model_tot = np.interp(x1_tot,xx_tmp,fm_tmp) + np.interp(x1_tot,x1_dust,model_dust) + model_tot_nl = np.interp(x1_tot,xx_tmp,fm_tmp_nl) + np.interp(x1_tot,x1_dust,model_dust) ytmp[kk,:] = model_tot[:] * c/np.square(x1_tot[:])/d + ytmp_nl[kk,:] = model_tot_nl[:] * c/np.square(x1_tot[:])/d else: x1_tot = xm_tmp ytmp[kk,:] = fm_tmp[:] * c / np.square(xm_tmp[:]) / d + ytmp_nl[kk,:] = fm_tmp_nl[:] * c / np.square(xm_tmp[:]) / d # # Grism plot + Fuv flux + LIR. @@ -745,6 +754,9 @@ def gaus(x,a,x0,sigma): ytmp16 = np.percentile(ytmp[:,:],16,axis=0) ytmp50 = np.percentile(ytmp[:,:],50,axis=0) ytmp84 = np.percentile(ytmp[:,:],84,axis=0) + ytmp16_nl = np.percentile(ytmp_nl[:,:],16,axis=0) + ytmp50_nl = np.percentile(ytmp_nl[:,:],50,axis=0) + ytmp84_nl = np.percentile(ytmp_nl[:,:],84,axis=0) if MB.f_dust: ytmp_dust50 = np.percentile(ytmp_dust[:,:],50, axis=0) @@ -769,7 +781,7 @@ def gaus(x,a,x0,sigma): plt.xlim(8000,20000) plt.show() ''' - if False: + if True: ax2t.plot(x1_tot[:], ytmp50, '-', lw=0.5, color='gray', zorder=3., alpha=1.0) else: ax2t.plot(x1_tot[:], spec_grsm50, '-', lw=0.5, color='gray', zorder=3., alpha=1.0) @@ -970,6 +982,12 @@ def func_tmp(xint,eobs,fmodel): col00.append(col3) col4 = fits.Column(name='f_model_84', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=ytmp84[:]) col00.append(col4) + col2 = fits.Column(name='f_model_noline_16', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=ytmp16_nl[:]) + col00.append(col2) + col3 = fits.Column(name='f_model_noline_50', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=ytmp50_nl[:]) + col00.append(col3) + col4 = fits.Column(name='f_model_noline_84', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=ytmp84_nl[:]) + col00.append(col4) # Each component # Stellar From 2f2c9579374324a997464e56284f7dc54c5fc0db Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Wed, 3 Aug 2022 20:54:45 -0700 Subject: [PATCH 12/14] plot scale removed addressing #31 . --- gsf/fitting.py | 1 + gsf/maketmp_z0.py | 4 ++-- gsf/plot_sed.py | 23 ++++++++++++----------- 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/gsf/fitting.py b/gsf/fitting.py index 2c7ec7c..2f30776 100644 --- a/gsf/fitting.py +++ b/gsf/fitting.py @@ -141,6 +141,7 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: self.zmcmin = self.zgal - float(self.fd_cat['ez_l'][iix]) self.zmcmax = self.zgal + float(self.fd_cat['ez_u'][iix]) except: + print('ZMCMIN and ZMCMAX cannot be found. z range is set to z \pm 1.0') self.zmcmin = None self.zmcmax = None diff --git a/gsf/maketmp_z0.py b/gsf/maketmp_z0.py index a9f40cc..3dd9523 100644 --- a/gsf/maketmp_z0.py +++ b/gsf/maketmp_z0.py @@ -349,7 +349,7 @@ def make_tmp_z0_bpass(MB, lammin=100, lammax=160000, Zforce=None, \ fd_sed = ascii.read(file_sed) fd_stm = ascii.read(file_stm) - wave0 = fd_sed['col1'] + wave0 = fd_sed['col1'] age_stm = fd_stm['col1'] mass_formed = fd_stm['col2'][0] @@ -429,7 +429,7 @@ def make_tmp_z0_bpass(MB, lammin=100, lammax=160000, Zforce=None, \ tautmp = ( 10**(lage_temp[iis]+0.05) - 10**(lage_temp[iis]-0.05) ) / 1e9 # Gyr agetmp = (age[ss]+age[ss-1])/2. - flux0 = fd_sed['col%d'%(iis+2)] #sp.get_spectrum(tage=age[ss], peraa=True) + flux0 = fd_sed['col%d'%(iis+2)] ms[ss] = fd_stm['col2'][iistm] mass_formed_tot += mass_formed diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py index 13287ff..8788d48 100644 --- a/gsf/plot_sed.py +++ b/gsf/plot_sed.py @@ -519,17 +519,18 @@ def gaus(x,a,x0,sigma): ax1.set_xticks(xticks) ax1.set_xticklabels(xlabels) - dely1 = 0.5 - while (ymax-0)/dely1<1: - dely1 /= 2. - while (ymax-0)/dely1>4: - dely1 *= 2. - - y1ticks = np.arange(0, ymax, dely1) - ax1.set_yticks(y1ticks) - ax1.set_yticklabels(np.arange(0, ymax, dely1), minor=False) - ax1.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) - ax1.yaxis.labelpad = 1.5 + if False: + dely1 = 0.5 + while (ymax-0)/dely1<1: + dely1 /= 2. + while (ymax-0)/dely1>4: + dely1 *= 2. + + y1ticks = np.arange(0, ymax, dely1) + ax1.set_yticks(y1ticks) + ax1.set_yticklabels(np.arange(0, ymax, dely1), minor=False) + ax1.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) + ax1.yaxis.labelpad = 1.5 xx = np.arange(100,400000) yy = xx * 0 From e9916d44fcca85e7b8c0ba808bbe5f9ae79a911f Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Sun, 7 Aug 2022 18:11:10 -0700 Subject: [PATCH 13/14] Update plot_sfh.py --- gsf/plot_sfh.py | 34 ++++++++++++++++++++++------------ 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/gsf/plot_sfh.py b/gsf/plot_sfh.py index 9311c29..cc3369f 100644 --- a/gsf/plot_sfh.py +++ b/gsf/plot_sfh.py @@ -164,14 +164,13 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4, delTu = np.zeros(len(age),dtype='float') if len(age) == 1: - #if tau0[0] < 0: # SSP; for aa in range(len(age)): try: tau_ssp = float(MB.inputs['TAU_SSP']) except: tau_ssp = tau_lim - delTl[aa] = tau_ssp/2 - delTu[aa] = tau_ssp/2 + delTl[aa] = tau_ssp/2. + delTu[aa] = tau_ssp/2. if age[aa] < tau_lim: # This is because fsps has the minimum tau = tau_lim delT[aa] = tau_lim @@ -182,25 +181,21 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4, if aa == 0: delTl[aa] = age[aa] delTu[aa] = (age[aa+1]-age[aa])/2. - delT[aa] = delTu[aa] + delTl[aa] - #print(age[aa],age[aa]-delTl[aa],age[aa]+delTu[aa]) + delT[aa] = delTu[aa] + delTl[aa] elif Tuni < age[aa]: delTl[aa] = (age[aa]-age[aa-1])/2. delTu[aa] = Tuni-age[aa] #delTl[aa] #10. delT[aa] = delTu[aa] + delTl[aa] - #print(age[aa],age[aa]-delTl[aa],age[aa]+delTu[aa]) elif aa == len(age)-1: delTl[aa] = (age[aa]-age[aa-1])/2. delTu[aa] = Tuni - age[aa] delT[aa] = delTu[aa] + delTl[aa] - #print(age[aa],age[aa]-delTl[aa],age[aa]+delTu[aa]) else: delTl[aa] = (age[aa]-age[aa-1])/2. delTu[aa] = (age[aa+1]-age[aa])/2. if age[aa]+delTu[aa]>Tuni: delTu[aa] = Tuni-age[aa] delT[aa] = delTu[aa] + delTl[aa] - #print(age[aa],age[aa]-delTl[aa],age[aa]+delTu[aa]) con_delt = (delT<=0) delT[con_delt] = 1e10 @@ -355,6 +350,10 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4, SFR_SED[mm] = -99 for aa in range(len(age)): + + if Tuni0[aa]<0: + continue + if np.sum(10**AM[aa:,mm])>0: AC[aa, mm] = np.log10(np.sum(10**AM[aa:,mm])) ZC[aa, mm] = np.log10(np.sum(10**ZMM[aa:,mm])/10**AC[aa, mm]) @@ -366,8 +365,13 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4, ACs = 0 ALs = 0 for bb in range(aa, len(age), 1): + + if Tuni0[bb]<0: + continue + tmpAA = 10**np.random.uniform(-eA[bb],eA[bb]) - tmpTT = np.random.uniform(-delT[bb]/1e9,delT[bb]/1e9) + tmpTT = np.random.uniform(-delT[bb]/1e9/2.,delT[bb]/1e9/2.) + TC[aa, mm] += (age[bb]+tmpTT) * 10**AAtmp[bb] * mslist[bb] * tmpAA TL[aa, mm] += (age[bb]+tmpTT) * 10**AAtmp[bb] * tmpAA ACs += 10**AAtmp[bb] * mslist[bb] * tmpAA @@ -375,10 +379,16 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4, TC[aa, mm] /= ACs TL[aa, mm] /= ALs + if TC[aa, mm]>0: TC[aa, mm] = np.log10(TC[aa, mm]) + else: + TC[aa, mm] = np.nan + if TL[aa, mm]>0: TL[aa, mm] = np.log10(TL[aa, mm]) + else: + TL[aa, mm] = np.nan # Do stuff... time.sleep(0.01) @@ -1593,10 +1603,10 @@ def get_evolv(MB, ID, Z=np.arange(-1.2,0.4249,0.05), age=[0.01, 0.1, 0.3, 0.7, 1 Av_tmp = samples['Av'][mtmp] - f0 = fits.open(DIR_TMP + 'ms_' + ID + '.fits') + f0 = fits.open(DIR_TMP + 'ms_' + ID + '.fits') sedpar = f0[1] - f1 = fits.open(DIR_TMP + 'ms.fits') - mloss = f1[1].data + f1 = fits.open(DIR_TMP + 'ms.fits') + mloss = f1[1].data Avrand = np.random.uniform(-eAv, eAv) if Av_tmp + Avrand<0: From b60a40f971df7dc3e6a87a3ccf22957a622c0207 Mon Sep 17 00:00:00 2001 From: Takahiro Morishita Date: Mon, 8 Aug 2022 12:33:37 -0700 Subject: [PATCH 14/14] Update function_class.py --- gsf/function_class.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/gsf/function_class.py b/gsf/function_class.py index 7cc7b91..916953e 100644 --- a/gsf/function_class.py +++ b/gsf/function_class.py @@ -241,7 +241,7 @@ def open_spec_fits_dir(self, nage:int, nz:int, kk, Av00:float, zgal:float, A00:f return A00 * yyd_sort, xxd_sort - def get_template(self, lib, Amp:float = 1.0, T:float = 1.0, Av:float = 0.0, Z:float = 0.0, zgal:float = 1.0, f_bb:bool = False): + def get_template(self, lib, Amp:float = 1.0, T:float = 1.0, Av:float = 0.0, Z:float = 0.0, zgal:float = 1.0, f_bb:bool = False, fneb=False): ''' Gets an element template given a set of parameters. Not necessarily the most efficient way, but easy to use. @@ -282,7 +282,10 @@ def get_template(self, lib, Amp:float = 1.0, T:float = 1.0, Av:float = 0.0, Z:fl if T - self.age[nmodel] != 0: print('T=%.2f is not found in age library. T=%.2f is used.'%(T,self.age[nmodel])) - coln = int(2 + pp*len(self.ZZ)*len(self.AA) + NZ*len(self.AA) + nmodel) + if fneb: + coln = int(2 + pp*len(self.ZZ)*1 + NZ*1 + 0) + else: + coln = int(2 + pp*len(self.ZZ)*len(self.AA) + NZ*len(self.AA) + nmodel) nr = lib[:, 0] xx = lib[:, 1] # This is OBSERVED wavelength range at z=zgal yy = lib[:, coln] @@ -780,7 +783,7 @@ def open_spec_fits(self, fall:int=0, orig:bool=False): NT = self.MB.ntau NA = self.MB.nage for zz,Z in enumerate(ZZ): - for tt,TT in enumerate(self.MB.tau): + for tt,TT in enumerate(self.MB.tau): for ss,TA in enumerate(self.MB.ageparam): if zz == 0 and tt == 0 and ss == 0: nr = hdu0['colnum']