Merge d5280ff into cc4b06f

.travis.yml

@@ -5,7 +5,7 @@ python:
     - "3.6"
     - "3.7"
 install:
-    - pip install -r requirements.txt
+    - pip install --upgrade -r requirements.txt
     - pip install -e .
 script:
     coverage run -m unittest discover -v

EXOSIMS/Completeness/BrownCompleteness.py

@@ -26,7 +26,7 @@ class BrownCompleteness(Completeness):
     Completeness Module calculations in exoplanet mission simulation.
     
     Args:
-        \*\*specs: 
+        specs: 
             user specified values
     
     Attributes:
@@ -51,7 +51,7 @@ 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.filename = self.PlanetPopulation.__class__.__name__ + self.PlanetPhysicalModel.__class__.__name__ + self.__class__.__name__
 
         # get path to dynamic completeness array in a pickled .dcomp file
         self.dfilename = self.PlanetPopulation.__class__.__name__ + \

EXOSIMS/Observatory/SotoStarshade.py

@@ -13,6 +13,7 @@
     import cPickle as pickle
 except:
     import pickle
+import sys
 
 EPS = np.finfo(float).eps
 
@@ -116,11 +117,17 @@ def generate_dVMap(self,TL,old_sInd,sInds,currentTime):
             self.vprint('Cached Starshade dV map file not found at "%s".' % dVpath)
             # looping over all target list and desired slew times to generate dV map
             self.vprint('Starting dV calculations for %s stars.' % TL.nStars)
-            tic = time.clock()
+            if sys.version_info[0] > 2:
+                tic = time.perf_counter()
+            else:
+                tic = time.clock()
             for i in range(len(dt)):
                 dVMap[i,:] = self.impulsiveSlew_dV(dt[i],TL,old_sInd,sInd_sorted,currentTime) #sorted
                 if not i % 5: self.vprint('   [%s / %s] completed.' % (i,len(dt)))
-            toc = time.clock()
+            if sys.version_info[0] > 2:
+                toc = time.perf_counter()
+            else:
+                toc = time.clock()
             B = {'dVMap':dVMap}
             with open(dVpath, 'wb') as ff:
                 pickle.dump(B, ff)

EXOSIMS/OpticalSystem/KasdinBraems.py

@@ -13,7 +13,7 @@ class KasdinBraems(OpticalSystem):
     the model from Kasdin & Braems 2006.
     
     Args:
-        \*\*specs:
+        specs:
             user specified values
     
     """

EXOSIMS/OpticalSystem/Nemati.py

@@ -14,7 +14,7 @@ class Nemati(OpticalSystem):
     the model from Nemati 2014.
     
     Args:
-        \*\*specs:
+        specs:
             user specified values
     
     """

EXOSIMS/PlanetPopulation/DulzPlavchan.py

@@ -20,7 +20,7 @@ class DulzPlavchan(PlanetPopulation):
     If occDataPath is not specified, the nominal Period-Radius table is loaded.
 
     Args:
-        \*\*specs:
+        specs:
             user specified values
 
     Attributes:
