diff --git a/EXOSIMS/Completeness/BrownCompleteness.py b/EXOSIMS/Completeness/BrownCompleteness.py
index 4a645fcfa..7a0dc0729 100644
--- a/EXOSIMS/Completeness/BrownCompleteness.py
+++ b/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
diff --git a/EXOSIMS/Completeness/SubtypeCompleteness.py b/EXOSIMS/Completeness/SubtypeCompleteness.py
index 38ecc4102..7cb047fbc 100644
--- a/EXOSIMS/Completeness/SubtypeCompleteness.py
+++ b/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))
diff --git a/EXOSIMS/Prototypes/PlanetPopulation.py b/EXOSIMS/Prototypes/PlanetPopulation.py
index af9310d59..a23fd526c 100644
--- a/EXOSIMS/Prototypes/PlanetPopulation.py
+++ b/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
diff --git a/EXOSIMS/Prototypes/SimulatedUniverse.py b/EXOSIMS/Prototypes/SimulatedUniverse.py
index 56ecec38d..157b678ab 100644
--- a/EXOSIMS/Prototypes/SimulatedUniverse.py
+++ b/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,7 +104,7 @@ 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 = {}
@@ -104,6 +112,9 @@ def __init__(self, fixedPlanPerStar=None, Min=None, cachedir=None,
         # 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
 
diff --git a/EXOSIMS/Prototypes/SurveySimulation.py b/EXOSIMS/Prototypes/SurveySimulation.py
index 6d87db7ce..c21329016 100644
--- a/EXOSIMS/Prototypes/SurveySimulation.py
+++ b/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)
diff --git a/EXOSIMS/Prototypes/TargetList.py b/EXOSIMS/Prototypes/TargetList.py
index 3417291cd..40d05be88 100644
--- a/EXOSIMS/Prototypes/TargetList.py
+++ b/EXOSIMS/Prototypes/TargetList.py
@@ -13,6 +13,15 @@
 import os.path
 import inspect
 import sys
+import json
+try:
+    import cPickle as pickle
+except:
+    import pickle
+try:
+    import urllib2
+except:
+    import urllib3
 
 
 class TargetList(object):
@@ -74,6 +83,15 @@ class TargetList(object):
             If not None, filters out any stars matching the names in the list.
         cachedir (str):
             Path to cache directory
+        filter_for_char (boolean):
+            TODO
+        earths_only (boolean):
+            TODO
+        getKnownPlanets (boolean):
+            a boolean indicating whether to grab the list of known planets from IPAC
+            and read the alias pkl file
+        I (numpy array):
+            array of star system inclinations
     
     """
 
@@ -82,7 +100,7 @@ class TargetList(object):
     def __init__(self, missionStart=60634, staticStars=True, 
         keepStarCatalog=False, fillPhotometry=False, explainFiltering=False, 
         filterBinaries=True, filterSubM=False, cachedir=None, filter_for_char=False,
-        earths_only=False, **specs):
+        earths_only=False, getKnownPlanets=False, **specs):
        
         #start the outspec
         self._outspec = {}
@@ -188,6 +206,15 @@ def __init__(self, missionStart=60634, staticStars=True,
                     c2=self.starprop(allInds, missionStart, eclip=True): \
                     c1[sInds] if eclip==False else c2[sInds]
 
+        #### Find Known Planets
+        self.getKnownPlanets = getKnownPlanets
+        self._outspec['getKnownPlanets'] = self.getKnownPlanets
+        if self.getKnownPlanets == True:
+            alias = self.loadAliasFile()
+            data = self.constructIPACurl()
+            starsWithPlanets = self.setOfStarsWithKnownPlanets(data)
+            knownPlanetBoolean = self.createKnownPlanetBoolean(alias,starsWithPlanets)
+
     def __str__(self):
         """String representation of the Target List object
         
@@ -278,6 +305,9 @@ def populate_target_list(self, popStars=None, **specs):
 
         # calculate 'true' and 'approximate' stellar masses
         self.stellar_mass()
