From 0aa8117b10dd155d64ff57f98790ec4536244d30 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <Sherelyn.Alejandro78@myhunter.cuny.edu>
Date: Tue, 21 Nov 2023 15:15:59 -0500
Subject: [PATCH 01/29] provide upper and lower errors when determining
 fundamental parameters.

---
 sedkit/isochrone.py |  95 +++++++++++++++++-------
 sedkit/sed.py       | 176 +++++++++++++++++++++++++++++++++-----------
 sedkit/utilities.py |   3 +-
 3 files changed, 205 insertions(+), 69 deletions(-)
 mode change 100755 => 100644 sedkit/isochrone.py
 mode change 100755 => 100644 sedkit/sed.py
 mode change 100755 => 100644 sedkit/utilities.py

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
old mode 100755
new mode 100644
index 25df79fa..75e0fe55
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -19,6 +19,7 @@
 
 from . import utilities as u
 
+
 # A dictionary of all supported moving group ages from Bell et al. (2015)
 NYMG_AGES = {'AB Dor': (149 * q.Myr, 51 * q.Myr, '2015MNRAS.454..593B'),
              'beta Pic': (24 * q.Myr, 3 * q.Myr, '2015MNRAS.454..593B'),
@@ -34,7 +35,7 @@
     EVO_MODELS = [os.path.basename(m).replace('.txt', '') for m in glob.glob(resource_filename('sedkit', 'data/models/evolutionary/*'))]
 # Fails RTD build for some reason
 except:
-    EVO_MODELS = ['COND03', 'dmestar_solar', 'DUSTY00', 'f2_solar_age', 'hybrid_solar_age', 'nc+03_age', 'nc-03_age', 'nc_solar_age', 'parsec12_solar']
+    EVO_MODELS = ['COND03', 'dmestar_solar', 'DUSTY00', 'f2_solar_age', 'hybrid_solar_age', 'nc+03_age','nc+0.0_age', 'nc-03_age', 'nc_solar_age','nc+0.5_age','nc-0.5_age', 'parsec12_solar','ATMO_NEQ_strong','ATMO_NEQ_strong_MIRI']
 
 class Isochrone:
     """A class to handle model isochrones"""
@@ -166,39 +167,77 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
             self.message('{}: age must be between {} and {} to infer {} from {} isochrones.'.format(*args))
             return None
 
-        # Get the lower, nominal, and upper values
-        lower = self.interpolate(xval[0] - xval[1], age[0]-age[1], xparam, yparam)
-        nominal = self.interpolate(xval[0], age[0], xparam, yparam)
-        upper = self.interpolate(xval[0] + xval[1], age[0]+age[1], xparam, yparam)
+        # Monte Carlo Approach
+        mu, sigma = xval[0], xval[1]  # mean and standard deviation for values on the x-axis
+        mu_a, sigma_a = age[0].value, age[1].value  # mean and standard deviation for the age range provided
 
-        if nominal is None:
-            return None
-        if lower is None:
-            lower = upper
-        if upper is None:
-            upper = lower
-        if upper is None:
-            return None
+        xsample = np.random.normal(mu, sigma, 1000)
+        ysample = np.random.uniform(mu_a-sigma_a, mu_a+sigma_a, 1000)
+        values_list = []
+        nan_counter = 0
+
+        for x, y in zip(xsample, ysample):
+            result = self.interpolate(x, y * age[0].unit, xparam, yparam)
+            if result is not None:
+                values_list.append(result.value if hasattr(result, 'unit') else result)
+            elif result is None:
+                print("Interpolate resulted in NaN")
+                nan_counter = nan_counter + 1
+
+        # Get final values
+        unit = self.data[yparam].unit or 1
+        # average = np.mean(np.array(values_list) * unit)
+        # # std = np.std(values_list)
+        # # error = std * unit
+        # min = np.min(np.array(values_list) * unit)
+        # max = np.max(np.array(values_list) * unit)
+        # error = abs(max-min)/2
+
+        #Error Propagation
+        dist_hist, hist_edges = np.histogram(values_list)
+
+        #Normalize Histogram
+        dist_hist = dist_hist/ np.sum(dist_hist)
+        # cumulative PDF
+        cum_PDF_hist_x = hist_edges[1:]  # right edges of the bins
+        cum_PDF_hist = np.zeros(hist_edges.size - 1)
+        for i in range(hist_edges.size - 1):
+            cum_PDF_hist[i] = np.sum(dist_hist[:i + 1])
+
+        # interquartile range of the PDF
+        error_lower_per = 0.16
+        error_upper_per = 0.84
+
+        Q1_hist = np.interp(error_lower_per, cum_PDF_hist, cum_PDF_hist_x)  # lower error
+        Q2_hist = np.interp(0.50, cum_PDF_hist, cum_PDF_hist_x )    # median
+        Q3_hist = np.interp(error_upper_per, cum_PDF_hist, cum_PDF_hist_x)  # upper error
+
+        # print("This is the lower error",Q1_hist)
+        # print("This is the median", Q2_hist)
+        # print("This is the upper error", Q3_hist)
+        lower_err = (Q2_hist - Q1_hist) * unit
+        upper_err = (Q3_hist - Q2_hist) * unit
+        average = Q2_hist * unit
 
-        # Calculate the symmetric error
-        error = max(abs(nominal - lower), abs(nominal - upper)) * 2
 
         # Plot the figure and evaluated point
         if plot:
-            val = nominal.value if hasattr(nominal, 'unit') else nominal
-            err = error.value if hasattr(error, 'unit') else error
+            # val = average.value if hasattr(average, 'unit') else average
+            # lower_err = lower_err.value if hasattr(lower_err, 'unit') else lower_err
+            # upper_err = upper_err.value if hasattr(upper_err, 'unit') else upper_err
             fig = self.plot(xparam, yparam)
-            legend = '{} = {:.3f} ({:.3f})'.format(yparam, val, err)
-            fig.circle(xval[0], val, color='red', legend=legend)
-            u.errorbars(fig, [xval[0]], [val], xerr=[xval[1]*2], yerr=[err], color='red')
-
+            legend = '{} = {:.3f} ({:-.3f} {:+.3f})'.format(yparam, average.value, lower_err.value,upper_err.value)
+            fig.circle(xval[0], average.value, color='purple', legend_label=legend)
+            u.errorbars(fig, [xval[0]], [average.value], xerr=[xval[1]*2], yerr=[lower_err.value*-1,upper_err.value], color='purple')
             show(fig)
 
-        # Balk at nans
-        if np.isnan(nominal):
-            raise ValueError("I got a nan from {} for some reason.".format(self.name))
-
-        return nominal, error, self.name
+        # # Balk at nans
+        if np.isnan(average):
+            print("I got a nan from {} for some reason.".format(self.name))
+            # raise ValueError("I got a nan from {} for some reason.".format(self.name))
+        return average, lower_err, upper_err, self.name
+        # return average, lower_error, upper_error, self.name
+        # return average, error, self.name
 
     def interpolate(self, xval, age, xparam, yparam):
         """Interpolate a value between two isochrones
@@ -250,6 +289,7 @@ def interpolate(self, xval, age, xparam, yparam):
         result = np.average([lower_val, upper_val], weights=weights)
         unit = self.data[yparam].unit or 1
 
+        # return result
         return result * unit
 
     def message(self, msg, pre='[sedkit]'):
@@ -301,6 +341,9 @@ def plot(self, xparam, yparam, draw=False, **kwargs):
         color_mapper = LinearColorMapper(palette='Viridis256',
                                          low=self.ages.min().value,
                                          high=self.ages.max().value)
+        # color_mapper = LinearColorMapper(palette='Viridis256',
+        #                                  low=.05,
+        #                                  high=.15)
         color_bar = ColorBar(color_mapper=color_mapper, ticker=BasicTicker(),
                              label_standoff=5, border_line_color=None,
                              title='Age [{}]'.format(self.age_units),
diff --git a/sedkit/sed.py b/sedkit/sed.py
old mode 100755
new mode 100644
index 4f7c39c4..d86dbe07
--- a/sedkit/sed.py
+++ b/sedkit/sed.py
@@ -8,13 +8,13 @@
 
 Author: Joe Filippazzo, jfilippazzo@stsci.edu
 """
-
 from copy import copy
 import os
 import shutil
 import time
 import warnings
 
+import itertools
 import astropy.table as at
 import astropy.units as q
 import astropy.io.ascii as ii
@@ -506,7 +506,7 @@ def age(self, age):
         age: sequence
             The age and uncertainty in distance units
         """
-        self._validate_and_set_param('age', age, q.Gyr, True, vmin=0 * q.Myr, vmax=13.8 * q.Gyr)
+        self._validate_and_set_param('age', age, q.Gyr, True, vmin=0 * q.Gyr, vmax=13.8 * q.Gyr)
 
     def _calculate_sed(self):
         """
@@ -589,6 +589,7 @@ def calculate_Lbol(self):
                 Lbol_unc = Lbol * np.sqrt((self.fbol[1] / self.fbol[0]).value**2 + (2 * self.distance[1] / self.distance[0]).value**2)
                 Lbol_sun_unc = round(abs(Lbol_unc / (Lbol * np.log(10))).value, 3)
 
+
             # Update the attributes
             self.Lbol = Lbol, Lbol_unc, 'This Work'
             self.Lbol_sun = Lbol_sun, Lbol_sun_unc, 'This Work'
@@ -704,9 +705,10 @@ def calculate_Teff(self):
                 Teff_unc = None
             else:
                 Teff_unc = (Teff * np.sqrt((self.Lbol[1] / self.Lbol[0]).value**2 + (2 * self.radius[1] / self.radius[0]).value**2) / 4.).astype(int)
+                Teff_unc_u = (Teff * np.sqrt((self.Lbol[1] / self.Lbol[0]).value**2 + (2 * self.radius[2] / self.radius[0]).value**2) / 4.).astype(int)
 
             # Update the attribute
-            self.Teff = Teff, Teff_unc, 'This Work'
+            self.Teff = Teff, Teff_unc,Teff_unc_u, 'This Work'
 
     def _calibrate_photometry(self, name='photometry'):
         """
@@ -745,13 +747,17 @@ def _calibrate_photometry(self, name='photometry'):
             if self.distance is not None:
 
                 # Calculate abs_mags
-                M, M_unc = u.flux_calibrate(row['app_magnitude'], self.distance[0], row['app_magnitude_unc'], self.distance[1])
-                table['abs_magnitude'][n] = M
-                table['abs_magnitude_unc'][n] = M_unc
+                for n,row in enumerate(table):
+
+                    M, M_unc = u.flux_calibrate(row['app_magnitude'], self.distance[0], row['app_magnitude_unc'], self.distance[1])
+                    table['abs_magnitude'][n] = M
+                    table['abs_magnitude_unc'][n] = M_unc
 
                 # Calculate abs_flux values
                 for n, row in enumerate(table):
+                    
                     abs_flux, abs_flux_unc = u.mag2flux(row['bandpass'], row['abs_magnitude'], sig_m=row['abs_magnitude_unc'], units=self.flux_units)
+
                     table['abs_flux'][n] = abs_flux
                     table['abs_flux_unc'][n] = abs_flux_unc
 
@@ -1166,6 +1172,7 @@ def export(self, parentdir='.', dirname=None, zipped=False):
         if self.photometry is not None:
             photpath = os.path.join(dirpath, '{}_photometry.txt'.format(name))
             phot_table = copy(self.photometry)
+            # print(phot_table)
             for colname in phot_table.colnames:
                 phot_table.rename_column(colname, colname.replace('/', '_'))
             phot_table.write(photpath, format='ipac')
@@ -1872,8 +1879,7 @@ def infer_logg(self, plot=False):
                 self.message("Could not calculate surface gravity.")
 
             # Store the value
-            self.logg = [logg[0].round(2), logg[1].round(2), logg[2]] if logg is not None else logg
-
+            self.logg = [logg[0].round(2), logg[1].round(2),logg[2].round(2),logg[3]] if logg is not None else logg
         # No dice
         else:
             self.message('Could not calculate logg without Lbol and age')
@@ -1908,7 +1914,7 @@ def infer_mass(self, mass_units=q.Msun, plot=False):
                 mass = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'mass', plot=plot)
 
             # Store the value
-            self.mass = [mass[0].round(3), mass[1].round(3), mass[2]] if mass is not None else mass
+            self.mass = [mass[0].round(0),  mass[1].round(0), mass[2].round(0), mass[3]] if mass is not None else mass
 
         else:
 
@@ -1935,7 +1941,7 @@ def infer_mass(self, mass_units=q.Msun, plot=False):
                 self.message('Could not calculate mass without Lbol, M_2MASS.J, or M_2MASS.Ks')
 
 
-    def infer_radius(self, radius_units=q.Rsun, infer_from=None, plot=False):
+    def infer_radius(self, radius_units=q.Rsun, infer_from='Lbol', plot=False):
         """
         Estimate the radius from model isochrones given an age and Lbol
 
@@ -1974,7 +1980,7 @@ def infer_radius(self, radius_units=q.Rsun, infer_from=None, plot=False):
                 radius = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'radius', plot=plot)
 
             # Store the value
