updating comments and which planet phase function behavior

EXOSIMS/Completeness/BrownCompleteness.py

@@ -51,7 +51,11 @@ def __init__(self, Nplanets=1e8, **specs):
         self.Nplanets = int(Nplanets)
 
         # get path to completeness interpolant stored in a pickled .comp file
-        self.filename = self.PlanetPopulation.__class__.__name__ + self.PlanetPhysicalModel.__class__.__name__ + self.__class__.__name__ + str(self.Nplanets) + self.PlanetPhysicalModel.whichPlanetPhaseFunction
+        if callable(self.PlanetPhysicalModel.whichPlanetPhaseFunction): 
+            pfName = self.PlanetPhysicalModel.whichPlanetPhaeFunction.__name__
+        else:
+            pfName = self.PlanetPhysicalModel.whichPlanetPhaseFunction
+        self.filename = self.PlanetPopulation.__class__.__name__ + self.PlanetPhysicalModel.__class__.__name__ + self.__class__.__name__ + str(self.Nplanets) + pfName
 
         # get path to dynamic completeness array in a pickled .dcomp file
         self.dfilename = self.PlanetPopulation.__class__.__name__ + \

EXOSIMS/Completeness/SubtypeCompleteness.py

@@ -98,7 +98,9 @@ def __init__(self, Nplanets=1e8, binTypes='kopparapuBins_extended', **specs):
             pmax = self.PlanetPopulation.prange[1]
             emax = self.PlanetPopulation.erange[1] #Maximum orbital Eccentricity
 
-            #THE RMAX AND RMIN HERE IS WRONG... NEED TO DO ANOTHER FORMULATION
+            #NEED TO DO ANOTHER FORMULATION for rmin and rmax
+            #rmax and rmin assume a fixed sma but should vary as a function of stellar luminosity
+            #the upper and lower dmaglimits of the planet subtype bin should be based on the brightest and faintest stellar luminosities of the bin
 
             #ravgt_rtLstar_lo = 1./np.sqrt(self.L_lo[ii,j]) #these are the classification bin edges being used
             #ravgt_rtLstar_hi = 1./np.sqrt(self.L_hi[ii,j]) #these are the classification bin edges being used

EXOSIMS/PlanetPopulation/Brown2005EarthLike.py

@@ -4,7 +4,7 @@
 
 class Brown2005EarthLike(PlanetPopulation):
     """
-    Population of Earth-Like Planets from Brown 2005 papernetary/factsheet/jupiterfact.html
+    Population of Earth-Like Planets from Brown 2005 paper
     
     This implementation is intended to enforce this population regardless
     of JSON inputs.  The only inputs that will not be disregarded are erange
@@ -41,14 +41,14 @@ def gen_plan_params(self, n):
         
         Returns:
             tuple:
-            a (astropy Quantity array):
-                Semi-major axis in units of AU
-            e (float ndarray):
-                Eccentricity
-            p (float ndarray):
-                Geometric albedo
-            Rp (astropy Quantity array):
-                Planetary radius in units of earthRad
+                a (astropy Quantity array):
+                    Semi-major axis in units of AU
+                e (float ndarray):
+                    Eccentricity
+                p (float ndarray):
+                    Geometric albedo
+                Rp (astropy Quantity array):
+                    Planetary radius in units of earthRad
         
         """
         n = self.gen_input_check(n)

EXOSIMS/Prototypes/Observatory.py

@@ -428,8 +428,8 @@ def keepout(self, TL, sInds, currentTime, koangles, returnExtra=False):
                     s x 11 x 2 element array of minimum and maximum keepouts used for each body.  
                     Same ordering as r_body.
         Returns:
-            kogood (boolean ndarray):
-                s x n x m array of boolean values. True is a target unobstructed and observable, 
+            boolean ndarray:
+                kogood s x n x m array of boolean values. True is a target unobstructed and observable, 
                 and False is a target unobservable due to obstructions in the keepout zone. 
         
         """

EXOSIMS/Prototypes/PlanetPhysicalModel.py

@@ -19,6 +19,8 @@ class PlanetPhysicalModel(object):
     Attributes:
         cachedir (str):
             Path to EXOSIMS cache directory
+        whichPlanetPhaseFunction (str or callable):
+            planet phase function to use
             
     """
 