+
+        # Calculate Star System Inclinations
+        self.I = self.gen_inclinations(self.PlanetPopulation.Irange)
         
         # include new attributes to the target list catalog attributes
         self.catalog_atts.append('comp0')
@@ -853,6 +883,37 @@ def stellarTeff(self, sInds):
         
         return Teff
 
+    def radiusFromMass(self,sInds):
+        """ Estimates the star radius based on its mass
+        Table 2, ZAMS models pg321 
+        STELLAR MASS-LUMINOSITY AND MASS-RADIUS RELATIONS OSMAN DEMIRCAN and GOKSEL KAHRAMAN 1991
+        Args:
+            sInds (list):
+                star indices
+        Return:
+            starRadius (numpy array):
+                star radius estimates
+        """
+        
+        M = self.MsTrue[sInds].value #Always use this??
+        a = -0.073
+        b = 0.668
+        starRadius = 10**(a+b*np.log(M))
+
+        return starRadius*u.R_sun
+
+    def gen_inclinations(self, Irange):
+        """Randomly Generate Inclination of Star System Orbital Plane
+        Args:
+            Irange (numpy array):
+                the range to generate inclinations over
+        Returns:
+            I (numpy array):
+                an array of star system inclinations
+        """
+        C = 0.5*(np.cos(Irange[0])-np.cos(Irange[1]))
+        return (np.arccos(np.cos(Irange[0]) - 2.*C*np.random.uniform(size=self.nStars))).to('deg')
+
     def dump_catalog(self):
         """Creates a dictionary of stellar properties for archiving use.
         
@@ -864,10 +925,146 @@ def dump_catalog(self):
                 Dictionary of star catalog properties
         
         """
-        atts = ['Name', 'Spec', 'parx', 'Umag', 'Bmag', 'Vmag', 'Rmag', 'Imag', 'Jmag', 'Hmag', 'Kmag', 'dist', 'BV', 'MV', 'BC', 'L', 'coords', 'pmra', 'pmdec', 'rv', 'Binary_Cut', 'MsEst', 'MsTrue', 'comp0', 'tint0']
-        # atts = ['Name', 'Spec', 'parx', 'Umag', 'Bmag', 'Vmag', 'Rmag', 'Imag', 'Jmag', 'Hmag', 'Kmag', 'dist', 'BV', 'MV', 'BC', 'L', 'coords', 'pmra', 'pmdec', 'rv', 'Binary_Cut']
+        atts = ['Name', 'Spec', 'parx', 'Umag', 'Bmag', 'Vmag', 'Rmag', 'Imag', 'Jmag', 'Hmag', 'Kmag', 'dist', 'BV', 'MV', 'BC', 'L', 'coords', 'pmra', 'pmdec', 'rv', 'Binary_Cut', 'MsEst', 'MsTrue', 'comp0', 'tint0', 'I']
         #Not sure if MsTrue and others can be dumped properly...
 
         catalog = {atts[i]: getattr(self,atts[i]) for i in np.arange(len(atts))}
 
         return catalog