-            self.radius = [radius[0].round(3), radius[1].round(3), radius[2]] if radius is not None else radius
+            self.radius = [radius[0].round(3), radius[1].round(3),radius[2].round(3),radius[3]] if radius is not None else radius
 
         else:
 
@@ -2033,7 +2039,7 @@ def infer_Teff(self, teff_units=q.K, plot=False):
                 teff = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'teff', plot=plot)
 
             # Store the value
-            self.Teff = [teff[0].round(0), teff[1].round(0), teff[2]] if teff is not None else teff
+            self.Teff = [teff[0].round(0), teff[1].round(0), teff[2].round(0),teff[3]] if teff is not None else teff
 
         else:
 
@@ -2487,9 +2493,9 @@ def photometry(self):
         self._photometry.sort('eff')
         return self._photometry
 
-    def plot(self, app=True, photometry=True, spectra=True, integral=True, synthetic_photometry=False,
+    def plot(self, app=False, photometry=True, spectra=True, integral=True, synthetic_photometry=False,
              best_fit=True, normalize=None, scale=['log', 'log'], output=False, fig=None,
-             color='#1f77b4', one_color=False, label=None, **kwargs):
+             color='#3a243b', one_color=False, label=None, **kwargs):
         """
         Plot the SED
 
@@ -2576,10 +2582,11 @@ def plot(self, app=True, photometry=True, spectra=True, integral=True, synthetic
 
         # ...or make a new plot
         else:
-            TOOLS = ['pan', 'reset', 'box_zoom', 'wheel_zoom', 'save']
+            # TOOLS = ['pan', 'reset', 'box_zoom', 'wheel_zoom', 'save']
+            TOOLS = ['pan', 'reset', 'box_zoom', 'save']
             xlab = 'Wavelength [{}]'.format(self.wave_units)
             ylab = 'Flux Density [{}]'.format(str(self.flux_units))
-            self.fig = figure(plot_width=900, plot_height=400, title=self.name,
+            self.fig = figure(width=900, height=400, title=self.name,
                               y_axis_type=scale[1], x_axis_type=scale[0],
                               x_axis_label=xlab, y_axis_label=ylab,
                               tools=TOOLS)
@@ -2587,10 +2594,13 @@ def plot(self, app=True, photometry=True, spectra=True, integral=True, synthetic
         # Plot spectra
         if spectra and len(self.spectra) > 0:
 
-            if spectra == 'all':
+            if (spectra == 'all' and app is False):
+                for n, spec in enumerate(self.spectra['spectrum']):
+                    self.fig = spec.plot(fig=self.fig, components=True, const=(spec.flux_calibrate(self.distance).flux/spec.flux))
+            elif (spectra == 'all' and app is True):
                 for n, spec in enumerate(self.spectra['spectrum']):
+                    print(const)
                     self.fig = spec.plot(fig=self.fig, components=True, const=const)
-
             else:
                 self.fig.line(spec_SED.wave, spec_SED.flux * const, color=color, alpha=0.8, legend_label='Spectrum')
 
@@ -2599,28 +2609,36 @@ def plot(self, app=True, photometry=True, spectra=True, integral=True, synthetic
 
             # Set up hover tool
             phot_tips = [('Band', '@desc'), ('Wave', '@x'), ('Flux', '@y'), ('Unc', '@z')]
-            hover = HoverTool(names=['photometry', 'nondetection'], tooltips=phot_tips)
-            self.fig.add_tools(hover)
+
 
             # Plot points with errors
             pts = np.array([(bnd, wav, flx * const, err * const) for bnd, wav, flx, err in np.array(self.photometry['band', 'eff', pre + 'flux', pre + 'flux_unc']) if not any([(np.isnan(i) or i <= 0) for i in [wav, flx, err]])], dtype=[('desc', 'S20'), ('x', float), ('y', float), ('z', float)])
             if len(pts) > 0:
                 source = ColumnDataSource(data=dict(x=pts['x'], y=pts['y'], z=pts['z'], desc=[b.decode("utf-8") for b in pts['desc']]))
-                self.fig.circle('x', 'y', source=source, legend_label='Photometry', name='photometry', color=color, fill_alpha=0.7, size=8)
+                c1 = self.fig.circle('x', 'y', source=source, legend_label='Photometry', name='photometry', color='#a64d79', fill_alpha=0.7, size=8)
+                self.fig.circle('x', 'y', source=source, legend_label='Photometry', name='photometry', color='#a64d79', fill_alpha=0.7, size=8)
                 self.fig = u.errorbars(self.fig, 'x', 'y', yerr='z', source=source, color=color)
+                hover = HoverTool(tooltips=phot_tips, renderers=[c1])
+                self.fig.add_tools(hover)
+
+
 
             # Plot points without errors
             pts = np.array([(bnd, wav, flx * const, err * const) for bnd, wav, flx, err in np.array(self.photometry['band', 'eff', pre + 'flux', pre + 'flux_unc']) if (np.isnan(err) or err <= 0) and not np.isnan(flx)], dtype=[('desc', 'S20'), ('x', float), ('y', float), ('z', float)])
             if len(pts) > 0:
                 source = ColumnDataSource(data=dict(x=pts['x'], y=pts['y'], z=pts['z'], desc=[b.decode("utf-8") for b in pts['desc']]))
+                c2 = self.fig.circle('x', 'y', source=source, legend_label='Limit', name='nondetection', color=color,fill_alpha=0, size=8)
                 self.fig.circle('x', 'y', source=source, legend_label='Limit', name='nondetection', color=color, fill_alpha=0, size=8)
+                hover1 = HoverTool(tooltips=phot_tips, renderers=[c2])
+                self.fig.add_tools(hover1)
+
 
         # Plot photometry
         if synthetic_photometry and len(self.synthetic_photometry) > 0:
 
             # Set up hover tool
-            phot_tips = [('Band', '@desc'), ('Wave', '@x'), ('Flux', '@y'), ('Unc', '@z')]
-            hover = HoverTool(names=['synthetic photometry'], tooltips=phot_tips, mode='vline')
+            phot_tips = [('Band', '@{desc}'), ('Wave', '@{x}'), ('Flux', '@{y}'), ('Unc', '@{z}')]
+            hover = HoverTool(name='synthetic photometry', tooltips=phot_tips, mode='vline')
             self.fig.add_tools(hover)
 
             # Plot points with errors
@@ -2642,9 +2660,9 @@ def plot(self, app=True, photometry=True, spectra=True, integral=True, synthetic
                 mod = mod_fit['full_model']
                 mod.wave_units = self.wave_units
                 if mod_fit['fit_to'] == 'phot':
-                    self.fig.square(mod.wave, mod.flux, alpha=0.3, color=color if one_color else next(col_list), legend_label=mod_fit['label'], size=12)
+                    self.fig.square(mod.wave, mod.flux, alpha=0.9, color=color if one_color else next(col_list), legend_label=mod_fit['label'], size=12)
                 else:
-                    self.fig.line(mod.wave, mod.flux, alpha=0.3, color=color if one_color else next(col_list), legend_label=mod_fit['label'], line_width=2)
+                    self.fig.line(mod.wave, mod.flux, alpha=0.9, color=color if one_color else next(col_list), legend_label=mod_fit['label'], line_width=2)
 
         self.fig.legend.location = "top_right"
         self.fig.legend.click_policy = "hide"
@@ -2722,25 +2740,42 @@ def results(self):
         # Get the params
         rows = []
         for param in params:
-
             # Get the values and format
             attr = getattr(self, param, None)
-
             if attr is None:
                 attr = '--'
-
+            # error = '$ ^{+0.004}_{-0.004}$'
             if isinstance(attr, (tuple, list)):
-                val, unc = attr[:2]
-                unit = val.unit if hasattr(val, 'unit') else '--'
-                val = val.value if hasattr(val, 'unit') else val
-                unc = unc.value if hasattr(unc, 'unit') else unc
-                if val < 1E-3 or val > 1e5:
-                    val = float('{:.2e}'.format(val))
-                    if unc is None:
-                        unc = '--'
+                if ((param == 'logg') or (param =='mass') or (param == 'radius') or (param == 'Teff')):
+                    val = attr[0]
+                    lower_err = attr[1]
+                    upper_err = attr[2]
+                    unit = val.unit if hasattr(val, 'unit') else '--'
+                    val = val.value if hasattr(val, 'unit') else val
+                    unc_l = lower_err.value * -1 if hasattr(lower_err, 'unit') else lower_err * -1
+                    unc_u = upper_err.value if hasattr(upper_err, 'unit') else upper_err
+                    if val < 1E-3 or val > 1e5:
+                        if param == 'mass':
+                            val = float('{:.0f}'.format(val))
+                            unc =  str('{:+.0f}'.format(unc_l) + ' ' + '{:+.0f}'.format(unc_u))
+                        else:
+                            val = float('{:.2e}'.format(val))
+                            unc = str('{:-.2e}'.format(unc_l) + ' ' + '{:+.2e}'.format(unc_u))
                     else:
-                        unc = float('{:.2e}'.format(unc))
-                rows.append([param, val, unc, unit])
+                        unc = str('{:-.2f}'.format(unc_l) + ' ' + '{:+.2f}'.format(unc_u))
+                    rows.append([param, val, unc, unit])
+                else:
+                    val, unc = attr[:2]
+                    unit = val.unit if hasattr(val, 'unit') else '--'
+                    val = val.value if hasattr(val, 'unit') else val
+                    unc = unc.value if hasattr(unc, 'unit') else unc
+                    if val < 1E-3 or val > 1e5:
+                        val = float('{:.2e}'.format(val))
+                        if unc is None:
+                            unc = '--'
+                        else:
+                            unc = float('{:.2e}'.format(unc))
+                    rows.append([param, val, unc, unit])
 
             elif isinstance(attr, (str, float, bytes, int)):
                 rows.append([param, attr, '--', '--'])
@@ -2755,6 +2790,62 @@ def results(self):
         self.calculated = True
 
         return at.Table(np.asarray(rows), names=('param', 'value', 'unc', 'units'))
+    # def results(self):
+    #     """
+    #     A property for displaying the results
+    #     """
+    #     # Make the SED to get the most recent results
+    #     if not self.calculated:
+    #         self.make_sed()
+    #
+    #     # Get the params to display
+    #     params = copy(self.params)
+    #     print(params)
+    #     # Add best fits
+    #     for name, fit in self.best_fit.items():
+    #         params.append(name)
+    #
+    #     # Get the params
+    #     rows = []
+    #     for param in params:
+    #
+    #         # Get the values and format
+    #         attr = getattr(self, param, None)
+    #
+    #         if attr is None:
+    #             attr = '--'
+    #
+    #         if isinstance(attr, (tuple, list)):
+    #             val,lower_unc,upper_unc = attr[:3]
+    #             unit = val.unit if hasattr(val, 'unit') else '--'
+    #             val = val.value if hasattr(val, 'unit') else val
+    #             lower_unc = lower_unc.value if hasattr(lower_unc, 'unit') else lower_unc
+    #             upper_unc = upper_unc.value if hasattr(upper_unc, 'unit') else upper_unc
+    #             if val < 1E-3 or val > 1e5:
+    #                 val = float('{:.2e}'.format(val))
+    #                 if lower_unc is None:
+    #                     lower_unc = '--'
+    #                 else:
+    #                     lower_unc = float('{:.2e}'.format(lower_unc))
+    #                 if upper_unc is None:
+    #                     upper_unc = '--'
+    #                 else:
+    #                     upper_unc = float('{:.2e}'.format(upper_unc))
+    #             rows.append([param, val, lower_unc,upper_unc, unit])
+    #
+    #         elif isinstance(attr, (str, float, bytes, int)):
+    #             rows.append([param, attr, '--', '--','--'])
+    #
+    #         elif hasattr(attr, 'deg'):
+    #             rows.append([param, attr.deg, '--', '--','deg'])
+    #
+    #         else:
+    #             pass
+    #
+    #     # Set as calculated
+    #     self.calculated = True
+    #
+    #     return at.Table(np.asarray(rows), names=('param', 'value', 'lower_unc','upper_unc', 'units'))
 
     def run_methods(self, method_list):
         """
@@ -2983,6 +3074,7 @@ def Teff(self, teff):
         """
         self._validate_and_set_param('Teff', teff, q.K, True, vmin=0 * q.K, vmax=50000 * q.K)
 
