Python astropy.units.d() Examples

def revisitFilter(self, sInds, tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering

        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead was added in MJD form
        Returns:
            sInds - indices of stars still in observation list
        """
        tovisit = np.zeros(self.TargetList.nStars, dtype=bool)#tovisit is a boolean array containing the 
        if len(sInds) > 0:#so long as there is at least 1 star left in sInds
            tovisit[sInds] = ((self.starVisits[sInds] == min(self.starVisits[sInds])) \
                    & (self.starVisits[sInds] < self.nVisitsMax))# Checks that no star has exceeded the number of revisits
            if self.starRevisit.size != 0:#There is at least one revisit planned in starRevisit
                dt_rev = self.starRevisit[:,1]*u.day - tmpCurrentTimeNorm#absolute temporal spacing between revisit and now.

                #return indices of all revisits within a threshold dt_max of revisit day and indices of all revisits with no detections past the revisit time
                ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds)]
                tovisit[ind_rev2] = (self.starVisits[ind_rev2] < self.nVisitsMax)
            sInds = np.where(tovisit)[0]

        return sInds 

def revisitFilter(self, sInds, tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering

        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead was added in MJD form
        Returns:
            sInds - indices of stars still in observation list
        """
        tovisit = np.zeros(self.TargetList.nStars, dtype=bool)#tovisit is a boolean array containing the 
        if len(sInds) > 0:#so long as there is at least 1 star left in sInds
            tovisit[sInds] = ((self.starVisits[sInds] == min(self.starVisits[sInds])) \
                    & (self.starVisits[sInds] < self.nVisitsMax))# Checks that no star has exceeded the number of revisits
            if self.starRevisit.size != 0:#There is at least one revisit planned in starRevisit
                dt_rev = self.starRevisit[:,1]*u.day - tmpCurrentTimeNorm#absolute temporal spacing between revisit and now.

                #return indices of all revisits within a threshold dt_max of revisit day and indices of all revisits with no detections past the revisit time
                ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds)]
                tovisit[ind_rev2] = (self.starVisits[ind_rev2] < self.nVisitsMax)
            sInds = np.where(tovisit)[0]

        return sInds 

def __init__(self, coeffs=[1,1,2,1], revisit_wait=91.25, **specs):
        
        SurveySimulation.__init__(self, **specs)
        
        #verify that coefficients input is iterable 4x1
        if not(isinstance(coeffs,(list,tuple,np.ndarray))) or (len(coeffs) != 4):
            raise TypeError("coeffs must be a 4 element iterable")


        #Add to outspec
        self._outspec['coeffs'] = coeffs
        self._outspec['revisit_wait'] = revisit_wait
        
        # normalize coefficients
        coeffs = np.array(coeffs)
        coeffs = coeffs/np.linalg.norm(coeffs)
        
        self.coeffs = coeffs

        self.revisit_wait = revisit_wait*u.d
        self.no_dets = np.ones(self.TargetList.nStars, dtype=bool) 

