new keepout method and koMap implementations

EXOSIMS/Prototypes/OpticalSystem.py

@@ -60,6 +60,14 @@ class OpticalSystem(object):
             Mission observing modes attributes
         cachedir (str):
             Path to EXOSIMS cache directory
+        koAngles_Sun (astropy Quantity):
+            Telescope minimum and maximum keepout angle in units of deg
+        koAngles_Earth (astropy Quantity):
+            Telescope minimum and maximum keepout angle in units of deg, for the Earth only
+        koAngles_Moon (astropy Quantity):
+            Telescope minimum and maximum keepout angle in units of deg, for the Moon only
+        koAngles_Small (astropy Quantity):
+            Telescope minimum and maximum keepout angle (for small bodies) in units of deg
 
     Common science instrument attributes:
         name (string):
@@ -196,6 +204,7 @@ def __init__(self, obscurFac=0.1, shapeFac=np.pi/4, pupilDiam=4, intCutoff=50,
             core_thruput=0.1, core_contrast=1e-10, core_platescale=None, 
             PSF=np.ones((3,3)), ohTime=1, observingModes=None, SNR=5, timeMultiplier=1., 
             IWA=None, OWA=None, ref_dMag=3, ref_Time=0, cachedir=None,
+            koAngles_Sun=[0,180], koAngles_Earth=[0,180], koAngles_Moon=[0,180], koAngles_Small=[0,180],
             use_char_minintTime=False, **specs):
 
         #start the outspec