+    # def _validate_and_set_param(self, param, values, units, set_uncalculated=True, trigger=[], vmin=None, vmax=None):
     def _validate_and_set_param(self, param, values, units, set_uncalculated=True, trigger=[], vmin=None, vmax=None):
         """
         Method to validate and set a calculated quantity
@@ -3014,7 +3106,6 @@ def _validate_and_set_param(self, param, values, units, set_uncalculated=True, t
             self.message("Setting {} to 'None'".format(param))
 
         else:
-
             # If the last value is string, it's the reference
             if isinstance(values[-1], str):
                 ref = values[-1]
@@ -3023,16 +3114,17 @@ def _validate_and_set_param(self, param, values, units, set_uncalculated=True, t
                 ref = None
 
             # Make sure it's a sequence
-            if not u.issequence(values, length=[2, 3]):
+            if not u.issequence(values, length=[2, 3, 4, 5]):
                 raise TypeError("{} must be a sequence of (value, error) or (value, lower_error, upper_error).".format(param))
 
             # Make sure it's in correct units
             if not all([u.equivalent(val, units) for val in values]):
+                print(values)
                 raise TypeError("{} values must be {}".format(param, 'unitless' if units is None else "astropy.units.quantity.Quantity of the appropriate units , e.g. '{}'".format(units)))
 
             # Ensure valid range but don't throw error
-            vmin = vmin or -np.inf * (units or 1)
-            vmax = vmax or np.inf * (units or 1)
+            vmin = vmin or (-np.inf * (units or 1))
+            vmax = vmax or (np.inf * (units or 1))
             if (values[0] < vmin) or (values[0] > vmax):
                 self.message("{}: {} value is not in valid range [{}, {}].".format(values, param, vmin, vmax))
 
diff --git a/sedkit/utilities.py b/sedkit/utilities.py
old mode 100755
new mode 100644
index 301fd7b5..f32cf033
--- a/sedkit/utilities.py
+++ b/sedkit/utilities.py
@@ -467,7 +467,8 @@ def flux2mag(flx, bandpass, photon=True):
     unit = flx.unit
 
     # Set uncertainty
-    unc = unc or np.nan * unit
+
+    # unc = (unc or (np.nan * unit))
 
     # Convert energy units to photon counts
     if photon:

From dd8dbadd0b5a83ecf8c8dad8df4ac62fcd81edf7 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <Sherelyn.Alejandro78@myhunter.cuny.edu>
Date: Tue, 28 Nov 2023 14:56:01 -0500
Subject: [PATCH 02/29] adding comment of absolute spectrum being at 10 pc in
 the header of the absolute sed file

---
 sedkit/sed.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sedkit/sed.py b/sedkit/sed.py
index d86dbe07..96d55f2c 100644
--- a/sedkit/sed.py
+++ b/sedkit/sed.py
@@ -1165,7 +1165,7 @@ def export(self, parentdir='.', dirname=None, zipped=False):
         # Absolute spectral SED
         if self.abs_spec_SED is not None:
             specpath = os.path.join(dirpath, '{}_absolute_SED.txt'.format(name))
-            header = '{} absolute spectrum (erg/s/cm2/A) as a function of wavelength (um)'.format(name)
+            header = '{} absolute spectrum (erg/s/cm2/A) at 10 pc as a function of wavelength (um)'.format(name)
             self.abs_spec_SED.export(specpath, header=header)
 
         # All photometry

From 7e2cfd142037d3626001c21a9f16353d4bcb918a Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <Sherelyn.Alejandro78@myhunter.cuny.edu>
Date: Mon, 19 Feb 2024 20:33:40 -0500
Subject: [PATCH 03/29] added chaubrier_2022 evolutionary models

---
 sedkit/isochrone.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
index 75e0fe55..44bffa80 100644
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -35,7 +35,7 @@
     EVO_MODELS = [os.path.basename(m).replace('.txt', '') for m in glob.glob(resource_filename('sedkit', 'data/models/evolutionary/*'))]
 # Fails RTD build for some reason
 except:
-    EVO_MODELS = ['COND03', 'dmestar_solar', 'DUSTY00', 'f2_solar_age', 'hybrid_solar_age', 'nc+03_age','nc+0.0_age', 'nc-03_age', 'nc_solar_age','nc+0.5_age','nc-0.5_age', 'parsec12_solar','ATMO_NEQ_strong','ATMO_NEQ_strong_MIRI']
+    EVO_MODELS = ['COND03', 'dmestar_solar', 'DUSTY00', 'f2_solar_age', 'hybrid_solar_age', 'nc+03_age','nc+0.0_age', 'nc-03_age', 'nc_solar_age','nc+0.5_age','nc-0.5_age', 'parsec12_solar','ATMO_NEQ_strong','ATMO_NEQ_strong_MIRI','chaubrier_2022.txt']
 
 class Isochrone:
     """A class to handle model isochrones"""

From 950fa6bf58480baed0e1bca84ed7af32f4503e41 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Mon, 15 Jul 2024 14:47:32 -0400
Subject: [PATCH 04/29] getting rid of comments

---
 sedkit/isochrone.py | 12 ++----------
 1 file changed, 2 insertions(+), 10 deletions(-)

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
index 44bffa80..4d773940 100644
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -186,16 +186,10 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
 
         # Get final values
         unit = self.data[yparam].unit or 1
-        # average = np.mean(np.array(values_list) * unit)
-        # # std = np.std(values_list)
-        # # error = std * unit
-        # min = np.min(np.array(values_list) * unit)
-        # max = np.max(np.array(values_list) * unit)
-        # error = abs(max-min)/2
 
         #Error Propagation
         dist_hist, hist_edges = np.histogram(values_list)
-
+        
         #Normalize Histogram
         dist_hist = dist_hist/ np.sum(dist_hist)
         # cumulative PDF
@@ -212,9 +206,7 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
         Q2_hist = np.interp(0.50, cum_PDF_hist, cum_PDF_hist_x )    # median
         Q3_hist = np.interp(error_upper_per, cum_PDF_hist, cum_PDF_hist_x)  # upper error
 
-        # print("This is the lower error",Q1_hist)
-        # print("This is the median", Q2_hist)
-        # print("This is the upper error", Q3_hist)
+
         lower_err = (Q2_hist - Q1_hist) * unit
         upper_err = (Q3_hist - Q2_hist) * unit
         average = Q2_hist * unit

From 69951398357bf7e9828f4633c1f8129a24410bd0 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Mon, 15 Jul 2024 14:49:06 -0400
Subject: [PATCH 05/29] increasing the number of samples

---
 sedkit/isochrone.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
index 4d773940..fd5b752e 100644
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -171,8 +171,8 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
         mu, sigma = xval[0], xval[1]  # mean and standard deviation for values on the x-axis
         mu_a, sigma_a = age[0].value, age[1].value  # mean and standard deviation for the age range provided
 
-        xsample = np.random.normal(mu, sigma, 1000)
-        ysample = np.random.uniform(mu_a-sigma_a, mu_a+sigma_a, 1000)
+        xsample = np.random.normal(mu, sigma, 5000)
+        ysample = np.random.uniform(mu_a-sigma_a, mu_a+sigma_a, 5000)
         values_list = []
         nan_counter = 0
 

From b4b376ca6201039bcb742654df9ef08415495ca6 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Mon, 15 Jul 2024 14:50:22 -0400
Subject: [PATCH 06/29] increasing the number of samples

---
 sedkit/isochrone.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
index fd5b752e..1bf3e8fe 100644
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -175,7 +175,7 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
         ysample = np.random.uniform(mu_a-sigma_a, mu_a+sigma_a, 5000)
         values_list = []
         nan_counter = 0
-
+        print('I AM USING NEW SAMPLE SIZE')
         for x, y in zip(xsample, ysample):
             result = self.interpolate(x, y * age[0].unit, xparam, yparam)
             if result is not None:

From 46cbcd9e700526bf6dd2d80e2b2f3d21a4a62f81 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Mon, 15 Jul 2024 14:52:09 -0400
Subject: [PATCH 07/29] getting rid of print statement

---
 sedkit/isochrone.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
index 1bf3e8fe..fd5b752e 100644
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -175,7 +175,7 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
         ysample = np.random.uniform(mu_a-sigma_a, mu_a+sigma_a, 5000)
         values_list = []
         nan_counter = 0
-        print('I AM USING NEW SAMPLE SIZE')
+
         for x, y in zip(xsample, ysample):
             result = self.interpolate(x, y * age[0].unit, xparam, yparam)
             if result is not None:

From 7de57bc9cbca5c8bb27244758b8b0c11a5aecb2c Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 25 Jul 2024 12:34:12 -0400
Subject: [PATCH 08/29] Increasing the number of iterations for monte-carlo
 approach

---
 sedkit/isochrone.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
index fd5b752e..fde13f2d 100644
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -171,8 +171,8 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
         mu, sigma = xval[0], xval[1]  # mean and standard deviation for values on the x-axis
         mu_a, sigma_a = age[0].value, age[1].value  # mean and standard deviation for the age range provided
 
-        xsample = np.random.normal(mu, sigma, 5000)
-        ysample = np.random.uniform(mu_a-sigma_a, mu_a+sigma_a, 5000)
+        xsample = np.random.normal(mu, sigma, 10000)
+        ysample = np.random.uniform(mu_a-sigma_a, mu_a+sigma_a, 10000)
         values_list = []
         nan_counter = 0
 

From 996a48235650066d36de014e1f16eca6a6fc106e Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 25 Jul 2024 12:39:40 -0400
Subject: [PATCH 09/29] Changing integrate function error functionality to take
 the median of the bootstrapped errors and not the max.Also, changing the
 number of samples in boostrap method to be 10^3.

---
 sedkit/spectrum.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/sedkit/spectrum.py b/sedkit/spectrum.py
index a4e900ad..5481699a 100755
--- a/sedkit/spectrum.py
+++ b/sedkit/spectrum.py
@@ -562,7 +562,7 @@ def flux_units(self, flux_units):
         self._flux_units = flux_units
         self._set_units()
 
-    def integrate(self, units=None, n_samples=10):
+    def integrate(self, units=None, n_samples=1000):
         """Calculate the area under the spectrum
 
         Parameters