+
+    def constructIPACurl(self, tableInput="exoplanets", columnsInputList=['pl_hostname','ra','dec','pl_discmethod','pl_pnum','pl_orbper','pl_orbsmax','pl_orbeccen',\
+        'pl_orbincl','pl_bmassj','pl_radj','st_dist','pl_tranflag','pl_rvflag','pl_imgflag',\
+        'pl_astflag','pl_omflag','pl_ttvflag', 'st_mass', 'pl_discmethod'],\
+        formatInput='json'):
+        """
+        Extracts Data from IPAC
+        Instructions for to interface with ipac using API
+        https://exoplanetarchive.ipac.caltech.edu/applications/DocSet/index.html?doctree=/docs/docmenu.xml&startdoc=item_1_01
+        Args:
+            tableInput (string):
+                describes which table to query
+            columnsInputList (list):
+                List of strings from https://exoplanetarchive.ipac.caltech.edu/docs/API_exoplanet_columns.html 
+            formatInput (string):
+                string describing output type. Only support JSON at this time
+        Returns:
+            data (dict):
+                a dictionary of IPAC data
+        """
+        baseURL = "https://exoplanetarchive.ipac.caltech.edu/cgi-bin/nstedAPI/nph-nstedAPI?"
+        tablebaseURL = "table="
+        # tableInput = "exoplanets" # exoplanets to query exoplanet table
+        columnsbaseURL = "&select=" # Each table input must be separated by a comma
+        # columnsInputList = ['pl_hostname','ra','dec','pl_discmethod','pl_pnum','pl_orbper','pl_orbsmax','pl_orbeccen',\
+        #                     'pl_orbincl','pl_bmassj','pl_radj','st_dist','pl_tranflag','pl_rvflag','pl_imgflag',\
+        #                     'pl_astflag','pl_omflag','pl_ttvflag', 'st_mass', 'pl_discmethod']
+                            #https://exoplanetarchive.ipac.caltech.edu/docs/API_exoplanet_columns.html for explanations
+
+        """
+        pl_hostname - Stellar name most commonly used in the literature.
+        ra - Right Ascension of the planetary system in decimal degrees.
+        dec - Declination of the planetary system in decimal degrees.
+        pl_discmethod - Method by which the planet was first identified.
+        pl_pnum - Number of planets in the planetary system.
+        pl_orbper - Time the planet takes to make a complete orbit around the host star or system.
+        pl_orbsmax - The longest radius of an elliptic orbit, or, for exoplanets detected via gravitational microlensing or direct imaging,\
+                    the projected separation in the plane of the sky. (AU)
+        pl_orbeccen - Amount by which the orbit of the planet deviates from a perfect circle.
+        pl_orbincl - Angular distance of the orbital plane from the line of sight.
+        pl_bmassj - Best planet mass estimate available, in order of preference: Mass, M*sin(i)/sin(i), or M*sin(i), depending on availability,\
+                    and measured in Jupiter masses. See Planet Mass M*sin(i) Provenance (pl_bmassprov) to determine which measure applies.
+        pl_radj - Length of a line segment from the center of the planet to its surface, measured in units of radius of Jupiter.
+        st_dist - Distance to the planetary system in units of parsecs. 
+        pl_tranflag - Flag indicating if the planet transits its host star (1=yes, 0=no)
+        pl_rvflag -     Flag indicating if the planet host star exhibits radial velocity variations due to the planet (1=yes, 0=no)
+        pl_imgflag - Flag indicating if the planet has been observed via imaging techniques (1=yes, 0=no)
+        pl_astflag - Flag indicating if the planet host star exhibits astrometrical variations due to the planet (1=yes, 0=no)
+        pl_omflag -     Flag indicating whether the planet exhibits orbital modulations on the phase curve (1=yes, 0=no)
+        pl_ttvflag -    Flag indicating if the planet orbit exhibits transit timing variations from another planet in the system (1=yes, 0=no).\
+                        Note: Non-transiting planets discovered via the transit timing variations of another planet in the system will not have\
+                         their TTV flag set, since they do not themselves demonstrate TTVs.
+        st_mass - Amount of matter contained in the star, measured in units of masses of the Sun.
+        pl_discmethod - Method by which the planet was first identified.
+        """
+
+        columnsInput = ','.join(columnsInputList)
+        formatbaseURL = '&format='
+        # formatInput = 'json' #https://exoplanetarchive.ipac.caltech.edu/docs/program_interfaces.html#format
+
+        # Different acceptable "Inputs" listed at https://exoplanetarchive.ipac.caltech.edu/applications/DocSet/index.html?doctree=/docs/docmenu.xml&startdoc=item_1_01
+
+        myURL = baseURL + tablebaseURL + tableInput + columnsbaseURL + columnsInput + formatbaseURL + formatInput
+        try:
+            response = urllib2.urlopen(myURL)
+            data = json.load(response)
+        except:
+            http = urllib3.PoolManager()
+            r = http.request('GET', myURL)
+            data = json.loads(r.data.decode('utf-8'))
+        return data
+
+    def setOfStarsWithKnownPlanets(self, data):
+        """ From the data dict created in this script, this method extracts the set of unique star names
+        Args:
+            data (dict):
+                dict containing the pl_hostname of each star
+        Returns:
+            list (list):
+                list of star names with a known planet
+
+        """
+        starNames = list()
+        for i in np.arange(len(data)):
+            starNames.append(data[i]['pl_hostname'])
+        return list(set(starNames))
+
+    def loadAliasFile(self):
+        """
+        Args:
+        Returns:
+            alias ():
+                list 
+        """
+        #OLD aliasname = 'alias_4_11_2019.pkl'
+        aliasname = 'alias_10_07_2019.pkl'
+        tmp1 = inspect.getfile(self.__class__).split('/')[:-2]
+        tmp1.append('util')
+        self.classpath = '/'.join(tmp1)
+        #self.classpath = os.path.split(inspect.getfile(self.__class__))[0]
+        #vprint(inspect.getfile(self.__class__))
+        self.alias_datapath = os.path.join(self.classpath, aliasname)
+        #Load pkl and outspec files
+        try:
+            with open(self.alias_datapath, 'rb') as f:#load from cache
+                alias = pickle.load(f, encoding='latin1')
+        except:
+            vprint('Failed to open fullPathPKL %s'%self.alias_datapath)
+            pass
+        return alias
+    ##########################################################
+
+    def createKnownPlanetBoolean(self, alias, starsWithPlanets):
+        """
+        Args:
+            alias ():
+
+            starsWithPlanets ():
+
+        Returns:
+            knownPlanetBoolean (numpy array):
+                boolean numpy array indicating whether the star has a planet (true)
+                or does not have a planet (false)
+
+        """
+        #Create List of Stars with Known Planets
+        knownPlanetBoolean = np.zeros(self.nStars, dtype=bool)
+        for i in np.arange(self.nStars):
+            #### Does the Star Have a Known Planet
+            starName = self.Name[i]#Get name of the current star
+            if starName in alias[:,1]:
+                indWhereStarName = np.where(alias[:,1] == starName)[0][0]# there should be only 1
+                starNum = alias[indWhereStarName,3]#this number is identical for all names of a target
+                aliases = [alias[j,1] for j in np.arange(len(alias)) if alias[j,3]==starNum] # creates a list of the known aliases
+                if np.any([True if aliases[j] in starsWithPlanets else False for j in np.arange(len(aliases))]):
+                    knownPlanetBoolean[i] = 1
+        return knownPlanetBoolean
diff --git a/EXOSIMS/Prototypes/ZodiacalLight.py b/EXOSIMS/Prototypes/ZodiacalLight.py
index a5331ac6f..ae3b686ea 100644
--- a/EXOSIMS/Prototypes/ZodiacalLight.py
+++ b/EXOSIMS/Prototypes/ZodiacalLight.py
@@ -43,7 +43,7 @@ class ZodiacalLight(object):
 
     _modtype = 'ZodiacalLight'
     