@@ -45,8 +45,7 @@ def __init__(self, starMass=1., occDataPath=None, esigma=0.175/np.sqrt(np.pi/2.)
         self.occDataPath = occDataPath
         self.esigma = float(esigma)
         er = self.erange
-        self.enorm = np.exp(-er[0]**2/(2.*self.esigma**2)) \
-                     - np.exp(-er[1]**2/(2.*self.esigma**2))
+        self.enorm = np.exp(-er[0]**2/(2.*self.esigma**2)) - np.exp(-er[1]**2/(2.*self.esigma**2))
         self.dist_sma_built = None
         self.dist_radius_built = None
         self.dist_albedo_built = None
@@ -372,7 +371,7 @@ def dist_sma(self, a):
         if self.dist_sma_built is None:
             agen, _ = self.gen_sma_radius(int(1e6))
             ar = self.arange.to('AU').value
-            ap, aedges = np.histogram(agen, bins=2000, range=(ar[0], ar[1]), density=True)
+            ap, aedges = np.histogram(agen.to('AU').value, bins=2000, range=(ar[0], ar[1]), density=True)
             aedges = 0.5*(aedges[1:] + aedges[:-1])
             aedges = np.hstack((ar[0], aedges, ar[1]))
             ap = np.hstack((0., ap, 0.))
@@ -401,7 +400,7 @@ def dist_radius(self, Rp):
         if self.dist_radius_built is None:
             _, Rgen = self.gen_sma_radius(int(1e6))
             Rpr = self.Rprange.to('earthRad').value
-            Rpp, Rpedges = np.histogram(Rgen, bins=2000, range=(Rpr[0], Rpr[1]), density=True)
+            Rpp, Rpedges = np.histogram(Rgen.to('earthRad').value, bins=2000, range=(Rpr[0], Rpr[1]), density=True)
             Rpedges = 0.5*(Rpedges[1:] + Rpedges[:-1])
             Rpedges = np.hstack((Rpr[0], Rpedges, Rpr[1]))
             Rpp = np.hstack((0., Rpp, 0.))
@@ -477,8 +476,7 @@ def dist_eccen_from_sma(self, e, a):
         if a.size not in [1, e.size]:
             elim, e = np.meshgrid(elim, e)
 
-        norm = np.exp(-self.erange[0]**2/(2.*self.esigma**2)) \
-               - np.exp(-elim**2/(2.*self.esigma**2))
+        norm = np.exp(-self.erange[0]**2/(2.*self.esigma**2)) - np.exp(-elim**2/(2.*self.esigma**2))
         ins = np.array((e >= self.erange[0])&(e <= elim), dtype=float, ndmin=1)
         f = np.zeros(e.shape)
         mask = (a >= arcon[0])&(a <= arcon[1])

EXOSIMS/Prototypes/BackgroundSources.py

@@ -11,7 +11,7 @@ class BackgroundSources(object):
     sources for a given target position and dark hole depth.
     
     Args:
-        \*\*specs:
+        specs:
             user specified values
     
     Attributes:

EXOSIMS/Prototypes/Completeness.py

@@ -12,7 +12,7 @@ class Completeness(object):
     Completeness Module calculations in exoplanet mission simulation.
     
     Args:
-        \*\*specs: 
+        specs: 
             user specified values
             
     Attributes:

EXOSIMS/Prototypes/Observatory.py

@@ -32,7 +32,7 @@ class Observatory(object):
     Observatory Definition Module calculations in exoplanet mission simulation.
     
     Args:
-        \*\*specs: 
+        specs: 
             user specified values
         spkpath (str):
             Path to SPK file on disk (Defaults to de432s.bsp). 
@@ -468,8 +468,8 @@ def keepout(self, TL, sInds, currentTime, koangles, returnExtra=False):
         r_targ = (r_targ - r_obs.reshape(nTimes,     1,3)  ).to('pc')  # (m  x n x 3)
         r_body = (r_body - r_obs.reshape(     1,nTimes,3)  ).to('AU')  # (11 x m x 3)
         # unit vectors wrt spacecraft
-        u_targ = (r_targ.value/np.linalg.norm(r_targ, axis=-1, keepdims=True)) 
-        u_body = (r_body.value/np.linalg.norm(r_body, axis=-1, keepdims=True))
+        u_targ = (r_targ/np.linalg.norm(r_targ, axis=-1, keepdims=True)).value
+        u_body = (r_body/np.linalg.norm(r_body, axis=-1, keepdims=True)).value
 
         # create array of koangles for all bodies, using minimum and maximum keepout
         # angles of each starlight suppression system in the telescope for
@@ -767,10 +767,10 @@ def star_angularSep(self,TL,old_sInd,sInds,currentTime):
         else:
             # position vector of previous target star
             r_old = TL.starprop(old_sInd, currentTime)[0]
-            u_old = r_old.value/np.linalg.norm(r_old)
+            u_old = r_old.to('AU').value/np.linalg.norm(r_old.to('AU').value)
             # position vector of new target stars
             r_new = TL.starprop(sInds, currentTime)
-            u_new = (r_new.value.T/np.linalg.norm(r_new, axis=1)).T
+            u_new = (r_new.to('AU').value.T/np.linalg.norm(r_new.to('AU').value, axis=1)).T
             # angle between old and new stars
             sd = np.arccos(np.clip(np.dot(u_old, u_new.T), -1, 1))*u.rad
 
@@ -1119,27 +1119,27 @@ def distForces(self, TL, sInd, currentTime):
         r_Es = self.solarSystem_body_position(currentTime, 'Earth')[0]
         # Telescope -> target vector and unit vector
         r_targ = TL.starprop(sInd, currentTime)[0] - r_obs
-        u_targ = r_targ.value/np.linalg.norm(r_targ)
+        u_targ = r_targ.to('AU').value/np.linalg.norm(r_targ.to('AU').value)
         # sun -> occulter vector
-        r_Os = r_obs + self.occulterSep*u_targ
+        r_Os = r_obs.to('AU') + self.occulterSep.to('AU')*u_targ
         # Earth-Moon barycenter -> spacecraft vectors
         r_TE = r_obs - r_Es
         r_OE = r_Os - r_Es
         # force on occulter
         Mfactor = -self.scMass*const.M_sun*const.G
-        F_sO = r_Os/(np.linalg.norm(r_Os)*r_Os.unit)**3. * Mfactor
-        F_EO = r_OE/(np.linalg.norm(r_OE)*r_OE.unit)**3. * Mfactor/328900.56
+        F_sO = r_Os/(np.linalg.norm(r_Os.to('AU').value)*r_Os.unit)**3. * Mfactor
+        F_EO = r_OE/(np.linalg.norm(r_OE.to('AU').value)*r_OE.unit)**3. * Mfactor/328900.56
         F_O = F_sO + F_EO
         # force on telescope
         Mfactor = -self.coMass*const.M_sun*const.G
-        F_sT = r_obs/(np.linalg.norm(r_obs)*r_obs.unit)**3. * Mfactor
-        F_ET = r_TE/(np.linalg.norm(r_TE)*r_TE.unit)**3. * Mfactor/328900.56
+        F_sT = r_obs/(np.linalg.norm(r_obs.to('AU').value)*r_obs.unit)**3. * Mfactor
+        F_ET = r_TE/(np.linalg.norm(r_TE.to('AU').value)*r_TE.unit)**3. * Mfactor/328900.56
         F_T = F_sT + F_ET
         # differential forces
         dF = F_O - F_T*self.scMass/self.coMass
-        dF_axial = dF.dot(u_targ).to('N')
+        dF_axial = (dF.dot(u_targ)).to('N')
         dF_lateral = (dF - dF_axial*u_targ).to('N')
-        dF_lateral = np.linalg.norm(dF_lateral)*dF_lateral.unit
+        dF_lateral = np.linalg.norm(dF_lateral.to('N').value)*dF_lateral.unit
         dF_axial = np.abs(dF_axial)
 
         return dF_lateral, dF_axial

EXOSIMS/Prototypes/OpticalSystem.py

@@ -22,7 +22,7 @@ class OpticalSystem(object):
     simulation.
     
     Args:
-        \*\*specs:
+        specs:
             User specified values.
             
     Attributes:

EXOSIMS/Prototypes/PlanetPhysicalModel.py

@@ -13,7 +13,7 @@ class PlanetPhysicalModel(object):
     Planet Physical Model Module calculations in exoplanet mission simulation.
     
     Args:
-        \*\*specs:
+        specs:
             user specified values
 
     Attributes:

EXOSIMS/Prototypes/PlanetPopulation.py

@@ -14,7 +14,7 @@ class PlanetPopulation(object):
     simulation.
     
     Args:
-        \*\*specs:
+        specs:
             user specified values
             
     Attributes:

EXOSIMS/Prototypes/PostProcessing.py

@@ -21,7 +21,7 @@ class PostProcessing(object):
     Post Processing Module calculations in exoplanet mission simulation.
     
     Args:
-        \*\*specs:
+        specs:
             user specified values
             
     Attributes:

EXOSIMS/Prototypes/SimulatedUniverse.py

@@ -7,6 +7,7 @@
 import numpy as np
 import astropy.units as u
 import astropy.constants as const
+import sys
 
 class SimulatedUniverse(object):
     """Simulated Universe class template
@@ -15,7 +16,7 @@ class SimulatedUniverse(object):
     Simulated Universe Module calculations in exoplanet mission simulation.
     
     Args:
-        \*\*specs:
+        specs:
             user specified values
     
     Attributes:
@@ -272,12 +273,45 @@ def init_systems(self):
 
         self.r = (A*r1 + B*r2).T.to('AU')                           # position
         self.v = (v1*(-A*r2 + B*v2)).T.to('AU/day')                 # velocity
-        self.d = np.linalg.norm(self.r, axis=1)*self.r.unit         # planet-star distance
-        self.s = np.linalg.norm(self.r[:,0:2], axis=1)*self.r.unit  # apparent separation
-        self.phi = PPMod.calc_Phi(np.arccos(self.r[:,2]/self.d))    # planet phase
+        self.s = np.linalg.norm(self.r[:,0:2], axis=1)              # apparent separation
+        self.d = np.linalg.norm(self.r, axis=1)                     # planet-star distance
+        try:
+            self.phi = PPMod.calc_Phi(np.arccos(self.r[:,2]/self.d))    # planet phase
+        except:
+            self.d = self.d*self.r.unit                    # planet-star distance
+            self.phi = PPMod.calc_Phi(np.arccos(self.r[:,2]/self.d))    # planet phase
+
+        # assert np.arccos(self.r[:,2]/self.d), "d and r do not have same unit >2.7"
+        # assert self.s.unit == TL.dist[0].to('AU').unit, "s and TL.dist do not have same unit >2.7"
+        # if sys.version_info[0] > 2:
+        #     self.s = np.linalg.norm(self.r[:,0:2], axis=1)#*self.r.unit          # apparent separation
+        #     self.d = np.linalg.norm(self.r, axis=1)#*self.r.unit                 # planet-star distance
+        #     assert self.d.unit == self.r.unit, "d and r do not have same unit >2.7"
+        #     assert self.s.unit == TL.dist[0].to('AU').unit, "s and TL.dist do not have same unit >2.7"
+        # else:
+        #     self.d = np.linalg.norm(self.r, axis=1)*self.r.unit         # planet-star distance #must have for 2.7
+        #     self.s = np.linalg.norm(self.r[:,0:2], axis=1)*self.r.unit  # apparent separation
+        #     assert self.d.unit == self.r.unit, "d and r do not have same unit 2.7"
+        #     assert self.s.unit == TL.dist[0].to('AU').unit, "s and TL.dist do not have same unit 2.7"
+        # try:
+        #     self.s.to('AU')
+        #     self.d.to('AU')
+        #     assert self.d.unit == self.r.unit, "d and r do not have same unit 1"
+        #     assert self.s.unit == TL.dist[0].to('AU').unit, "s and TL.dist do not have same unit 1"
+        # except:
+        #     self.s = np.linalg.norm(self.r[:,0:2], axis=1)*self.r.unit              # apparent separation
+        #     self.d = np.linalg.norm(self.r, axis=1)*self.r.unit                     # planet-star distance
+        #     assert self.d.unit == self.r.unit, "d and r do not have same unit 2"
+        #     assert self.s.unit == TL.dist[0].to('AU').unit, "s and TL.dist do not have same unit 2"
+
+        #self.phi = PPMod.calc_Phi(np.arccos(self.r[:,2]/self.d))    # planet phase
         self.fEZ = ZL.fEZ(TL.MV[self.plan2star], self.I, self.d)    # exozodi brightness
         self.dMag = deltaMag(self.p, self.Rp, self.d, self.phi)     # delta magnitude
-        self.WA = np.arctan(self.s/TL.dist[self.plan2star]).to('arcsec')# working angle
+        try:
+            self.WA = np.arctan(self.s/TL.dist[self.plan2star]).to('arcsec')# working angle
+        except:
+            self.s = self.s*self.r.unit
+            self.WA = np.arctan(self.s/TL.dist[self.plan2star]).to('arcsec')# working angle
 
     def propag_system(self, sInd, dt):
         """Propagates planet time-dependant parameters: position, velocity, 
@@ -345,13 +379,29 @@ def propag_system(self, sInd, dt):
         # phase function, exozodi surface brightness, delta magnitude and working angle
         self.r[pInds] = x1[rind]*u.AU
         self.v[pInds] = x1[vind]*u.AU/u.day
-        self.d[pInds] = np.linalg.norm(self.r[pInds], axis=1)*self.r.unit
-        self.s[pInds] = np.linalg.norm(self.r[pInds,0:2], axis=1)*self.r.unit
-        self.phi[pInds] = PPMod.calc_Phi(np.arccos(self.r[pInds,2]/self.d[pInds]))
+        # if sys.version_info[0] > 2:
+        #     self.d[pInds] = np.linalg.norm(self.r[pInds], axis=1)
+        #     self.s[pInds] = np.linalg.norm(self.r[pInds,0:2], axis=1)
+        # else:
+        #     self.d[pInds] = np.linalg.norm(self.r[pInds], axis=1)*self.r.unit
+        #     self.s[pInds] = np.linalg.norm(self.r[pInds,0:2], axis=1)*self.r.unit
+
+        try:
+            self.d[pInds] = np.linalg.norm(self.r[pInds], axis=1)
+            self.phi[pInds] = PPMod.calc_Phi(np.arccos(self.r[pInds,2]/self.d[pInds]))
+        except:
+            self.d[pInds] = np.linalg.norm(self.r[pInds], axis=1)*self.r.unit
+            self.phi[pInds] = PPMod.calc_Phi(np.arccos(self.r[pInds,2]/self.d[pInds]))
+
         # self.fEZ[pInds] = ZL.fEZ(TL.MV[sInd], self.I[pInds], self.d[pInds])
         self.dMag[pInds] = deltaMag(self.p[pInds], self.Rp[pInds], self.d[pInds],
                 self.phi[pInds])
-        self.WA[pInds] = np.arctan(self.s[pInds]/TL.dist[sInd]).to('arcsec')
+        try:
+            self.s[pInds] = np.linalg.norm(self.r[pInds,0:2], axis=1)
+            self.WA[pInds] = np.arctan(self.s[pInds]/TL.dist[sInd]).to('arcsec')
+        except:
+            self.s[pInds] = np.linalg.norm(self.r[pInds,0:2], axis=1)*self.r.unit
+            self.WA[pInds] = np.arctan(self.s[pInds]/TL.dist[sInd]).to('arcsec')
 
     def set_planet_phase(self, beta=np.pi/2):
         """Positions all planets at input star-planet-observer phase angle
@@ -412,12 +462,29 @@ def set_planet_phase(self, beta=np.pi/2):
 
         self.r = (A*r*np.cos(nu) + B*r*np.sin(nu)).T*u.AU           # position
         self.v = (A*v1 + B*v2).T.to('AU/day')                       # velocity
-        self.d = np.linalg.norm(self.r, axis=1)*self.r.unit         # planet-star distance
-        self.s = np.linalg.norm(self.r[:,0:2], axis=1)*self.r.unit  # apparent separation
-        self.phi = PPMod.calc_Phi(np.arccos(self.r[:,2]/self.d))    # planet phase
+        # if sys.version_info[0] > 2:
+        #     self.d = np.linalg.norm(self.r, axis=1)    # planet-star distance
+        #     self.s = np.linalg.norm(self.r[:,0:2], axis=1)  # apparent separation
+        # else:
+        #     self.d = np.linalg.norm(self.r, axis=1)*self.r.unit  # planet-star distance
+        #     self.s = np.linalg.norm(self.r[:,0:2], axis=1)*self.r.unit  # apparent separation   
+
+        try:
+            self.d = np.linalg.norm(self.r, axis=1)    # planet-star distance       
+            self.phi = PPMod.calc_Phi(np.arccos(self.r[:,2].to('AU').value/self.d.to('AU').value)*u.rad)    # planet phase
+        except:
+            self.d = np.linalg.norm(self.r, axis=1)*self.r.unit    # planet-star distance       
+            self.phi = PPMod.calc_Phi(np.arccos(self.r[:,2].to('AU').value/self.d.to('AU').value)*u.rad)    # planet phase
+
         self.fEZ = ZL.fEZ(TL.MV[self.plan2star], self.I, self.d)    # exozodi brightness
         self.dMag = deltaMag(self.p, self.Rp, self.d, self.phi)     # delta magnitude
-        self.WA = np.arctan(self.s/TL.dist[self.plan2star]).to('arcsec')# working angle
+
+        try:
+            self.s = np.linalg.norm(self.r[:,0:2], axis=1)  # apparent separation
+            self.WA = np.arctan(self.s/TL.dist[self.plan2star]).to('arcsec')# working angle
+        except:
+            self.s = np.linalg.norm(self.r[:,0:2], axis=1)*self.r.unit  # apparent separation
+            self.WA = np.arctan(self.s/TL.dist[self.plan2star]).to('arcsec')# working angle
 
     def dump_systems(self):
         """Create a dictionary of planetary properties for archiving use.

EXOSIMS/Prototypes/SurveyEnsemble.py

@@ -6,7 +6,7 @@ class SurveyEnsemble(object):
     """Survey Ensemble prototype
     
     Args:
-        \*\*specs:
+        specs:
             user specified values
             
     Attributes:

EXOSIMS/Prototypes/SurveySimulation.py

@@ -29,7 +29,7 @@ class SurveySimulation(object):
     Simulated Universe, Observatory, TimeKeeping, PostProcessing
     
     Args:
-        \*\*specs:
+        specs:
             user specified values
         scriptfile (string):
             JSON script file.  If not set, assumes that dictionary has been