From 246a63a0087e39e8fe99df938d0b866dcfbc20a6 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Tue, 14 Jul 2020 11:36:29 -0400
Subject: [PATCH 1/9] version

---
 gsf/__init__.py | 2 +-
 setup.py        | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/gsf/__init__.py b/gsf/__init__.py
index 1f6f5a3..20c96a6 100644
--- a/gsf/__init__.py
+++ b/gsf/__init__.py
@@ -1,4 +1,4 @@
 __author__ = 'Takahiro Morishita'
 __email__ = 'tmorishita@stsci.edu'
-__version__ = '1.2'
+__version__ = '1.3'
 __credits__ = 'STScI'
diff --git a/setup.py b/setup.py
index d439e34..4a9e214 100644
--- a/setup.py
+++ b/setup.py
@@ -51,7 +51,7 @@ def read(fname):
 
 setup(
     name = "gsf",
-    version = "1.2",
+    version = "1.3",
     author = "Takahiro Morishita",
     author_email = "tmorishita@stsci.edu",
     description = "SED Fitting Code for HST Grism",

From bfff28094406adc209fecc3862a4866dcdf56c0a Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Tue, 14 Jul 2020 11:36:49 -0400
Subject: [PATCH 2/9] all plot minor update

---
 gsf/plot_sed.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py
index f8e4ed0..0ef73a3 100644
--- a/gsf/plot_sed.py
+++ b/gsf/plot_sed.py
@@ -2061,7 +2061,7 @@ def density_estimation(m1, m2):
     ax0.set_xlabel('Observed wavelength ($\mathrm{\mu m}$)', fontsize=14)
     ax0.set_ylabel('Flux ($\mathrm{erg}/\mathrm{s}/\mathrm{cm}^{2}/\mathrm{\AA}$)', fontsize=13)
     ax1.set_xlabel('$t$ (Gyr)', fontsize=12)
-    ax1.set_ylabel('$\log \dot{M_*}/M_\odot$yr$^{-1}$', fontsize=12)
+    ax1.set_ylabel('$\dot{M_*}/M_\odot$yr$^{-1}$', fontsize=12)
     ax1.set_xlim(0.008, Txmax)
     ax1.set_ylim(0, SFmax)
     ax1.set_xscale('log')

From 992e337c725be61cca6af595fc8608daf533147d Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Tue, 14 Jul 2020 11:37:19 -0400
Subject: [PATCH 3/9] Catalog function update for those bright

---
 gsf/maketmp_filt.py | 29 ++++++++++++++++-------------
 1 file changed, 16 insertions(+), 13 deletions(-)

diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py
index 1b9557d..9a8d6d1 100644
--- a/gsf/maketmp_filt.py
+++ b/gsf/maketmp_filt.py
@@ -18,14 +18,16 @@
 
 def sim_spec(lmin, fin, sn):
     '''
-    ###################################################
-    ### SIMULATION of SPECTRA.
-    ###################################################
-    #
-    # wave_obs, wave_temp, flux_temp, sn_obs
-    # Return: frand, erand
-    #
+    Purpose:
+    ========
+    SIMULATION of SPECTRA.
+    
+    Input:
+    ======
+    wave_obs, wave_temp, flux_temp, sn_obs
+    Return: frand, erand
     '''
+
     frand = fin * 0
     erand = fin * 0
     for ii in range(len(lmin)):
@@ -40,10 +42,12 @@ def sim_spec(lmin, fin, sn):
 
 def maketemp(MB):
     '''
-    ###################################################
-    # Make SPECTRA at given z and filter set.
-    ###################################################
-    #
+    Purpose:
+    ========
+    Make SPECTRA at given z and filter set.
+    
+    Input:
+    ======
     # inputs      : Configuration file.
     # zbest(float): Best redshift at this iteration. Templates are generated based on this reshift.
     # Z (array)   : Stellar phase metallicity in logZsun.
@@ -61,7 +65,6 @@ def maketemp(MB):
     DIR_TMP = MB.DIR_TMP# './templates/'
     zbest = MB.zgal
     tau0 = MB.tau0
-    #tau0 = [float(x.strip()) for x in tau0.split(',')]
 
     fnc  = MB.fnc #Func(ID, PA, Z, nage) # Set up the number of Age/ZZ
     bfnc = MB.bfnc #Basic(Z)
@@ -566,7 +569,7 @@ def maketemp(MB):
     ##########################################
     lamliml = 0.
     lamlimu = 20000.
-    ebblim  = 1e5
+    ebblim  = 1e10
     ncolbb  = 10000
     fw = open(DIR_TMP + 'spec_obs_' + ID + '_PA' + PA + '.cat', 'w')
     fw.write('# BB data (>%d) in this file are not used in fitting.\n'%(ncolbb))

From 8f661abf46262487e40b645753c478b415a1e8b4 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Tue, 14 Jul 2020 15:47:36 -0400
Subject: [PATCH 4/9] plot functions have assined SED/SFH

---
 gsf/fitting.py  | 104 +++++++++-----------------------
 gsf/plot_sed.py | 155 ++++++++++++++++++++----------------------------
 gsf/plot_sfh.py |  91 +++++++---------------------
 3 files changed, 112 insertions(+), 238 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index 27c1973..4b6afe4 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -513,11 +513,11 @@ def residual_z(pars,z):
     def fit_redshift(self, dict, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, zlimu=6., snlim=0):
         '''
         Purpose:
-        ==========
+        ========
         Find an optimal redshift, before going into a big fit, by using several templates.
 
         Input:
-        ==========
+        ======
         delzz  : Delta z in redshift search space
         zliml  : Lower limit range for redshift
         zlimu  : Upper limit range for redshift
@@ -525,7 +525,7 @@ def fit_redshift(self, dict, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
         snlim  : SN limit for data points. Those below the number will be cut from the fit.
 
         Note:
-        ==========
+        =====
         Can be used for any SFH.
 
         '''
@@ -538,7 +538,6 @@ def fit_redshift(self, dict, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
         sn = dict['fy']/dict['ey']
         # Only spec data?
         con_cz = (dict['NR']<10000) & (sn>snlim)
-        #con_cz = (dict['NR']<100000) & (sn>snlim)
         fy_cz  = dict['fy'][con_cz] # Already scaled by self.Cz0
         ey_cz  = dict['ey'][con_cz]
         x_cz   = dict['x'][con_cz] # Observed range
@@ -800,19 +799,24 @@ def add_param(self, fit_params, sigz=1.0):
 
 
     #def main(self, zgal, flag_m, Cz0, Cz1, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, f_move=False):
-    def main(self, flag_m, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, f_move=False, verbose=False):
+    def main(self, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, f_move=False, verbose=False, skip_fitz=False):
         '''
         Input:
-        ========
-        flag_m : related to redshift error in redshift check func.
+        ======
         ferr   : For error parameter
         zgal   : Input redshift.
-        #
-        #
-        # sigz (float): confidence interval for redshift fit.
-        # ezmin (float): minimum error in redshift.
-        #
+        skip_fitz (bool): Skip redshift fit.
+        sigz (float): confidence interval for redshift fit.
+        ezmin (float): minimum error in redshift.p
+
         '''
+        import emcee
+        try:
+            import multiprocess
+        except:
+            import multiprocessing as multiprocess
+
+        from .posterior_flexible import Post
 
         print('########################')
         print('### Fitting Function ###')
@@ -886,7 +890,6 @@ def main(self, flag_m, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0
         self.dict = dict
 
         # Call likelihood/prior/posterior function;
-        from .posterior_flexible import Post
         class_post = Post(self)
 
         ###############################
@@ -986,22 +989,20 @@ def main(self, flag_m, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0
         AA_tmp = np.zeros(len(self.age), dtype='float64')
         ZZ_tmp = np.zeros(len(self.age), dtype='float64')
         fm_tmp, xm_tmp = fnc.tmp04_val(out, self.zgal, self.lib)
-        #print(self.lib[:,0])
 
         ########################
         # Check redshift
         ########################
-        flag_z = self.fit_redshift(dict, xm_tmp, fm_tmp)
+        if skip_fitz:
+            flag_z = 'y'
+        else:
+            flag_z = self.fit_redshift(dict, xm_tmp, fm_tmp)
 
         #################################################
         # Gor for mcmc phase
         #################################################
         if flag_z == 'y' or flag_z == '':
-            #zrecom  = self.zprev
-            #Czrec0  = self.Cz0
-            #Czrec1  = self.Cz1
 
-            # plot z-distribution
             self.get_zdist()
 
             #######################
@@ -1040,17 +1041,11 @@ def main(self, flag_m, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0
 
             ################################
             print('\nMinimizer Defined\n')
-            import emcee
             mini = Minimizer(class_post.lnprob, out.params, fcn_args=[dict['fy'],dict['ey'],dict['wht2'],self.f_dust], f_disp=self.f_disp, moves=emcee.moves.DEMove(sigma=1e-05, gamma0=None)) #, f_move=f_move)
             print('######################')
             print('### Starting emcee ###')
             print('######################')
 
-            try:
-                import multiprocess
-            except:
-                import multiprocessing as multiprocess
-
             ncpu0 = int(multiprocess.cpu_count()/2)
             try:
                 ncpu = int(inputs['NCPU'])
@@ -1086,44 +1081,6 @@ def main(self, flag_m, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0
                 print('### Saving chain took %.1f sec'%(tcalc_mc))
                 print('#################################')
 
-            '''
-            Avmc = np.percentile(res.flatchain['Av'], [16,50,84])
-            Avpar = np.zeros((1,3), dtype='float64')
-            Avpar[0,:] = Avmc
-
-            ####################
-            # Best parameters
-            ####################
-            Amc  = np.zeros((len(self.age),3), dtype='float64')
-            Ab   = np.zeros(len(self.age), dtype='float64')
-            Zmc  = np.zeros((len(self.age),3), dtype='float64')
-            Zb   = np.zeros(len(self.age), dtype='float64')
-            NZbest = np.zeros(len(self.age), dtype='int')
-            f0     = fits.open(self.DIR_TMP + 'ms_' + self.ID + '_PA' + self.PA + '.fits')
-            sedpar = f0[1]
-            ms     = np.zeros(len(self.age), dtype='float64')
-            for aa in range(len(self.age)):
-                Ab[aa]    = res.params['A'+str(aa)].value
-                Amc[aa,:] = np.percentile(res.flatchain['A'+str(aa)], [16,50,84])
-                try:
-                    Zb[aa]    = res.params['Z'+str(aa)].value
-                    Zmc[aa,:] = np.percentile(res.flatchain['Z'+str(aa)], [16,50,84])
-                except:
-                    Zb[aa]    = res.params['Z0'].value
-                    Zmc[aa,:] = np.percentile(res.flatchain['Z0'], [16,50,84])
-                NZbest[aa]= bfnc.Z2NZ(Zb[aa])
-                ms[aa]    = sedpar.data['ML_' +  str(NZbest[aa])][aa]
-
-            Avb = res.params['Av'].value
-
-            if self.f_dust:
-                Mdust_mc = np.zeros(3, dtype='float64')
-                Tdust_mc = np.zeros(3, dtype='float64')
-                Mdust_mc[:] = np.percentile(res.flatchain['MDUST'], [16,50,84])
-                Tdust_mc[:] = np.percentile(res.flatchain['TDUST'], [16,50,84])
-                print(Mdust_mc)
-                print(Tdust_mc)
-            '''
 
             ####################
             # MCMC corner plot.
@@ -1136,7 +1093,7 @@ def main(self, flag_m, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0
                 plot_density=False, levels=[0.68, 0.95, 0.997], truth_color='gray', color='#4682b4')
                 fig1.savefig('SPEC_' + self.ID + '_PA' + self.PA + '_corner.pdf')
 
-            # Do analysis on MCMC results.
+            # Analyze MCMC results.
             # Write to file.
             stop  = timeit.default_timer()
             tcalc = stop - start
@@ -1147,7 +1104,7 @@ def main(self, flag_m, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0
             stop_mc  = timeit.default_timer()
             tcalc_mc = stop_mc - start_mc
 
-            return False #, self.zgal, self.Cz0, self.Cz1
+            return False
 
 
         elif flag_z == 'm':
@@ -1176,7 +1133,7 @@ def main(self, flag_m, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0
             print('\n\n')
             print('Generate model templates with input redshift and Scale.')
             print('\n\n')
-            return True #, self.zgal, self.Cz0, self.Cz1
+            return True
 
         else:
             print('\n\n')
@@ -1192,18 +1149,18 @@ def main(self, flag_m, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0
                 self.Cz0   = self.Czrec0
                 self.Cz1   = self.Czrec1
 
-                return True #, self.zgal, self.Cz0, self.Cz1
+                return True
 
             else:
                 print('\n\n')
                 print('There is nothing to do.')
                 print('Terminating process.')
                 print('\n\n')
-                return -1 #, self.zgal, self.Czrec0, self.Czrec1
+                return -1
 
 
 
-    def quick_fit(self, zgal, flag_m, Cz0, Cz1, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, f_move=False):
+    def quick_fit(self, zgal, Cz0, Cz1, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, f_move=False):
         '''
         Purpose:
         ==========
@@ -1211,14 +1168,11 @@ def quick_fit(self, zgal, flag_m, Cz0, Cz1, specplot=1, sigz=1.0, ezmin=0.01, fe
 
         Input:
         ==========
-        flag_m : related to redshift error in redshift check func.
         ferr   : For error parameter
         zgal   : Input redshift.
-        #
-        #
-        # sigz (float): confidence interval for redshift fit.
-        # ezmin (float): minimum error in redshift.
-        #
+        sigz (float): confidence interval for redshift fit.
+        ezmin (float): minimum error in redshift.
+
         '''
         from .posterior_flexible import Post
 
diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py
index 0ef73a3..86fd71b 100644
--- a/gsf/plot_sed.py
+++ b/gsf/plot_sed.py
@@ -19,7 +19,7 @@
 col = ['violet', 'indigo', 'b', 'lightblue', 'lightgreen', 'g', 'orange', 'coral', 'r', 'darkred']#, 'k']
 #col = ['darkred', 'r', 'coral','orange','g','lightgreen', 'lightblue', 'b','indigo','violet','k']
 
-def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, inputs=False, nmc_rand=300, dust_model=0, DIR_TMP='./templates/', f_label=False, f_bbbox=False, verbose=False, f_silence=True, f_fill=False):
+def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, inputs=False, nmc_rand=300, dust_model=0, DIR_TMP='./templates/', f_label=False, f_bbbox=False, verbose=False, f_silence=True, f_fill=False, f_fancyplot=False):
     '''
     Input:
     ============
@@ -282,27 +282,24 @@ def gaus(x,a,x0,sigma):
             pass
         try:
             conebb_ls = (fybb/eybb<=SNlim) & (eybb>0)
-            #ax1.errorbar(xbb[conebb_ls], eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) / d * SNlim, \
-            #yerr=fybb[conebb_ls]*0+np.max(fybb[conebb_ls]*c/np.square(xbb[conebb_ls])/d)*0.05, \
-            #uplims=eybb[conebb_ls]*c/np.square(xbb[conebb_ls])/d*SNlim, color='r', linestyle='', linewidth=0.5, zorder=4, ms=1, capsize=0.)
             ax1.errorbar(xbb[conebb_ls], eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) / d, \
             color=col_dat, linestyle='', linewidth=0.5, zorder=4, ms=4, marker='v')
 
         except:
             pass
 
-    # If any IRAC excess:
-    #if True:
+
+    # For any data removed fron fit (i.e. IRAC excess):
     try:
         col_ex = 'lawngreen'
         #col_ex = 'limegreen'
         #col_ex = 'r'
-        # Currently, this file is made manually;
+        # Currently, this file is made after FILTER_SKIP;
         data_ex = ascii.read(DIR_TMP + 'bb_obs_' + ID + '_PA' + PA + '_removed.cat')
-        x_ex  = data_ex['col2']
-        fy_ex = data_ex['col3']
-        ey_ex = data_ex['col4']
-        ex_ex  = data_ex['col5']
+        x_ex    = data_ex['col2']
+        fy_ex   = data_ex['col3']
+        ey_ex   = data_ex['col4']
+        ex_ex   = data_ex['col5']
 
         ax1.errorbar(x_ex, fy_ex * c / np.square(x_ex) / d, \
         xerr=ex_ex, yerr=ey_ex*c/np.square(x_ex)/d, color='k', linestyle='', linewidth=0.5, zorder=5)
@@ -353,10 +350,6 @@ def gaus(x,a,x0,sigma):
             if f_fill:
                 ax1.fill_between(x0[::nstep_plot], ((ysum - y0_r) * c/ np.square(x0) / d)[::nstep_plot], (ysum * c/ np.square(x0) / d)[::nstep_plot], linestyle='None', lw=0.5, color=col[ii], alpha=alp, zorder=-1, label='')
 
-        #if jj == len(II0) - 1:
-        #    ax1.plot(x0[::nstep_plot], (ysum * c/ np.square(x0) / d)[::nstep_plot], linestyle='-', lw=0.5, color='k', alpha=1., zorder=-1, label='')
-
-
         try:
             ysum_wid = ysum * 0
             for kk in range(0,ii+1,1):
@@ -392,24 +385,6 @@ def gaus(x,a,x0,sigma):
         ax1.plot(x0d, y0d * c/ np.square(x0d) / d, '--', lw=0.5, color='purple', zorder=-1, label='')
         ax3t.plot(x0d, y0d * c/ np.square(x0d) / d, '--', lw=0.5, color='purple', zorder=-1, label='')
 
-        '''
-        print(MD16,MD50,MD84)
-        print(nTD16,nTD50,nTD84)
-        # 16;
-        par['MDUST'].value=MD16
-        par['TDUST'].value=nTD16
-        y0d, x0d = fnc.tmp04_dust(par.valuesdict(), zbes, lib_dust_all, tau0=tau0)
-        ax3t.plot(x0d, y0d * c/ np.square(x0d) / d, '--', lw=0.5, color='purple', zorder=-1, label='')
-        # 84
-        par['MDUST'].value=MD84
-        par['TDUST'].value=nTD84
-        y0d, x0d = fnc.tmp04_dust(par.valuesdict(), zbes, lib_dust_all, tau0=tau0)
-        ax3t.plot(x0d, y0d * c/ np.square(x0d) / d, '--', lw=0.5, color='purple', zorder=-1, label='')
-        '''
-
-        #ax1.plot(x0d, y0d, '--', lw=0.5, color='purple', zorder=-1, label='')
-        #ax3t.plot(x0d, y0d, '--', lw=0.5, color='purple', zorder=-1, label='')
-
         # data;
         ddat  = np.loadtxt(DIR_TMP + 'bb_dust_obs_' + ID + '_PA' + PA + '.cat', comments='#')
         NRbbd = ddat[:, 0]
@@ -418,6 +393,7 @@ def gaus(x,a,x0,sigma):
         eybbd = ddat[:, 3]
         exbbd = ddat[:, 4]
         snbbd = fybbd/eybbd
+
         try:
             conbbd_hs = (fybbd/eybbd>SNlim)
             ax1.errorbar(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) / d, \
@@ -428,6 +404,7 @@ def gaus(x,a,x0,sigma):
             '.r', linestyle='', linewidth=0, zorder=4)#, label='Obs.(BB)')
         except:
             pass
+
         try:
             conebbd_ls = (fybbd/eybbd<=SNlim)
             ax1.errorbar(xbbd[conebbd_ls], eybbd[conebbd_ls] * c / np.square(xbbd[conebbd_ls]) / d, \
@@ -556,9 +533,9 @@ def gaus(x,a,x0,sigma):
     #
     # From MCMC chain
     #
-    file = 'chain_' + ID + '_PA' + PA + '_corner.cpkl'
+    file  = 'chain_' + ID + '_PA' + PA + '_corner.cpkl'
     niter = 0
-    data = loadcpkl(os.path.join('./'+file))
+    data  = loadcpkl(os.path.join('./'+file))
     try:
         ndim   = data['ndim']     # By default, use ndim and burnin values contained in the cpkl file, if present.
         burnin = data['burnin']
@@ -629,7 +606,6 @@ def gaus(x,a,x0,sigma):
 
             # Each;
             ytmp_each[kk,:,ss] = ferr_tmp * mod0_tmp[:] * c / np.square(xm_tmp[:]) / d
-            #print(ferr_tmp, AA_tmp, Av_tmp, ss, ZZ_tmp, zbes)
 
         #
         # Dust component;
@@ -693,7 +669,17 @@ def gaus(x,a,x0,sigma):
         ax1.fill_between(x1_tot[::nstep_plot], ytmp16[::nstep_plot], ytmp84[::nstep_plot], ls='-', lw=.5, color='gray', zorder=-2, alpha=0.5)
     ax1.plot(x1_tot[::nstep_plot], ytmp50[::nstep_plot], '-', lw=.5, color='gray', zorder=-1, alpha=1.)
 
-    f_fancyplot=False
+    # Attach the data point in MB;
+    MB.sed_wave_obs   = xbb
+    MB.sed_flux_obs   = fybb * c / np.square(xbb) / d
+    MB.sed_eflux_obs  = eybb * c / np.square(xbb) / d
+    # Attach the best SED to MB;
+    MB.sed_wave   = x1_tot
+    MB.sed_flux16 = ytmp16
+    MB.sed_flux50 = ytmp50
+    MB.sed_flux84 = ytmp84
+
+
     if f_fancyplot:
         # For each age;
         ytmp_each50       = np.zeros(len(xm_tmp), dtype='float64')
@@ -705,12 +691,6 @@ def gaus(x,a,x0,sigma):
             ax1.fill_between(x1_tot[::nstep_plot], ytmp_each50_prior[::nstep_plot], ytmp_each50[::nstep_plot], linestyle='None', lw=0.5, color=col[ii])
             ytmp_each50_prior[:] += ytmp_each50[:]
 
-    '''
-    if f_dust:
-        l50_bb, ytmp50_bb, l50_fwhm = filconv(ALLFILT, x1_tot, ytmp50, DIR_FILT, fw=True)
-    else:
-        l50_bb, ytmp50_bb, l50_fwhm = filconv(SFILT, x1_tot, ytmp50, DIR_FILT, fw=True)
-    '''
 
     #########################
     # Calculate non-det chi2
@@ -729,7 +709,6 @@ def func_tmp(xint,eobs,fmodel):
 
 
     chi2 = sum((np.square(fy-ysump)*np.sqrt(wht3))[conw])
-    #chi2 = sum((np.square(fy-ysump))[conw])
 
     # Effective ndim;
     ndim_eff = MB.ndim
@@ -785,7 +764,7 @@ def func_tmp(xint,eobs,fmodel):
         zorder=2, alpha=1.0, marker='d', s=50)
 
     else:
-        lbb, fbb, lfwhm = filconv(SFILT, x1_tot, ytmp50, DIR_FILT, fw=True)
+        lbb, fbb, lfwhm   = filconv(SFILT, x1_tot, ytmp50, DIR_FILT, fw=True)
         lbb, fbb16, lfwhm = filconv(SFILT, x1_tot, ytmp16, DIR_FILT, fw=True)
         lbb, fbb84, lfwhm = filconv(SFILT, x1_tot, ytmp84, DIR_FILT, fw=True)
 
@@ -1043,6 +1022,8 @@ def func_tmp(xint,eobs,fmodel):
 def plot_corner_TZ(ID, PA, Zall=np.arange(-1.2,0.4249,0.05), age=[0.01, 0.1, 0.3, 0.7, 1.0, 3.0]):
     '''
     '''
+    import matplotlib
+    import matplotlib.cm as cm
     col = ['violet', 'indigo', 'b', 'lightblue', 'lightgreen', 'g', 'orange', 'coral', 'r', 'darkred']#, 'k']
     nage = np.arange(0,len(age),1)
     fnc  = Func(ID, PA, Zall, age, dust_model=dust_model) # Set up the number of Age/ZZ
@@ -1074,8 +1055,6 @@ def plot_corner_TZ(ID, PA, Zall=np.arange(-1.2,0.4249,0.05), age=[0.01, 0.1, 0.3
     f0     = fits.open(DIR_TMP + 'ms_' + ID + '_PA' + PA + '.fits')
     sedpar = f0[1]
 
-    import matplotlib
-    import matplotlib.cm as cm
     getcmap   = matplotlib.cm.get_cmap('jet')
     nc        = np.arange(0, nmc, 1)
     col = getcmap((nc-0)/(nmc-0))
@@ -1583,9 +1562,7 @@ def write_lines(ID, PA, zbes, R_grs=45, dw=4, umag=1.0):
 def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1000):
     '''
     Purpose:
-    ==========
-    Creat temporal png for gif image.
-
+    ========
     For summary. In the same format as plot_corner_physparam_frame.
 
     '''
@@ -1605,9 +1582,13 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1
         age = MB.age_fix
     except:
         age  = MB.age
+
     tau0 = MB.tau0 #[0.1,0.2,0.3]
     dust_model = MB.dust_model
-    DIR_TMP = MB.DIR_TMP #'./templates/'
+    DIR_TMP = MB.DIR_TMP
+
+    #Txmax = 4 # Max x value
+    Txmax = np.max(age) + 1.0 # Max x value
 
     ###########################
     # Open result file
@@ -1636,14 +1617,17 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1
         A16[aa] = hdul[1].data['A'+str(aa)][0]
         A84[aa] = hdul[1].data['A'+str(aa)][2]
 
-    Asum  = np.sum(A50)
+    Asum = np.sum(A50)
     aa   = 0
-    Av50 = hdul[1].data['Av'+str(aa)][1]
+
     Av16 = hdul[1].data['Av'+str(aa)][0]
+    Av50 = hdul[1].data['Av'+str(aa)][1]
     Av84 = hdul[1].data['Av'+str(aa)][2]
-    Z50 = np.zeros(len(age), dtype='float64')
-    Z16 = np.zeros(len(age), dtype='float64')
-    Z84 = np.zeros(len(age), dtype='float64')
+
+    Z50  = np.zeros(len(age), dtype='float64')
+    Z16  = np.zeros(len(age), dtype='float64')
+    Z84  = np.zeros(len(age), dtype='float64')
+
     NZbest = np.zeros(len(age), dtype='int')
     for aa in range(len(age)):
         Z50[aa] = hdul[1].data['Z'+str(aa)][1]
@@ -1651,16 +1635,14 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1
         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.
+    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)
-
+    zscl  = (1.+zbes)
 
     try:
         os.makedirs(DIR_OUT)
@@ -1668,14 +1650,14 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1
         pass
 
     # plot Configuration
-    K = 4 # No of params.
     Par = ['$\log M_*/M_\odot$', '$\log T_*$/Gyr', '$A_V$/mag', '$\log Z_* / Z_\odot$']
+    K   = len(Par) # No of params.
     factor = 2.0           # size of one side of one panel
-    lbdim = 0.5 * factor   # size of left/bottom margin
-    trdim = 0.2 * factor   # size of top/right margin
-    whspace = 0.02         # w/hspace size
-    plotdim = factor * K + factor * (K - 1.) * whspace
-    dim = lbdim + plotdim + trdim
+    lbdim  = 0.5 * factor   # size of left/bottom margin
+    trdim  = 0.2 * factor   # size of top/right margin
+    whspace= 0.02         # w/hspace size
+    plotdim= factor * K + factor * (K - 1.) * whspace
+    dim    = lbdim + plotdim + trdim
     sclfig = 0.7
 
 
@@ -1698,21 +1680,18 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1
     else:
         MB.dict = MB.read_data(MB.Cz0, MB.Cz1, MB.zgal)
 
-    #dat  = np.loadtxt(DIR_TMP + 'spec_obs_' + ID + '_PA' + PA + '.cat', comments='#')
-    NRbb   = MB.dict['NRbb'] #dat[:, 0]
-    xbb    = MB.dict['xbb'] #dat[:, 1]
-    fybb   = MB.dict['fybb'] #dat[:, 2]
-    eybb   = MB.dict['eybb'] #dat[:, 3]
+    # 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
+    snbb   = fybb/eybb
 
-    # Create a new figure if one wasn't provided.
-    ###############################
-    # Data taken from
-    ###############################
+    # Get spec data points;
     dat = np.loadtxt(DIR_TMP + 'spec_obs_' + ID + '_PA' + PA + '.cat', comments='#')
-    NR = dat[:, 0]
-    x  = dat[:, 1]
+    NR  = dat[:, 0]
+    x   = dat[:, 1]
     fy00 = dat[:, 2]
     ey00 = dat[:, 3]
 
@@ -1729,27 +1708,23 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1
     fg2  = fy00[con2]
     eg2  = ey00[con2]
     fy01 = np.append(fg0,fg1)
-    fy   = np.append(fy01,fg2)
     ey01 = np.append(eg0,eg1)
+
+    fy   = np.append(fy01,fg2)
     ey   = np.append(ey01,eg2)
     wht=1./np.square(ey)
 
-
     conspec = (NR<10000) #& (fy/ey>1)
-    #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>1)
+    conbb   = (fybb/eybb>1)
     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)
     ax0.plot(xbb[conbb], fybb[conbb] * c / np.square(xbb[conbb]) / d, '.r', ms=10, linestyle='', linewidth=0, zorder=4)#, label='Obs.(BB)')