@@ -334,6 +343,12 @@ def __init__(self, obscurFac=0.1, shapeFac=np.pi/4, pupilDiam=4, intCutoff=50,
             if nsyst == 0:
                 lam, BW = syst.get('lam').value, syst.get('BW')
 
+            # get keepout angles for specific instrument
+            syst['koAngles_Sun']   = [float(x) for x in syst.get('koAngles_Sun',  koAngles_Sun)]*u.deg
+            syst['koAngles_Earth'] = [float(x) for x in syst.get('koAngles_Earth',koAngles_Earth)]*u.deg
+            syst['koAngles_Moon']  = [float(x) for x in syst.get('koAngles_Moon', koAngles_Moon)]*u.deg
+            syst['koAngles_Small'] = [float(x) for x in syst.get('koAngles_Small',koAngles_Small)]*u.deg
+
             # get coronagraph input parameters
             syst = self.get_coro_param(syst, 'occ_trans')
             syst = self.get_coro_param(syst, 'core_thruput')

EXOSIMS/Prototypes/SurveySimulation.py

@@ -290,16 +290,34 @@ def __init__(self, scriptfile=None, ntFlux=1, nVisitsMax=5, charMargin=0.15,
         #Generate File Hashnames and loction
         self.cachefname = self.generateHashfName(specs)
 
-        # getting keepout map for entire mission
-        startTime = self.TimeKeeping.missionStart.copy()
-        endTime   = self.TimeKeeping.missionFinishAbs.copy()
-        if not(nokoMap):
-            self.koMap,self.koTimes = self.Observatory.generate_koMap(TL,startTime,endTime)
-
         # choose observing modes selected for detection (default marked with a flag)
         allModes = OS.observingModes
         det_mode = list(filter(lambda mode: mode['detectionMode'] == True, allModes))[0]
         self.mode = det_mode
+
+        # getting keepout map for entire mission
+        startTime = self.TimeKeeping.missionStart.copy()
+        endTime   = self.TimeKeeping.missionFinishAbs.copy()
+
+        nSystems  = len(allModes)
+        systNames = np.unique([allModes[x]['syst']['name'] for x in np.arange(nSystems)]).tolist()
+        koStr     = list(filter(lambda syst: syst.startswith('koAngles_') , allModes[0]['syst'].keys()))
+        koangles  = np.zeros([len(systNames),4,2])
+        tmpNames  = list(systNames)
+        cnt = 0
+
+        for x in np.arange(nSystems):
+            name = allModes[x]['syst']['name']
+            if name in tmpNames:
+                koangles[cnt] = np.asarray([allModes[x]['syst'][k] for k in koStr])
+                cnt += 1
+                tmpNames.remove(name)
+
+        if not(nokoMap):
+            koMaps,self.koTimes = self.Observatory.generate_koMap(TL,startTime,endTime,koangles)
+            self.koMaps = {}
+            for x,n in enumerate(systNames):
+                self.koMaps[n] = koMaps[x,:,:]
 
         # Precalculating intTimeFilter
         sInds = np.arange(TL.nStars) #Initialize some sInds array
@@ -471,7 +489,7 @@ def run_sim(self):
                     self.vprint('waitTime is not None')
                 else:
                     startTimes = TK.currentTimeAbs.copy() + np.zeros(TL.nStars)*u.d # Start Times of Observations
-                    observableTimes = Obs.calculate_observableTimes(TL,np.arange(TL.nStars),startTimes,self.koMap,self.koTimes,self.mode)[0]
+                    observableTimes = Obs.calculate_observableTimes(TL,np.arange(TL.nStars),startTimes,self.koMaps,self.koTimes,self.mode)[0]
                     #CASE 2 If There are no observable targets for the rest of the mission
                     if((observableTimes[(TK.missionFinishAbs.copy().value*u.d > observableTimes.value*u.d)*(observableTimes.value*u.d >= TK.currentTimeAbs.copy().value*u.d)].shape[0]) == 0):#Are there any stars coming out of keepout before end of mission
                         self.vprint('No Observable Targets for Remainder of mission at currentTimeNorm= ' + str(TK.currentTimeNorm.copy()))
@@ -550,7 +568,9 @@ def next_target(self, old_sInd, mode):
         # allocate settling time + overhead time
         tmpCurrentTimeAbs = TK.currentTimeAbs.copy() + Obs.settlingTime + mode['syst']['ohTime']
         tmpCurrentTimeNorm = TK.currentTimeNorm.copy() + Obs.settlingTime + mode['syst']['ohTime']
-
+
+        #create appropriate koMap
+        koMap = self.koMaps[mode['syst']['name']]
 
         # look for available targets
         # 1. initialize arrays
@@ -566,7 +586,7 @@ def next_target(self, old_sInd, mode):
         sd = None
         if OS.haveOcculter == True:
             sd        = Obs.star_angularSep(TL, old_sInd, sInds, tmpCurrentTimeAbs)
-            obsTimes  = Obs.calculate_observableTimes(TL,sInds,tmpCurrentTimeAbs,self.koMap,self.koTimes,mode)
+            obsTimes  = Obs.calculate_observableTimes(TL,sInds,tmpCurrentTimeAbs,self.koMaps,self.koTimes,mode)
             slewTimes = Obs.calculate_slewTimes(TL, old_sInd, sInds, sd, obsTimes, tmpCurrentTimeAbs)  
 
         # 2.1 filter out totTimes > integration cutoff
@@ -582,7 +602,7 @@ def next_target(self, old_sInd, mode):
             tmpIndsbool = list()
             for i in np.arange(len(sInds)):
                 koTimeInd = np.where(np.round(startTimes[sInds[i]].value)-self.koTimes.value==0)[0][0] # find indice where koTime is startTime[0]
-                tmpIndsbool.append(self.koMap[sInds[i]][koTimeInd].astype(bool)) #Is star observable at time ind
+                tmpIndsbool.append(koMap[sInds[i]][koTimeInd].astype(bool)) #Is star observable at time ind
             sInds = sInds[tmpIndsbool]
             del tmpIndsbool
         except:#If there are no target stars to observe 
@@ -603,7 +623,7 @@ def next_target(self, old_sInd, mode):
             else:                
                 intTimes[sInds] = self.calc_targ_intTime(sInds, startTimes[sInds], mode)
                 sInds = sInds[np.where(intTimes[sInds] <= maxIntTime)]  # Filters targets exceeding end of OB
-                endTimes = startTimes + intTimes
+                endTimes = tmpCurrentTimeAbs.copy() + intTimes
 
                 if maxIntTime.value <= 0:
                     sInds = np.asarray([],dtype=int)
@@ -617,7 +637,7 @@ def next_target(self, old_sInd, mode):
                 tmpIndsbool = list()
                 for i in np.arange(len(sInds)):
                     koTimeInd = np.where(np.round(endTimes[sInds[i]].value)-self.koTimes.value==0)[0][0] # find indice where koTime is endTime[0]
-                    tmpIndsbool.append(self.koMap[sInds[i]][koTimeInd].astype(bool)) #Is star observable at time ind
+                    tmpIndsbool.append(koMap[sInds[i]][koTimeInd].astype(bool)) #Is star observable at time ind
                 sInds = sInds[tmpIndsbool]
                 del tmpIndsbool
             except:
@@ -1420,8 +1440,11 @@ def observation_characterization(self, sInd, mode):
             startTime = TK.currentTimeAbs.copy() + mode['syst']['ohTime'] + Obs.settlingTime
             startTimeNorm = TK.currentTimeNorm.copy() + mode['syst']['ohTime'] + Obs.settlingTime
             # planets to characterize
-            tochar[tochar] = Obs.keepout(TL, sInd, startTime)
-
+            koTimeInd = np.where(np.round(startTime.value)-self.koTimes.value==0)[0][0]  # find indice where koTime is startTime[0]
+            #wherever koMap is 1, the target is observable
+            koMap = self.koMaps[mode['syst']['name']]
+            tochar[tochar] = koMap[sInd][koTimeInd]
+
         # 2/ if any planet to characterize, find the characterization times
         # at the detected fEZ, dMag, and WA
         if np.any(tochar):
@@ -1441,9 +1464,16 @@ def observation_characterization(self, sInd, mode):
             # planets to characterize
             tochar = ((totTimes > 0) & (totTimes <= OS.intCutoff) & 
                     (endTimesNorm <= TK.OBendTimes[TK.OBnumber]))
+
         # 3/ is target still observable at the end of any char time?
         if np.any(tochar) and Obs.checkKeepoutEnd:
-            tochar[tochar] = Obs.keepout(TL, sInd, endTimes[tochar])
+            koTimeInds = np.zeros(len(endTimes.value),dtype=int)
+            for t,endTime in enumerate(endTimes.value):
+                if endTime > self.koTimes.value[-1]:
+                    koTimeInds[t] = np.where(np.floor(endTime)-self.koTimes.value==0)[0][0]
+                else:
+                    koTimeInds[t] = np.where(np.round(endTime)-self.koTimes.value==0)[0][0]  # find indice where koTime is endTimes[0]
+            tochar[tochar] = koMap[sInd][koTimeInds]
 
         # 4/ if yes, allocate the overhead time, and perform the characterization 
         # for the maximum char time

EXOSIMS/Prototypes/TargetList.py

@@ -666,27 +666,36 @@ def starprop(self, sInds, currentTime, eclip=False):
                 False, corresponding to heliocentric equatorial frame.
 
         Returns:
-            astropy Quantity nx3 array: 
+            r_targ (astropy Quantity array): 
                 Target star positions vector in heliocentric equatorial (default)
-                or ecliptic frame in units of pc
+                or ecliptic frame in units of pc. Will return an m x n x 3 array 
+                where m is size of currentTime, n is size of sInds. If either m or 
+                n is 1, will return n x 3 or m x 3. 
 
         Note: Use eclip=True to get ecliptic coordinates.
 
         """
 
+        # if multiple time values, check they are different otherwise reduce to scalar
+        if currentTime.size > 1:
+            if np.all(currentTime == currentTime[0]):
+                currentTime = currentTime[0]
+
         # cast sInds to array
         sInds = np.array(sInds, ndmin=1, copy=False)
 
-        # if the starprop_static method was created (staticStars is True), then use it
-        if self.starprop_static is not None:
-            return self.starprop_static(sInds, currentTime, eclip)
-
         # get all array sizes
         nStars = sInds.size
         nTimes = currentTime.size
-        assert nStars==1 or nTimes==1 or nTimes==nStars, \
-                "If multiple times and targets, currentTime and sInds sizes must match"
-
+
+        # if the starprop_static method was created (staticStars is True), then use it
+        if self.starprop_static is not None:
+            r_targ = self.starprop_static(sInds, currentTime, eclip)
+            if (nTimes == 1 or nStars == 1 or nTimes == nStars):
+                return r_targ
+            else:
+                return np.tile(r_targ, (nTimes, 1, 1))
+
         # target star ICRS coordinates
         coord_old = self.coords[sInds]
         # right ascension and declination
@@ -702,17 +711,34 @@ def starprop(self, sInds, currentTime, eclip=False):
         v = mu0/self.parx[sInds]*u.AU + r0*self.rv[sInds]
         # set J2000 epoch
         j2000 = Time(2000., format='jyear')
-        # target star positions vector in heliocentric equatorial frame
-        dr = v*(currentTime.mjd - j2000.mjd)*u.day
-        r_targ = (coord_old.cartesian.xyz + dr).T.to('pc')
-
-        if eclip:
-            # transform to heliocentric true ecliptic frame
-            coord_new = SkyCoord(r_targ[:,0], r_targ[:,1], r_targ[:,2], 
-                    representation='cartesian')
-            r_targ = coord_new.heliocentrictrueecliptic.cartesian.xyz.T.to('pc')
-
-        return r_targ
+
+        # if only 1 time in currentTime
+        if (nTimes == 1 or nStars == 1 or nTimes == nStars):
+            # target star positions vector in heliocentric equatorial frame
+            dr = v*(currentTime.mjd - j2000.mjd)*u.day
+            r_targ = (coord_old.cartesian.xyz + dr).T.to('pc')
+
+            if eclip:
+                # transform to heliocentric true ecliptic frame
+                coord_new = SkyCoord(r_targ[:,0], r_targ[:,1], r_targ[:,2], 
+                            representation='cartesian')
+                r_targ = coord_new.heliocentrictrueecliptic.cartesian.xyz.T.to('pc')
+            return r_targ
+
+        # create multi-dimensional array for r_targ
+        else:
+            # target star positions vector in heliocentric equatorial frame
+            r_targ = np.zeros([nTimes,nStars,3])*u.pc
+            for i,m in enumerate(currentTime):
+                 dr = v*(m.mjd - j2000.mjd)*u.day
+                 r_targ[i,:,:] = (coord_old.cartesian.xyz + dr).T.to('pc')
+
+            if eclip:
+                # transform to heliocentric true ecliptic frame
+                coord_new = SkyCoord(r_targ[i,:,0], r_targ[i,:,1], r_targ[i,:,2], 
+                            representation='cartesian')
+                r_targ[i,:,:] = coord_new.heliocentrictrueecliptic.cartesian.xyz.T.to('pc')
+            return r_targ
 
     def starMag(self, sInds, lam):
         """Calculates star visual magnitudes with B-V color using empirical fit 

EXOSIMS/SurveySimulation/ExoC_Scheduler.py

@@ -173,7 +173,7 @@ def run_sim(self):
                     self.vprint('waitTime is not None')
                 else:
                     startTimes = TK.currentTimeAbs.copy() + np.zeros(TL.nStars)*u.d # Start Times of Observations
-                    observableTimes = Obs.calculate_observableTimes(TL,np.arange(TL.nStars),startTimes,self.koMap,self.koTimes,self.mode)[0]
+                    observableTimes = Obs.calculate_observableTimes(TL,np.arange(TL.nStars),startTimes,self.koMaps,self.koTimes,self.mode)[0]
                     #CASE 2 If There are no observable targets for the rest of the mission
                     if((observableTimes[(TK.missionFinishAbs.copy().value*u.d > observableTimes.value*u.d)*(observableTimes.value*u.d >= TK.currentTimeAbs.copy().value*u.d)].shape[0]) == 0):#Are there any stars coming out of keepout before end of mission
                         self.vprint('No Observable Targets for Remainder of mission at currentTimeNorm= ' + str(TK.currentTimeNorm.copy()))
@@ -271,7 +271,10 @@ def observation_characterization(self, sInd, mode, mode_index):
             startTime = TK.currentTimeAbs.copy() + mode['syst']['ohTime'] + Obs.settlingTime
             startTimeNorm = TK.currentTimeNorm.copy() + mode['syst']['ohTime'] + Obs.settlingTime
             # planets to characterize
-            tochar[tochar] = Obs.keepout(TL, sInd, startTime)
+            koTimeInd = np.where(np.round(startTime.value)-self.koTimes.value==0)[0][0]  # find indice where koTime is startTime[0]
+            #wherever koMap is 1, the target is observable
+            koMap = self.koMaps[mode['syst']['name']]
+            tochar[tochar] = koMap[sInd][koTimeInd]
 
         # 2/ if any planet to characterize, find the characterization times
         # at the detected fEZ, dMag, and WA
@@ -294,7 +297,13 @@ def observation_characterization(self, sInd, mode, mode_index):
                     (endTimesNorm <= TK.OBendTimes[TK.OBnumber]))
         # 3/ is target still observable at the end of any char time?
         if np.any(tochar) and Obs.checkKeepoutEnd:
-            tochar[tochar] = Obs.keepout(TL, sInd, endTimes[tochar])
+            koTimeInds = np.zeros(len(endTimes.value),dtype=int)
+            for t,endTime in enumerate(endTimes.value):
+                if endTime > self.koTimes.value[-1]:
+                    koTimeInds[t] = np.where(np.floor(endTime)-self.koTimes.value==0)[0][0]
+                else:
+                    koTimeInds[t] = np.where(np.round(endTime)-self.koTimes.value==0)[0][0]  # find indice where koTime is endTimes[0]
+            tochar[tochar] = koMap[sInd][koTimeInds]
 
         # 4/ if yes, allocate the overhead time, and perform the characterization 
         # for the maximum char time

EXOSIMS/SurveySimulation/SS_char_only.py

@@ -286,6 +286,9 @@ def observation_characterization(self, sInd, mode):
         Obs = self.Observatory
         TK = self.TimeKeeping
 
+        # selecting appropriate koMap
+        koMap = self.koMaps[mode['syst']['name']]
+
         # find indices of planets around the target
         pInds = np.where(SU.plan2star == sInd)[0]
         # get the last detected planets, and check if there was a FA
@@ -332,7 +335,10 @@ def observation_characterization(self, sInd, mode):
             startTime = TK.currentTimeAbs
             startTimeNorm = TK.currentTimeNorm
             # planets to characterize
-            tochar[tochar] = Obs.keepout(TL, sInd, startTime, mode)
+            koTimeInd = np.where(np.round(startTime.value)-self.koTimes.value==0)[0][0]  # find indice where koTime is startTime[0]
+            #wherever koMap is 1, the target is observable
+            tochar[tochar] = koMap[sInd][koTimeInd]
+
 
         # 2/ if any planet to characterize, find the characterization times
         if np.any(tochar):
@@ -361,7 +367,13 @@ def observation_characterization(self, sInd, mode):
 
         # 3/ is target still observable at the end of any char time?
         if np.any(tochar) and Obs.checkKeepoutEnd:
-            tochar[tochar] = Obs.keepout(TL, sInd, endTimes[tochar], mode)
+            koTimeInds = np.zeros(len(endTimes.value),dtype=int)
+            for t,endTime in enumerate(endTimes.value):
+                if endTime > self.koTimes.value[-1]:
+                    koTimeInds[t] = np.where(np.floor(endTime)-self.koTimes.value==0)[0][0]
+                else:
+                    koTimeInds[t] = np.where(np.round(endTime)-self.koTimes.value==0)[0][0]  # find indice where koTime is endTimes[0]
+            tochar[tochar] = koMap[sInd][koTimeInds]
 
         # 4/ if yes, perform the characterization for the maximum char time
         if np.any(tochar):