-    def __init__(self, magZ=23, magEZ=22, varEZ=0, cachedir=None, **specs):
+    def __init__(self, magZ=23, magEZ=22, varEZ=0, cachedir=None, commonSystemfEZ=False, **specs):
 
         #start the outspec
         self._outspec = {}
@@ -64,6 +64,10 @@ def __init__(self, magZ=23, magEZ=22, varEZ=0, cachedir=None, **specs):
         
         assert self.varEZ >= 0, "Exozodi variation must be >= 0"
         
+        #### Common Star System Number of Exo-zodi
+        self.commonSystemfEZ = commonSystemfEZ #ZL.nEZ must be calculated in SU
+        self._outspec['commonSystemfEZ'] = self.commonSystemfEZ
+
         # populate outspec
         for att in self.__dict__:
             if att not in ['vprint','_outspec']:
@@ -139,13 +143,11 @@ def fEZ(self, MV, I, d):
         # apparent magnitude of the Sun (in the V band)
         MVsun = 4.83
         
-        # assume log-normal distribution of variance
-        nEZ = np.ones(len(MV))
-        if self.varEZ != 0:
-            mu = np.log(nEZ) - 0.5*np.log(1. + self.varEZ/nEZ**2)
-            v = np.sqrt(np.log(self.varEZ/nEZ**2 + 1.))
-            nEZ = np.random.lognormal(mean=mu, sigma=v, size=len(MV))
-        
+        if self.commonSystemfEZ:
+            nEZ = self.nEZ
+        else:
+            nEZ = self.gen_systemnEZ(len(MV))
+
         # supplementary angle for inclination > 90 degrees
         beta = I.to('deg').value
         mask = np.where(beta > 90)[0]
