Merge 17b90fe into 43ba5bc

EXOSIMS/Completeness/BrownCompleteness.py

@@ -66,6 +66,7 @@ def __init__(self, Nplanets=1e8, **specs):
                 self.extstr += '%s: ' % att + str(getattr(self.PlanetPopulation, att)) + ' '
         ext = hashlib.md5(self.extstr.encode("utf-8")).hexdigest()
         self.filename += ext
+        self.filename.replace(" ","") #Remove spaces from string (in the case of prototype use)
 
     def target_completeness(self, TL, calc_char_comp0=False):
         """Generates completeness values for target stars

EXOSIMS/Completeness/SubtypeCompleteness.py

@@ -30,7 +30,7 @@ class SubtypeCompleteness(BrownCompleteness):
     Completeness Module calculations in exoplanet mission simulation.
     
     Args:
-        \*\*specs: 
+        specs: 
             user specified values
     
     Attributes:
@@ -475,12 +475,26 @@ def SubtypeHist(self, nplan, xedges, yedges, TL):
         Returns:
             float ndarray:
                 2D numpy ndarray containing completeness frequencies
-            hs
-            binInds
-            h_earthLike
-            h
-            xedges
-            yedges
+            hs (dict):
+                dict with index [bini,binj] containing arrays of counts per Array bin
+            bini (float):
+                planet size type index
+            binj (float):
+                planet incident stellar flux index
+            h_earthLike (float):
+                number of earth like exoplanets
+            h (numpy array):
+                2D numpy array of bin counts over all dmag vs s
+            xedges (numpy array):
+                array of bin edges originally input and used in histograms
+            yedges (numpy array):
+                array of bin edges originally input and used in histograms
+            counts (dict):
+                dict with index [bini,binj] containing total number of planets in bini,binj
+            count (float):
+                total number of planets in the whole populatiom
+            count_earthLike (float):
+                total number of earth-like planets in the whole population
         
         """
         tStartSTH = time.time()
@@ -525,8 +539,7 @@ def genplans(self, nplan, TL):
         """Generates planet data needed for Monte Carlo simulation
         
         Args:
-            nplan (integer):
-                Number of planets
+            nplan (integer) - Number of planets
             TL (target list object)
                 
         Returns:
@@ -535,9 +548,13 @@ def genplans(self, nplan, TL):
                 Planet apparent separations in units of AU
             dMag (ndarray):
                 Difference in brightness
-            bini (int) - planet size-type: 0-rocky, 1- Super-Earths, 2- sub-Neptunes, 3- sub-Jovians, 4- Jovians
-            binj (int) - planet incident stellar-flux: 0- hot, 1- warm, 2- cold
-            earthLike (bool) - boolean indicating whether the planet is earthLike or not earthLike
+            bini (int):
+                planet size-type: 0-rocky, 1- Super-Earths, 2- sub-Neptunes, 3- sub-Jovians, 4- Jovians
+            binj (int):
+                planet incident stellar-flux: 0- hot, 1- warm, 2- cold
+            earthLike (bool):
+                boolean indicating whether the planet is earthLike or not earthLike
+
         """
 
         PPop = self.PlanetPopulation