-    #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='')
 
     ####################
     # MCMC corner plot.
     ####################
-    file = 'chain_' + ID + '_PA' + PA + '_corner.cpkl'
+    file  = 'chain_' + ID + '_PA' + PA + '_corner.cpkl'
     niter = 0
-    data = loadcpkl(os.path.join('./'+file))
+    data  = loadcpkl(os.path.join('./'+file))
 
     try:
         ndim   = data['ndim']     # By default, use ndim and burnin values contained in the cpkl file, if present.
@@ -1777,7 +1752,6 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1
     ACtmp= np.zeros(nplot, dtype='float64')
 
     # Time bin
-    Txmax = 4 # Max x value
     Tuni = MB.cosmo.age(zbes).value
     Tuni0 = (Tuni - age[:])
     delT  = np.zeros(len(age),dtype='float64')
@@ -1817,7 +1791,7 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1
     delTu[:] *= 1e9 # Gyr to yr
 
     ######
-    files = [] # For movie
+    files = [] # For gif animation
     SFmax = 0
     Tsmin = 0
     Tsmax = 0
@@ -1833,8 +1807,6 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1
         AA_tmp = np.max(samples['A'+str(ii)][:])
         AA_tmp84 = np.percentile(samples['A'+str(ii)][:],95)
         AA_tmp16 = np.percentile(samples['A'+str(ii)][:],5)
