Fixing Infer Radius Method

sedkit/isochrone.py

@@ -121,7 +121,7 @@ def age_units(self, unit):
         # Get the min and max ages
         self.ages = np.array(np.unique(self.data['age']))*self.age_units
 
-    def evaluate(self, xval, age, xparam, yparam, plot=False):
+    def evaluate(self, xval, age, xparam, yparam, distribution, plot=False):
         """Interpolate the value and uncertainty of *yparam* given an
         x-value and age range
 
@@ -135,6 +135,8 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
             The name of the parameter on the x-axis
         yparam: str
             The name of the parameter on the y-axis
+        distribution: str
+            The type of distribution used for the Monte Carlo approach, either 'uniform' or 'normal'
         plot: bool
             Plot all isochrones and the interpolated value
 
@@ -168,8 +170,17 @@ 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
 
+        # Generate random sample for values on the x-axis
         xsample = np.random.normal(mu, sigma, 10000)
-        ysample = np.random.uniform(mu_a-sigma_a, mu_a+sigma_a, 10000)
+
+        # Generate random sample for the age range provided
+        if distribution == 'uniform':
+            ysample = np.random.uniform(mu_a-sigma_a, mu_a+sigma_a, 10000)
+        elif distribution == 'normal':
+            ysample = np.random.normal(mu_a, sigma_a, 10000)
+        else:
+            raise ValueError("Distribution must be set to 'uniform' or 'normal'.")
+
         values_list = []
         nan_counter = 0
 

sedkit/relations.py

@@ -236,9 +236,9 @@ def evaluate(self, rel_name, x_val, xunits=None, yunits=None, fit_local=False, p
 
                 # Get the relation
                 full_rel = self.relations[rel_name]
-                out_xunits = full_rel['xunit'].to(xunits) * xunits if xunits is not None else full_rel['xunit'] or 1
-                out_yunits = full_rel['yunit'].to(yunits) * yunits if yunits is not None else full_rel['yunit'] or 1
 
+                out_xunits = '{:.{p}}'.format(full_rel['xunit'].to(xunits), p=3) * xunits if xunits is not None else full_rel['xunit'] or 1
+                out_yunits = '{:.{p}}'.format(full_rel['yunit'].to(yunits), p=3) * yunits if yunits is not None else full_rel['yunit'] or 1
                 # Use local points for relation
                 if isinstance(fit_local, int) and fit_local is not False:
 
@@ -479,7 +479,7 @@ def get_radius(self, spt, plot=False):
             fig.triangle([spt], [radius.value], color='red', size=15, legend_label=u.specType(spt))
             show(fig)
 
-        return radius.round(3), radius_unc.round(3)
+        return radius, radius_unc
 
     def generate(self, orders):
         """

sedkit/sed.py

@@ -157,6 +157,7 @@ def __init__(self, name='My Target', verbose=True, method_list=None, substellar=
         self._membership = None
         self._sky_coords = None
         self._evo_model = None
+        self._age_distribution = None
 
         # Static attributes
         self.evo_model = None
@@ -166,6 +167,7 @@ def __init__(self, name='My Target', verbose=True, method_list=None, substellar=
         self.multiple = False
         self.HRsequence = rel.DwarfSequence()
         self.substellar = substellar
+        self.age_distribution = None
 
         # Dictionary to keep track of references
         self._refs = {}
@@ -1074,6 +1076,33 @@ def evo_model(self, model):
         # Set as uncalculated
         self.calculated = False
 
+    @property
+    def age_distribution(self):
+        """
+        A getter for the age distribution
+        """
+        return self._age_distribution
+
+    @age_distribution.setter
+    def age_distribution(self, distribution):
+        """
+        A setter for the age distribution
+
+        Parameters
+        ----------
+        distribution: str
+            The chosen age distribution
+        """
+        if distribution is None:
+            self._age_distribution = 'uniform'
+            self.message("Setting age distribution to uniform")
+        else:
+            if distribution == 'normal':
+                self._age_distribution = 'normal'
+                self.message("Setting age distribution to normal")
+            else:
+                raise ValueError("Please use either 'uniform' or 'normal' for age_distribution parameter")
+
     def export(self, parentdir='.', dirname=None, zipped=False):
         """
         Exports the photometry and results tables and a file of the
@@ -1813,6 +1842,9 @@ def infer_logg(self, plot=False):
         """
         self._logg = None
 
+        # Checking the age distribution
+        age_distribution = self.age_distribution
+
         # Try model isochrones
         if self.evo_model is not None and self.age is not None and self.Lbol_sun is not None:
 
@@ -1823,7 +1855,8 @@ def infer_logg(self, plot=False):
             if self.Lbol_sun[1] is None:
                 self.message('Lbol={0.Lbol}. Uncertainties are needed to calculate the surface gravity.'.format(self))
             else:
-                logg = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'logg', plot=plot)
+                logg = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'logg', distribution=age_distribution,
+                                               plot=plot)
 
             # Print a message if None
             if logg is None:
@@ -1848,6 +1881,9 @@ def infer_mass(self, mass_units=q.Msun, plot=False):
         """
         self._mass = None
 
+        # Checking the age distribution
+        age_distribution = self.age_distribution
+
         if self.substellar:
             mass_units = q.Mjup
 
@@ -1862,7 +1898,8 @@ def infer_mass(self, mass_units=q.Msun, plot=False):
                 self.message('Lbol={0.Lbol}. Uncertainties are needed to calculate the mass.'.format(self))
             else:
                 self.evo_model.mass_units = mass_units
-                mass = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'mass', plot=plot)
+                mass = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'mass', distribution=age_distribution,
+                                               plot=plot)
 
             # Store the value
             self.mass = [mass[0].round(0),  mass[1].round(0), mass[2].round(0), mass[3]] if mass is not None else mass
@@ -1907,6 +1944,9 @@ def infer_radius(self, radius_units=q.Rsun, infer_from='Lbol', plot=False):
         """
         self._radius = None
 
+        # Checking the age distribution
+        age_distribution = self.age_distribution
+
         # Change units to Jupiter radii for substellar objects
         if self.substellar:
             radius_units = q.Rjup
@@ -1946,7 +1986,8 @@ def infer_radius(self, radius_units=q.Rsun, infer_from='Lbol', plot=False):
                 self.message('Lbol={0.Lbol}. Uncertainties are needed to calculate the radius.'.format(self))
             else:
                 self.evo_model.radius_units = radius_units
-                radius = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'radius', plot=plot)
+                radius = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'radius',
+                                                 distribution=age_distribution, 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
@@ -2038,6 +2079,9 @@ def infer_Teff(self, teff_units=q.K, plot=False):
         """
         self._Teff = None
 
+        # Checking the age distribution
+        age_distribution = self.age_distribution
+
         # Try model isochrones
         if self.evo_model is not None and self.age is not None and self.Lbol_sun is not None:
 
@@ -2049,7 +2093,7 @@ def infer_Teff(self, teff_units=q.K, plot=False):
                 self.message('Lbol={0.Lbol}. Uncertainties are needed to calculate the teff.'.format(self))
             else:
                 self.evo_model.teff_units = teff_units
-                teff = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'teff', plot=plot)
+                teff = self.evo_model.evaluate(self.Lbol_sun, self.age, 'Lbol', 'teff', distribution=age_distribution, 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