@@ -771,6 +788,7 @@ def comps_input_reshape(self, intTimes, TL, sInds, fZ, fEZ, WA, mode, C_b=None,
                 Maximum projected separations in AU
             dMag (ndarray):
                 Difference in brightness magnitude
+
         """
 
         # cast inputs to arrays and check
@@ -857,16 +875,20 @@ def calc_fdmag(self, dMag, smin, smax, subpop=-2):
     def putPlanetsInBoxes(self,out,TL):
         """ Classifies planets in a gen_summary out file by their hot/warm/cold and rocky/superearth/subneptune/subjovian/jovian bins
         Args:
-            out () - a gen_summary output list
-            TL () - 
+            out (dict):
+                a gen_summary output dict
+            TL (object):
+                a target list object
         Returns:
-            aggbins (list) - dims [# simulations, 5x3 numpy array]
-            earthLikeBins (list) - dims [# simulations]
+            aggbins (list):
+                dims [# simulations, 5x3 numpy array]
+            earthLikeBins (list):
+                dims [# simulations]
+
         """
         aggbins = list()
         earthLikeBins = list()
         bins = np.zeros((self.L_bins.shape[0]-1,self.L_bins.shape[1]-1)) # planet type, planet temperature
-        #DELETE bins = np.zeros((5,3)) # planet type, planet temperature
         #planet types: rockey, super-Earths, sub-Neptunes, sub-Jovians, Jovians
         #planet temperatures: cold, warm, hot
         for i in np.arange(len(out['starinds'])): # iterate over simulations
@@ -898,14 +920,22 @@ def putPlanetsInBoxes(self,out,TL):
     def classifyPlanets(self, Rp, TL, starind, sma, ej):
         """ Determine Kopparapu bin of an individual planet. Verified with Kopparapu Extended
         Args:
-            Rp (float) - planet radius in Earth Radii
-            TL (object) - EXOSIMS target list object
-            sma (float) - planet semi-major axis in AU
-            ej (float) - planet eccentricity
+            Rp (float):
+                planet radius in Earth Radii
+            TL (object):
+                EXOSIMS target list object
+            sma (float):
+                planet semi-major axis in AU
+            ej (float):
+                planet eccentricity
         Returns:
-            bini (int) - planet size-type: 0-rocky, 1- Super-Earths, 2- sub-Neptunes, 3- sub-Jovians, 4- Jovians
-            binj (int) - planet incident stellar-flux: 0- hot, 1- warm, 2- cold
-            earthLike (bool) - boolean indicating whether the planet is earthLike or not earthLike
+            bini (int):
+                planet size-type: 0-rocky, 1- Super-Earths, 2- sub-Neptunes, 3- sub-Jovians, 4- Jovians
+            binj (int):
+                planet incident stellar-flux: 0- hot, 1- warm, 2- cold
+            earthLike (bool):
+                boolean indicating whether the planet is earthLike or not earthLike
+
         """
         Rp = Rp.to('earthRad').value
         sma = sma.to('AU').value
@@ -958,14 +988,22 @@ def classifyPlanets(self, Rp, TL, starind, sma, ej):
     def classifyPlanet(self, Rp, TL, starind, sma, ej):
         """ Determine Kopparapu bin of an individual planet
         Args:
-            Rp (float) - planet radius in Earth Radii
-            TL (object) - EXOSIMS target list object
-            sma (float) - planet semi-major axis in AU
-            ej (float) - planet eccentricity
+            Rp (float):
+                planet radius in Earth Radii
+            TL (object):
+                EXOSIMS target list object
+            sma (float):
+                planet semi-major axis in AU
+            ej (float):
+                planet eccentricity
         Returns:
-            bini (int) - planet size-type: 0-rocky, 1- Super-Earths, 2- sub-Neptunes, 3- sub-Jovians, 4- Jovians
-            binj (int) - planet incident stellar-flux: 0- hot, 1- warm, 2- cold
-            earthLike (bool) - boolean indicating whether the planet is earthLike or not earthLike
+            bini (int):
+                planet size-type: 0-rocky, 1- Super-Earths, 2- sub-Neptunes, 3- sub-Jovians, 4- Jovians
+            binj (int):
+                planet incident stellar-flux: 0- hot, 1- warm, 2- cold
+            earthLike (bool):
+                boolean indicating whether the planet is earthLike or not earthLike
+
         """
         #Find Planet Rp range
         bini = np.where((self.Rp_lo < Rp)*(Rp < self.Rp_hi))[0] # index of planet size, rocky,...,jovian
@@ -1010,7 +1048,13 @@ def classifyPlanet(self, Rp, TL, starind, sma, ej):
 
 
     def kopparapuBins_old(self):
-        """
+        """ A function containing the Inner 12 Kopparapu bins
+        Updates the Rp_bins, Rp_lo, Rp_hi, L_bins, L_lo, and L_hi attributes
+        Args:
+
+        Returns:
+            None
+
         """
         #1: planet-radius bin-edges  [units = Earth radii]
         self.Rp_bins = np.array([0.5, 1.4, 4.0, 14.3]) # Old early 2018 had 3 bins
@@ -1032,7 +1076,13 @@ def kopparapuBins_old(self):
         return None
 
     def kopparapuBins(self):
-        """
+        """A function containing the Center 15 Kopparapu bins
+        Updates the Rp_bins, Rp_lo, Rp_hi, L_bins, L_lo, and L_hi attributes
+        Args:
+
+        Returns:
+            None
+
         """
         #1: planet-radius bin-edges  [units = Earth radii]
         # New (May 2018, 5 bins x 3 bins, see Kopparapu et al, arxiv:1802.09602v1,
@@ -1059,7 +1109,13 @@ def kopparapuBins(self):
         return None
 
     def kopparapuBins_extended(self):
-        """
+        """A function containing the Full 35 Kopparapu bins
+        Updates the Rp_bins, Rp_lo, Rp_hi, L_bins, L_lo, L_hi, and type_names attributes
+        Args:
+
+        Returns:
+            None
+
         """
         #1: planet-radius bin-edges  [units = Earth radii]
         self.Rp_bins = np.array([0., 0.5, 1.0, 1.75, 3.5, 6.0, 14.3, 11.2*4.6])
@@ -1139,20 +1195,29 @@ def dmag_limits(self,rmin,rmax,pmax,pmin,Rmax,Rmin,phaseFunc):
         Limits on dmag vs s JPDF from Garrett 2016
         #See https://github.com/dgarrett622/FuncComp/blob/master/FuncComp/util.py
         Args:
-            rmin (float) - minimum planet-star distance possible in AU
-            rmax (float) - maximum planet-star distance possible in AU
-            pmax (float) - maximum planet albedo
-            Rmax (float) - maximum planet radius in earthRad
+            rmin (float):
+                minimum planet-star distance possible in AU
+            rmax (float):
+                maximum planet-star distance possible in AU
+            pmax (float):
+                maximum planet albedo
+            Rmax (float):
+                maximum planet radius in earthRad
             phaseFunc (function) - with input in units of rad
         Returns:
-            dmag_limit_functions (list) - list of lambda functions taking in 's' with units of AU
-            lower_limits (list) - list of floats representing lower bounds on 's'
-            upper_limits (list) - list of floats representing upper bounds on 's'
+            dmag_limit_functions (list):
+                list of lambda functions taking in 's' with units of AU
+            lower_limits (list):
+                list of floats representing lower bounds on 's'
+            upper_limits (list):
+                list of floats representing upper bounds on 's'
+                
         """
         def betaStarFinder(beta,phaseFunc):
             """From Garrett 2016
             Args:
                 beta (float) in radians
+            
             """
             return -phaseFunc(beta*u.rad)*np.sin(beta)**2.
         res = minimize_scalar(betaStarFinder,args=(self.PlanetPhysicalModel.calc_Phi,),method='golden',tol=1e-4, bracket=(0.,np.pi/3.,np.pi))

EXOSIMS/Prototypes/PlanetPopulation.py

@@ -45,10 +45,6 @@ class PlanetPopulation(object):
         eta (float):
             Global occurrence rate defined as expected number of planets 
             per star in a given universe
-        uniform (float, callable):
-            Uniform distribution over a given range
-        logunif (float, callable):
-            Log-uniform distribution over a given range
         cachedir (str):
             Path to cache directory
         
@@ -186,7 +182,7 @@ def gen_mass(self, n):
 
         return Mp
 
-    def gen_angles(self, n):
+    def gen_angles(self, n, commonSystemInclinations=None):
         """Generate inclination, longitude of the ascending node, and argument
         of periapse in degrees
         
@@ -197,11 +193,14 @@ def gen_angles(self, n):
         Args:
             n (integer):
                 Number of samples to generate
+            commonSystemInclinations (None or tuple):
+                None if inclinations are to be generated for each planet individually
+                (mean, standard deviation )
                 
         Returns:
             tuple:
             I (astropy Quantity array):
-                Inclination in units of degrees
+                Inclination in units of degrees OR deviation in inclination in deg
             O (astropy Quantity array):
                 Longitude of the ascending node in units of degrees
             w (astropy Quantity array):
@@ -211,7 +210,10 @@ def gen_angles(self, n):
         n = self.gen_input_check(n)
         # inclination
         C = 0.5*(np.cos(self.Irange[0])-np.cos(self.Irange[1]))
-        I = (np.arccos(np.cos(self.Irange[0]) - 2.*C*np.random.uniform(size=n))).to('deg')
+        if commonSystemInclinations == None:
+            I = (np.arccos(np.cos(self.Irange[0]) - 2.*C*np.random.uniform(size=n))).to('deg')
+        else:
+            I = np.random.normal(loc=[0],scale=commonSystemInclinations[1],size=n)*u.deg
         # longitude of the ascending node
         Or = self.Orange.to('deg').value
         O = np.random.uniform(low=Or[0], high=Or[1], size=n)*u.deg

EXOSIMS/Prototypes/SimulatedUniverse.py

@@ -82,10 +82,18 @@ class SimulatedUniverse(object):
         WA (astropy Quantity array)
             Working angles of the planets of interest in units of arcsec
         fixedPlanPerStar (int or None):
-            Fixed number of planets to generate for each star
+            Fixed number of planets to generate for each star\
+        Min (float):
+            Initial constant Mean Anomaly
         cachedir (str):
             Path to cache directory
-    
+        lucky_planets (boolean):
+            TODO
+        commonSystemInclinations (boolean):
+            If True, all planets in the same star systems will share a common
+            system inclination with deviations specified to individual planets.
+            Also used for Exozodi
+
     Notes:
         PlanetPopulation.eta is treated as the rate parameter of a Poisson distribution.
         Each target's number of planets is a Poisson random variable sampled with \lambda=\eta.
@@ -96,14 +104,17 @@ class SimulatedUniverse(object):
     _modtype = 'SimulatedUniverse'
 
     def __init__(self, fixedPlanPerStar=None, Min=None, cachedir=None,
-                 lucky_planets=False, **specs):
+                 lucky_planets=False, commonSystemInclinations=None, **specs):
 
         #start the outspec
         self._outspec = {}
 
         # load the vprint function (same line in all prototype module constructors)
         self.vprint = vprint(specs.get('verbose', True))
         self.lucky_planets = lucky_planets
+        self._outspec['lucky_planets'] = lucky_planets
+        self.commonSystemInclinations = commonSystemInclinations
+        self._outspec['commonSystemInclinations'] = commonSystemInclinations
 
         # save fixed number of planets to generate
         self.fixedPlanPerStar = fixedPlanPerStar
@@ -196,13 +207,19 @@ def gen_physical_properties(self, **specs):
         self.nPlans = len(self.plan2star)
 
         # sample all of the orbital and physical parameters
-        self.I, self.O, self.w = PPop.gen_angles(self.nPlans)
+        self.I, self.O, self.w = PPop.gen_angles(self.nPlans, self.commonSystemInclinations)
+        if not self.commonSystemInclinations == None: #OVERWRITE I with TL.I+dI
+            self.I = self.I.copy() + TL.I[self.plan2star]
         self.a, self.e, self.p, self.Rp = PPop.gen_plan_params(self.nPlans)
         if PPop.scaleOrbits:
             self.a *= np.sqrt(TL.L[self.plan2star])
         self.gen_M0()                           # initial mean anomaly
         self.Mp = PPop.gen_mass(self.nPlans)    # mass
 
+        if self.ZodiacalLight.commonSystemfEZ == True:
+            if not hasattr(self.ZodiacalLight, 'nEZ'):#NEED TO CHECK IF ZL.nEZ WAS LOADED FROM OUTSPEC
+                self.ZodiacalLight.nEZ = self.ZodiacalLight.gen_systemnEZ(TL.nStars)
+
         # The prototype StarCatalog module is made of one single G star at 1pc. 
         # In that case, the SimulatedUniverse prototype generates one Jupiter 
         # at 5 AU to allow for characterization testing.
@@ -510,6 +527,10 @@ def dump_systems(self):
                'p':self.p,
                'plan2star':self.plan2star,
                'star':self.TargetList.Name[self.plan2star]}
+        if not self.commonSystemInclinations == None:
+            systems['starI'] = self.TargetList.I
+        if self.ZodiacalLight.commonSystemfEZ:
+            systems['starnEZ'] = self.ZodiacalLight.nEZ
 
         return systems
 

EXOSIMS/Prototypes/SurveySimulation.py

@@ -95,6 +95,10 @@ class SurveySimulation(object):
         scaleWAdMag (bool):
             If True, rescale dMagint and WAint for all stars based on luminosity and 
             to ensure that WA is within the IWA/OWA. Defaults False.
+        record_counts_path (TODO):
+            TODO
+        nokoMap (bool):
+            TODO
         cachedir (str):
             Path to cache directory
         defaultAddExoplanetObsTime (boolean):
@@ -1816,6 +1820,7 @@ def reset_sim(self, genNewPlanets=True, rewindPlanets=True, seed=None):
         # generate new planets if requested (default)
         if genNewPlanets:
             TL.stellar_mass()
+            TL.I = TL.gen_inclinations(TL.PlanetPopulation.Irange)
             SU.gen_physical_properties(**SU._outspec)
             rewindPlanets = True
         # re-initialize systems if requested (default)