@@ -159,6 +161,24 @@ def fEZ(self, MV, I, d):
         
         return fEZ
 
+    def gen_systemnEZ(self, nStars):
+        """ Ranomly generates the number of Exo-Zodi
+        Args:
+            nStars (int):
+                number of exo-zodi to generate
+        Returns:
+            nEZ (numpy array):
+                numpy array of exo-zodi randomly selected from fitsdata
+        """
+        # assume log-normal distribution of variance
+        nEZ = np.ones(nStars)
+        if self.varEZ != 0:
+            mu = np.log(nEZ) - 0.5*np.log(1. + self.varEZ/nEZ**2)
+            v = np.sqrt(np.log(self.varEZ/nEZ**2 + 1.))
+            nEZ = np.random.lognormal(mean=mu, sigma=v, size=nStars)
+
+        return nEZ
+
     def generate_fZ(self, Obs, TL, TK, mode, hashname):
         """Calculates fZ values for all stars over an entire orbit of the sun
             
diff --git a/EXOSIMS/ZodiacalLight/Mennesson.py b/EXOSIMS/ZodiacalLight/Mennesson.py
index 01413d82f..276b09db6 100644
--- a/EXOSIMS/ZodiacalLight/Mennesson.py
+++ b/EXOSIMS/ZodiacalLight/Mennesson.py
@@ -54,16 +54,11 @@ def fEZ(self, MV, I, d):
         # apparent magnitude of the Sun (in the V band)
         MVsun = 4.83
         
-        # assume log-normal distribution of variance
-        # nEZ = np.ones(len(MV))
-        # if self.varEZ != 0:
-        #     mu = np.log(nEZ) - 0.5*np.log(1. + self.varEZ/nEZ**2)
-        #     v = np.sqrt(np.log(self.varEZ/nEZ**2 + 1.))
-        #     nEZ = np.random.lognormal(mean=mu, sigma=v, size=len(MV))
+        if self.commonSystemfEZ:
+            nEZ = self.nEZ
+        else:
+            nEZ = self.gen_systemnEZ(len(MV))
 
-        nEZ_seed = np.random.randint(len(self.fitsdata) - len(MV))
-        nEZ = self.fitsdata[nEZ_seed:(nEZ_seed + len(MV))]
-        
         # supplementary angle for inclination > 90 degrees
         beta = I.to('deg').value
         mask = np.where(beta > 90)[0]
@@ -77,3 +72,16 @@ def fEZ(self, MV, I, d):
                 MVsun))*2*fbeta/d.to('AU').value**2/u.arcsec**2
         
         return fEZ
+
+    def gen_systemnEZ(self, nStars):
+        """ Ranomly generates the number of Exo-Zodi
+        Args:
+            nStars (int):
+                number of exo-zodi to generate
+        Returns:
+            nEZ (numpy array):
+                numpy array of exo-zodi randomly selected from fitsdata
+        """
+        nEZ_seed = np.random.randint(len(self.fitsdata) - nStars)
+        nEZ = self.fitsdata[nEZ_seed:(nEZ_seed + nStars)]
+        return nEZ
diff --git a/EXOSIMS/ZodiacalLight/Stark.py b/EXOSIMS/ZodiacalLight/Stark.py
index ceea916e3..2564396d8 100644
--- a/EXOSIMS/ZodiacalLight/Stark.py
+++ b/EXOSIMS/ZodiacalLight/Stark.py
@@ -125,8 +125,14 @@ def calcfbetaInput(self):
         return  points, values
 
     def calclogf(self):
+        """
         # wavelength dependence, from Table 19 in Leinert et al 1998
         # interpolated w/ a quadratic in log-log space