-        #AA_tmp84 = A84[ii]
-        #AA_tmp16 = A16[ii]
         nZtmp  = bfnc.Z2NZ(ZZ_tmp)
         mslist = sedpar.data['ML_'+str(nZtmp)][ii]
         AMtmp16 += mslist*AA_tmp16
@@ -1857,7 +1829,6 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', nplot=1
     # For redshift
     if zbes<2:
         zred  = [zbes, 2, 3, 6]
-        #zredl = ['$z_\mathrm{obs.}$', 2, 3, 6]
         zredl = ['$z_\mathrm{obs.}$', 2, 3, 6]
     elif zbes<2.5:
         zred  = [zbes, 2.5, 3, 6]
diff --git a/gsf/plot_sfh.py b/gsf/plot_sfh.py
index 0e39e32..3431f7d 100644
--- a/gsf/plot_sfh.py
+++ b/gsf/plot_sfh.py
@@ -23,14 +23,12 @@
 def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, fil_path = './FILT/', inputs=None, dust_model=0, DIR_TMP='./templates/',f_SFMS=False, verbose=False):
     '''
     Purpose:
-    ==========
-
+    ========
     Star formation history plot.
 
 
     Input:
-    ==========
-
+    ======
     flim  : Lower limit for plotting an age bin.
     lsfrl : Lower limit for SFR, in logMsun/yr
 
@@ -90,16 +88,6 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     zbes = hdul[0].header['z']
     chinu= hdul[1].data['chi']
 
-    '''
-    rek  = 0
-    erekl= 0
-    ereku= 0
-    mu = 1.0
-    nn = 0
-    qq = 0
-    enn = 0
-    eqq = 0
-    '''
     try:
         RA   = hdul[0].header['RA']
         DEC  = hdul[0].header['DEC']
@@ -125,7 +113,6 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
         else:
             SN = 1
 
-
     Asum = 0
     A50 = np.arange(len(age), dtype='float32')
     for aa in range(len(A50)):
@@ -144,6 +131,7 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     delT  = np.zeros(len(age),dtype='float32')
     delTl = np.zeros(len(age),dtype='float32')
     delTu = np.zeros(len(age),dtype='float32')
+
     if len(age) == 1:
         for aa in range(len(age)):
             try:
@@ -153,16 +141,7 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
             delTl[aa] = tau_ssp/2
             delTu[aa] = tau_ssp/2
             delT[aa]  = delTu[aa] + delTl[aa]
-        """
-        elif MB.f_bpass == 1:
-            file_all = MB.DIR_TMP + 'spec_all.fits'
-            hd = fits.open(file_all)[1].header
-            for aa in range(len(age)):
-                tau_ssp   = hd['realtau%d(Gyr)'%(aa)]
-                delT[aa]  = tau_ssp
-                delTl[aa] = tau_ssp/2
-                delTu[aa] = tau_ssp/2
-        """
+
     else: # This is only true when CSP...
         for aa in range(len(age)):
             if aa == 0:
@@ -186,6 +165,7 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     delT[:]  *= 1e9 # Gyr to yr
     delTl[:] *= 1e9 # Gyr to yr
     delTu[:] *= 1e9 # Gyr to yr
+
     ##############################
     # Load Pickle
     ##############################
@@ -195,7 +175,6 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     niter = 0
     data = loadcpkl(os.path.join(samplepath+'/'+pfile))
     try:
-    #if 1>0:
         ndim   = data['ndim']     # By default, use ndim and burnin values contained in the cpkl file, if present.
         burnin = data['burnin']
         nmc    = data['niter']
@@ -211,8 +190,6 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     ######################
     # Mass-to-Light ratio.
     ######################
-    #ms     = np.zeros(len(age), dtype='float32')
-    # Wht do you want from MCMC sampler?
     AM = np.zeros((len(age), mmax), dtype='float32') # Mass in each bin.
     AC = np.zeros((len(age), mmax), dtype='float32') # Cumulative mass in each bin.
     AL = np.zeros((len(age), mmax), dtype='float32') # Cumulative light in each bin.
@@ -366,7 +343,6 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     SFR_SED_med = np.percentile(SFR_SED[:],[16,50,84])
     f_SFRSED_plot = False
     if f_SFRSED_plot:
-        #ax1.plot(delt_tot/2./1e9, np.log10(SFR_SED_med[1]), color='k', marker='*',ms=10)
         ax1.errorbar(delt_tot/2./1e9, np.log10(SFR_SED_med[1]), xerr=[[delt_tot/2./1e9],[delt_tot/2./1e9]], \
         yerr=[[np.log10(SFR_SED_med[1])-np.log10(SFR_SED_med[0])],[np.log10(SFR_SED_med[2])-np.log10(SFR_SED_med[1])]], \
         linestyle='', color='orange', lw=1., marker='*',ms=8,zorder=-2)
@@ -382,6 +358,7 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     ax1.scatter(age[conA], np.log10(SFp[:,1])[conA], marker='.', c='k', s=msize[conA])
     ax1.errorbar(age[conA], np.log10(SFp[:,1])[conA], xerr=[delTl[:][conA]/1e9,delTu[:][conA]/1e9], yerr=[np.log10(SFp[:,1])[conA]-np.log10(SFp[:,0])[conA], np.log10(SFp[:,2])[conA]-np.log10(SFp[:,1])[conA]], linestyle='-', color='k', lw=0.5, marker='')
 
+
     #############
     # Get SFMS in log10;
     #############
@@ -487,13 +464,10 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
         %(t0-delTl[ii]/1e9, t0+delTl[ii]/1e9, SFp[ii,1], SFp[ii,0], SFp[ii,2], ACp[ii,0], ACp[ii,1], ACp[ii,2], ZCp[ii,1], ZCp[ii,0], ZCp[ii,2]))
     fw_sfr.close()
 
-
     #############
     # Axis
     #############
-    #ax1.set_xlabel('$t$ (Gyr)', fontsize=12)
     ax1.set_ylabel('$\log \dot{M}_*/M_\odot$yr$^{-1}$', fontsize=12)
-    #ax1.set_ylabel('$\log M_*/M_\odot$', fontsize=12)
 
     lsfru = 2.8
     if np.max(np.log10(SFp[:,2]))>2.8:
@@ -502,12 +476,10 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     if f_comp == 1:
         lsfru = np.max([lsfru, np.log10(np.max(sfr_exp*C_exp))])
 
-
     ax1.set_xlim(0.008, Txmax)
     ax1.set_ylim(lsfrl, lsfru)
     ax1.set_xscale('log')
 
-    #ax2.set_xlabel('$t$ (Gyr)', fontsize=12)
     ax2.set_ylabel('$\log M_*/M_\odot$', fontsize=12)
 
     ax2.set_xlim(0.008, Txmax)
@@ -527,14 +499,6 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     prihdr['z']      = zbes
     prihdr['RA']     = RA
     prihdr['DEC']    = DEC
-    '''
-    prihdr['Re_kpc'] = rek
-    prihdr['Ser_n']  = nn
-    prihdr['Axis_q'] = qq
-    prihdr['e_Re']   = (erekl+ereku)/2.
-    prihdr['e_n']    = enn
-    prihdr['e_q']    = eqq
-    '''
     # Add rejuv properties;
     prihdr['f_rejuv']= f_rejuv
     prihdr['t_quen'] = t_quench
@@ -584,42 +548,30 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     col50 = fits.Column(name='AV', format='E', unit='mag', array=para[:])
     col01.append(col50)
 
-    '''
-    Ths is in gsf_sed_*.fits
-    # U-V
-    para = [uv[0], uv[1], uv[2]]
-    col50 = fits.Column(name='U-V', format='E', unit='mag', array=para[:])
-    col01.append(col50)
-
-    # V-J
-    para = [vj[0], vj[1], vj[2]]
-    col50 = fits.Column(name='V-J', format='E', unit='mag', array=para[:])
-    col01.append(col50)
-    '''
-
+    # 
     colms  = fits.ColDefs(col01)
     dathdu = fits.BinTableHDU.from_columns(colms)
     hdu    = fits.HDUList([prihdu, dathdu])
     hdu.writeto('SFH_' + ID + '_PA' + PA + '_param.fits', overwrite=True)
 
+    # Attach to MB;
+    MB.sfh_tlook = age
+    MB.sfh_tlookl= delTl[:][conA]/1e9
+    MB.sfh_tlooku= delTu[:][conA]/1e9
+    MB.sfh_sfr16 = SFp[:,0]
+    MB.sfh_sfr50 = SFp[:,1]
+    MB.sfh_sfr84 = SFp[:,2]
+    MB.sfh_mfr16 = ACp[:,0]
+    MB.sfh_mfr50 = ACp[:,1]
+    MB.sfh_mfr84 = ACp[:,2]
+    MB.sfh_zfr16 = ZCp[:,0]
+    MB.sfh_zfr50 = ZCp[:,1]
+    MB.sfh_zfr84 = ZCp[:,2]
 
-    #
     # SFH
-    #
     zzall = np.arange(1.,12,0.01)
     Tall  = MB.cosmo.age(zzall).value # , use_flat=True, **cosmo)
 
-    '''
-    fw = open('SFH_' + ID + '_PA' + PA + '_sfh.cat', 'w')
-    fw.write('%s'%(ID))
-    for mm in range(len(age)):
-        mmtmp = np.argmin(np.abs(Tall - Tuni0[mm]))
-        zztmp = zzall[mmtmp]
-        fw.write(' %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f'%(zztmp, Tuni0[mm], np.log10(ACp[mm,1]), (np.log10(ACp[mm,1])-np.log10(ACp[mm,0])), (np.log10(ACp[mm,2])-np.log10(ACp[mm,1])), ZCp[mm,1], ZCp[mm,1]-ZCp[mm,0], ZCp[mm,2]-ZCp[mm,1], SFp[mm,1], SFp[mm,1]-SFp[mm,0], SFp[mm,2]-SFp[mm,1]))
-    fw.write('\n')
-    fw.close()
-    '''
-
     dely2 = 0.1
     while (y2max-y2min)/dely2>7:
         dely2 *= 2.
@@ -646,7 +598,6 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     if zbes<4:
         if zbes<2:
             zred  = [zbes, 2, 3, 6]
-            #zredl = ['$z_\mathrm{obs.}$', 2, 3, 6]
             zredl = ['$z_\mathrm{obs.}$', 2, 3, 6]
         elif zbes<2.5:
             zred  = [zbes, 2.5, 3, 6]
@@ -670,7 +621,6 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
         if Tzz[zz] < 0.01:
             Tzz[zz] = 0.01
 
-    #print(zred, Tzz)
     #ax3t.set_xscale('log')
     #ax3t.set_xlim(0.008, Txmax)
 
@@ -693,7 +643,6 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     ax4.set_yticks([-0.8, -0.4, 0., 0.4])
     ax4.set_yticklabels(['-0.8', '-0.4', '0', '0.4'])
     #ax4.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
-
     #ax3.yaxis.labelpad = -2
     ax4.yaxis.labelpad = -2
 

From efaed2b360e3c87462b8c92e9326050a772120ed Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Tue, 14 Jul 2020 15:48:03 -0400
Subject: [PATCH 5/9] minor

---
 gsf/maketmp_filt.py | 19 ++++++++-----------
 gsf/writing.py      | 32 ++++++++++++++++++++++++--------
 gsf/zfit.py         | 10 ++++------
 3 files changed, 36 insertions(+), 25 deletions(-)

diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py
index 9a8d6d1..575e20a 100644
--- a/gsf/maketmp_filt.py
+++ b/gsf/maketmp_filt.py
@@ -40,7 +40,7 @@ def sim_spec(lmin, fin, sn):
     return frand, erand
 
 
-def maketemp(MB):
+def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=20000., ncolbb=10000):
     '''
     Purpose:
     ========