def revisitFilter(self, sInds, tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering

        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead was added in MJD form
        Returns:
            sInds - indices of stars still in observation list
        """
        tovisit = np.zeros(self.TargetList.nStars, dtype=bool)#tovisit is a boolean array containing the 
        if len(sInds) > 0:#so long as there is at least 1 star left in sInds
            tovisit[sInds] = ((self.starVisits[sInds] == min(self.starVisits[sInds])) \
                    & (self.starVisits[sInds] < self.nVisitsMax))# Checks that no star has exceeded the number of revisits
            if self.starRevisit.size != 0:#There is at least one revisit planned in starRevisit
                dt_rev = self.starRevisit[:,1]*u.day - tmpCurrentTimeNorm#absolute temporal spacing between revisit and now.

                #return indices of all revisits within a threshold dt_max of revisit day and indices of all revisits with no detections past the revisit time
                # ind_rev = [int(x) for x in self.starRevisit[np.abs(dt_rev) < self.dt_max, 0] if (x in sInds and self.no_dets[int(x)] == False)]
                # ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds and self.no_dets[int(x)] == True)]
                # tovisit[ind_rev] = (self.starVisits[ind_rev] < self.nVisitsMax)#IF duplicates exist in ind_rev, the second occurence takes priority
                ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds)]
                tovisit[ind_rev2] = (self.starVisits[ind_rev2] < self.nVisitsMax)
            sInds = np.where(tovisit)[0]

        return sInds 

def revisitFilter(self, sInds, tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering

        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead was added in MJD form
        Returns:
            sInds - indices of stars still in observation list
        """
        tovisit = np.zeros(self.TargetList.nStars, dtype=bool)#tovisit is a boolean array containing the 
        if len(sInds) > 0:#so long as there is at least 1 star left in sInds
            tovisit[sInds] = ((self.starVisits[sInds] == min(self.starVisits[sInds])) \
                    & (self.starVisits[sInds] < self.nVisitsMax))# Checks that no star has exceeded the number of revisits
            if self.starRevisit.size != 0:#There is at least one revisit planned in starRevisit
                dt_rev = self.starRevisit[:,1]*u.day - tmpCurrentTimeNorm#absolute temporal spacing between revisit and now.

                #return indices of all revisits within a threshold dt_max of revisit day and indices of all revisits with no detections past the revisit time
                ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds)]
                tovisit[ind_rev2] = (self.starVisits[ind_rev2] < self.nVisitsMax)
            sInds = np.where(tovisit)[0]

        return sInds 

def revisitFilter(self, sInds, tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering

        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead was added in MJD form
        Returns:
            sInds - indices of stars still in observation list
        """
        tovisit = np.zeros(self.TargetList.nStars, dtype=bool)#tovisit is a boolean array containing the 
        if len(sInds) > 0:#so long as there is at least 1 star left in sInds
            tovisit[sInds] = ((self.starVisits[sInds] == min(self.starVisits[sInds])) \
                    & (self.starVisits[sInds] < self.nVisitsMax))# Checks that no star has exceeded the number of revisits
            if self.starRevisit.size != 0:#There is at least one revisit planned in starRevisit
                dt_rev = self.starRevisit[:,1]*u.day - tmpCurrentTimeNorm#absolute temporal spacing between revisit and now.

                ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds)]
                tovisit[ind_rev2] = (self.starVisits[ind_rev2] < self.nVisitsMax)
            sInds = np.where(tovisit)[0]
        return sInds 

def __init__(self, coeffs=[1,1,2,1], revisit_wait=91.25, **specs):
        
        SurveySimulation.__init__(self, **specs)
        
        #verify that coefficients input is iterable 4x1
        if not(isinstance(coeffs,(list,tuple,np.ndarray))) or (len(coeffs) != 4):
            raise TypeError("coeffs must be a 4 element iterable")


        #Add to outspec
        self._outspec['coeffs'] = coeffs
        self._outspec['revisit_wait'] = revisit_wait
        
        # normalize coefficients
        coeffs = np.array(coeffs)
        coeffs = coeffs/np.linalg.norm(coeffs)
        
        self.coeffs = coeffs

        self.revisit_wait = revisit_wait*u.d
        self.no_dets = np.ones(self.TargetList.nStars, dtype=bool) 

def revisitFilter(self, sInds, tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering

        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead was added in MJD form
        Returns:
            sInds - indices of stars still in observation list
        """
        tovisit = np.zeros(self.TargetList.nStars, dtype=bool)#tovisit is a boolean array containing the 
        if len(sInds) > 0:#so long as there is at least 1 star left in sInds
            tovisit[sInds] = ((self.starVisits[sInds] == min(self.starVisits[sInds])) \
                    & (self.starVisits[sInds] < self.nVisitsMax))# Checks that no star has exceeded the number of revisits
            if self.starRevisit.size != 0:#There is at least one revisit planned in starRevisit
                dt_rev = self.starRevisit[:,1]*u.day - tmpCurrentTimeNorm#absolute temporal spacing between revisit and now.

                ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds)]
                tovisit[ind_rev2] = (self.starVisits[ind_rev2] < self.nVisitsMax)
            sInds = np.where(tovisit)[0]
        return sInds 