+        Returns:
+            interpolant (object):
+                a 1D quadratic interpolant of intensity vs wavelength
+
+        """
         self.zodi_lam = np.array([0.2, 0.3, 0.4, 0.5, 0.7, 0.9, 1.0, 1.2, 2.2, 3.5,
                 4.8, 12, 25, 60, 100, 140]) # um
         self.zodi_Blam = np.array([2.5e-8, 5.3e-7, 2.2e-6, 2.6e-6, 2.0e-6, 1.3e-6,
@@ -134,7 +140,7 @@ def calclogf(self):
                 3.2e-9, 6.9e-10]) # W/m2/sr/um
         x = np.log10(self.zodi_lam)
         y = np.log10(self.zodi_Blam)
-        return interp1d(x, y, kind='quadratic')# logf = 
+        return interp1d(x, y, kind='quadratic')
 
     def calcfZmax(self, sInds, Obs, TL, TK, mode, hashname):
         """Finds the maximum zodiacal light values for each star over an entire orbit of the sun not including keeoput angles
@@ -156,6 +162,7 @@ def calcfZmax(self, sInds, Obs, TL, TK, mode, hashname):
                 the maximum fZ
             absTimefZmax[sInds] (astropy Time array):
                 returns the absolute Time the maximum fZ occurs (for the prototype, these all have the same value)
+
         """
         #Generate cache Name########################################################################
         cachefname = hashname + 'fZmax'
diff --git a/EXOSIMS/util/plotKeepoutMap.py b/EXOSIMS/util/plotKeepoutMap.py
index 97a63bec8..1aded4133 100644
--- a/EXOSIMS/util/plotKeepoutMap.py
+++ b/EXOSIMS/util/plotKeepoutMap.py
@@ -158,7 +158,7 @@ def singleRunPostProcessing(self, PPoutpath=None, folder=None):
 
         outline=PathEffects.withStroke(linewidth=5, foreground='black')
         plt.plot([-1.,-1.],[-1.,-1.],color=cmap.colors[0],label='Visible',path_effects=[outline])
-        plt.plot([-1.,-1.],[-1.,-1.],color=cmap.colors[1],label=ur'$\u2609$')
+        plt.plot([-1.,-1.],[-1.,-1.],color=cmap.colors[1],label=ur"$\u2609$")
         plt.plot([-1.,-1.],[-1.,-1.],color=cmap.colors[2],label=ur'$\oplus$')##\u2641$')
         plt.plot([-1.,-1.],[-1.,-1.],color=cmap.colors[3],label=ur'$\u263D$')
         plt.plot([-1.,-1.],[-1.,-1.],color=cmap.colors[4],label=ur'$\u2642\u263F$')
diff --git a/EXOSIMS/util/waypoint.py b/EXOSIMS/util/waypoint.py
index 3d15ff981..c91f7b485 100644
--- a/EXOSIMS/util/waypoint.py
+++ b/EXOSIMS/util/waypoint.py
@@ -4,7 +4,8 @@
 import inspect
 
 def waypoint(comps, intTimes, duration, mpath, tofile):
-    """generates waypoint dictionary for MissionSim
+    """ Generates waypoint dictionary for MissionSim
+
     Args:
         comps (array):
             An array of completeness values for all stars
@@ -20,6 +21,7 @@ def waypoint(comps, intTimes, duration, mpath, tofile):
         out (dictionary):
             Output dictionary containing the number of stars visited, the total completness
             achieved, and the amount of time spent integrating.
+            
     """
 
     CbT = comps/intTimes