@@ -48,13 +48,14 @@ def maketemp(MB):
     
     Input:
     ======
-    # inputs      : Configuration file.
-    # zbest(float): Best redshift at this iteration. Templates are generated based on this reshift.
-    # Z (array)   : Stellar phase metallicity in logZsun.
-    # age (array) : Age, in Gyr.
-    # fneb (int)  : flag for adding nebular emissionself.
-    #
+    inputs      : Configuration file.
+    zbest(float): Best redshift at this iteration. Templates are generated based on this reshift.
+    Z (array)   : Stellar phase metallicity in logZsun.
+    age (array) : Age, in Gyr.
+    fneb (int)  : flag for adding nebular emissionself.
+    
     '''
+    
     inputs = MB.inputs
     ID = MB.ID #inputs['ID']
     PA = MB.PA #inputs['PA']
@@ -567,10 +568,6 @@ def maketemp(MB):
     # For observation.
     # Write out for the Multi-component fitting.
     ##########################################
-    lamliml = 0.
-    lamlimu = 20000.
-    ebblim  = 1e10
-    ncolbb  = 10000
     fw = open(DIR_TMP + 'spec_obs_' + ID + '_PA' + PA + '.cat', 'w')
     fw.write('# BB data (>%d) in this file are not used in fitting.\n'%(ncolbb))
     for ii in range(len(lm)):
diff --git a/gsf/writing.py b/gsf/writing.py
index c6408c5..301f960 100644
--- a/gsf/writing.py
+++ b/gsf/writing.py
@@ -4,21 +4,37 @@
 
 from .function import filconv, calc_Dn4
 
+
 def get_param(self, res, fitc, tcalc=1., burnin=-1):
     '''