@@ -619,7 +619,8 @@ def integrate(self, units=None, n_samples=10):
                 uvals.append(abs(err.value - val.value))
 
             # Get 1-sigma of distribution
-            vunc = np.max(abs(np.asarray(uvals))) * units
+            # vunc = np.max(abs(np.asarray(uvals))) * units
+            vunc = np.median(abs(np.asarray(uvals))) * units
 
         return val, vunc
 

From 82b7d1a85418b61a330b50159ece526e2adc14c1 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 01:13:54 -0400
Subject: [PATCH 10/29] changed the set Vega radius to be in the format of
 upper and lower error.

---
 sedkit/sed.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/sedkit/sed.py b/sedkit/sed.py
index c21ebb44..31842c5b 100644
--- a/sedkit/sed.py
+++ b/sedkit/sed.py
@@ -1943,8 +1943,8 @@ def infer_radius(self, radius_units=q.Rsun, infer_from='Lbol', plot=False):
                 radius = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'radius', plot=plot)
 
             # Store the value
-            self.radius = [radius[0].round(3), radius[1].round(3),radius[2].round(3),radius[3]] if radius is not None else radius
-
+            self.radius = [radius[0].round(3), radius[1].round(3), radius[2].round(3), radius[3]] if radius is not None else radius
+            print(self.radius)
         # Try radius(spt) relation
         elif infer_from == 'spt':
 