diff --git a/documentation/EXOSIMS.Completeness.rst b/documentation/EXOSIMS.Completeness.rst
index 283766174..37344a2ce 100644
--- a/documentation/EXOSIMS.Completeness.rst
+++ b/documentation/EXOSIMS.Completeness.rst
@@ -25,3 +25,11 @@ EXOSIMS.Completeness.GarrettCompleteness module
    :undoc-members:
    :show-inheritance:
 
+EXOSIMS.Completeness.SubtypeCompleteness module
+-----------------------------------------------
+
+.. automodule:: EXOSIMS.Completeness.SubtypeCompleteness
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
diff --git a/documentation/EXOSIMS.util.KeplerSTM_C.rst b/documentation/EXOSIMS.util.KeplerSTM_C.rst
index 534323875..4276cfd97 100644
--- a/documentation/EXOSIMS.util.KeplerSTM_C.rst
+++ b/documentation/EXOSIMS.util.KeplerSTM_C.rst
@@ -9,14 +9,6 @@ EXOSIMS.util.KeplerSTM\_C package
 Submodules
 ----------
 
-EXOSIMS.util.KeplerSTM\_C.CyKeplerSTM.cpython\-37m\-darwin module
------------------------------------------------------------------
-
-.. automodule:: EXOSIMS.util.KeplerSTM_C.CyKeplerSTM.cpython-37m-darwin
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
 EXOSIMS.util.KeplerSTM\_C.CyKeplerSTM module
 --------------------------------------------
 
diff --git a/documentation/EXOSIMS.util.rst b/documentation/EXOSIMS.util.rst
index 012fed645..c3ad95476 100644
--- a/documentation/EXOSIMS.util.rst
+++ b/documentation/EXOSIMS.util.rst
@@ -64,14 +64,6 @@ EXOSIMS.util.deltaMag module
    :undoc-members:
    :show-inheritance:
 
-EXOSIMS.util.depthOfSearch module
----------------------------------
-
-.. automodule:: EXOSIMS.util.depthOfSearch
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
 EXOSIMS.util.eccanom module
 ---------------------------
 
@@ -176,14 +168,6 @@ EXOSIMS.util.plotCompletenessJointPDFs module
    :undoc-members:
    :show-inheritance:
 
-EXOSIMS.util.plotKeepoutMap module
-----------------------------------
-
-.. automodule:: EXOSIMS.util.plotKeepoutMap
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
 EXOSIMS.util.plotPlanetPopRvsAandDetectedRvsA module
 ----------------------------------------------------
 
diff --git a/documentation/_static/icd.pdf b/documentation/_static/icd.pdf
index a76ef0370..cfbf17e93 100644
Binary files a/documentation/_static/icd.pdf and b/documentation/_static/icd.pdf differ
diff --git a/documentation/builddocs.sh b/documentation/builddocs.sh
index 4f700dc2a..5922b804b 100755
--- a/documentation/builddocs.sh
+++ b/documentation/builddocs.sh
@@ -8,7 +8,7 @@ if [ ! -d "../EXOSIMS/Prototypes" ] || [ `basename $PWD` != "documentation" ] ;
 fi
 
 #sphinx-apidoc -f -o . ../EXOSIMS/
-sphinx-apidoc -M -f -o . ../EXOSIMS/ ../EXOSIMS/util/runPostProcessing.py ../EXOSIMS/util/plotConvergencevsNumberofRuns.py ../EXOSIMS/util/plotTimeline.py ../EXOSIMS/util/evenlyDistributePointsOnSphere.py ../EXOSIMS/util/KeplerSTM_C/CyKeplerSTM_setup.py
+sphinx-apidoc -M -f -o . ../EXOSIMS/ ../EXOSIMS/util/runPostProcessing.py ../EXOSIMS/util/plotConvergencevsNumberofRuns.py ../EXOSIMS/util/plotTimeline.py ../EXOSIMS/util/evenlyDistributePointsOnSphere.py ../EXOSIMS/util/KeplerSTM_C/CyKeplerSTM_setup.py ../EXOSIMS/util/plotKeepoutMap.py ../EXOSIMS/util/depthOfSearch.py
 
 rm modules.rst