def sacrificeStarCbyT(self, sInds, t_dets, fZ, fEZ, WA, overheadTime):
        """Sacrifice the worst performing CbyT star
        Args:
            sInds[nStars] - indicies of stars in the list
            t_dets[nStars] - time to observe each star (in days)
            fZ[nStars] - zodiacal light for each target
            fEZ - 0 
            WA - inner working angle of the instrument
            overheadTime - overheadTime added to each observation
        Return:
            sInds[nStars] - indicies of stars in the list
            t_dets[nStars] - time to observe each star (in days)
            sacrificedStarTime - time to distribute in days       
        """
        CbyT = self.Completeness.comp_per_intTime(t_dets*u.d, self.TargetList, sInds, self.valfZmin[sInds], fEZ, WA, self.mode, self.Cb[sInds], self.Csp[sInds])/t_dets#takes 5 seconds to do 1 time for all stars

        sacrificeIndex = np.argmin(CbyT)#finds index of star to sacrifice

        #Need index of sacrificed star by this point
        sacrificedStarTime = t_dets[sacrificeIndex] + overheadTime#saves time being sacrificed
        sInds = np.delete(sInds,sacrificeIndex)
        t_dets = np.delete(t_dets,sacrificeIndex)
        return sInds, t_dets, sacrificedStarTime 

def objfun(self,t,sInds,fZ):
        """
        Objective Function for SLSQP minimization. Purpose is to maximize summed completeness

        Args:
            t (ndarray):
                Integration times in days. NB: NOT an astropy quantity.
            sInds (ndarray):
                Target star indices (of same size as t)
            fZ (astropy Quantity):
                Surface brightness of local zodiacal light in units of 1/arcsec2
                Same size as t

        """
        good = t*u.d >= 0.1*u.s # inds that were not downselected by initial MIP

        comp = self.Completeness.comp_per_intTime(t[good]*u.d, self.TargetList, sInds[good], fZ[good], 
                self.ZodiacalLight.fEZ0, self.WAint[sInds][good], self.detmode)
        #self.vprint(-comp.sum()) # for verifying SLSQP output
        return -comp.sum() 

def objfun_deriv(self,t,sInds,fZ):
        """
        Jacobian of objective Function for SLSQP minimization. 

        Args:
            t (astropy Quantity):
                Integration times in days. NB: NOT an astropy quantity.
            sInds (ndarray):
                Target star indices (of same size as t)
            fZ (astropy Quantity):
                Surface brightness of local zodiacal light in units of 1/arcsec2
                Same size as t

        """
        good = t*u.d >= 0.1*u.s # inds that were not downselected by initial MIP

        tmp = self.Completeness.dcomp_dt(t[good]*u.d, self.TargetList, sInds[good], fZ[good], 
                self.ZodiacalLight.fEZ0, self.WAint[sInds][good], self.detmode).to("1/d").value

        jac = np.zeros(len(t))
        jac[good] = tmp
        return -jac 

def revisitFilter(self, sInds, tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering

        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead was added in MJD form
        Returns:
            sInds - indices of stars still in observation list
        """
        tovisit = np.zeros(self.TargetList.nStars, dtype=bool)#tovisit is a boolean array containing the 
        if len(sInds) > 0:#so long as there is at least 1 star left in sInds
            tovisit[sInds] = ((self.starVisits[sInds] == min(self.starVisits[sInds])) \
                    & (self.starVisits[sInds] < self.nVisitsMax))# Checks that no star has exceeded the number of revisits
            if self.starRevisit.size != 0:#There is at least one revisit planned in starRevisit
                dt_rev = self.starRevisit[:,1]*u.day - tmpCurrentTimeNorm#absolute temporal spacing between revisit and now.

                #return indices of all revisits within a threshold dt_max of revisit day and indices of all revisits with no detections past the revisit time
                ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds)]
                tovisit[ind_rev2] = (self.starVisits[ind_rev2] < self.nVisitsMax)
            sInds = np.where(tovisit)[0]

        return sInds 

def revisitFilter(self, sInds, tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering

        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead was added in MJD form
        Returns:
            sInds - indices of stars still in observation list
        """
        tovisit = np.zeros(self.TargetList.nStars, dtype=bool)#tovisit is a boolean array containing the 
        if len(sInds) > 0:#so long as there is at least 1 star left in sInds
            tovisit[sInds] = ((self.starVisits[sInds] == min(self.starVisits[sInds])) \
                    & (self.starVisits[sInds] < self.nVisitsMax))# Checks that no star has exceeded the number of revisits
            if self.starRevisit.size != 0:#There is at least one revisit planned in starRevisit
                dt_rev = self.starRevisit[:,1]*u.day - tmpCurrentTimeNorm#absolute temporal spacing between revisit and now.

                #return indices of all revisits within a threshold dt_max of revisit day and indices of all revisits with no detections past the revisit time
                ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds)]
                tovisit[ind_rev2] = (self.starVisits[ind_rev2] < self.nVisitsMax)
            sInds = np.where(tovisit)[0]

        return sInds 