@@ -38,13 +40,16 @@ def __init__(self, cachedir=None, whichPlanetPhaseFunction='lambert', **specs):
         self.vprint = vprint(specs.get('verbose', True))
 
         #Select which Phase Function to use
+        assert isinstance(whichPlanetPhaseFunction, str), "whichPlanetPhaseFunction is not a string"
         self.whichPlanetPhaseFunction = whichPlanetPhaseFunction
         if whichPlanetPhaseFunction == 'quasiLambertPhaseFunction':
             from EXOSIMS.util.phaseFunctions import quasiLambertPhaseFunction
             self.calc_Phi = quasiLambertPhaseFunction
         elif whichPlanetPhaseFunction == 'hyperbolicTangentPhaseFunc':
             from EXOSIMS.util.phaseFunctions import hyperbolicTangentPhaseFunc
             self.calc_Phi = hyperbolicTangentPhaseFunc
+        if callable(whichPlanetPhaseFunction):
+            self.calc_Phi = whichPlanetPhaseFunction
         #else: if whichPlanetPhaseFunction == 'lambert': Default, Do nothing
         self._outspec['whichPlanetPhaseFunction'] = whichPlanetPhaseFunction
 

EXOSIMS/util/phaseFunctions.py

@@ -5,13 +5,13 @@
 import numpy as np
 
 def phi_lambert(alpha):
-    """ Lambert phase function as presented in Garrett2016
+    """ Lambert phase function most easily found in Garrett2016 and initially presented in Sobolev 1975
     Args:
         alpha (float):
             phase angle in radians
     Returns:
-        Phi (float):
-            phase function value between 0 and 1
+        float:
+            Phi, phase function value between 0 and 1
     """
     phi = (np.sin(alpha) + (np.pi-alpha)*np.cos(alpha))/np.pi
     return phi
@@ -24,8 +24,8 @@ def transitionStart(x,a,b):
         a (float):
             transition midpoint in deg
     Returns:
-        s (float):
-            Transition value from 0 to 1
+        float:
+            s, Transition value from 0 to 1
     """
     s = 0.5+0.5*np.tanh((x-a)/b)
     return s
@@ -40,23 +40,23 @@ def transitionEnd(x,a,b):
         a (float):
             transition midpoint in deg
     Returns:
-        s (float):
-            Transition value from 1 to 0
+        float:
+            s, transition value from 1 to 0
     """
     s = 0.5-0.5*np.tanh((x-a)/b)
     return s
 
 
 def quasiLambertPhaseFunction(beta):
     """ Quasi Lambert Phase Function as presented
-    Analyticall Invertible Phase function from Agol 2007, 'Rounding up the wanderers: optimizing
+    Analytically Invertible Phase function from Agol 2007, 'Rounding up the wanderers: optimizing
     coronagraphic searches for extrasolar planets'
     Args:
         beta (numpy array):
             planet phase angles in radians
     Returns:
-        Phi (numpy array):
-            phase function value
+        ndarray:
+            Phi, phase function value
     """
     Phi = np.cos(beta/2.)**4
     return Phi
@@ -68,8 +68,8 @@ def quasiLambertPhaseFunctionInverse(Phi):
             phase function value
         
     Returns:
-        beta (numpy array):
-            planet phase angles
+        ndarray:
+            beta, planet phase angles
     """
     beta = 2.*np.arccos((Phi)**(1./4.))
     return beta
@@ -95,8 +95,8 @@ def hyperbolicTangentPhaseFunc(beta,A,B,C,D,planetName=None):
             planet name string all lower case for one of 8 solar system planets
 
     Returns:
-        Phi (float):
-            phase angle in degrees
+        float:
+            Phi, phase angle in degrees
     """
     if planetName == None:
         None #do nothing
@@ -141,8 +141,8 @@ def hyperbolicTangentPhaseFuncInverse(Phi,A,B,C,D,planetName=None):
             planet name string all lower case for one of 8 solar system planets
 
     Returns:
-        beta (float):
-            Phase Angle  in degrees
+        float:
+            beta, Phase Angle  in degrees
     """
     if planetName == None:
         None #do nothing
@@ -175,8 +175,8 @@ def betaFunc(inc,v,w):
         w (numpy array):
             planet argument of periapsis
     Returns:
-        beta (numpy array):
-            planet phase angle
+        ndarray:
+            beta, planet phase angle
     """
     beta = np.arccos(np.sin(inc)*np.sin(v+w))
     return beta