+    Purpose:
+    ========
+    Write a parameter file.
+
     '''
     print('##########################')
     print('### Writing parameters ###')
     print('##########################')
     lib_all = self.lib_all
-    zrecom = self.zrecom
-    Czrec0 = self.Czrec0
-    Czrec1 = self.Czrec1
-    z_cz   = self.z_cz
-    scl_cz0= self.scl_cz0
-    scl_cz1= self.scl_cz1
-    tau0   = self.tau0
 
+    # Those are from redshiftfit;
+    #zrecom = self.zrecom
+    #Czrec0 = self.Czrec0
+    #Czrec1 = self.Czrec1
+    zrecom = self.zgal
+    Czrec0 = self.Cz0
+    Czrec1 = self.Cz1
+
+    try:
+        z_cz   = self.z_cz
+        scl_cz0= self.scl_cz0
+        scl_cz1= self.scl_cz1
+    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])
+
+    tau0= self.tau0
     ID0 = self.ID
     PA0 = self.PA
     try:
@@ -218,7 +234,7 @@ def get_param(self, res, fitc, tcalc=1., burnin=-1):
 def get_index(mmax=300):
     '''
     Purpose:
-    ==========
+    ========
     Retrieve spectral indices from each realization.
     '''
 
diff --git a/gsf/zfit.py b/gsf/zfit.py
index 35c3585..b80916e 100644
--- a/gsf/zfit.py
+++ b/gsf/zfit.py
@@ -3,15 +3,14 @@
 from .function import check_line_cz_man
 
 
-def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, prior, NR, zliml, zlimu, \
-delzz=0.01, nmc_cz=100, nwalk_cz=10):
+def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, prior, NR, zliml, zlimu, delzz=0.01, nmc_cz=100, nwalk_cz=10):
     '''
     Purpose:
