diff --git a/CONAN3/VERSION.dat b/CONAN3/VERSION.dat
index a4f52a5..0444f32 100755
--- a/CONAN3/VERSION.dat
+++ b/CONAN3/VERSION.dat
@@ -1 +1 @@
-3.2.0
\ No newline at end of file
+3.2.1
\ No newline at end of file
diff --git a/CONAN3/_classes.py b/CONAN3/_classes.py
index 07afe38..ce91e75 100755
--- a/CONAN3/_classes.py
+++ b/CONAN3/_classes.py
@@ -16,6 +16,7 @@
 from .utils import phase_fold, supersampling, convert_LD, get_transit_time, bin_data_with_gaps
 from copy import deepcopy
 from scipy.interpolate import LSQUnivariateSpline,LSQBivariateSpline
+from uncertainties import ufloat
 
 __all__ = ["load_lightcurves", "load_rvs", "fit_setup", "load_result", "__default_backend__"]
 
@@ -900,7 +901,7 @@ def __init__(self, file_list=None, data_filepath=None, filters=None, wl=None, np
                 print(f"writing ones to the missing columns of file: {f}")
                 new_cols = np.ones((nrow,9-ncol))
                 fdata = np.hstack((fdata,new_cols))
-                np.savetxt(self._fpath+f,fdata,fmt='%.8f')
+                # np.savetxt(self._fpath+f,fdata,fmt='%.8f')   # no need to replace the file since its stored in lc_obj
             #remove nan rows from fdata and print number of removed rows
             n_nan = np.sum(np.isnan(fdata).any(axis=1))
             if n_nan > 0: print(f"removed {n_nan} row(s) with NaN values from file: {f}")            
@@ -1270,7 +1271,7 @@ def get_decorr(self, T_0=None, Period=None, rho_star=None, Duration=None, D_occ=
         return self._decorr_result
     
     
-    def clip_outliers(self, lc_list="all", clip=5, width=15, show_plot=False, verbose=True):
+    def clip_outliers(self, lc_list="all", clip=5, width=15, select_column=["col1"], niter=1, show_plot=False, verbose=True):
 
         """
         Remove outliers using a running median method. Points > clip*M.A.D are removed
@@ -1280,16 +1281,16 @@ def clip_outliers(self, lc_list="all", clip=5, width=15, show_plot=False, verbos
         ------------
         lc_list: list of string, None, 'all';
             list of lightcurve filenames on which perform outlier clipping. Default is 'all' which clips all lightcurves in the object.
-        
         clip: int/float,list;
-            cut off value above the median. Default is 5
-
+            cut off value above the median deviation. Default is 5. If list, must be same length as lc_list.
         width: int,list;
-            Number of points in window to use when computing the running median. Must be odd. Default is 15
-
+            Number of points in window to use when computing the running median. Must be odd. Default is 15. If list, must be same length as lc_list.
+        select_column: list of str;
+            list of column names on which to perform clipping. Default is only ["col1"] which is the flux column.
+        niter: int;
+            Number of iterations to perform clipping. Default is 1
         show_plot: bool;
             set True to plot the data and show clipped points.
-        
         verbose: bool;
             Prints number of points that have been cut. Default is True
 
@@ -1302,6 +1303,12 @@ def clip_outliers(self, lc_list="all", clip=5, width=15, show_plot=False, verbos
             print("lc_list is None: No lightcurve to clip outliers.")
             return None
         
+        if isinstance(select_column,str): select_column = [select_column]
+        if isinstance(select_column, list):
+            for col in select_column: 
+                    assert col in ["col1","col3","col4","col5","col6","col7","col8"],\
+                            f'clip_outliers(): elements of select_column must be in ["col1","col3","col4","col5","col6","col7","col8"] but "{col}" given.'
+        
         if isinstance(lc_list, str) and (lc_list != 'all'): lc_list = [lc_list]
         if lc_list == "all": lc_list = self._names
 
@@ -1333,23 +1340,31 @@ def clip_outliers(self, lc_list="all", clip=5, width=15, show_plot=False, verbos
             
             self._clipped_data.config[self._names.index(file)] = f"W{width[i]}C{clip[i]}"
 
-            thisLCdata = self._input_lc[file]
-            _,_,clpd_mask = outlier_clipping(x=thisLCdata["col0"],y=thisLCdata["col1"],clip=clip[i],width=width[i],
-                                                verbose=False, return_clipped_indices=True)   #returns mask of the clipped points
-            ok = ~clpd_mask     #invert mask to get indices of points that are not clipped
-            if verbose and (not show_plot): print(f'\n{file}: Rejected {sum(~ok)} pts > {clip[i]:0.1f}MAD from the median')
+            thisLCdata = deepcopy(self._input_lc[file])
+            ok      = np.ones(len(thisLCdata["col0"]), dtype=bool)  #initialize mask to all True, used to store indices of points that are not clipped
+            ok_iter = np.ones(len(thisLCdata["col0"]), dtype=bool)  #initialize mask to all True, ok points for each iteration
+
+            for col in select_column:
+                for _ in range(niter):
+                    thisLCdata  = {k:v[ok_iter] for k,v in thisLCdata.items()}   #remove clipped points from previous iteration
+                    _,_,clpd_mask = outlier_clipping(x=thisLCdata["col0"],y=thisLCdata[col],clip=clip[i],width=width[i],
+                                                        verbose=False, return_clipped_indices=True)   #returns mask of the clipped points
+                    ok_iter = ~clpd_mask     #invert mask to get indices of points that are not clipped
+                    ok[ok] &= ok_iter        #update points in ok mask with the new iteration clipping
+            if verbose and (not show_plot): print(f'\n{file}: Rejected {sum(~ok)}pts > {clip[i]:0.1f}MAD from the median of columns {select_column}')
 
             if show_plot:
-                ax[i].set_title(f'{file}: Rejected {sum(~ok)} pts > {clip[i]:0.1f}MAD')
-                ax[i].plot(thisLCdata["col0"][ok],  thisLCdata["col1"][ok], '.C0', ms=5)
-                ax[i].plot(thisLCdata["col0"][~ok], thisLCdata["col1"][~ok], '.r', ms=5)
+                ax[i].set_title(f'{file}: Rejected {sum(~ok)}pts>{clip[i]:0.1f}MAD')
+                ax[i].plot(self._input_lc[file]["col0"][ok],  self._input_lc[file]["col1"][ok], '.C0', ms=5)
+                ax[i].plot(self._input_lc[file]["col0"][~ok], self._input_lc[file]["col1"][~ok], '.r', ms=5)
+
             
             #replace all columns of input file with the clipped data
-            self._input_lc[file] = {k:v[ok] for k,v in thisLCdata.items()}
+            self._input_lc[file] = {k:v[ok] for k,v in self._input_lc[file].items()}
 
             #recompute rms estimate and multiplicative jitter
-            self._rms_estimate[self._names.index(file)]  = np.std(np.diff(thisLCdata["col1"][ok]))/np.sqrt(2)
-            self._jitt_estimate[self._names.index(file)] = np.sqrt(self._rms_estimate[self._names.index(file)]**2 - np.mean(thisLCdata["col2"][ok]**2))
+            self._rms_estimate[self._names.index(file)]  = np.std(np.diff(self._input_lc[file]["col1"]))/np.sqrt(2)
+            self._jitt_estimate[self._names.index(file)] = np.sqrt(self._rms_estimate[self._names.index(file)]**2 - np.mean(self._input_lc[file]["col2"]**2))
             if np.isnan(self._jitt_estimate[self._names.index(file)]): self._jitt_estimate[self._names.index(file)] = 1e-20
         
         self._clipped_data.flag = True # SimpleNamespace(flag=True, width=width, clip=clip, lc_list=lc_list, config=conf)
@@ -1635,7 +1650,7 @@ def add_GP(self ,lc_list=None, par=["col0"], kernel=["mat32"], operation=[""],
         celerite_allowed = dict(kernels = ["real","mat32","sho"], columns= ["col0","col3","col4","col5","col6","col7","col8"])
 
         self._GP_dict  = {}
-        self._sameLCgp  = SimpleNamespace(flag = False, first_index =None)
+        self._sameLCgp  = SimpleNamespace(flag = False, first_index =None) #flag to indicate if same GP is to be used for all lcs
 
         if lc_list is None or lc_list == []:
             if self._nphot>0:
@@ -1744,10 +1759,10 @@ def add_GP(self ,lc_list=None, par=["col0"], kernel=["mat32"], operation=[""],
                                                                         user_data = [this_kern, this_par])
                     elif isinstance(v, tuple):
                         if len(v)==2:
-                            steps = 0 if (self._sameLCgp.flag and i!=0) else 0.1*v[1]   #if sameRVgp is set, only first pars will jump and be used for all rvs
+                            steps = 0 if (self._sameLCgp.flag and i!=0) else 0.1*v[1]   #if sameLCgp is set, only first pars will jump and be used for all rvs
                             self._GP_dict[lc][p+str(j)] = _param_obj(to_fit="y", start_value=v[0],step_size=steps,prior="p", 
                                                                         prior_mean=v[0], prior_width_lo=v[1], prior_width_hi=v[1], 
-                                                                        bounds_lo=v[0]-10*v[1], bounds_hi=v[0]+10*v[1],
+                                                                        bounds_lo=v[0]-10*v[1], bounds_hi=v[0]+10*v[1],     #10sigma cutoff
                                                                         user_data=[this_kern, this_par])
                         elif len(v)==3:
                             steps = 0 if (self._sameLCgp.flag and i!=0) else min(0.001,0.001*np.ptp(v))
@@ -3333,6 +3348,30 @@ class load_result:
         Examples:
         ---------
         >>> result = CONAN3.load_result(folder="output")
+        # different plots from the result object
+        >>> fig    = result.plot_corner()                       # corner plot
+        >>> fig    = result.plot_burnin_chains()                # burn-in chains
+        >>> fig    = result.plot_chains()                       # posterior chains
+        >>> fig    = result.lc.plot_bestfit(detrend=True)        # model of the light curves
+        >>> fig    = result.rv.plot_bestfit(detrend=True)        # model of the RV curves
+
+        #get the best-fit parameters
+        >>> med_pars = result.params.median()      # median values of the fitted parameters
+        >>> stdev    = result.params.stdev()       # standard deviation of the fitted parameters
+        >>> pars_dict= result.get_all_params_dict(stat="med")              # get all parameters (fitted, derived, and fixed) as a dictionary
+        
+        # load files
+        >>> out_lc  = result.lc.out_data()            # output data of the light curves i.e *_lcout.dat files
+        >>> out_rv  = result.rv.out_data()            # output data of the RV curves i.e *_rvout.dat files
+        >>> in_lc   = result.lc.in_data()             # input light curves
+        >>> in_rv   = result.rv.in_data()             # input RV data 
+
+        # evaluate model (lc or rv) at user-defined times
+        >>> t = np.linspace(0,1,1000)
+        >>> model  = result.lc.evaluate(file="lc1.dat", time=t, params= result.params.median, return_std=True)                # model of the light curve "lc1.dat" at user time t
+        >>> lc_mod = model.planet_model      # model of the planet
+        >>> comps  = model.components        # for multiplanet fit, this will be a dict with lc_mod for each planet. i.e. comps["pl_1"] for planet 1 
+        >>> sigma_low, sigma_hi = model.sigma_low, model.sigma_hi    # lower and upper 1-sigma model uncertainties that can be plotted along with lc_mod
     """
 
     def __init__(self, folder="output",chain_file = "chains_dict.pkl", burnin_chain_file="burnin_chains_dict.pkl",verbose=True):
@@ -3380,6 +3419,7 @@ def __init__(self, folder="output",chain_file = "chains_dict.pkl", burnin_chain_
             self.params     = SimpleNamespace(  names   = list(self._par_names),
                                                 median  = self._stat_vals["med"],
                                                 max     = self._stat_vals["max"],
+                                                stdev   = self._stat_vals["stdev"] if "stdev" in self._stat_vals.keys() else np.zeros_like(self._stat_vals["med"]),
                                                 bestfit = self._stat_vals["bf"],
                                                 T0      = self._stat_vals["T0"],
                                                 P       = self._stat_vals["P"],
@@ -3387,10 +3427,11 @@ def __init__(self, folder="output",chain_file = "chains_dict.pkl", burnin_chain_
             assert len(self.params.median)==len(self.params.names), "load_result(): number of parameter names and values do not match."
         except:
             self.params     = SimpleNamespace(  names   = list(self._par_names),
-                                                median  = np.median(self.flat_posterior,axis=0))
+                                                median  = np.median(self.flat_posterior,axis=0),
+                                                stdev   = np.std(self.flat_posterior,axis=0))
             
         self.evidence   = self._stat_vals["evidence"] if hasattr(self,"_stat_vals") and "evidence" in self._stat_vals.keys() else None
-        self.params_dict  = {k:v for k,v in zip(self.params.names, self.params.median)}
+        self.params_dict  = {k:ufloat(v,e) for k,v,e in zip(self.params.names, self.params.median,self.params.stdev)}
 
         if hasattr(self,"_stat_vals"):     #summary statistics are available only if fit completed
             # evaluate model of each lc at a smooth time grid
@@ -3421,7 +3462,7 @@ def __init__(self, folder="output",chain_file = "chains_dict.pkl", burnin_chain_
                 self._rv_smooth_time_mod[rv] = SimpleNamespace()
                 tmin, tmax = input_rvs[rv]["col0"].min(), input_rvs[rv]["col0"].max()
                 t_sm  = np.linspace(tmin,tmax,max(2000, len(input_rvs[rv]["col0"])))
-                gam = self.params_dict[f"rv{i+1}_gamma"]
+                gam = self.params_dict[f"rv{i+1}_gamma"].n
                 self._rv_smooth_time_mod[rv].time    = t_sm
                 self._rv_smooth_time_mod[rv].model   = self._evaluate_rv(file=rv, time=self._rv_smooth_time_mod[rv].time).planet_model + gam
                 self._rv_smooth_time_mod[rv].gamma   = gam
@@ -3816,6 +3857,58 @@ def _plot_bestfit_rv(self, plot_cols=(0,1,2), detrend=False, col_labels=None, nr
 
             if return_fig: return fig
 
+    def get_all_params_dict(self, stat="med",uncertainty="1sigma", return_type="ufloat"):
+        """
+            Get all parameters(jumping,derived,fixed) from the result_**.dat and load in a dictionary with uncertainties.
+
+            Parameters:
+            -----------
+            stat : str;
+                summary statistic to load for model calculation. Must be one of ["med","max","bf"] for median, maximum and bestfit respectively.
+                Default is "med" to load the 'result_med.dat' file.
+            uncertainty : str;
+                uncertainty to load from file. Must be one of ["1sigma","3sigma"] for 1sigma or 3sigma uncertainties.
+                Default is "1sigma".
+            return_type : str;
+                return type of the values. Must be one of ["ufloat","array"] to return each parameter as ufloat(val,+/-sigma) or array of [val,lo_sigma,hi_sigma] .
+                Default is "ufloat".
+            Returns:
+            --------
+            results_dict : dict;
+                dictionary of all parameters with uncertainties for jumping and derived parameters
+
+        """
+        results_dict = {}
+        assert uncertainty in  ["1sigma","3sigma"], "get_all_params_dict(): uncertainty must be one of ['1sigma','3sigma']"
+        assert stat in ["med","max","bf"], "get_all_params_dict(): stat must be one of ['med','max','bf']"
+        assert return_type in ["ufloat","array"], "get_all_params_dict(): return_type must be one of ['ufloat','array']"
+        
+        print(f"Loading file: {self._folder}/results_{stat}.dat ...\n")
+
+        with open(f"{self._folder}/results_{stat}.dat", 'r') as file:
+            for line in file:
+                # Ignore lines that start with '#'
+                if not line.startswith('#'):
+                    # Split the line into words
+                    words = line.split()
+                    # Use the first word as the key
+                    val = float(words[1])
+                    
+                    if len(words) > 2: # jumping and derived parameters
+                        if uncertainty == "1sigma":
+                            lo, up = abs(float(words[2])), float(words[3])
+                            sigma  = np.median([lo, up]) if return_type == "ufloat" else [lo,up]
+                        else:
+                            lo, up = abs(float(words[4])), float(words[5])
+                            sigma  = np.median([lo, up]) if return_type == "ufloat" else [lo,up]
+                        
+                        results_dict[words[0]] = ufloat(val, sigma) if return_type=="ufloat" else np.array([val]+sigma)
+                    
+                    else:
+                        results_dict[words[0]] = val
+                        
+        return results_dict
+    
     def _load_result_array(self, data=["lc","rv"],verbose=True):
         """
             Load result array from CONAN3 fit allowing for customised plots.
diff --git a/CONAN3/conf.py b/CONAN3/conf.py
index 5b46aa5..a7c0e1d 100755
--- a/CONAN3/conf.py
+++ b/CONAN3/conf.py
@@ -1,5 +1,4 @@
 from ._classes import load_lightcurves, load_rvs, fit_setup,_print_output
-from .fit_data import run_fit
 from .utils import ecc_om_par
 import numpy as np 
 
@@ -27,6 +26,7 @@ def fit_configfile(config_file = "input_config.dat", out_folder = "output",
         verbose: bool;
             show print statements
     """
+    from .fit_data import run_fit
 
     lc_obj, rv_obj, fit_obj = load_configfile(config_file, init_decorr=init_decorr, verbose=verbose)
     result = run_fit(lc_obj, rv_obj, fit_obj,out_folder=out_folder,
@@ -64,6 +64,9 @@ def create_configfile(lc_obj=None, rv_obj=None, fit_obj=None, filename="input_co
         
         fit_obj : object;
             Instance of CONAN.fit_setup() object and its attributes.
+
+        filename : str;
+            name of the configuration file to be saved.
     """
     if lc_obj is None:
         lc_obj = load_lightcurves()
diff --git a/CONAN3/fit_data.py b/CONAN3/fit_data.py
index 1e7b308..69a6431 100755
--- a/CONAN3/fit_data.py
+++ b/CONAN3/fit_data.py
@@ -24,6 +24,7 @@
 from celerite import GP as cGP
 from copy import deepcopy
 from .utils import gp_params_convert
+from .conf import create_configfile
 from scipy.stats import norm, uniform, lognorm, loguniform,truncnorm
 
 
@@ -112,6 +113,9 @@ def run_fit(lc_obj=None, rv_obj=None, fit_obj=None, statistic = "median", out_fo
     if not os.path.exists(out_folder):
         print(f"Creating output folder...{out_folder}")
         os.mkdir(out_folder)
+        # print(f"Saving configuration to file: {out_folder}/config_save.dat")
+        create_configfile(lc_obj, rv_obj, fit_obj, f"{out_folder}/config_save.dat")  #create config file
+
 
     if os.path.exists(f'{out_folder}/chains_dict.pkl'):
         if not rerun_result:
@@ -1453,10 +1457,10 @@ def run_fit(lc_obj=None, rv_obj=None, fit_obj=None, statistic = "median", out_fo
     bpfull[[ijnames][0]] = bp
 
     try:
-        medvals,maxvals = mcmc_outputs(posterior,jnames, ijnames, njnames, nijnames, bpfull, uwl, Rs_in, Ms_in, Rs_PDF, Ms_PDF, 
+        medvals,maxvals,sigma1s = mcmc_outputs(posterior,jnames, ijnames, njnames, nijnames, bpfull, uwl, Rs_in, Ms_in, Rs_PDF, Ms_PDF, 
                                         nfilt, filnames, howstellar, extinpars, RVunit, extind_PDF,npl,out_folder)
     except:
-        medvals,maxvals = mcmc_outputs(posterior,jnames, ijnames, njnames, nijnames, bpfull, uwl, Rs_in, Ms_in, Rs_PDF, Ms_PDF, 
+        medvals,maxvals,sigma1s = mcmc_outputs(posterior,jnames, ijnames, njnames, nijnames, bpfull, uwl, Rs_in, Ms_in, Rs_PDF, Ms_PDF, 
                                         nfilt, filnames, howstellar, extinpars, RVunit, extind_PDF,npl,out_folder)
 
     npar=len(jnames)
@@ -1468,6 +1472,8 @@ def run_fit(lc_obj=None, rv_obj=None, fit_obj=None, statistic = "median", out_fo
     medp[[ijnames][0]]=medvals
     maxp=initial
     maxp[[ijnames][0]]=maxvals
+    stdev = np.zeros_like(initial)
+    stdev[[ijnames][0]] = sigma1s
 
     #============================== PLOTTING ===========================================
     print('\nPlotting output figures')
@@ -1481,7 +1487,7 @@ def run_fit(lc_obj=None, rv_obj=None, fit_obj=None, statistic = "median", out_fo
 
     #save summary_stats and as a hidden files. 
     #can be used to run logprob_multi() to generate out_full.dat files for median posterior, max posterior and best fit values
-    stat_vals = dict(med = medp[jumping], max = maxp[jumping], bf  = bpfull[jumping], 
+    stat_vals = dict(med = medp[jumping], max = maxp[jumping], bf  = bpfull[jumping], stdev=stdev[jumping],
                         T0 = T0_post,  P = p_post, dur = Dur_post, evidence=evidence)
     pickle.dump(stat_vals, open(out_folder+"/.stat_vals.pkl","wb"))
 
diff --git a/CONAN3/get_files.py b/CONAN3/get_files.py
index 7345d82..dc97d5c 100755
--- a/CONAN3/get_files.py
+++ b/CONAN3/get_files.py
@@ -100,6 +100,7 @@ def download(self,sectors=None,author=None,exptime=None, select_flux="pdcsap_flu
             print(f"downloaded lightcurve for sector {s}")
 
         if hasattr(self,"_ok_mask"): del self._ok_mask
+        self.data_splitted = False
 
     def discard_ramp(self,length=0.25, gap_size=1, start=True, end=True):
         """
@@ -156,7 +157,75 @@ def scatter(self):
         """
         assert self.lc != {}, "No light curves downloaded yet. Run `download()` first."
         for s in self.sectors:
-            self.lc[s].scatter()
+            ax = self.lc[s].scatter()
+            ax.set_title(f"sector {s}")
+
+    def split_data(self, gap_size=None, split_times=None, show_plot=True):
+        """
+        Split the light curves into sections based on (1)gaps in the data, (2)given time intervals (3)given split time(s).
+
+        Parameters
+        ----------
+        gap_size : float
+            minimum gap size at which data will be splitted. same unit as the time axis. Default is None.
+        split_times : float,list
+            Time to split the data. Default is None.
+        show_plot : bool
+            Show plot of the split times. Default is True.
+        """
+
+        assert self.lc != {}, "No light curves downloaded yet. Run `download()` first."
+        if isinstance(split_times,(int,float)): split_times = [split_times] #convert to list
+        if self.data_splitted: 
+            print("data has already being split. Run `.download()` again to restart.")
+            return
+
+        remove_sectors = []
+        curr_sector = self.sectors.copy()
+        for s in curr_sector:
+            t = self.lc[s].time.value
+    
+            if gap_size is not None:
+                self._find_gaps(t,gap_size)
+                if show_plot:
+                    ax = self.lc[s].scatter()
+                    [ax.axvline(t[i], color='r', linestyle='--') for i in self.gap_ind[1:-1]]
+                nsplit = len(self.gap_ind)-1
+                for i in range(len(self.gap_ind)-1):
+                    self.lc[f"{s}_{i+1}"] = self.lc[s][(t >= t[self.gap_ind[i]]) & (t < t[self.gap_ind[i+1]])]
+            
+            if split_times is not None:
+                split_times = np.array([0] +split_times +[t[-1]])
+                if show_plot:
+                    ax = self.lc[s].scatter()
+                    [ax.axvline(ti, color='r', linestyle='--') for ti in split_times[1:-1]]
+                nsplit = len(split_times)-1
+                for i in range(len(split_times)-1):
+                    self.lc[f"{s}_{i+1}"] = self.lc[s][(t >= split_times[i]) & (t < split_times[i+1])]
+            
+            if split_times is None and gap_size is None:
+                raise ValueError("either gap_size or split_time must be given")
+            
+            remove_sectors.append(s)
+            self.sectors += [f"{s}_{i+1}" for i in range(nsplit)]
+        
+        print(f"Sector {s} data splitted into {nsplit} chunks: {list(self.sectors)}")
+        #remove sectors that have been split        
+        _ = [self.sectors.pop(self.sectors.index(s)) for s in remove_sectors]
+        self.data_splitted = True
+        
+
+    def _find_gaps(self,t,gap_size):
+        """
+        find gaps in the data
+        """
+        gap = np.diff(t) 
+        gap = np.insert(gap,len(gap),0) #insert diff of 0 at the beginning
+        self.gap_ind = np.where(gap > gap_size)[0]
+        self.gap_ind = np.array([0]+list(self.gap_ind)+[len(t)-1])
+        # self.gap_ind = np.append(0,self.gap_ind)
+        # self.gap_ind = np.append(self.gap_ind, len(t)-1)
+
     
     def save_CONAN_lcfile(self,bjd_ref = 2450000, folder="data", out_filename=None):
         """
@@ -227,7 +296,8 @@ def search(self,filters=None):
         print(pd.DataFrame(data, columns=["obj_id_catname","file_key", "pi_name", "date_mjd_start",
                                             "obs_total_exptime", "data_arch_rev",'status_published']))
 
-    def download(self, file_keys="all", aperture="DEFAULT", decontaminate=True, reject_highpoints=True, bg_MAD_reject=3,verbose=True):
+    def download(self, file_keys="all", aperture="DEFAULT", decontaminate=True, reject_highpoints=True, bg_MAD_reject=3,
+                 outlier_clip=4, outlier_window_width=11,outlier_clip_niter=1, verbose=True):
         """   
         Download CHEOPS light curves from the cheops database using the dace_query package
 
@@ -235,18 +305,20 @@ def download(self, file_keys="all", aperture="DEFAULT", decontaminate=True, reje
         ----------
         file_keys : list of str
             List of file keys to download. if "all" is given, all file_keys shown in `.search()` result will be downloaded.
-
         aperture : str
             Aperture to use for the light curve. Default is "DEFAULT". The options are: "R15","R16",...,"R40","RINF","RSUP".
-
         decontaminate : bool
             Decontaminate the light curve. Default is True.
-
         reject_highpoints : bool
             Reject high points in the light curve. Default is True. 
-        
         bg_MAD_reject : float
             Number of median absolute deviations to reject high background points. Default is 3.
+        outlier_clip : float
+            Number of standard deviations to clip outliers. Default is 4.
+        outlier_window_width : int
+            Window width to use for clipping outliers. Default is 11.
+        outlier_clip_niter : int
+            Number of iterations to clip outliers. Default is 3.
         """
         
         
@@ -268,6 +340,9 @@ def download(self, file_keys="all", aperture="DEFAULT", decontaminate=True, reje
             mask  = d.lc["bg"] > bg_MAD_reject*np.nanmedian(d.lc["bg"])
             _,_,_ = d.mask_data(mask, verbose=False)
 
+            #iterative sigma clipping
+            for _ in range(outlier_clip_niter): t_,f_,e_ = d.clip_outliers(outlier_clip,outlier_window_width, verbose=verbose)
+
             self.lc.append(d)
 
     def scatter(self, figsize=(10,4)):
@@ -551,7 +626,7 @@ def save_CONAN_lcfile(self,bjd_ref = 2450000, folder="data"):
 
 
 
-def get_parameters(planet_name, database="exoplanetarchive", table="pscomppars", overwrite_cache=False):
+def get_parameters(planet_name, database="exoplanetarchive", table="pscomppars", ps_select_index=None,overwrite_cache=False):
     """
     get stellar and planet parameters from nasa exoplanet archive or exoplanet.eu
 
@@ -561,6 +636,9 @@ def get_parameters(planet_name, database="exoplanetarchive", table="pscomppars",
         Name of the planet.
     database : str
         Database to use. Default is "exoplanetarchive". Options are "exoplanetarchive" or "exoplanet.eu".
+    table : str
+        Table to use. Default is "pscomppars". Options are "pscomppars" or "ps".
+
     overwrite_cache : bool
         Overwrite the cached parameters. Default is True.
     """
@@ -580,14 +658,23 @@ def get_parameters(planet_name, database="exoplanetarchive", table="pscomppars",
         if table=="pscomppars":
             df = NasaExoplanetArchive.query_object(planet_name,table=table)
         elif table=="ps":
-            df = NasaExoplanetArchive.query_object(planet_name,table=table) 
-            df = df[df['default_flag']==1]
+            ref_table = NasaExoplanetArchive.query_object(planet_name,table=table)
+            ref_table["pl_refname"] = [ref.split('refstr=')[1].split(' href')[0] for ref in ref_table["pl_refname"]]
+            ref_table["sy_refname"] = [ref.split('refstr=')[1].split(' href')[0] for ref in ref_table["sy_refname"]]
+            print("\nAvailable parameters:\n---------------------")
+            print(ref_table[["pl_name","default_flag","pl_refname","sy_refname",]].to_pandas())
+            if ps_select_index is None:
+                df = ref_table[ref_table['default_flag']==1]
+                print(f"\nselecting default parameters: {df['pl_refname'][0]}. Set `ps_select_index` [between 0 – {len(ref_table)-1}] to select another")
+            else:
+                ref_table["index"] = np.arange(len(ref_table))
+                df = ref_table[ref_table["index"]==ps_select_index]
+                print(f"\nselecting parameters: {df['pl_refname'][0]}")
         else: 
             raise ValueError("table must be either 'pscomppars' or 'ps'")
         
         df = df.to_pandas()
         
-
         params["star"]["Teff"]     = (df["st_teff"][0],df["st_tefferr1"][0])
         params["star"]["logg"]     = (df["st_logg"][0],df["st_loggerr1"][0])
         params["star"]["FeH"]      = (df["st_met"][0],df["st_meterr1"][0])
diff --git a/CONAN3/logprob_multi.py b/CONAN3/logprob_multi.py
index a154e80..c15b185 100755
--- a/CONAN3/logprob_multi.py
+++ b/CONAN3/logprob_multi.py
@@ -5,6 +5,7 @@
 from .utils import rho_to_tdur, gp_params_convert
 import matplotlib
 from ._classes import __default_backend__
+from os.path import splitext
 
 
 def logprob_multi(p, args,t=None,make_outfile=False,verbose=False,debug=False,get_model=False,out_folder=""):
@@ -324,9 +325,9 @@ def logprob_multi(p, args,t=None,make_outfile=False,verbose=False,debug=False,ge
 
                 out_data = np.hstack((out_data,phases))
                 if make_outfile:
-                    outfile=out_folder+"/"+name.split(".")[0]+'_lcout.dat' 
+                    outfile=out_folder+"/"+splitext(name)[0]+'_lcout.dat' 
                     if verbose: print(f"Writing LC output to file: {outfile}")
-                    np.savetxt(outfile,out_data,header=header_fmt.format(*header),fmt='%-16.6f')
+                    np.savetxt(outfile,out_data,header=header_fmt.format(*header),fmt='%-16.6f',delimiter="\t")
         
         elif useGPphot[j] in ['y','ce']: 
             gppars   = params_all[len(params):len(params)+len(GPparams)]   # the GP parameters for all lcs
@@ -384,9 +385,9 @@ def logprob_multi(p, args,t=None,make_outfile=False,verbose=False,debug=False,ge
 
                 out_data = np.hstack((out_data,phases))
                 if make_outfile:
-                    outfile=out_folder+"/"+name.split(".")[0]+'_lcout.dat'
+                    outfile=out_folder+"/"+splitext(name)[0]+'_lcout.dat'
                     if verbose: print(f"Writing LC output with GP({which_GP}) to file: {outfile}")
-                    np.savetxt(outfile,out_data,header=header_fmt.format(*header),fmt='%-16.6f')
+                    np.savetxt(outfile,out_data,header=header_fmt.format(*header),fmt='%-16.6f',delimiter="\t")
 
 
     # now do the RVs and add their probabilities to the model
@@ -534,10 +535,10 @@ def logprob_multi(p, args,t=None,make_outfile=False,verbose=False,debug=False,ge
                     header_fmt += "{:<16s}\t"
                     header     += [f"phase_{n+1}"] if npl>1 else ["phase"]
                 if make_outfile:   
-                    outfile  = out_folder+"/"+RVnames[j].split(".")[0]+'_rvout.dat'
+                    outfile  = out_folder+"/"+splitext(RVnames[j])[0]+'_rvout.dat'
                     out_data = np.hstack((out_data,phases))
                     if verbose: print(f"Writing RV output to file: {outfile}")
-                    np.savetxt(outfile,out_data,header=header_fmt.format(*header),fmt='%-16.6f')
+                    np.savetxt(outfile,out_data,header=header_fmt.format(*header),fmt='%-16.6f',delimiter="\t")
                     # pd.DataFrame(out_data,columns=header).to_csv(outfile,sep="\t",index=False, float_format='%-16.6f')
 
 
@@ -593,9 +594,9 @@ def logprob_multi(p, args,t=None,make_outfile=False,verbose=False,debug=False,ge
                 
                 if make_outfile:   
                     out_data = np.hstack((out_data,phases))
-                    outfile  = out_folder+"/"+RVnames[j].split(".")[0]+'_rvout.dat'
+                    outfile  = out_folder+"/"+splitext(RVnames[j])[0]+'_rvout.dat'
                     if verbose: print(f"Writing RV output with GP({which_GP}) to file: {outfile}")
-                    np.savetxt(outfile,out_data,header=header_fmt.format(*header),fmt='%-16.6f')
+                    np.savetxt(outfile,out_data,header=header_fmt.format(*header),fmt='%-16.6f',delimiter="\t")
 
 
     if get_model:
diff --git a/CONAN3/models.py b/CONAN3/models.py
index c9aa6e0..c85b11c 100755
--- a/CONAN3/models.py
+++ b/CONAN3/models.py
@@ -339,33 +339,6 @@ def get_value(self, tarr, ss=None,grprs=0, vcont=0):
         return mm, model_components
 
 
-        # #==== set up for LC and baseline creation
-        # col5 = np.copy(col5_in)
-        # col3 = np.copy(col3_in)
-        # col4 = np.copy(col4_in)
-        # col6 = np.copy(col6_in)
-        # col7 = np.copy(col7_in)
-        # col8 = np.copy(col8_in)
-        # ts   = tt-np.median(tt)  
-
-        # # MONIKA: added least square optimisation for baselines
-        # if (baseLSQ == 'y'):
-        #     #bvar contains the indices of the non-fixed baseline variables        
-        #     coeffstart  = np.copy(bases[bvar])   
-        #     icoeff,dump = scipy.optimize.leastsq(para_minfunc, coeffstart, args=(bvar, mm, ft, ts, col5, col3, col4, col6, col7,col8))
-        #     coeff = np.copy(bases)
-        #     coeff[bvar] = np.copy(icoeff)   
-        # else:        
-        #     coeff = np.copy(bases)  # the input coefficients 
-
-        # bfunc,spl_comp = basefunc_noCNM(coeff, ts, col5, col3, col4, col6, col7,col8,ft/mm,useSpline)
-        # mod=mm*bfunc
-        
-        # lc_result = SimpleNamespace(planet_LC=mm, full_LCmod=mod, LC_bl=bfunc, spline=spl_comp )
-
-        # return lc_result
-
-
 def basefunc_noCNM(coeff, LCdata,res,useSpline):
     # the full baseline function calculated with the coefficients given; of which some are not jumping and set to 0
     t, flux, err, col3, col4, col5, col6, col7, col8 = LCdata.values()
@@ -548,33 +521,6 @@ def RadialVelocity_Model(tt,T0,per,K,sesinw=0,secosw=0,Gamma=0,npl=None):
 
     return mod_RV, model_components
 
-    # bfstartRV= 1+7*npl + nttv + nddf + nocc*5 + nfilt*2 + nphot+ 2*nRV + nphot*22 +j*12  #the first index in the param array that refers to a baseline function
-    # incoeff = list(range(bfstartRV,bfstartRV+12))  # the indices for the coefficients for the base function        
-
-    # ts = tt-np.median(tt)
-
-    # if (baseLSQ == 'y'):
-    #     #get the indices of the variable baseline parameters via bvar (0 for fixed, 1 for variable)
-    #     ivars = np.copy(bvarsRV[j][0])
-
-    #     incoeff=np.array(incoeff)
-    #     coeffstart = np.copy(params[incoeff[ivars]]) 
-    #     if len(ivars) > 0:
-    #         icoeff,dump = scipy.optimize.leastsq(para_minfuncRV, coeffstart, args=(ivars, mod_RV, RVmes, ts, bis, fwhm, contra))
-    #         coeff = np.copy(params[incoeff])   # the full coefficients -- set to those defined in params (in case any are fixed non-zero)
-    #         coeff[ivars] = np.copy(icoeff)     # and the variable ones are set to the result from the minimization
-    #     else:
-    #         coeff = np.copy(params[incoeff])
-    # else:        
-    #     coeff = np.copy(params[incoeff])   # the coefficients for the base function
-    
-    # bfuncRV,spl_comp=basefuncRV(coeff, ts, bis, fwhm, contra,RVmes-mod_RV ,useSpline)
-    # mod_RVbl = mod_RV + bfuncRV
-
-    # rv_result = SimpleNamespace(planet_RV=mod_RV-Gamma_in, full_RVmod=mod_RVbl, RV_bl=bfuncRV, spline=spl_comp ) 
-    # return rv_result
-
-
 def para_minfuncRV(icoeff, ivars, mod_RV, RVdata):
     """
     least square fit to the baseline after subtracting model RV mod_RV
diff --git a/CONAN3/outputs.py b/CONAN3/outputs.py
index 4bfd28e..7b15cc1 100755
--- a/CONAN3/outputs.py
+++ b/CONAN3/outputs.py
@@ -334,17 +334,17 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
     sig1md =  np.zeros([nderived,2]) 
     sig3md =  np.zeros([nderived,2]) 
     
-    of.write('====================================================================================================\n')
-    of.write(f'{"Jump parameters":25s} {"median":14s} {"-1sigma":14s} {"+1sigma":14s} {"-3sigma":14s} {"+3sigma":14s}\n')
-    of.write('====================================================================================================\n')
+    of.write('#====================================================================================================\n')
+    of.write(f'#{"Jump parameters":25s} {"median":14s} {"-1sigma":14s} {"+1sigma":14s} {"-3sigma":14s} {"+3sigma":14s}\n')
+    of.write('#====================================================================================================\n')
     
-    of2.write('====================================================================================================\n')
-    of2.write(f'{"Jump parameters":25s} {"median":14s} {"-1sigma":14s} {"+1sigma":14s} {"-3sigma":14s} {"+3sigma":14s}\n')
-    of2.write('====================================================================================================\n')
+    of2.write('#====================================================================================================\n')
+    of2.write(f'#{"Jump parameters":25s} {"median":14s} {"-1sigma":14s} {"+1sigma":14s} {"-3sigma":14s} {"+3sigma":14s}\n')
+    of2.write('#====================================================================================================\n')
 
-    of3.write('====================================================================================================\n')
-    of3.write(f'{"Jump parameters":25s} {"median":14s} {"-1sigma":14s} {"+1sigma":14s} {"-3sigma":14s} {"+3sigma":14s}\n')
-    of3.write('====================================================================================================\n')
+    of3.write('#====================================================================================================\n')
+    of3.write(f'#{"Jump parameters":25s} {"median":14s} {"-1sigma":14s} {"+1sigma":14s} {"-3sigma":14s} {"+3sigma":14s}\n')
+    of3.write('#====================================================================================================\n')
 
     for i in range(npara):
         vals=posterior[:,i]
@@ -357,9 +357,9 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         #print jnames[i], medvals[i],sig1[i,0],sig1[i,1], sig3[i,0], sig3[i,1]
         of.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % (jnames[i], medvals[i],sig1[i,0],sig1[i,1], sig3[i,0], sig3[i,1]))       
     
-    of.write('====================================================================================================\n')
-    of.write('Stellar input parameters: \n')
-    of.write('====================================================================================================\n')
+    of.write('#====================================================================================================\n')
+    of.write('#Stellar input parameters: \n')
+    of.write('#====================================================================================================\n')
     if howstellar == 'Rrho':
         vals=Rs_PDF
         # calculate median
@@ -368,7 +368,7 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         sval   = np.sort(dval)
         sig1s  = sval[i1sig] # the 1-sigma intervals (the left side is naturally negative) 
         sig3s  = sval[i3sig] # the 1-sigma intervals (the left side is naturally negative) 
-        of.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('Rstar', medval,sig1s[0],sig1s[1], sig3s[0], sig3s[1]))  
+        of.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('#Rstar', medval,sig1s[0],sig1s[1], sig3s[0], sig3s[1]))  
     if howstellar == 'Mrho':
         vals   = Ms_PDF
         # calculate median
@@ -377,7 +377,7 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         sval   = np.sort(dval)
         sig1s  = sval[i1sig] # the 1-sigma intervals (the left side is naturally negative) 
         sig3s  = sval[i3sig] # the 1-sigma intervals (the left side is naturally negative) 
-        of.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('Mstar', medval,sig1s[0],sig1s[1], sig3s[0], sig3s[1]))  
+        of.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('#Mstar', medval,sig1s[0],sig1s[1], sig3s[0], sig3s[1]))  
 
     for i in range(len(extinpars)):
         vals   = extind_PDF[:,i]
@@ -389,9 +389,9 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         of.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % (extinpars[i], medval,sig1s[0],sig1s[1], sig3s[0], sig3s[1]))  
 
 
-    of.write('====================================================================================================\n')
-    of.write('Derived parameters: ('+starstring+') \n')
-    of.write('====================================================================================================\n')
+    of.write('#====================================================================================================\n')
+    of.write('#Derived parameters: ('+starstring+') \n')
+    of.write('#====================================================================================================\n')
         
     for i in range(nderived):
         vals=derived_PDFs[i]
@@ -403,9 +403,9 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         sig3d[i]    = sval[i3sig] # the 1-sigma intervals (the left side is naturally negative) 
         of.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % (derived_pnames[i], medvalsd[i],sig1d[i,0],sig1d[i,1], sig3d[i,0], sig3d[i,1])) 
 
-    of.write('====================================================================================================\n')     
-    of.write('Fixed parameters: \n')
-    of.write('====================================================================================================\n')
+    of.write('#====================================================================================================\n')     
+    of.write('#Fixed parameters: \n')
+    of.write('#====================================================================================================\n')
 
     nfix = len(njnames)
     for i in range(nfix):
@@ -414,7 +414,7 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         else:
             of.write('%-25s %14.8f \n' % (njnames[i], bp[nijnames[0][i]])) 
         
-    of.write('====================================================================================================\n')
+    of.write('#====================================================================================================\n')
     of.close()
     
     #now write out outfile2: the peak of the posterior and the area containing 68% of points
@@ -431,9 +431,9 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         of2.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % (jnames[i], maxvals[i],sig1m[i,0],sig1m[i,1], sig3m[i,0], sig3m[i,1])) 
         param_histbp(vals,jnames[i],medvals[i],sig1[i],sig3[i],maxvals[i],sig1m[i],sig3m[i],bp[ijnames[0][i]],s1bps,out_folder)
 
-    of2.write('====================================================================================================\n')
-    of2.write('Stellar input parameters: \n')
-    of2.write('====================================================================================================\n')
+    of2.write('#====================================================================================================\n')
+    of2.write('#Stellar input parameters: \n')
+    of2.write('#====================================================================================================\n')
     if howstellar == 'Rrho':
         vals=Rs_PDF
         pdf, xpdf, HPDmin, iHDP = credregionML(vals)
@@ -443,7 +443,7 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         pdf, xpdf, HPDmin, iHDP = credregionML(vals,pdf=pdf, xpdf=xpdf, percentile=0.9973)
         sig3ms[0] = np.amin(xpdf[pdf>HPDmin]) - maxval
         sig3ms[1] = np.amax(xpdf[pdf>HPDmin]) - maxval      
-        of2.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('Rstar', maxval,sig1ms[0],sig1ms[1], sig3ms[0], sig3ms[1])) 
+        of2.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('#Rstar', maxval,sig1ms[0],sig1ms[1], sig3ms[0], sig3ms[1])) 
         param_hist(vals,'Rstar',medval,sig1s,sig3s,maxval,sig1ms,sig3ms,out_folder=out_folder)
     if howstellar == 'Mrho':
         vals=Ms_PDF
@@ -454,7 +454,7 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         pdf, xpdf, HPDmin, iHDP = credregionML(vals,pdf=pdf, xpdf=xpdf, percentile=0.9973)
         sig3ms[0] = np.amin(xpdf[pdf>HPDmin]) - maxval
         sig3ms[1] = np.amax(xpdf[pdf>HPDmin]) - maxval      
-        of2.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('Mstar', maxval,sig1ms[0],sig1ms[1], sig3ms[0], sig3ms[1])) 
+        of2.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('#Mstar', maxval,sig1ms[0],sig1ms[1], sig3ms[0], sig3ms[1])) 
         param_hist(vals,'Mstar',medval,sig1s,sig3s,maxval,sig1ms,sig3ms,out_folder=out_folder)
 
     for i in range(len(extinpars)):
@@ -468,9 +468,9 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         sig3ms[1] = np.amax(xpdf[pdf>HPDmin]) - maxval      
         of2.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % (extinpars[i], maxval,sig1ms[0],sig1ms[1], sig3ms[0], sig3ms[1]))  
     
-    of2.write('====================================================================================================\n')
-    of2.write('Derived parameters: ('+starstring+') \n')
-    of2.write('====================================================================================================\n')    
+    of2.write('#====================================================================================================\n')
+    of2.write('#Derived parameters: ('+starstring+') \n')
+    of2.write('#====================================================================================================\n')    
     
     for i in range(nderived):     
         vals=derived_PDFs[i]
@@ -501,9 +501,9 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         of2.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % (derived_pnames[i], maxvalsd[i], sig1md[i,0],sig1md[i,1], sig3md[i,0], sig3md[i,1])) 
         param_hist(vals,derived_pnames[i],medvalsd[i],sig1d[i],sig3d[i],maxvalsd[i],sig1md[i],sig3md[i],out_folder=out_folder)
     
-    of2.write('====================================================================================================\n')
-    of2.write('Fixed parameters: \n')
-    of2.write('====================================================================================================\n')
+    of2.write('#====================================================================================================\n')
+    of2.write('#Fixed parameters: \n')
+    of2.write('#====================================================================================================\n')
 
     nfix = len(njnames)
     for i in range(nfix):
@@ -513,7 +513,7 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
             of2.write('%-25s %14.8f \n' % (njnames[i], bp[nijnames[0][i]])) 
 
     
-    of2.write('====================================================================================================\n')    
+    of2.write('#====================================================================================================\n')    
     of2.close()
         
         
@@ -533,9 +533,9 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         
         of3.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % (jnames[i],bp[ijnames[0][i]],s1bps[0],s1bps[1],s3bps[0],s3bps[1]))     
         
-    of3.write('====================================================================================================\n')
-    of3.write('Stellar input parameters: \n')
-    of3.write('====================================================================================================\n')
+    of3.write('#====================================================================================================\n')
+    of3.write('#Stellar input parameters: \n')
+    of3.write('#====================================================================================================\n')
     if howstellar == 'Rrho':
         vals=Rs_PDF
         pdf, xpdf, HPDmin, iHDP = credregionML(vals)
@@ -545,7 +545,7 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         pdf, xpdf, HPDmin, iHDP = credregionML(vals,pdf=pdf, xpdf=xpdf, percentile=0.9973)
         sig3ms[0] = np.amin(xpdf[pdf>HPDmin]) - maxval
         sig3ms[1] = np.amax(xpdf[pdf>HPDmin]) - maxval      
-        of3.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('Rstar', maxval,sig1ms[0],sig1ms[1], sig3ms[0], sig3ms[1])) 
+        of3.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('#Rstar', maxval,sig1ms[0],sig1ms[1], sig3ms[0], sig3ms[1])) 
         param_hist(vals,'Rstar',medval,sig1s,sig3s,maxval,sig1ms,sig3ms,out_folder=out_folder)
     if howstellar == 'Mrho':
         vals=Ms_PDF
@@ -556,7 +556,7 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         pdf, xpdf, HPDmin, iHDP = credregionML(vals,pdf=pdf, xpdf=xpdf, percentile=0.9973)
         sig3ms[0] = np.amin(xpdf[pdf>HPDmin]) - maxval
         sig3ms[1] = np.amax(xpdf[pdf>HPDmin]) - maxval      
-        of3.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('Mstar', maxval,sig1ms[0],sig1ms[1], sig3ms[0], sig3ms[1])) 
+        of3.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % ('#Mstar', maxval,sig1ms[0],sig1ms[1], sig3ms[0], sig3ms[1])) 
         param_hist(vals,'Mstar',medval,sig1s,sig3s,maxval,sig1ms,sig3ms,out_folder=out_folder)
     for i in range(len(extinpars)):
         vals = extind_PDF[:,i]
@@ -577,9 +577,9 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
     
     # FIRST: derive the bp values, call them der_bp
 
-    of3.write('====================================================================================================\n')
-    of3.write('Derived parameters: ('+starstring+') \n')
-    of3.write('====================================================================================================\n')    
+    of3.write('#====================================================================================================\n')
+    of3.write('#Derived parameters: ('+starstring+') \n')
+    of3.write('#====================================================================================================\n')    
     
     for i in range(nderived):     
         
@@ -588,9 +588,9 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         
         of3.write('%-25s %14.8f %14.8f %14.8f %14.8f %14.8f\n' % (derived_pnames[i], derived_bp[i], sig1bp[0],sig1bp[1], sig3bp[0], sig3bp[1])) 
     
-    of3.write('====================================================================================================\n')
-    of3.write('Fixed parameters: \n')
-    of3.write('====================================================================================================\n')
+    of3.write('#====================================================================================================\n')
+    of3.write('#Fixed parameters: \n')
+    of3.write('#====================================================================================================\n')
 
     nfix = len(njnames)
     for i in range(nfix):
@@ -599,7 +599,7 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         else:
             of3.write('%-25s %14.8f \n' % (njnames[i], bp[nijnames[0][i]])) 
     
-    of3.write('====================================================================================================\n')    
+    of3.write('#====================================================================================================\n')    
     of3.close()        
         
         
@@ -616,7 +616,7 @@ def mcmc_outputs(posterior, jnames, ijnames, njnames, nijnames, bp, uwl, Rs_in,
         plot_traspec(dRpRsres, edRpRsres, uwl,out_folder)
 
     
-    return medvals, maxvals
+    return medvals, maxvals,np.median(abs(sig1),axis=1)
 
 
 def gr_print(jnames,GRvals, out_folder):
@@ -731,7 +731,7 @@ def derive_parameters(filnames, nm, Rs_PDF, Ms_PDF, RpRs_PDF, Period_PDF, b_PDF,
     a_PDF = aRs_PDF * Rs_PDF * Rsolar / au
     Rsa_PDF = 1./aRs_PDF
 
-    rhoS_PDF = rhoS_PDF / (rhoSolar/1000)   #solar units
+    rhoS_PDF = rhoS_PDF #/ (rhoSolar/1000)   #solar units
 
     ome_PDF = ome_PDF * 180. / np.pi
 
@@ -796,7 +796,7 @@ def derive_parameters(filnames, nm, Rs_PDF, Ms_PDF, RpRs_PDF, Period_PDF, b_PDF,
             Fn_PDFs.append(Fn)
 
 
-    derived_pnames = [f"Rp{nm}_[Rjup]",f"Mp{nm}_[Mjup]", f"rho{nm}_p_[rhoJup]", f"g_p{nm}_[SI]", f"dF{nm}", f"aRs{nm}", f"a{nm}_[au]", f"rhoS{nm}_[rhoSun]", "Ms_[Msun]", "Rs_[Rsun]",
+    derived_pnames = [f"Rp{nm}_[Rjup]",f"Mp{nm}_[Mjup]", f"rho{nm}_p_[rhoJup]", f"g_p{nm}_[SI]", f"dF{nm}", f"aRs{nm}", f"a{nm}_[au]", f"rho_star{nm}_[g_cm3]", "Ms_[Msun]", "Rs_[Rsun]",
                         f"inclination{nm}_[deg]", f"eccentricity{nm}", f"omega{nm}_[deg]", f"Occult_dur{nm}", f"Rs_a{nm}", "MF_PDF_[Msun]",f"Dur{nm}_[d]"]
     
     derived_pnames =  derived_pnames + pnames_LD + pnames_Fn
diff --git a/CONAN3/plotting.py b/CONAN3/plotting.py
index 09b85f6..ebc3e59 100755
--- a/CONAN3/plotting.py
+++ b/CONAN3/plotting.py
@@ -6,6 +6,7 @@
 import matplotlib.pyplot as plt    
 from CONAN3.logprob_multi import logprob_multi
 import pickle
+from os.path import splitext
 
 
 def fit_plots(nttv, nphot, nRV, filters,names,RVnames,out_folder,prefix="/",RVunit="km/s",params=None,T0=None,period=None,Dur=None):
@@ -20,7 +21,7 @@ def fit_plots(nttv, nphot, nRV, filters,names,RVnames,out_folder,prefix="/",RVun
     
     #model plot for each LC
     for j in range(nphot):
-        infile=out_folder + "/" + names[j].split(".")[0]+'_lcout.dat'
+        infile=out_folder + "/" + splitext(names[j])[0]+'_lcout.dat'
         tt, flux, err, full_mod, bfunc, mm, det_flux = np.loadtxt(infile, usecols=(0,1,2,3,7,8,9), unpack = True)  # reading in the lightcurve data
         
         #evaluate rv model on smooth time grid
@@ -37,7 +38,7 @@ def fit_plots(nttv, nphot, nRV, filters,names,RVnames,out_folder,prefix="/",RVun
         _,    resbin          = bin_data_with_gaps(tt, flux-full_mod,binsize=binsize)    #residuals
 
         ########## Plot and save lightcurve with fit ########
-        outname=out_folder+prefix+names[j].split(".")[0]+'_fit.png'
+        outname=out_folder+prefix+splitext(names[j])[0]+'_fit.png'
 
         fig,ax = plt.subplots(3,1, figsize=(12,12), sharex=True,gridspec_kw={"height_ratios":(3,3,1)})
         ax[0].set_title('Fit for lightcurve '+names[j][:-4])
@@ -136,8 +137,8 @@ def fit_plots(nttv, nphot, nRV, filters,names,RVnames,out_folder,prefix="/",RVun
 
     ############ RVs#####################
     for j in range(nRV):
-        infile  = out_folder + "/" + RVnames[j].split(".")[0]+'_rvout.dat'
-        outname = out_folder+prefix+RVnames[j].split(".")[0]+'_fit.png'
+        infile  = out_folder + "/" + splitext(RVnames[j])[0]+'_rvout.dat'
+        outname = out_folder+prefix+splitext(RVnames[j])[0]+'_fit.png'
         tt, y_rv , e_rv, full_mod, base, rv_mod, det_RV = np.loadtxt(infile, usecols=(0,1,2,3,7,8,9),unpack = True)  # reading in the rvcurve data
         rv_resid = y_rv-full_mod
 
@@ -188,7 +189,7 @@ def fit_plots(nttv, nphot, nRV, filters,names,RVnames,out_folder,prefix="/",RVun
             phase_all = []
 
             for j in range(nRV):
-                infile  = out_folder + "/" + RVnames[j].split(".")[0]+'_rvout.dat'
+                infile  = out_folder + "/" + splitext(RVnames[j])[0]+'_rvout.dat'
                 tt, y_rv , e_rv, full_mod, base, rv_mod, det_RV = np.loadtxt(infile, usecols=(0,1,2,3,7,8,9), unpack = True)
                 rv_resid = y_rv-full_mod
 
diff --git a/CONAN3/utils.py b/CONAN3/utils.py
index bc392a6..00aabc0 100755
--- a/CONAN3/utils.py
+++ b/CONAN3/utils.py
@@ -417,7 +417,7 @@ def aR_to_Tdur(aR, b, Rp, P,e=0,w=90, tra_occ="tra"):
     factr =  ((1+Rp)**2 - b**2)/(aR**2-b**2)
     ecc_fac = np.sqrt(1-e**2)/(1+e*np.sin(np.deg2rad(w)))  if tra_occ=="tra" else np.sqrt(1-e**2)/(1-e*np.sin(np.deg2rad(w)))
     Tdur = (P/np.pi)*np.arcsin( np.sqrt(factr) ) * ecc_fac
-    return Tdur
+    return np.round(Tdur,8)
 
 
 def Tdur_to_aR(Tdur, b, Rp, P,e=0,w=90, tra_occ = "tra"):
diff --git a/change_log.rst b/change_log.rst
index c7a97ca..dbb3e30 100755
--- a/change_log.rst
+++ b/change_log.rst
@@ -1,3 +1,14 @@
+3-Jun-2024: version 3.2.1
+~~~~~~~~~~~~~
+* added function to read the parameters values and errors from the result_**.dat file --> result.get_all_params_dict()
+* automatically save fit config file in the output folder as config_save.dat
+* input files no longer overwritten when modified within CONAN. CONAN injests the file and works with a copy of it.
+* fixed bug in splitting filenames with multiple dots
+* selection of system parameter reference when using 'ps' table of the NASA exoplanet archive
+* changed reported stellar density unit from rho_sun to g/cm3
+* outlier rejection can be performed on selected columns of the data. e.g. lc_obj.clip_outliers(clip=4,width=11,selected_column=["col1","col3","col4"])
+* when Downloading TESS data, user can now split the sector into orbits to create different file that can be detrended differently
+* fixed bug in writing result files when duration is a fixed parameter
 
 7-May-2024: version 3.2.0
 ~~~~~~~~~~~~~