@@ -2046,7 +2046,7 @@ def infer_Teff(self, teff_units=q.K, plot=False):
                 teff = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'teff', plot=plot)
 
             # Store the value
-            self.Teff = [teff[0].round(0), teff[1].round(0), teff[2].round(0),teff[3]] if teff is not None else teff
+            self.Teff = [teff[0].round(0), teff[1].round(0), teff[2].round(0), teff[3]] if teff is not None else teff
 
         else:
 
@@ -3131,8 +3131,10 @@ def _validate_and_set_param(self, param, values, units, set_uncalculated=True, t
                 raise TypeError("{} values must be {}".format(param, 'unitless' if units is None else "astropy.units.quantity.Quantity of the appropriate units , e.g. '{}'".format(units)))
 
             # Ensure valid range but don't throw error
-            vmin = vmin or (-np.inf * (units or 1))
-            vmax = vmax or (np.inf * (units or 1))
+
+            vmin = vmin if vmin is not None else -np.inf * (1 if units is None else units)
+            vmax = vmax if vmax is not None else np.inf * (1 if units is None else units)
+
 
             if (values[0] < vmin) or (values[0] > vmax):
                 self.message("{}: {} value is not in valid range [{}, {}].".format(values, param, vmin, vmax))
@@ -3211,7 +3213,7 @@ def __init__(self, **kwargs):
         self.find_SDSS()
         self.find_2MASS()
         self.find_WISE()
-        self.radius = 2.818 * q.Rsun, 0.008 * q.Rsun, '2010ApJ...708...71Y'
+        self.radius = 2.818 * q.Rsun, 0.008 * q.Rsun, 0.008 * q.Rsun, '2010ApJ...708...71Y'
         self.age = 455 * q.Myr, 13 * q.Myr, '2010ApJ...708...71Y'
         self.logg = 4.1, 0.1, '2006ApJ...645..664A'
 

From 3ba7a1b670ace751a3a9e283b6f8da750aa3d3dd Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 01:42:41 -0400
Subject: [PATCH 11/29] With the change in uncertainty calculation, we still
 get a result that is not zero regardless of giving ucertainties for xparam or
 yparam.

---
 sedkit/tests/test_isochrone.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/sedkit/tests/test_isochrone.py b/sedkit/tests/test_isochrone.py
index 2b9cd3ee..5a51dcee 100644
--- a/sedkit/tests/test_isochrone.py
+++ b/sedkit/tests/test_isochrone.py
@@ -27,9 +27,9 @@ def test_evaluate(self):
         result = self.hsa.evaluate((-4, 0.1), 4*q.Gyr, 'Lbol', 'mass')
         self.assertTrue(isinstance(result, tuple))
 
-        # No xparam or yparam uncertainties
+        # No xparam and yparam uncertainties
         result = self.hsa.evaluate(-4, 4*q.Gyr, 'Lbol', 'mass')
-        self.assertTrue(isinstance(result, tuple) and result[1] == 0)
+        self.assertTrue(isinstance(result, tuple))
 
     def test_interp(self):
         """Test that the model isochrone can be interpolated"""

From 1e75b16a08eff19eb7449f9a2d54f48abc2bf6ba Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 01:43:04 -0400
Subject: [PATCH 12/29] getting rid of random bool statements

---
 sedkit/isochrone.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
index 79e2eec5..9c10a16a 100644
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -155,10 +155,8 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
         age = age[0], age[1].to(age[0].unit)
 
         # Check if the xval has an uncertainty
-        no_xerr = False
         if not isinstance(xval, (tuple, list)):
             xval = (xval, xval * 0)
-            no_xerr = True
 
         # Test the age range is inbounds
         if age[0] < self.ages.min() or age[0] > self.ages.max():

From e2998197e3208661052c26a383b536c1653882a2 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 03:00:30 -0400
Subject: [PATCH 13/29] getting rid of random print statements

---
 sedkit/spectrum.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/sedkit/spectrum.py b/sedkit/spectrum.py
index b09e6e4a..daa3241e 100644
--- a/sedkit/spectrum.py
+++ b/sedkit/spectrum.py
@@ -1058,9 +1058,9 @@ def resamp(self, wave=None, resolution=None, name=None):
         wave = wave.value
 
         # Bin the spectrum
-        print(wave, self.wave)
+
         binned = u.spectres(wave, self.wave, self.flux, self.unc)
-        print(binned[0])
+
 
         # Update the spectrum
         spectrum = [i * Q for i, Q in zip(binned, self.units)]

From e89ab4de0bc10fe3ec61052c266abb10fa589ce2 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 03:20:04 -0400
Subject: [PATCH 14/29] Changed the plotting to keep in mind whether a distance
 is provided or not. Also, got rid of a few unecessary print statements.

---
 sedkit/sed.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/sedkit/sed.py b/sedkit/sed.py
index 31842c5b..b2111222 100644
--- a/sedkit/sed.py
+++ b/sedkit/sed.py
@@ -1944,7 +1944,7 @@ def infer_radius(self, radius_units=q.Rsun, infer_from='Lbol', plot=False):
 
             # Store the value
             self.radius = [radius[0].round(3), radius[1].round(3), radius[2].round(3), radius[3]] if radius is not None else radius
-            print(self.radius)
+
         # Try radius(spt) relation
         elif infer_from == 'spt':
 
@@ -2546,7 +2546,10 @@ def plot(self, app=False, photometry=True, spectra=True, integral=True, syntheti
             self.make_sed()
 
         # Distinguish between apparent and absolute magnitude
-        pre = 'app_' if app else 'abs_'
+        if self.distance is None:
+            pre = 'app_'
+        else:
+            pre = 'app_' if app else 'abs_'
 
         # Calculate reasonable axis limits
         full_SED = getattr(self, pre + 'SED')
@@ -2601,13 +2604,11 @@ def plot(self, app=False, photometry=True, spectra=True, integral=True, syntheti
         # Plot spectra
         if spectra and len(self.spectra) > 0:
 
-
             if (spectra == 'all' and app is False):
                 for n, spec in enumerate(self.spectra['spectrum']):
                     self.fig = spec.plot(fig=self.fig, components=True, const=(spec.flux_calibrate(self.distance).flux/spec.flux))
             elif (spectra == 'all' and app is True):
                 for n, spec in enumerate(self.spectra['spectrum']):
-                    print(const)
                     self.fig = spec.plot(fig=self.fig, components=True, const=const)
             else:
                 self.fig.line(spec_SED.wave, spec_SED.flux * const, color=color, alpha=0.8, name='spectrum', legend_label='Spectrum')
@@ -3127,7 +3128,6 @@ def _validate_and_set_param(self, param, values, units, set_uncalculated=True, t
 
             # Make sure it's in correct units
             if not all([u.equivalent(val, units) for val in values]):
-                print(values)
                 raise TypeError("{} values must be {}".format(param, 'unitless' if units is None else "astropy.units.quantity.Quantity of the appropriate units , e.g. '{}'".format(units)))
 
             # Ensure valid range but don't throw error
@@ -3215,7 +3215,7 @@ def __init__(self, **kwargs):
         self.find_WISE()
         self.radius = 2.818 * q.Rsun, 0.008 * q.Rsun, 0.008 * q.Rsun, '2010ApJ...708...71Y'
         self.age = 455 * q.Myr, 13 * q.Myr, '2010ApJ...708...71Y'
-        self.logg = 4.1, 0.1, '2006ApJ...645..664A'
+        self.logg = 4.1, 0.1, 0.1, '2006ApJ...645..664A'
 
         # Get the spectrum
         self.add_spectrum(sp.Vega(snr=100))

From 5550483fbb0da291bceb2efac2fe10de1639cc72 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 03:20:47 -0400
Subject: [PATCH 15/29] Changed the formatting of the unceratainty to have and
 upper and a lower unc.

---
 sedkit/tests/test_sed.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/sedkit/tests/test_sed.py b/sedkit/tests/test_sed.py
index de51af49..c08fddeb 100644
--- a/sedkit/tests/test_sed.py
+++ b/sedkit/tests/test_sed.py
@@ -148,7 +148,7 @@ def test_no_spectra(self):
         """Test that a purely photometric SED can be creted"""
         s = copy.copy(self.sed)
         s.age = 455*q.Myr, 13*q.Myr
-        s.radius = 2.362*q.Rsun, 0.02*q.Rjup
+        s.radius = 2.362*q.Rsun, 0.02*q.Rjup, 0.02*q.Rjup
         s.parallax = 130.23*q.mas, 0.36*q.mas
         s.spectral_type = 'A0V'
         s.add_photometry('2MASS.J', -0.177, 0.206)
@@ -185,10 +185,10 @@ def test_plot(self):
         fig = v.plot(integral=True, synthetic_photometry=True, best_fit=True)
 
     def test_no_photometry(self):
-        """Test that a purely photometric SED can be creted"""
+        """Test that a purely photometric SED can be created"""
         s = copy.copy(self.sed)
         s.age = 455*q.Myr, 13*q.Myr
-        s.radius = 2.362*q.Rsun, 0.02*q.Rjup
+        s.radius = 2.362*q.Rsun, 0.02*q.Rjup,0.02*q.Rjup
         s.parallax = 130.23*q.mas, 0.36*q.mas
         s.spectral_type = 'A0V'
         s.add_spectrum(self.spec1)

From 08c8adf3663d47d061508a9c29c842d1c316a420 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 11:18:28 -0400
Subject: [PATCH 16/29] It seems I already have them fixed here!

---
 sedkit/sed.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/sedkit/sed.py b/sedkit/sed.py
index b2111222..41048303 100644
--- a/sedkit/sed.py
+++ b/sedkit/sed.py
@@ -762,7 +762,6 @@ def _calibrate_photometry(self, name='photometry'):
                   # Calculate abs_flux values
                   for n, row in enumerate(table):  
                       abs_flux, abs_flux_unc = u.mag2flux(row['bandpass'], row['abs_magnitude'], sig_m=row['abs_magnitude_unc'], units=self.flux_units)
-
                       table['abs_flux'][n] = abs_flux
                       table['abs_flux_unc'][n] = abs_flux_unc
                        

From f4d5e24daf68c7edd9df2c35af31d0f2bbbaa3a2 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 12:40:45 -0400
Subject: [PATCH 17/29] deleting comment

Co-authored-by: Kelle Cruz <kellecruz@gmail.com>
---
 sedkit/spectrum.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/sedkit/spectrum.py b/sedkit/spectrum.py
index daa3241e..6b0421e7 100644
--- a/sedkit/spectrum.py
+++ b/sedkit/spectrum.py
@@ -619,7 +619,6 @@ def integrate(self, units=None, n_samples=1000):
                 uvals.append(abs(err.value - val.value))
 
             # Get 1-sigma of distribution
-            # vunc = np.max(abs(np.asarray(uvals))) * units
             vunc = np.median(abs(np.asarray(uvals))) * units
 
         return val, vunc

From 1a596a72d7ca0ffac89e40bb9e4015fc95c12170 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 12:41:06 -0400
Subject: [PATCH 18/29] deleting commented function

Co-authored-by: Kelle Cruz <kellecruz@gmail.com>
---
 sedkit/sed.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/sedkit/sed.py b/sedkit/sed.py
index 41048303..9798772a 100644
--- a/sedkit/sed.py
+++ b/sedkit/sed.py
@@ -3082,7 +3082,6 @@ def Teff(self, teff):
         """
         self._validate_and_set_param('Teff', teff, q.K, True, vmin=0 * q.K, vmax=50000 * q.K)
 
-    # def _validate_and_set_param(self, param, values, units, set_uncalculated=True, trigger=[], vmin=None, vmax=None):
     def _validate_and_set_param(self, param, values, units, set_uncalculated=True, trigger=[], vmin=None, vmax=None):
         """
         Method to validate and set a calculated quantity

From 96e041ed99338d2937303fb0eec2e94f4137137c Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 12:41:43 -0400
Subject: [PATCH 19/29] deleting uncertainty comment

Co-authored-by: Kelle Cruz <kellecruz@gmail.com>
---
 sedkit/utilities.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/sedkit/utilities.py b/sedkit/utilities.py
index 8c94d0e8..a3e181af 100644
--- a/sedkit/utilities.py
+++ b/sedkit/utilities.py
@@ -467,8 +467,6 @@ def flux2mag(flx, bandpass, photon=True):
     unit = flx.unit
 
     # Set uncertainty
-    # unc = (unc or (np.nan * unit))
-
 
     # Convert energy units to photon counts
     if photon:

From d0c6dca973117400f3b1ad1f0edff52cc73e9356 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 12:42:14 -0400
Subject: [PATCH 20/29] deleting commented color mapper

Co-authored-by: Kelle Cruz <kellecruz@gmail.com>
---
 sedkit/isochrone.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
index 9c10a16a..67241b2c 100644
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -330,9 +330,7 @@ def plot(self, xparam, yparam, draw=False, labels=True, **kwargs):
         color_mapper = LinearColorMapper(palette='Viridis256',
                                          low=self.ages.min().value,
                                          high=self.ages.max().value)
-        # color_mapper = LinearColorMapper(palette='Viridis256',
-        #                                  low=.05,
-        #                                  high=.15)
+     
         color_bar = ColorBar(color_mapper=color_mapper, ticker=BasicTicker(),
                              label_standoff=5, border_line_color=None,
                              title='Age [{}]'.format(self.age_units),

From 86726d0fa52a3b372345268bb127c9880c7c2631 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 12:45:24 -0400
Subject: [PATCH 21/29] deleting commented result function

---
 sedkit/sed.py | 56 ---------------------------------------------------
 1 file changed, 56 deletions(-)

diff --git a/sedkit/sed.py b/sedkit/sed.py
index 9798772a..b6c33d3b 100644
--- a/sedkit/sed.py
+++ b/sedkit/sed.py
@@ -2798,62 +2798,6 @@ def results(self):
         self.calculated = True
 
         return at.Table(np.asarray(rows), names=('param', 'value', 'unc', 'units'))
-    # def results(self):
-    #     """
-    #     A property for displaying the results
-    #     """
-    #     # Make the SED to get the most recent results
-    #     if not self.calculated:
-    #         self.make_sed()
-    #
-    #     # Get the params to display
-    #     params = copy(self.params)
-    #     print(params)
-    #     # Add best fits
-    #     for name, fit in self.best_fit.items():
-    #         params.append(name)
-    #
-    #     # Get the params
-    #     rows = []
-    #     for param in params:
-    #
-    #         # Get the values and format
-    #         attr = getattr(self, param, None)
-    #
-    #         if attr is None:
-    #             attr = '--'
-    #
-    #         if isinstance(attr, (tuple, list)):
-    #             val,lower_unc,upper_unc = attr[:3]
-    #             unit = val.unit if hasattr(val, 'unit') else '--'
-    #             val = val.value if hasattr(val, 'unit') else val
-    #             lower_unc = lower_unc.value if hasattr(lower_unc, 'unit') else lower_unc
-    #             upper_unc = upper_unc.value if hasattr(upper_unc, 'unit') else upper_unc
-    #             if val < 1E-3 or val > 1e5:
-    #                 val = float('{:.2e}'.format(val))
-    #                 if lower_unc is None:
-    #                     lower_unc = '--'
-    #                 else:
-    #                     lower_unc = float('{:.2e}'.format(lower_unc))
-    #                 if upper_unc is None:
-    #                     upper_unc = '--'
-    #                 else:
-    #                     upper_unc = float('{:.2e}'.format(upper_unc))
-    #             rows.append([param, val, lower_unc,upper_unc, unit])
-    #
-    #         elif isinstance(attr, (str, float, bytes, int)):
-    #             rows.append([param, attr, '--', '--','--'])
-    #
-    #         elif hasattr(attr, 'deg'):
-    #             rows.append([param, attr.deg, '--', '--','deg'])
-    #
-    #         else:
-    #             pass
-    #
-    #     # Set as calculated
-    #     self.calculated = True
-    #
-    #     return at.Table(np.asarray(rows), names=('param', 'value', 'lower_unc','upper_unc', 'units'))
 
     def run_methods(self, method_list):
         """

From ccc0eae5a10f50a9a1149befffa1acdab8ca54b6 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Thu, 31 Oct 2024 15:47:35 -0400
Subject: [PATCH 22/29] changed test evaluate to use mark parametrize

---
 sedkit/tests/test_isochrone.py | 40 ++++++++++++++++++----------------
 1 file changed, 21 insertions(+), 19 deletions(-)

diff --git a/sedkit/tests/test_isochrone.py b/sedkit/tests/test_isochrone.py
index 5a51dcee..d2e4b222 100644
--- a/sedkit/tests/test_isochrone.py
+++ b/sedkit/tests/test_isochrone.py
@@ -1,36 +1,38 @@
 """Series of unit tests for the isochrone.py module"""
 import unittest
-
+import pytest
 import astropy.units as q
 
 from .. import isochrone as iso
 from .. import utilities as u
 
 
+@pytest.mark.parametrize('xval,age,xparam,yparam', [
+    ((-4, 0.1), (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass'),  # With uncertainties
+    (-4, (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass'),  # No xparam uncertainty
+    ((-4, 0.1), 4 * q.Gyr, 'Lbol', 'mass'),  # No yparam uncertainty
+    (-4, 4 * q.Gyr, 'Lbol', 'mass')  # No xparam and yparam uncertainties
+])
+def test_evaluate( xval, age, xparam, yparam):
+    # average, lower, upper
+    """Test the evaluate method"""
+    hsa = iso.Isochrone('hybrid_solar_age')
+    result = hsa.evaluate(xval, age, xparam, yparam)
+    assert (isinstance(result, tuple)) is True
+    # assert (np.isclose(result[0],value,))
+    # try == first but if it can't happen then use isclose
+    # test the value of the y param uncertainties (the three values - lower,average and upper)
+    # test lbol to radius
+    # test lbol to logg
+    # test for different lbols and age/age ranges (cold, mid-temp, warm object)
+    # Add the three results for the three different values
+
 class TestIsochrone(unittest.TestCase):
     """Tests for the hybrid_solar_age model isochrones"""
     def setUp(self):
         # Make Spectrum class for testing
         self.hsa = iso.Isochrone('hybrid_solar_age')
 
-    def test_evaluate(self):
-        """Test the evaluate method"""
-        # With uncertainties
-        result = self.hsa.evaluate((-4, 0.1), (4*q.Gyr, 0.1*q.Gyr), 'Lbol', 'mass')
-        self.assertTrue(isinstance(result, tuple))
-
-        # No xparam uncertainty
-        result = self.hsa.evaluate(-4, (4*q.Gyr, 0.1*q.Gyr), 'Lbol', 'mass')
-        self.assertTrue(isinstance(result, tuple))
-
-        # No yparam uncertainty
-        result = self.hsa.evaluate((-4, 0.1), 4*q.Gyr, 'Lbol', 'mass')
-        self.assertTrue(isinstance(result, tuple))
-
-        # No xparam and yparam uncertainties
-        result = self.hsa.evaluate(-4, 4*q.Gyr, 'Lbol', 'mass')
-        self.assertTrue(isinstance(result, tuple))
-
     def test_interp(self):
         """Test that the model isochrone can be interpolated"""
         # Successful interpolation

From f376a261c8ef19a524c0010bd3e666107fb493c9 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Mon, 4 Nov 2024 17:00:45 -0500
Subject: [PATCH 23/29] Adding More Tests

Co-authored-by: Kelle Cruz <kellecruz@gmail.com>
---
 sedkit/tests/test_isochrone.py | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/sedkit/tests/test_isochrone.py b/sedkit/tests/test_isochrone.py
index d2e4b222..039c0899 100644
--- a/sedkit/tests/test_isochrone.py
+++ b/sedkit/tests/test_isochrone.py
@@ -7,18 +7,19 @@
 from .. import utilities as u
 
 
-@pytest.mark.parametrize('xval,age,xparam,yparam', [
-    ((-4, 0.1), (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass'),  # With uncertainties
-    (-4, (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass'),  # No xparam uncertainty
-    ((-4, 0.1), 4 * q.Gyr, 'Lbol', 'mass'),  # No yparam uncertainty
-    (-4, 4 * q.Gyr, 'Lbol', 'mass')  # No xparam and yparam uncertainties
+@pytest.mark.parametrize('xval,age,xparam,yparam,expected_result', [
+    ((-4, 0.1), (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass',0.072),  # With uncertainties
+    (-4, (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass',0.072),  # No xparam uncertainty
+    ((-4, 0.1), 4 * q.Gyr, 'Lbol', 'mass',0.072),  # No yparam uncertainty
+    (-4, 4 * q.Gyr, 'Lbol', 'mass',0.020)  # No xparam and yparam uncertainties
 ])
-def test_evaluate( xval, age, xparam, yparam):
+def test_evaluate( xval, age, xparam, yparam,expected_result):
     # average, lower, upper
     """Test the evaluate method"""
     hsa = iso.Isochrone('hybrid_solar_age')
     result = hsa.evaluate(xval, age, xparam, yparam)
     assert (isinstance(result, tuple)) is True
+    assert (np.isclose(result[0].value, expected_result, atol=0.005))
     # assert (np.isclose(result[0],value,))
     # try == first but if it can't happen then use isclose
     # test the value of the y param uncertainties (the three values - lower,average and upper)

From 72cc44266ef583d47f7da283cab4537a318bcc64 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Mon, 4 Nov 2024 17:13:16 -0500
Subject: [PATCH 24/29] added numpy import

---
 sedkit/tests/test_isochrone.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/sedkit/tests/test_isochrone.py b/sedkit/tests/test_isochrone.py
index 039c0899..fad39f4a 100644
--- a/sedkit/tests/test_isochrone.py
+++ b/sedkit/tests/test_isochrone.py
@@ -2,18 +2,19 @@
 import unittest
 import pytest
 import astropy.units as q
+import numpy as np
 
 from .. import isochrone as iso
 from .. import utilities as u
 
 
 @pytest.mark.parametrize('xval,age,xparam,yparam,expected_result', [
-    ((-4, 0.1), (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass',0.072),  # With uncertainties
-    (-4, (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass',0.072),  # No xparam uncertainty
-    ((-4, 0.1), 4 * q.Gyr, 'Lbol', 'mass',0.072),  # No yparam uncertainty
-    (-4, 4 * q.Gyr, 'Lbol', 'mass',0.020)  # No xparam and yparam uncertainties
+    ((-4, 0.1), (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass', 0.072),  # With uncertainties
+    (-4, (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass', 0.072),  # No xparam uncertainty
+    ((-4, 0.1), 4 * q.Gyr, 'Lbol', 'mass', 0.072),  # No yparam uncertainty
+    (-4, 4 * q.Gyr, 'Lbol', 'mass', 0.020)  # No xparam and yparam uncertainties
 ])
-def test_evaluate( xval, age, xparam, yparam,expected_result):
+def test_evaluate(xval, age, xparam, yparam, expected_result):
     # average, lower, upper
     """Test the evaluate method"""
     hsa = iso.Isochrone('hybrid_solar_age')

From 1730b4294ee5b1c04df51aed14cc06e43c06dc69 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Mon, 4 Nov 2024 17:25:55 -0500
Subject: [PATCH 25/29] removing deprecated pkg file call

---
 sedkit/tests/test_utilities.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/sedkit/tests/test_utilities.py b/sedkit/tests/test_utilities.py
index 2f438e6d..3a54072b 100644
--- a/sedkit/tests/test_utilities.py
+++ b/sedkit/tests/test_utilities.py
@@ -1,6 +1,6 @@
 """A suite of tests for the utilities.py module"""
 import copy
-from pkg_resources import resource_filename
+import importlib_resources
 import pytest
 import unittest
 
@@ -294,7 +294,7 @@ def test_group_spectra():
 def test_spectrum_from_fits():
     """Test spectrum_from_fits function"""
     # Get the file
-    f = resource_filename('sedkit', '/data/Gl752B_NIR.fits')
+    f = str(importlib_resources.files('sedkit')/'/data/Gl752B_NIR.fits')
 
     # Get the spectrum
     spec = u.spectrum_from_fits(f)

From 037a91dfa479b0fe042b9bda7ecc854f5beb3e3b Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Mon, 4 Nov 2024 19:05:02 -0500
Subject: [PATCH 26/29] changes the lower and upper errors to not use a range
 and directly use the interpolated lower and upper errors

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

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
index 67241b2c..709ce86e 100644
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -203,9 +203,9 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
         Q2_hist = np.interp(0.50, cum_PDF_hist, cum_PDF_hist_x )    # median
         Q3_hist = np.interp(error_upper_per, cum_PDF_hist, cum_PDF_hist_x)  # upper error
 
-
-        lower_err = (Q2_hist - Q1_hist) * unit
-        upper_err = (Q3_hist - Q2_hist) * unit
+        # assigning units to values
+        lower_err = Q1_hist * unit
+        upper_err = Q3_hist * unit
         average = Q2_hist * unit
 
 

From e1477bc8e9b8bcc2128a07a8f9b0dd4fcb2fdce0 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Mon, 4 Nov 2024 19:13:52 -0500
Subject: [PATCH 27/29] Revert "changes the lower and upper errors to not use a
 range and directly use the interpolated lower and upper errors"

This reverts commit 037a91dfa479b0fe042b9bda7ecc854f5beb3e3b.
---
 sedkit/isochrone.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/sedkit/isochrone.py b/sedkit/isochrone.py
index 709ce86e..67241b2c 100644
--- a/sedkit/isochrone.py
+++ b/sedkit/isochrone.py
@@ -203,9 +203,9 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
         Q2_hist = np.interp(0.50, cum_PDF_hist, cum_PDF_hist_x )    # median
         Q3_hist = np.interp(error_upper_per, cum_PDF_hist, cum_PDF_hist_x)  # upper error
 
-        # assigning units to values
-        lower_err = Q1_hist * unit
-        upper_err = Q3_hist * unit
+
+        lower_err = (Q2_hist - Q1_hist) * unit
+        upper_err = (Q3_hist - Q2_hist) * unit
         average = Q2_hist * unit
 
 

From 4ffe0b8d2b0d71631e8abcc30830ab2605fbf780 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Mon, 4 Nov 2024 20:01:46 -0500
Subject: [PATCH 28/29] Added more tests for mass,radius and logg

---
 sedkit/tests/test_isochrone.py | 41 ++++++++++++++++++++++++----------
 1 file changed, 29 insertions(+), 12 deletions(-)

diff --git a/sedkit/tests/test_isochrone.py b/sedkit/tests/test_isochrone.py
index fad39f4a..c1aa1a75 100644
--- a/sedkit/tests/test_isochrone.py
+++ b/sedkit/tests/test_isochrone.py
@@ -8,24 +8,41 @@
 from .. import utilities as u
 
 
-@pytest.mark.parametrize('xval,age,xparam,yparam,expected_result', [
-    ((-4, 0.1), (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass', 0.072),  # With uncertainties
-    (-4, (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass', 0.072),  # No xparam uncertainty
-    ((-4, 0.1), 4 * q.Gyr, 'Lbol', 'mass', 0.072),  # No yparam uncertainty
-    (-4, 4 * q.Gyr, 'Lbol', 'mass', 0.020)  # No xparam and yparam uncertainties
+@pytest.mark.parametrize('xval,age,xparam,yparam,expected_result,expected_result_low,expected_result_up', [
+    # Mass
+    ((-4, 0.1), (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass', 0.072, 0.072, 0.072),  # With uncertainties
+    (-4, (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass', 0.072, 0.072, 0.072),  # No xparam uncertainty
+    ((-4, 0.1), 4 * q.Gyr, 'Lbol', 'mass', 0.072, 0.072, 0.072),  # No yparam uncertainty
+    (-4, 4 * q.Gyr, 'Lbol', 'mass', 0.020, 0, 0.058),  # No xparam and yparam uncertainties
+    # Radius
+    ((-4, 0.1), (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'radius', 0.095, 0.095, 0.095),  # With uncertainties
+    (-4, (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'radius', 0.095, 0.095, 0.095),  # No xparam uncertainty
+    ((-4, 0.1), 4 * q.Gyr, 'Lbol', 'radius', 0.095, 0.095, 0.095),  # No yparam uncertainty
+    (-4, 4 * q.Gyr, 'Lbol', 'radius', 0.045, 0.01, 0.080),  # No xparam and yparam uncertainties
+    # Logg
+    ((-4, 0.1), (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'logg', 5.345, 5.34, 5.35),  # With uncertainties
+    (-4, (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'logg', 5.345, 5.34, 5.35),  # No xparam uncertainty
+    ((-4, 0.1), 4 * q.Gyr, 'Lbol', 'logg', 5.345, 5.34, 5.35),  # No yparam uncertainty
+    (-4, 4 * q.Gyr, 'Lbol', 'logg', 5.395, 5.36, 5.43)  # No xparam and yparam uncertainties
 ])
-def test_evaluate(xval, age, xparam, yparam, expected_result):
+def test_evaluate(xval, age, xparam, yparam, expected_result, expected_result_low, expected_result_up):
     # average, lower, upper
     """Test the evaluate method"""
     hsa = iso.Isochrone('hybrid_solar_age')
     result = hsa.evaluate(xval, age, xparam, yparam)
+    average = result[0]  # Average yparam value
+    lower = result[0] - result[1]   # Lower yparam value
+    upper = result[0] + result[2]   # Upper yparam value
     assert (isinstance(result, tuple)) is True
-    assert (np.isclose(result[0].value, expected_result, atol=0.005))
-    # assert (np.isclose(result[0],value,))
-    # try == first but if it can't happen then use isclose
-    # test the value of the y param uncertainties (the three values - lower,average and upper)
-    # test lbol to radius
-    # test lbol to logg
+    if yparam == 'logg':
+        assert (np.isclose(average, expected_result, atol=0.005))
+        assert (np.isclose(lower, expected_result_low, atol=0.01))
+        assert (np.isclose(upper, expected_result_up, atol=0.01))
+    else:
+        assert (np.isclose(average.value, expected_result, atol=0.005))
+        assert (np.isclose(lower.value, expected_result_low, atol=0.01))
+        assert (np.isclose(upper.value, expected_result_up, atol=0.01))
+
     # test for different lbols and age/age ranges (cold, mid-temp, warm object)
     # Add the three results for the three different values
 

From 588d0affcc1e582f3e630373340fa570b126d821 Mon Sep 17 00:00:00 2001
From: Sherelyn Alejandro <80668767+SherelynA@users.noreply.github.com>
Date: Wed, 6 Nov 2024 22:36:07 -0500
Subject: [PATCH 29/29] Added test for young and old ages.

---
 sedkit/tests/test_isochrone.py | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/sedkit/tests/test_isochrone.py b/sedkit/tests/test_isochrone.py
index c1aa1a75..aec4c105 100644
--- a/sedkit/tests/test_isochrone.py
+++ b/sedkit/tests/test_isochrone.py
@@ -23,10 +23,17 @@
     ((-4, 0.1), (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'logg', 5.345, 5.34, 5.35),  # With uncertainties
     (-4, (4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'logg', 5.345, 5.34, 5.35),  # No xparam uncertainty
     ((-4, 0.1), 4 * q.Gyr, 'Lbol', 'logg', 5.345, 5.34, 5.35),  # No yparam uncertainty
-    (-4, 4 * q.Gyr, 'Lbol', 'logg', 5.395, 5.36, 5.43)  # No xparam and yparam uncertainties
+    (-4, 4 * q.Gyr, 'Lbol', 'logg', 5.395, 5.36, 5.43),  # No xparam and yparam uncertaintiesd
+    # Young age
+    ((-4, 0.1), (0.4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass', 0.0515, 0.045, 0.055),   # mass with uncertainties
+    ((-4, 0.1), (0.4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'logg', 5.08, 5.01, 5.15),   # logg with uncertainties
+    ((-4, 0.1), (0.4 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'radius', 0.105, 0.10, 0.11),  # radius with uncertainties
+    # Old age
+    ((-4, 0.1), (9 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'mass', 0.074, 0.070, 0.080),  # mass with uncertainties
+    ((-4, 0.1), (9 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'logg', 5.345, 5.34, 5.35),  # logg with uncertainties
+    ((-4, 0.1), (9 * q.Gyr, 0.1 * q.Gyr), 'Lbol', 'radius', 0.095, 0.09, 0.10)  # radius with uncertainties
 ])
 def test_evaluate(xval, age, xparam, yparam, expected_result, expected_result_low, expected_result_up):
-    # average, lower, upper
     """Test the evaluate method"""
     hsa = iso.Isochrone('hybrid_solar_age')
     result = hsa.evaluate(xval, age, xparam, yparam)
@@ -43,8 +50,6 @@ def test_evaluate(xval, age, xparam, yparam, expected_result, expected_result_lo
         assert (np.isclose(lower.value, expected_result_low, atol=0.01))
         assert (np.isclose(upper.value, expected_result_up, atol=0.01))
 
-    # test for different lbols and age/age ranges (cold, mid-temp, warm object)
-    # Add the three results for the three different values
 
 class TestIsochrone(unittest.TestCase):
     """Tests for the hybrid_solar_age model isochrones"""