-    =========
+    ========
     Fit observed flux with a template to get redshift probability.
 
     Input:
-    =========
+    ======
 
     zbest : Initial value for redshift.
     prior : Prior for redshift determination. E.g., Eazy z-probability.
@@ -22,7 +21,7 @@ def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, prior, NR, zliml, zl
     fobs, eobs, xobs: Observed spectrum. (Already scaled with Cz0prev.)
 
     Return:
-    =========
+    =======
 
     res_cz  :
     fitc_cz :
@@ -43,7 +42,6 @@ def residual_z(pars):
         Cz1s  = vals['Cz1']
 
         xm_s = xm_tmp * (1+z)
-        #fm_s = np.interp(xobs, xm_s, fm_tmp)
         fint = interpolate.interp1d(xm_s, fm_tmp, kind='nearest', fill_value="extrapolate")
         fm_s = fint(xobs)
 

From a0b2e7471320c457e25e30a5e3f61a2b2af14317 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Fri, 17 Jul 2020 15:44:51 -0400
Subject: [PATCH 6/9] Update in README

---
 README.rst          | 31 ++++++++++++++++++-------------
 gsf_environment.yml |  5 ++---
 requirements.txt    |  8 ++++++++
 3 files changed, 28 insertions(+), 16 deletions(-)
 create mode 100755 requirements.txt

diff --git a/README.rst b/README.rst
index 42d2b2e..cd33f1f 100644
--- a/README.rst
+++ b/README.rst
@@ -10,48 +10,53 @@ version 1.2
 - Far IR data set can be fitted simultaneously with a simple gray body spectrum (to be implemented in future version).
 - Due to the complex nature of grism data and code, feel free to contact the author.
 - From v1.2, BPASS templates can also be used. Those who wish to try the functionality, please contact the author.
+
 ========================================================================================
 
 
 Demonstration
-~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~
 .. image:: ./sample.png
-- Fitting movie can be found `here <https://youtu.be/pdkA9Judd-M>`__.
-
 
+- Fitting movie can be found `here <https://youtu.be/pdkA9Judd-M>`__.
 
 Pre-requirement
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~
 
 - astropy
 - Pandas
-- multiprocessing
+- multiprocess
 - lmfit (no older than ver1.0.0)
 - emcee (no older than ver3.0.0)
 - corner
 - python-fsps (ver0.3.0)
 - fsps
 
-These packages will be installed if users create a new conda environment, by executing;
-
-.. code-block:: bash
 
-    conda create --name gsf_environment.yml
+To run stellar population model, users need to install `fsps <https://github.com/cconroy20/fsps>`__ by following their instruction.
 
 
 Installation & Documentation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+Required packages will be installed by;
+
+.. code-block:: bash
+
+    pip install -r requirements.txt 
+
+
+Or if users create a new conda environment, by executing;
+
 .. code-block:: bash
 
-    git clone https://github.com/mtakahiro/gsf
-    cd gsf
-    python setup.py install
+    conda create --name gsf_environment.yml
+
 
 
 Examples
 ~~~~~~~~
-Example data set is stored in ./gsf/examples/
+An example data set is stored in ./gsf/examples/
 
 .. code-block:: bash
 
diff --git a/gsf_environment.yml b/gsf_environment.yml
index 6ecf56a..b43ea6f 100644
--- a/gsf_environment.yml
+++ b/gsf_environment.yml
@@ -55,8 +55,7 @@ dependencies:
   - xz=5.2.5=h0b31af3_0
   - zlib=1.2.11=h0b31af3_1006
   - pip:
-    - cosmolopy==0.4
-    - fsps==0.3.0
-    - gsf==1.1.0
+    - fsps
+    - git+https://github.com/mtakahiro/gsf#1.2
 prefix: /Users/tmorishita/anaconda3/envs/gsf
 
diff --git a/requirements.txt b/requirements.txt
new file mode 100755
index 0000000..f1af177
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,8 @@
+astropy==4.0.1.post1
+pandas==1.0.3
+multiprocess
+lmfit==1.0.0
+emcee==3.0.2
+corner==2.0.1
+git+https://github.com/dfm/python-fsps
+git+https://github.com/mtakahiro/gsf#1.2

From 814220d35e94f96d88f8c1d28e3dab9354acedeb Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Wed, 22 Jul 2020 09:18:06 -0400
Subject: [PATCH 7/9] format

---
 gsf/function.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/gsf/function.py b/gsf/function.py
index b8632ab..e3659a4 100644
--- a/gsf/function.py
+++ b/gsf/function.py
@@ -34,13 +34,13 @@ def get_input():
 def read_input(parfile):
     '''
     Purpose:
-    ==========
+    ========
     #
     # Get info from param file.
     #
 
     Return:
-    ===========
+    =======
     Input dictionary.
 
     '''
@@ -818,7 +818,7 @@ def filconv_fast(filts, band, l0, f0, fw=False):
 def filconv(band0, l0, f0, DIR, fw=False):
     '''
     Input:
-    ============
+    ======
     f0: Flux for spectrum, in fnu
     l0: Wavelength for spectrum, in AA (that matches filter response curve's.)
 

From 11778f65d57c74c307c218b907b5ef58c97b989f Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Fri, 24 Jul 2020 23:30:04 -0400
Subject: [PATCH 8/9] progress bar, just for fun

---
 gsf/function.py | 35 +++++++++++++++++++++++++++++++++++
 gsf/plot_sed.py | 12 +++++++++---
 gsf/plot_sfh.py |  8 +++++++-
 3 files changed, 51 insertions(+), 4 deletions(-)

diff --git a/gsf/function.py b/gsf/function.py
index e3659a4..c3de68f 100644
--- a/gsf/function.py
+++ b/gsf/function.py
@@ -17,6 +17,41 @@
 LW0 = [2800, 3347, 3727, 3799, 3836, 3869, 4102, 4341, 4861, 4960, 5008, 5175, 6563, 6717, 6731]
 fLW = np.zeros(len(LW0), dtype='int') # flag.
 
+
+def printProgressBar (iteration, total, prefix = '', suffix = '', decimals = 1, length = 100, fill = None, printEnd = "\r", emojis = ['🥚','🐣','🐥','🦆']):
+    '''
+    Call in a loop to create terminal progress bar
+    @params:
+        iteration   - Required  : current iteration (Int)
+        total       - Required  : total iterations (Int)
+        prefix      - Optional  : prefix string (Str)
+        suffix      - Optional  : suffix string (Str)
+        decimals    - Optional  : positive number of decimals in percent complete (Int)
+        length      - Optional  : character length of bar (Int)
+        fill        - Optional  : bar fill character (Str)
+        printEnd    - Optional  : end character (e.g. "\r", "\r\n") (Str)
+    '''
+    percent = ("{0:." + str(decimals) + "f}").format(100 * (iteration / float(total)))
+    if fill == None:
+        if float(percent) < 33:
+            fill = emojis[0]
+        elif float(percent) < 33:
+            fill = emojis[1]
+        elif float(percent) < 99:
+            fill = emojis[2]
+        else:
+            fill = emojis[3]
+
+
+    filledLength = int(length * iteration // total)
+    bar = fill * filledLength + '-' * (length - filledLength)
+    string = '(%d/%d)'%(iteration,total)
+    print(f'\r{prefix} |{bar}| {percent}% {suffix} {string}', end = printEnd)
+    # Print New Line on Complete
+    if iteration == total: 
+        print()
+
+
 def get_input():
     '''
     This returns somewhat a common default input dictionary.
diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py
index 86fd71b..16dbb01 100644
--- a/gsf/plot_sed.py
+++ b/gsf/plot_sed.py
@@ -22,13 +22,13 @@
 def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, inputs=False, nmc_rand=300, dust_model=0, DIR_TMP='./templates/', f_label=False, f_bbbox=False, verbose=False, f_silence=True, f_fill=False, f_fancyplot=False):
     '''
     Input:
-    ============
+    ======
 
     SNlim   : SN limit to show flux or up lim in SED.
     f_chind : If include non-detection in chi2 calculation, using Sawicki12.
 
     Returns:
-    ============
+    ========
     plots
 
     '''
@@ -43,6 +43,7 @@ def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=Fal
     import scipy.special as special
     import os.path
     from astropy.io import ascii
+    import time
 
     if f_silence:
         import matplotlib
@@ -567,7 +568,6 @@ def gaus(x,a,x0,sigma):
 
     # From random chain;
     alp=0.02
-    #nmc_rand = 10
     for kk in range(0,nmc_rand,1):
         nr = np.random.randint(len(samples['A0']))
         #nr = np.random.randint(int(len(samples['A0'])*0.99), len(samples['A0']))