def revisitFilter(self, sInds, tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering
        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead was added in MJD form
        Returns:
            sInds - indices of stars still in observation list
        """
        tovisit = np.zeros(self.TargetList.nStars, dtype=bool)#tovisit is a boolean array containing the 
        if len(sInds) > 0:#so long as there is at least 1 star left in sInds
            tovisit[sInds] = ((self.starVisits[sInds] == min(self.starVisits[sInds])) \
                    & (self.starVisits[sInds] < self.nVisitsMax))# Checks that no star has exceeded the number of revisits
            if self.starRevisit.size != 0:#There is at least one revisit planned in starRevisit
                dt_rev = self.starRevisit[:,1]*u.day - tmpCurrentTimeNorm#absolute temporal spacing between revisit and now.

                #return indices of all revisits within a threshold dt_max of revisit day and indices of all revisits with no detections past the revisit time
                # ind_rev = [int(x) for x in self.starRevisit[np.abs(dt_rev) < self.dt_max, 0] if (x in sInds and self.no_dets[int(x)] == False)]
                # ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds and self.no_dets[int(x)] == True)]
                # tovisit[ind_rev] = (self.starVisits[ind_rev] < self.nVisitsMax)#IF duplicates exist in ind_rev, the second occurence takes priority
                ind_rev2 = [int(x) for x in self.starRevisit[dt_rev < 0*u.d, 0] if (x in sInds)]
                tovisit[ind_rev2] = (self.starVisits[ind_rev2] < self.nVisitsMax)
            sInds = np.where(tovisit)[0]

        return sInds 

def semi_major_axis(P, Mtotal):
    """Semi-major axis

    Kepler's third law

    Args:
        P (float): Orbital period [days]
        Mtotal (float): Mass [Msun]

    Returns:
        float or array: semi-major axis in AU
    """

    # convert inputs to array so they work with units
    P = np.array(P)
    Mtotal = np.array(Mtotal)

    Mtotal = Mtotal*c.M_sun.value
    P = (P * u.d).to(u.second).value
    G = c.G.value
    
    a = ((P**2)*G*Mtotal/(4*(np.pi)**2))**(1/3.)
    a = a/c.au.value

    return a 

def test_timedelta_setitem():
    t = TimeDelta([1, 2, 3] * u.d, format='jd')

    t[0] = 0.5
    assert allclose_jd(t.value, [0.5, 2, 3])

    t[1:] = 4.5
    assert allclose_jd(t.value, [0.5, 4.5, 4.5])

    t[:] = 86400 * u.s
    assert allclose_jd(t.value, [1, 1, 1])

    t[1] = TimeDelta(2, format='jd')
    assert allclose_jd(t.value, [1, 2, 1])

    with pytest.raises(ValueError) as err:
        t[1] = 1 * u.m
    assert 'cannot convert value to a compatible TimeDelta' in str(err.value) 

def _read_leap_seconds(cls, file, **kwargs):
        """Read a file, identifying expiration by matching 'File expires'"""
        expires = None
        # Find expiration date.
        with get_readable_fileobj(file) as fh:
            lines = fh.readlines()
            for line in lines:
                match = cls._re_expires.match(line)
                if match:
                    expires = Time.strptime(match.groups()[0], '%d %B %Y',
                                            scale='tai', format='iso',
                                            out_subfmt='date')
                    break
            else:
                raise ValueError(f'did not find expiration date in {file}')

        self = cls.read(lines, format='ascii.no_header', **kwargs)
        self._expires = expires
        return self 

def tests_comp_per_intTime(self):

        comp = self.fixture.comp_per_intTime(np.ones(self.nStars)*u.d,self,np.arange(self.nStars),0/(u.arcsec**2),0/(u.arcsec**2),0*u.arcsec,{})
        
        self.assertTrue(len(comp),self.nStars)

        with self.assertRaises(AssertionError):
            comp = self.fixture.comp_per_intTime(np.ones(self.nStars-1)*u.d,self,np.arange(self.nStars),0/(u.arcsec**2),0/(u.arcsec**2),0*u.arcsec,{})

        with self.assertRaises(AssertionError):
            comp = self.fixture.comp_per_intTime(np.ones(self.nStars)*u.d,self,np.arange(self.nStars),np.zeros(2)/(u.arcsec**2),0/(u.arcsec**2),0*u.arcsec,{})

        with self.assertRaises(AssertionError):
            comp = self.fixture.comp_per_intTime(np.ones(self.nStars)*u.d,self,np.arange(self.nStars),0/(u.arcsec**2),np.zeros(2)/(u.arcsec**2),0*u.arcsec,{})

        with self.assertRaises(AssertionError):
            comp = self.fixture.comp_per_intTime(np.ones(self.nStars-1)*u.d,self,np.arange(self.nStars),0/(u.arcsec**2),0/(u.arcsec**2),np.zeros(2)*u.arcsec,{}) 

def test_reset_sim(self):
        r"""Test reset_sim method.

        Approach: Ensures the simulation resets completely
        """

        with RedirectStreams(stdout=self.dev_null):
            sim = self.fixture(SimpleScript)
            sim.run_sim()
        
        self.assertGreater(len(sim.DRM), 0)
        self.assertGreater(sim.TimeKeeping.currentTimeNorm,0.0*u.d)

        sim.reset_sim()
        self.assertEqual(sim.TimeKeeping.currentTimeNorm,0.0*u.d)
        self.assertEqual(len(sim.DRM), 0) 

def test_advancetToStartOfNextOB(self):
        r""" Test advancetToStartOfNextOB method
                Strategy is to call the method once and ensure it advances the Observing Block
        """  
        tk = self.fixture()

        tk.OBstartTimes = [0,10,20,30]*u.d
        tk.OBendTimes = [5,15,25,35]*u.d
        tk.OBnumber = 0

        # 1) Set current Time to End of OB
        tk.currentTimeNorm = tk.OBendTimes[0]
        tk.currentTimeAbs = tk.missionStart + tk.currentTimeNorm
        tmpOBnumber = tk.OBnumber
        tmpcurrentTimeNorm = tk.currentTimeNorm.copy()
        tmpcurrentTimeAbs = tk.currentTimeAbs.copy()
        tmpexoplanetObsTime = tk.exoplanetObsTime.copy()
        tmpOBstartTimes = tk.OBstartTimes.copy()
        tmpOBendTimes = tk.OBendTimes.copy()
        tk.advancetToStartOfNextOB()
        self.assertTrue(tmpOBnumber + 1 == tk.OBnumber)
        self.assertTrue(tk.currentTimeNorm == tk.OBstartTimes[tk.OBnumber])
        self.assertTrue(tk.currentTimeAbs == tk.OBstartTimes[tk.OBnumber] + tk.missionStart) 

def test_zeroFAP(self):
        """
        Test case where false alarm probability is 0 and missed det prob is 0
        All observations above SNR limit should be detected
        """

        SNRin = np.array([1.0,2.0,3.0,4.0,5.0,5.1,6.0])
        expected_MDresult = [True,True,True,True,False,False,False] 

        for mod in self.allmods:
            with RedirectStreams(stdout=self.dev_null):
                obj = mod(FAP=0.0,MDP=0.0,**self.specs)

            FA, MD = obj.det_occur(SNRin,self.TL.OpticalSystem.observingModes[0],self.TL,0,1*u.d)     
            for x,y in zip(MD,expected_MDresult):
                self.assertEqual( x,y )
            self.assertEqual( FA, False)

    #another limiting case 

def test_oneFAP(self):
        """
        Test case where false alarm probability is 1 and missed det prob is 0
        """

        SNRin = np.array([1.0,2.0,3.0,4.0,5.0,5.1,6.0])
        expected_MDresult = [True,True,True,True,False,False,False] 

        for mod in self.allmods:
            with RedirectStreams(stdout=self.dev_null):
                obj = mod(FAP=1.0,MDP=0.0,**self.specs)

            FA, MD = obj.det_occur(SNRin,self.TL.OpticalSystem.observingModes[0],self.TL,0,1*u.d)     
            for x,y in zip(MD,expected_MDresult):
                self.assertEqual( x,y )
            self.assertEqual( FA, True) 

def test_observation_detection(self):
        r"""Test observation_detection method.

        Approach: Ensure that all outputs are set as expected
        """

        exclude_mods = []
        exclude_mod_type = 'sotoSS'

        for mod in self.allmods:
            if mod.__name__ in exclude_mods:
                continue
            if 'observation_detection' in mod.__dict__ or exclude_mod_type in mod.__name__:

                spec = copy.deepcopy(self.spec)
                if 'tieredScheduler' in mod.__name__:
                    self.script = resource_path('test-scripts/simplest_occ.json')
                    with open(self.script) as f:
                        spec = json.loads(f.read())
                    spec['occHIPs'] = resource_path('SurveySimulation/top100stars.txt')

                with RedirectStreams(stdout=self.dev_null):
                    sim = mod(**spec)

                    #default settings should create dummy planet around first star
                    sInd = 0
                    pInds = np.where(sim.SimulatedUniverse.plan2star == sInd)[0]
                    detected, fZ, systemParams, SNR, FA = \
                            sim.observation_detection(sInd,1.0*u.d,\
                            sim.OpticalSystem.observingModes[0])
                
                self.assertEqual(len(detected),len(pInds),\
                        'len(detected) != len(pInds) for %s'%mod.__name__)
                self.assertIsInstance(detected[0],(int,np.int32,np.int64,np.int_),\
                        'detected elements not ints for %s'%mod.__name__)
                for s in SNR[detected == 1]:
                    self.assertGreaterEqual(s,sim.OpticalSystem.observingModes[0]['SNR'],\
                        'detection SNR < mode requirement for %s'%mod.__name__)
                self.assertIsInstance(FA, bool,\
                        'False Alarm not boolean for %s'%mod.__name__) 

def initializeStorageArrays(self):
        """
        Initialize all storage arrays based on # of stars and targets
        """

        self.DRM = []
        self.fullSpectra = np.zeros(self.SimulatedUniverse.nPlans, dtype=int)
        self.partialSpectra = np.zeros(self.SimulatedUniverse.nPlans, dtype=int)
        self.propagTimes = np.zeros(self.TargetList.nStars)*u.d
        self.lastObsTimes = np.zeros(self.TargetList.nStars)*u.d
        self.starVisits = np.zeros(self.TargetList.nStars, dtype=int)#contains the number of times each star was visited
        self.starRevisit = np.array([])
        self.starExtended = np.array([], dtype=int)
        self.lastDetected = np.empty((self.TargetList.nStars, 4), dtype=object) 

def get_ObsDetectionMaxIntTime(self,Obs,mode,currentTimeNorm=None,OBnumber=None):
        """Tells you the maximum Detection Observation Integration Time you can pass into observation_detection(X,intTime,X)
        Args:
            Obs (Observatory Object):
                Observatory module for Obs.settlingTime
            mode (dict):
                Selected observing mode for detection

        Returns:
            tuple:
            maxIntTimeOBendTime (astropy Quantity):
                The maximum integration time bounded by Observation Block end Time
            maxIntTimeExoplanetObsTime (astropy Quantity):
                The maximum integration time bounded by exoplanetObsTime
            maxIntTimeMissionLife (astropy Quantity):
                The maximum integration time bounded by MissionLife
        """
        if OBnumber == None:
            OBnumber = self.OBnumber
        if currentTimeNorm == None:
            currentTimeNorm = self.currentTimeNorm.copy()
            
        maxTimeOBendTime = self.OBendTimes[OBnumber] - currentTimeNorm
        maxIntTimeOBendTime = (maxTimeOBendTime - Obs.settlingTime - mode['syst']['ohTime'])/(1. + mode['timeMultiplier'] -1.)

        maxTimeExoplanetObsTime = self.missionLife.to('day')*self.missionPortion - self.exoplanetObsTime
        maxIntTimeExoplanetObsTime = (maxTimeExoplanetObsTime - Obs.settlingTime - mode['syst']['ohTime'])/(1. + mode['timeMultiplier'] -1.)

        maxTimeMissionLife = self.missionLife.to('day') - currentTimeNorm
        maxIntTimeMissionLife = (maxTimeMissionLife - Obs.settlingTime - mode['syst']['ohTime'])/(1. + mode['timeMultiplier'] -1.)

        #Ensure all are positive or zero
        if maxIntTimeOBendTime < 0.:
            maxIntTimeOBendTime = 0.*u.d
        if maxIntTimeExoplanetObsTime < 0.:
            maxIntTimeExoplanetObsTime = 0.*u.d
        if maxIntTimeMissionLife < 0.:
            maxIntTimeMissionLife = 0.*u.d

        return maxIntTimeOBendTime, maxIntTimeExoplanetObsTime, maxIntTimeMissionLife 

def mission_is_over(self, OS, Obs, mode):
        r"""Is the time allocated for the mission used up?
        
        This supplies an abstraction around the test: ::
            
            (currentTimeNorm > missionFinishNorm)
            
        so that users of the class do not have to perform arithmetic
        on class variables.

        Args:
            OS (Optical System object):
                Optical System module for OS.haveOcculter
            Obs (Observatory Object):
                Observatory module for Obs.settlingTime
            mode (dict):
                Selected observing mode for detection (uses only overhead time)

        Returns:
            boolean:
                True if the mission time is used up, else False.
        """
        
        is_over = ((self.currentTimeNorm + Obs.settlingTime + mode['syst']['ohTime'] >= self.missionLife.to('day')) \
            or (self.exoplanetObsTime.to('day') + Obs.settlingTime + mode['syst']['ohTime'] >= self.missionLife.to('day')*self.missionPortion) \
            or (self.currentTimeNorm + Obs.settlingTime + mode['syst']['ohTime'] >= self.OBendTimes[-1]) \
            or (OS.haveOcculter and Obs.scMass < Obs.dryMass))

        if (OS.haveOcculter and Obs.scMass < Obs.dryMass):
            self.vprint('Total fuel mass (Obs.dryMass %.2f) less than (Obs.scMass %.2f) at currentTimeNorm %sd'%(Obs.dryMass.value, Obs.scMass.value, self.currentTimeNorm.to('day').round(2)))
        if (self.currentTimeNorm + Obs.settlingTime + mode['syst']['ohTime'] >= self.OBendTimes[-1]):
            self.vprint('Last Observing Block (OBnum %d, OBendTime[-1] %.2f) would be exceeded at currentTimeNorm %sd'%(self.OBnumber, self.OBendTimes[-1].value, self.currentTimeNorm.to('day').round(2)))
        if (self.exoplanetObsTime.to('day') + Obs.settlingTime + mode['syst']['ohTime'] >= self.missionLife.to('day')*self.missionPortion):
            self.vprint('exoplanetObstime (%.2f) would exceed (missionPortion*missionLife= %.2f) at currentTimeNorm %sd'%(self.exoplanetObsTime.value, self.missionPortion*self.missionLife.to('day').value, self.currentTimeNorm.to('day').round(2)))
        if (self.currentTimeNorm + Obs.settlingTime + mode['syst']['ohTime'] >= self.missionLife.to('day')):
            self.vprint('missionLife would be exceeded at %s'%self.currentTimeNorm.to('day').round(2))
        
        return is_over 

def calculate_slewTimes(self,TL,old_sInd,sInds,sd,obsTimes,currentTime):
        """Finds slew times and separation angles between target stars
        
        This method determines the slew times of an occulter spacecraft needed
        to transfer from one star's line of sight to all others in a given 
        target list.
        
        Args:
            TL (TargetList module):
                TargetList class object
            old_sInd (integer):
                Integer index of the most recently observed star
            sInds (integer ndarray):
                Integer indeces of the star of interest
            sd (astropy Quantity):
                Angular separation between stars in rad
            currentTime (astropy Time):
                Current absolute mission time in MJD
                
        Returns:
            astropy Quantity:
                Time to transfer to new star line of sight in units of days
        """
    
        self.ao = self.thrust/self.scMass
        slewTime_fac = (2.*self.occulterSep/np.abs(self.ao)/(self.defburnPortion/2. - 
            self.defburnPortion**2./4.)).decompose().to('d2')

        if old_sInd is None:
            slewTimes = np.zeros(TL.nStars)*u.d
        else:
            # calculate slew time
            slewTimes = np.sqrt(slewTime_fac*np.sin(abs(sd)/2.)) #an issue exists if sd is negative
            
            #The following are debugging 
            assert np.where(np.isnan(slewTimes))[0].shape[0] == 0, 'At least one slewTime is nan'
        
        return slewTimes 

def fEZ(self, MV, I, d):
        """Returns surface brightness of exo-zodiacal light
        
        Args:
            MV (integer ndarray):
                Apparent magnitude of the star (in the V band)
            I (astropy Quantity array):
                Inclination of the planets of interest in units of deg
            d (astropy Quantity nx3 array):
                Distance to star of the planets of interest in units of AU
        
        Returns:
            astropy Quantity array:
                Surface brightness of exo-zodiacal light in units of 1/arcsec2
        
        """
        
        # apparent magnitude of the star (in the V band)
        MV = np.array(MV, ndmin=1, copy=False)
        # apparent magnitude of the Sun (in the V band)
        MVsun = 4.83
        
        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]
        beta[mask] = 180.0 - beta[mask]
        beta = 90.0 - beta

        fbeta = 2.44 - 0.0403*beta + 0.000269*beta**2
        
        fEZ = nEZ*10**(-0.4*self.magEZ)*10.**(-0.4*(MV - 
                MVsun))*2*fbeta/d.to('AU').value**2/u.arcsec**2
        
        return fEZ 

def eclip2rot(self,TL,sInd,currentTime):
        """Rotates star position vectors from ecliptic to rotating frame in CRTBP
        
        This method returns a star's position vector in the rotating frame of 
        the Circular Restricted Three Body Problem.  
        
        Args:
            TL (TargetList module):
                TargetList class object
            sInd (integer):
                Integer index of the star of interest
            currentTime (astropy Time):
                Current absolute mission time in MJD

        Returns:
            star_rot (astropy Quantity 1x3 array):
                Star position vector in rotating frame in units of AU
        """
        
        star_pos = TL.starprop(sInd,currentTime).to('au')
        theta    = (np.mod(currentTime.value,self.equinox.value[0])*u.d).to('yr') / u.yr * (2.*np.pi) * u.rad
        
        if currentTime.size == 1:
            star_rot = np.array([np.dot(self.rot(theta, 3), 
                star_pos[x,:].to('AU').value) for x in range(len(star_pos))])[0]*u.AU
        else:
            star_rot = np.array([np.dot(self.rot(theta[x], 3), 
                star_pos[x,:].to('AU').value) for x in range(len(star_pos))])*u.AU

        return star_rot 

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.
        
        Note:
            Prototype includes keepout angles because the values are all the same

        Args:
            sInds[sInds] (integer array):
                the star indicies we would like fZmax and fZmaxInds returned for
            Obs (module):
                Observatory module
            TL (module):
                Target List Module
            TK (TimeKeeping object):
                TimeKeeping object
            mode (dict):
                Selected observing mode
            hashname (string):
                hashname describing the files specific to the current json script
                
        Returns:
            tuple:
            valfZmax[sInds] (astropy Quantity array):
                the maximum fZ (for the prototype, these all have the same value) with units 1/arcsec**2
            absTimefZmax[sInds] (astropy Time array):
                returns the absolute Time the maximum fZ occurs (for the prototype, these all have the same value)
        """
        # cast sInds to array
        sInds = np.array(sInds, ndmin=1, copy=False)
        # get all array sizes
        nStars = sInds.size

        nZ = np.ones(nStars)
        valfZmax = nZ*10**(-0.4*self.magZ)/u.arcsec**2

        absTimefZmax = nZ*u.d + TK.currentTimeAbs

        return valfZmax[sInds], absTimefZmax[sInds]