@@ -652,6 +652,12 @@ def gaus(x,a,x0,sigma):
         UVJ[kk,2] = -2.5*np.log10(fconv[2]/fconv[3])
         UVJ[kk,3] = -2.5*np.log10(fconv[4]/fconv[3])
 
+        # Do stuff...
+        time.sleep(0.01)
+        # Update Progress Bar
+        printProgressBar(kk, nmc_rand, prefix = 'Progress:', suffix = 'Complete', length = 40)
+
+
     #
     # Plot Median SED;
     #
diff --git a/gsf/plot_sfh.py b/gsf/plot_sfh.py
index 3431f7d..bfd05ea 100644
--- a/gsf/plot_sfh.py
+++ b/gsf/plot_sfh.py
@@ -26,7 +26,6 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     ========
     Star formation history plot.
 
-
     Input:
     ======
     flim  : Lower limit for plotting an age bin.
@@ -34,6 +33,7 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
 
     '''
     import os.path
+    import time
 
     fnc  = MB.fnc #Func(ID, PA, Z, nage, dust_model=dust_model, DIR_TMP=DIR_TMP) # Set up the number of Age/ZZ
     bfnc = MB.bfnc #Basic(Z)
@@ -323,6 +323,12 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
             TC[aa, mm] /= ACs
             TL[aa, mm] /= ALs
 
+        # Do stuff...
+        time.sleep(0.01)
+        # Update Progress Bar
+        printProgressBar(mm, mmax, prefix = 'Progress:', suffix = 'Complete', length = 40)
+
+
     Avtmp  = np.percentile(Av[:],[16,50,84])
 
     #############

From cf3e93f10e1d65237efb48bd16fc95b64b4b9e33 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Sat, 25 Jul 2020 21:03:24 -0400
Subject: [PATCH 9/9] minor.

---
 gsf/fitting.py            | 17 +++++++++++------
 gsf/function.py           |  2 +-
 gsf/gsf.py                | 12 +++++++-----
 gsf/maketmp_filt.py       |  1 +
 gsf/plot_sed.py           | 23 +++++++++++------------
 gsf/plot_sfh.py           | 22 ++++++++++++----------
 gsf/posterior_flexible.py | 13 +++++++------
 7 files changed, 50 insertions(+), 40 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index 4b6afe4..861e5e6 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -46,7 +46,7 @@ def __init__(self, inputs, c=3e18, Mpc_cm=3.08568025e+24, m0set=25.0, pixelscale
     def update_input(self, inputs, c=3e18, Mpc_cm=3.08568025e+24, m0set=25.0, pixelscale=0.06, Lsun=3.839*1e33, cosmo=None):
         '''
         INPUT:
-        ==========
+        ======
         parfile: Ascii file that lists parameters for everything.
 
         Mpc_cm : cm/Mpc
@@ -292,11 +292,13 @@ def get_lines(self, LW0):
 
     def read_data(self, Cz0, Cz1, zgal, add_fir=False):
         '''
+        Input:
+        ======
         Cz0, Cz1 : Normalization coeffs for grism spectra.
         zgal     : Current redshift estimate.
 
         Note:
-        =======
+        =====
         Can be used for any SFH
 
         '''
@@ -408,12 +410,12 @@ def read_data(self, Cz0, Cz1, zgal, add_fir=False):
     def search_redshift(self, dict, xm_tmp, fm_tmp, zliml=0.01, zlimu=6.0, delzz=0.01, lines=False, prior=None, method='powell'):
         '''
         Purpose:
-        =========
+        ========
         Search redshift space to find the best redshift and probability distribution.
 
 
         Input:
-        =========
+        ======
 
         fm_tmp : a library for various templates. Should be in [ n * len(wavelength)].
         xm_tmp : a wavelength array, common for the templates above, at z=0. Should be in [len(wavelength)].
@@ -425,7 +427,7 @@ def search_redshift(self, dict, xm_tmp, fm_tmp, zliml=0.01, zlimu=6.0, delzz=0.0
         method : powell is more accurate. nelder is faster.
 
         Return:
-        =========
+        =======
 
         zspace :
         chi2s  :
@@ -1041,7 +1043,10 @@ def main(self, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, f_move
 
             ################################
             print('\nMinimizer Defined\n')
-            mini = Minimizer(class_post.lnprob, out.params, fcn_args=[dict['fy'],dict['ey'],dict['wht2'],self.f_dust], f_disp=self.f_disp, moves=emcee.moves.DEMove(sigma=1e-05, gamma0=None)) #, f_move=f_move)
+            mini = Minimizer(class_post.lnprob, out.params, fcn_args=[dict['fy'],dict['ey'],dict['wht2'],self.f_dust], f_disp=self.f_disp, \
+                moves=[(emcee.moves.DEMove(), 0.8), (emcee.moves.DESnookerMove(), 0.2),]
+                )
+                #moves=emcee.moves.DEMove(sigma=1e-05, gamma0=None))
             print('######################')
             print('### Starting emcee ###')
             print('######################')
diff --git a/gsf/function.py b/gsf/function.py
index c3de68f..2efcf7f 100644
--- a/gsf/function.py
+++ b/gsf/function.py
@@ -35,7 +35,7 @@ def printProgressBar (iteration, total, prefix = '', suffix = '', decimals = 1,
     if fill == None:
         if float(percent) < 33:
             fill = emojis[0]
-        elif float(percent) < 33:
+        elif float(percent) < 66:
             fill = emojis[1]
         elif float(percent) < 99:
             fill = emojis[2]
diff --git a/gsf/gsf.py b/gsf/gsf.py
index dc78f24..b2828bf 100644
--- a/gsf/gsf.py
+++ b/gsf/gsf.py
@@ -19,7 +19,7 @@
 def run_gsf_template(inputs, fplt=0):
     '''
     Purpose:
-    ==========
+    ========
     This is only for 0 and 1, to get templates.
     Not for fitting, nor plotting.
 
@@ -130,7 +130,7 @@ def run_gsf_all(parfile, fplt, cornerplot=True):
         #
         MB.zprev = MB.zgal #zrecom # redshift from previous run
 
-        flag_suc = MB.main(0, cornerplot=cornerplot)
+        flag_suc = MB.main(cornerplot=cornerplot)
 
         while (flag_suc and flag_suc!=-1):
 
@@ -142,7 +142,7 @@ def run_gsf_all(parfile, fplt, cornerplot=True):
             print('Going into another trial with updated templates and redshift.')
             print('\n\n')
 
-            flag_suc = MB.main(1, cornerplot=cornerplot)
+            flag_suc = MB.main(cornerplot=cornerplot)
 
         # Total calculation time
         stop = timeit.default_timer()
@@ -152,10 +152,12 @@ def run_gsf_all(parfile, fplt, cornerplot=True):
     if fplt <= 3 and flag_suc != -1:
         from .plot_sfh import plot_sfh
         from .plot_sed import plot_sed
-        plot_sfh(MB, f_comp=MB.ftaucomp, fil_path=MB.DIR_FILT,
+        
+        plot_sfh(MB, f_comp=MB.ftaucomp, fil_path=MB.DIR_FILT, mmax=100,
         inputs=MB.inputs, dust_model=MB.dust_model, DIR_TMP=MB.DIR_TMP, f_SFMS=True)
+
         plot_sed(MB, fil_path=MB.DIR_FILT,
-        figpdf=False, save_sed=True, inputs=MB.inputs, nmc_rand=1000,
+        figpdf=False, save_sed=True, inputs=MB.inputs, mmax=30,
         dust_model=MB.dust_model, DIR_TMP=MB.DIR_TMP, f_label=True)
 
 
diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py
index 575e20a..e4b3e8b 100644
--- a/gsf/maketmp_filt.py
+++ b/gsf/maketmp_filt.py
@@ -585,6 +585,7 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=20000., ncolbb=10000):
     for ii in range(len(ltmpbb[0,:])):
         if SFILT[ii] in SKIPFILT:# data point to be skiped;
             fw.write('%d %.5f %.5e %.5e\n'%(ii+ncolbb, ltmpbb[0,ii], 0.0, fbb[ii]))
+            #fw.write('%d %.5f %.5e %.5e\n'%(ii+ncolbb, ltmpbb[0,ii], 0.0, 1000))
         elif  ebb[ii]>ebblim:
             fw.write('%d %.5f 0 1000\n'%(ii+ncolbb, ltmpbb[0,ii]))
         else:
diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py
index 16dbb01..e500a98 100644
--- a/gsf/plot_sed.py
+++ b/gsf/plot_sed.py
@@ -19,7 +19,7 @@
 col = ['violet', 'indigo', 'b', 'lightblue', 'lightgreen', 'g', 'orange', 'coral', 'r', 'darkred']#, 'k']
 #col = ['darkred', 'r', 'coral','orange','g','lightgreen', 'lightblue', 'b','indigo','violet','k']
 
-def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, inputs=False, nmc_rand=300, dust_model=0, DIR_TMP='./templates/', f_label=False, f_bbbox=False, verbose=False, f_silence=True, f_fill=False, f_fancyplot=False):
+def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, inputs=False, 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):
     '''
     Input:
     ======
@@ -550,8 +550,8 @@ def gaus(x,a,x0,sigma):
     samples  = res
 
     # Saved template;
-    ytmp = np.zeros((nmc_rand,len(ysum)), dtype='float64')
-    ytmp_each = np.zeros((nmc_rand,len(ysum),len(age)), dtype='float64')
+    ytmp = np.zeros((mmax,len(ysum)), dtype='float64')
+    ytmp_each = np.zeros((mmax,len(ysum),len(age)), dtype='float64')
 
     ytmpmax = np.zeros(len(ysum), dtype='float64')
     ytmpmin = np.zeros(len(ysum), dtype='float64')
@@ -559,18 +559,17 @@ def gaus(x,a,x0,sigma):
     # MUV;
     DL      = MB.cosmo.luminosity_distance(zbes).value * Mpc_cm # Luminositydistance in cm
     DL10    = Mpc_cm/1e6 * 10 # 10pc in cm
-    Fuv     = np.zeros(nmc_rand, dtype='float64') # For Muv
-    Fuv28   = np.zeros(nmc_rand, dtype='float64') # For Fuv(1500-2800)
-    Lir     = np.zeros(nmc_rand, dtype='float64') # For L(8-1000um)
-    UVJ     = np.zeros((nmc_rand,4), dtype='float64') # For UVJ color;
+    Fuv     = np.zeros(mmax, dtype='float64') # For Muv
+    Fuv28   = np.zeros(mmax, dtype='float64') # For Fuv(1500-2800)
+    Lir     = np.zeros(mmax, dtype='float64') # For L(8-1000um)
+    UVJ     = np.zeros((mmax,4), dtype='float64') # For UVJ color;
 
     Cmznu   = 10**((48.6+m0set)/(-2.5)) # Conversion from m0_25 to fnu
 
     # From random chain;
     alp=0.02
-    for kk in range(0,nmc_rand,1):
+    for kk in range(0,mmax,1):
         nr = np.random.randint(len(samples['A0']))
-        #nr = np.random.randint(int(len(samples['A0'])*0.99), len(samples['A0']))
         try:
             Av_tmp = samples['Av'][nr]
         except:
@@ -655,7 +654,7 @@ def gaus(x,a,x0,sigma):
         # Do stuff...
         time.sleep(0.01)
         # Update Progress Bar
-        printProgressBar(kk, nmc_rand, prefix = 'Progress:', suffix = 'Complete', length = 40)
+        printProgressBar(kk, mmax, prefix = 'Progress:', suffix = 'Complete', length = 40)
 
 
     #
@@ -726,13 +725,13 @@ def func_tmp(xint,eobs,fmodel):
             ndim_eff -= 1
     '''
 
-    con_up = (ey>0) & (fy/ey<=SNlim)
+    con_up = (fy==0) & (ey>0) & (fy/ey<=SNlim)
     chi_nd = 0.0
     if f_chind:
         # Chi2 for non detection;
         for nn in range(len(ey[con_up])):
             #result  = integrate.quad(lambda xint: func_tmp(xint, ey[con_up][nn]/SNlim, ysump[con_up][nn]), -ey[con_up][nn]/SNlim, ey[con_up][nn]/SNlim, limit=100)
-            result  = integrate.quad(lambda xint: func_tmp(xint, ey[con_up][nn], ysump[con_up][nn]), -ey[con_up][nn], ey[con_up][nn], limit=100)
+            result  = integrate.quad(lambda xint: func_tmp(xint, ey[con_up][nn], ysump[con_up][nn]), -ey[con_up][nn]*100, ey[con_up][nn], limit=100)
             chi_nd += np.log(result[0])
 
     # Number of degree;
diff --git a/gsf/plot_sfh.py b/gsf/plot_sfh.py
index bfd05ea..8150713 100644
--- a/gsf/plot_sfh.py
+++ b/gsf/plot_sfh.py
@@ -20,7 +20,7 @@
 
 lcb   = '#4682b4' # line color, blue
 
-def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, fil_path = './FILT/', inputs=None, dust_model=0, DIR_TMP='./templates/',f_SFMS=False, verbose=False):
+def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmin=0.08, Txmax=4, lmmin=9.5, fil_path = './FILT/', inputs=None, dust_model=0, DIR_TMP='./templates/',f_SFMS=False, verbose=False):
     '''
     Purpose:
     ========
@@ -47,9 +47,11 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
         age  = MB.age
     nage = MB.nage #np.arange(0,len(age),1)
     tau0 = MB.tau0 #[0.1,0.2,0.3]
-
     age = np.asarray(age)
 
+    if Txmin > np.min(age):
+        Txmin = np.min(age) * 0.8
+
     ################
     # RF colors.
     home = os.path.expanduser('~')
@@ -482,13 +484,13 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     if f_comp == 1:
         lsfru = np.max([lsfru, np.log10(np.max(sfr_exp*C_exp))])
 
-    ax1.set_xlim(0.008, Txmax)
+    ax1.set_xlim(Txmin, Txmax)
     ax1.set_ylim(lsfrl, lsfru)
     ax1.set_xscale('log')
 
     ax2.set_ylabel('$\log M_*/M_\odot$', fontsize=12)
 
-    ax2.set_xlim(0.008, Txmax)
+    ax2.set_xlim(Txmin, Txmax)
     ax2.set_ylim(y2min, y2max)
     ax2.set_xscale('log')
 
@@ -591,7 +593,7 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     #ax3.set_xlabel('$t$ (Gyr)', fontsize=12)
     #ax3.set_ylabel('$\log Z_*/Z_\odot$', fontsize=12)
     y3min, y3max = np.min(Z), np.max(Z)
-    #ax3.set_xlim(0.008, Txmax)
+    #ax3.set_xlim(Txmin, Txmax)
     #ax3.set_ylim(y3min, y3max)
     #ax3.set_xscale('log')
     #ax3.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
@@ -628,7 +630,7 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
             Tzz[zz] = 0.01
 
     #ax3t.set_xscale('log')
-    #ax3t.set_xlim(0.008, Txmax)
+    #ax3t.set_xlim(Txmin, Txmax)
 
     ax1.set_xlabel('$t_\mathrm{lookback}$/Gyr', fontsize=12)
     ax2.set_xlabel('$t_\mathrm{lookback}$/Gyr', fontsize=12)
@@ -636,13 +638,13 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
     ax4.set_ylabel('$\log Z_*/Z_\odot$', fontsize=12)
 
     ax1t.set_xscale('log')
-    ax1t.set_xlim(0.008, Txmax)
+    ax1t.set_xlim(Txmin, Txmax)
     ax2t.set_xscale('log')
-    ax2t.set_xlim(0.008, Txmax)
+    ax2t.set_xlim(Txmin, Txmax)
     ax4t.set_xscale('log')
-    ax4t.set_xlim(0.008, Txmax)
+    ax4t.set_xlim(Txmin, Txmax)
 
-    ax4.set_xlim(0.008, Txmax)
+    ax4.set_xlim(Txmin, Txmax)
     ax4.set_ylim(y3min-0.05, y3max)
     ax4.set_xscale('log')
 
diff --git a/gsf/posterior_flexible.py b/gsf/posterior_flexible.py
index 6252f0e..66f4cb6 100644
--- a/gsf/posterior_flexible.py
+++ b/gsf/posterior_flexible.py
@@ -18,12 +18,12 @@ def __init__(self, mainbody):
     def residual(self, pars, fy, ey, wht, f_fir=False, out=False):
         '''
         Input:
-        ==========
+        ======
         out   : model as second output. For lnprob func.
         f_fir : Bool. If dust component is on or off.
 
         Returns:
-        ==========
+        ========
         residual of model and data.
 
         '''
@@ -76,7 +76,7 @@ def lnprob(self, pars, fy, ey, wht, f_fir, f_chind=True, SNlim=1.0):
         '''
 
         Returns:
-        =========
+        ========
         log posterior
 
         '''
@@ -89,16 +89,17 @@ def lnprob(self, pars, fy, ey, wht, f_fir, f_chind=True, SNlim=1.0):
 
         resid, model = self.residual(pars, fy, ey, wht, f_fir, out=True)
         sig     = np.sqrt(1./wht+f**2*model**2)
+        chi_nd  = 0.0
 
-        chi_nd = 0
-        con_up = (fy==0) & (ey>0)
         if f_chind:
-            # This does not improve, but costs time;
+            con_up = (fy==0) & (fy/ey<=SNlim) & (ey>0)
+            # This may be a bit cost of time;
             for nn in range(len(ey[con_up])):
                 result  = integrate.quad(lambda xint: self.func_tmp(xint, ey[con_up][nn]/SNlim, model[con_up][nn]), -ey[con_up][nn], ey[con_up][nn], limit=100)
                 chi_nd += np.log(result[0])
 
             con_res = (model>=0) & (wht>0) & (fy>0) # Instead of model>0, model>=0 is for Lyman limit where flux=0.
+            #print(len(fy[con_up]),len(fy[con_res]))
             lnlike  = -0.5 * (np.sum(resid[con_res]**2 + np.log(2 * 3.14 * sig[con_res]**2)) - 2 * chi_nd)
 
         else: