Python astropy.units.AU Examples

def test_fEZ(self):
        """
        Test that fEZ returns correct shape and units.
        """
        exclude_mods=[]

        for mod in self.allmods:
            if mod.__name__ in exclude_mods:
                continue
            if 'fEZ' in mod.__dict__:
                obj = mod()
                # use 3 planets
                d = 10.*np.random.rand(3)*u.AU
                I = np.random.uniform(0.0, 180.0, 3)*u.deg
                fEZs = obj.fEZ(self.TL.MV[0], I, d)
                self.assertEqual(len(fEZs), 3, 'fEZ does not return same number of values as planets tested for {}'.format(mod.__name__))
                self.assertEqual(fEZs.unit, self.unit, 'fEZ does not return 1/arcsec**2 for {}'.format(mod.__name__)) 

def get_gcrs_pos(self, location):
        """ Transforms ICRS frame to GCRS frame
        Parameters
        ----------
        loc : numpy.ndarray shape (3,1) units: AU
              Cartesian Coordinates of the location
              in ICRS frame
        Returns
        ----------
        loc : numpy.ndarray shape (3,1) units: meters
              GCRS coordinates in cartesian system
        """
        loc = location
        loc = coordinates.SkyCoord(x=loc[0], y=loc[1], z=loc[2], unit=units.AU,
                frame='icrs', representation_type='cartesian').transform_to('gcrs')
        loc.representation_type = 'cartesian'
        conv = np.float32(((loc.x.unit/units.m).decompose()).to_string())
        loc = np.array([np.float32(loc.x), np.float32(loc.y),
                        np.float32(loc.z)])*conv
        return loc 

def hist(self, nplan, xedges, yedges):
        """Returns completeness histogram for Monte Carlo simulation
        
        This function uses the inherited Planet Population module.
        
        Args:
            nplan (float):
                number of planets used
            xedges (float ndarray):
                x edge of 2d histogram (separation)
            yedges (float ndarray):
                y edge of 2d histogram (dMag)
        
        Returns:
            h (ndarray):
                2D numpy ndarray containing completeness histogram
        
        """
        
        s, dMag = self.genplans(nplan)
        # get histogram
        h, yedges, xedges = np.histogram2d(dMag, s.to('AU').value, bins=1000,
                range=[[yedges.min(), yedges.max()], [xedges.min(), xedges.max()]])
        
        return h, xedges, yedges 

def test_calc_Teff(self):
        """
        Tests that calc_Teff returns values within appropriate ranges.
        """

        for mod in self.allmods:
            if 'calc_Teff' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    obj = mod()
                starL = np.random.uniform(0.7,1.3,100)
                p = np.random.uniform(0.0, 1.0, 100)
                d = np.random.uniform(0.2, 10.0, 100)*u.AU
                Teff = obj.calc_Teff(starL,d,p)

                self.assertTrue(len(Teff) == len(starL),'length of Teff returned does not match inputs for %s'%mod.__name__)
                self.assertTrue(np.all(np.isfinite(Teff)),'calc_Teff returned infinite value for %s'%mod.__name__)
                self.assertTrue(np.all(Teff >= 0.0),'calc_Teff returned negative value for %s'%mod.__name__) 

def test_equivalent_units2():
    units = set(u.Hz.find_equivalent_units(u.spectral()))
    match = set(
        [u.AU, u.Angstrom, u.Hz, u.J, u.Ry, u.cm, u.eV, u.erg, u.lyr,
         u.m, u.micron, u.pc, u.solRad, u.Bq, u.Ci, u.k, u.earthRad,
         u.jupiterRad])
    assert units == match

    from astropy.units import imperial
    with u.add_enabled_units(imperial):
        units = set(u.Hz.find_equivalent_units(u.spectral()))
        match = set(
            [u.AU, u.Angstrom, imperial.BTU, u.Hz, u.J, u.Ry,
             imperial.cal, u.cm, u.eV, u.erg, imperial.ft, imperial.fur,
             imperial.inch, imperial.kcal, u.lyr, u.m, imperial.mi,
             imperial.mil, u.micron, u.pc, u.solRad, imperial.yd, u.Bq, u.Ci,
             imperial.nmi, u.k, u.earthRad, u.jupiterRad])
        assert units == match

    units = set(u.Hz.find_equivalent_units(u.spectral()))
    match = set(
        [u.AU, u.Angstrom, u.Hz, u.J, u.Ry, u.cm, u.eV, u.erg, u.lyr,
         u.m, u.micron, u.pc, u.solRad, u.Bq, u.Ci, u.k, u.earthRad,
         u.jupiterRad])
    assert units == match 

def test_regression_4210():
    """
    Issue: https://github.com/astropy/astropy/issues/4210
    Related PR with actual change: https://github.com/astropy/astropy/pull/4211
    """
    crd = SkyCoord(0*u.deg, 0*u.deg, distance=1*u.AU)
    ecl = crd.geocentricmeanecliptic
    # bug was that "lambda", which at the time was the name of the geocentric
    # ecliptic longitude, is a reserved keyword. So this just makes sure the
    # new name is are all valid
    ecl.lon

    # and for good measure, check the other ecliptic systems are all the same
    # names for their attributes
    from astropy.coordinates.builtin_frames import ecliptic
    for frame_name in ecliptic.__all__:
        eclcls = getattr(ecliptic, frame_name)
        eclobj = eclcls(1*u.deg, 2*u.deg, 3*u.AU)

        eclobj.lat
        eclobj.lon
        eclobj.distance 

def gen_sma(self, n):
        """Generate semi-major axis values in AU
        
        Samples a power law distribution with exponential turn-off 
        determined by class attribute smaknee
        
        Args:
            n (integer):
                Number of samples to generate
                
        Returns:
            astropy Quantity array:
                Semi-major axis values in units of AU
        
        """
        n = self.gen_input_check(n)
        ar = self.arange.to('AU').value
        a = statsFun.simpSample(self.dist_sma, n, ar[0], ar[1])*u.AU
        
        return a 

def test_component_names_repr():
    from astropy.coordinates.baseframe import BaseCoordinateFrame, RepresentationMapping

    # Frame class with new component names that includes a name swap
    class NameChangeFrame(BaseCoordinateFrame):
        default_representation = r.PhysicsSphericalRepresentation

        frame_specific_representation_info = {
            r.PhysicsSphericalRepresentation: [
                RepresentationMapping('phi', 'theta', u.deg),
                RepresentationMapping('theta', 'phi', u.arcsec),
                RepresentationMapping('r', 'JUSTONCE', u.AU)]
        }
    frame = NameChangeFrame(0*u.deg, 0*u.arcsec, 0*u.AU)

    # Check for the new names in the Frame repr
    assert "(theta, phi, JUSTONCE)" in repr(frame)

    # Check that the letter "r" has not been replaced more than once in the Frame repr
    assert repr(frame).count("JUSTONCE") == 1 

def test_period_semi():
    # Check that an error is raised if neither a nor porb is given
    with pytest.raises(ValueError) as e:
        body = starry.Secondary(starry.Map())
    assert "Must provide a value for either `porb` or `a`" in str(e.value)

    # Check that the semi --> period conversion works
    pri = starry.Primary(starry.Map(), m=1.0, mass_unit=u.Msun)
    sec = starry.Secondary(
        starry.Map(), a=10.0, m=1.0, length_unit=u.AU, mass_unit=u.Mearth
    )
    sys = starry.System(pri, sec)
    period = sys._get_periods()[0]
    true_period = (
        (
            (2 * np.pi)
            * (sec.a * sec.length_unit) ** (3 / 2)
            / (np.sqrt(G * (pri.m * pri.mass_unit + sec.m * sec.mass_unit)))
        )
        .to(u.day)
        .value
    )
    assert np.allclose(period, true_period) 

def __init__(self, a, e, I, O, w, lM):
            
            # store list of semimajor axis values (convert from AU to km)
            self.a = (a*u.AU).to('km').value
            if not isinstance(self.a, float):
                self.a = self.a.tolist()
            # store list of dimensionless eccentricity values
            self.e = e
            # store list of inclination values (degrees)
            self.I = I
            # store list of right ascension of ascending node values (degrees)
            self.O = O 
            # store list of longitude of periapsis values (degrees)
            self.w = w 
            # store list of mean longitude values (degrees)
            self.lM = lM 

def eclip2equat(self, r_eclip, currentTime):
        """Rotates heliocentric coordinates from ecliptic to equatorial frame.
        
        Args:
            r_eclip (astropy Quantity nx3 array):
                Positions vector in heliocentric ecliptic frame in units of AU
            currentTime (astropy Time array):
                Current absolute mission time in MJD
        
        Returns:
            r_equat (astropy Quantity nx3 array):
                Positions vector in heliocentric equatorial frame in units of AU
        
        """
        
        r_equat = self.equat2eclip(r_eclip, currentTime, rotsign=-1)
        
        return r_equat 

def dist_sma(self, a):
        """Probability density function for semi-major axis in AU
        
        The prototype provides a log-uniform distribution between the minimum
        and maximum values.
        
        Args:
            a (float ndarray):
                Semi-major axis value(s) in AU. Not an astropy quantity.
                
        Returns:
            float ndarray:
                Semi-major axis probability density
        
        """
        
        return self.logunif(a, self.arange.to('AU').value) 

def gen_sma(self, n):
        """Generate semi-major axis values in AU
        
        The Earth-twin population assumes a uniform distribution between the minimum and 
        maximum values.
        
        Args:
            n (integer):
                Number of samples to generate
                
        Returns:
            a (astropy Quantity array):
                Semi-major axis values in units of AU
        
        """
        n = self.gen_input_check(n)
        ar = self.arange.to('AU').value
        a = np.random.uniform(low=ar[0], high=ar[1], size=n)*u.AU
        
        return a 

def gen_sma(self, n):
        """Generate semi-major axis values in AU
        
        Samples a power law distribution with exponential turn-off 
        determined by class attribute smaknee
        
        Args:
            n (integer):
                Number of samples to generate
                
        Returns:
            a (astropy Quantity array):
                Semi-major axis in units of AU
        
        """
        n = self.gen_input_check(n)
        a = self.sma_sampler(n)*u.AU
        
        return a 

def calc_fdmag(self, dMag, smin, smax):
        """Calculates probability density of dMag by integrating over projected
        separation
        
        Args:
            dMag (float ndarray):
                Planet delta magnitude(s)
            smin (float ndarray):
                Value of minimum projected separation (AU) from instrument
            smax (float ndarray):
                Value of maximum projected separation (AU) from instrument
        
        Returns:
            float:
                Value of probability density
        
        """
        
        f = np.zeros(len(smin))
        for k, dm in enumerate(dMag):
            f[k] = interpolate.InterpolatedUnivariateSpline(self.xnew,self.EVPOCpdf(self.xnew,dm),ext=1).integral(smin[k],smax[k])
            
        return f 

def hist(self, nplan, xedges, yedges):
        """Returns completeness histogram for Monte Carlo simulation
        
        This function uses the inherited Planet Population module.
        
        Args:
            nplan (float):
                number of planets used
            xedges (float ndarray):
                x edge of 2d histogram (separation)
            yedges (float ndarray):
                y edge of 2d histogram (dMag)
        
        Returns:
            float ndarray:
                2D numpy ndarray containing completeness frequencies
        
        """
        
        s, dMag = self.genplans(nplan)
        # get histogram
        h, yedges, xedges = np.histogram2d(dMag, s.to('AU').value, bins=1000,
                range=[[yedges.min(), yedges.max()], [xedges.min(), xedges.max()]])
        
        return h, xedges, yedges 

def test_earth_barycentric_velocity_multi_d():
    # Might as well test it with a multidimensional array too.
    t = Time('2016-03-20T12:30:00') + np.arange(8.).reshape(2, 2, 2) * u.yr / 2.
    ep, ev = get_body_barycentric_posvel('earth', t)
    # note: assert_quantity_allclose doesn't like the shape mismatch.
    # this is a problem with np.testing.assert_allclose.
    assert quantity_allclose(ep.get_xyz(xyz_axis=-1),
                             [[-1., 0., 0.], [+1., 0., 0.]]*u.AU,
                             atol=0.06*u.AU)
    expected = u.Quantity([0.*u.one,
                           np.cos(23.5*u.deg),
                           np.sin(23.5*u.deg)]) * ([[-30.], [30.]] * u.km / u.s)
    assert quantity_allclose(ev.get_xyz(xyz_axis=-1), expected,
                             atol=2.*u.km/u.s) 

def test_load_systems(self):
        """
        Test that load systems loads correctly
        """

        with open(self.script) as f:
            spec = json.loads(f.read())
        spec['modules']['StarCatalog'] = 'EXOCAT1'
        with RedirectStreams(stdout=self.dev_null):
            SU = SimulatedUniverse(**spec)

        # dictionary of planetary parameter keys
        param_keys = ['a', 'e', 'I', 'O', 'w', 'M0', 'Mp', 'Rp', 'p', 'plan2star']
        systems = {'a': np.array([5.])*u.AU,
                   'e': np.array([0.]),
                   'I': np.array([0.])*u.deg,
                   'O': np.array([0.])*u.deg,
                   'w': np.array([0.])*u.deg,
                   'M0': np.array([0.])*u.deg,
                   'Mp': np.array([300.])*u.earthMass,
                   'Rp': np.array([10.])*u.earthRad,
                   'p': np.array([0.6]),
                   'plan2star': np.array([0], dtype=int)}
        SU.load_systems(systems)
        for key in param_keys:
            self.assertTrue(systems[key] == getattr(SU, key), 'Value for %s not assigned with load_systems'%key) 

def __init__(self, coords):
        valid_coords = ['KRTP', 'KSM', 'KSO', 'RTN']
        if coords not in valid_coords:
            raise ValueError('coords must be one of {}'.format(valid_coords))
        self.coords = coords

        Rs = u.def_unit('saturnRad', 60268 * u.km)
        if (coords == 'KRTP'):
            self.units = OrderedDict([('Bx', u.nT), ('By', u.nT), ('Bz', u.nT),
                                      ('X', Rs), ('|B|', u.nT),
                                      ('Y', u.deg),
                                      ('Z', u.deg),
                                      ('Local hour', u.dimensionless_unscaled),
                                      ('n points', u.dimensionless_unscaled)])
        if (coords == 'RTN'):
            self.units = OrderedDict([('Bx', u.nT), ('By', u.nT), ('Bz', u.nT),
                                      ('X', u.AU), ('Y', u.AU), ('Z', u.AU),
                                      ('|B|', u.nT),
                                      ('Local hour', u.dimensionless_unscaled),
                                      ('n points', u.dimensionless_unscaled)])
        if (coords == 'KSM' or coords == 'KSO'):
            self.units = OrderedDict([('Bx', u.nT), ('By', u.nT), ('Bz', u.nT),
                                      ('X', Rs), ('Y', Rs), ('Z', Rs),
                                      ('|B|', u.nT),
                                      ('Local hour', u.dimensionless_unscaled),
                                      ('n points', u.dimensionless_unscaled)]) 

def solarSystem_body_position(self, currentTime, bodyname, eclip=False):
        """Finds solar system body positions vector in heliocentric equatorial (default)
        or ecliptic frame for current time (MJD).
        
        This passes all arguments to one of spk_body or keplerplanet, depending
        on the value of self.havejplephem.
        
        Args:
            currentTime (astropy Time array):
                Current absolute mission time in MJD
            bodyname (string):
                Solar system object name
            eclip (boolean):
                Boolean used to switch to heliocentric ecliptic frame. Defaults to 
                False, corresponding to heliocentric equatorial frame.
        
        Returns:
            astropy Quantity nx3 array:
                Solar system body positions in heliocentric equatorial (default)
                or ecliptic frame in units of AU
        
        Note: 
            Use eclip=True to get ecliptic coordinates.
        
        """
        
        # heliocentric
        if bodyname == 'Sun':
            return np.zeros((currentTime.size, 3))*u.AU
        
        # choose JPL or static ephemerides
        if self.havejplephem:
            r_body = self.spk_body(currentTime, bodyname, eclip=eclip).to('AU')
        else:
            r_body = self.keplerplanet(currentTime, bodyname, eclip=eclip).to('AU')
        
        return r_body 

def __init__(self, probe):
        self.probe = _check_probe(probe)
        self.units = OrderedDict([
            ('B instrument', u.dimensionless_unscaled),
            ('Bx', u.nT), ('By', u.nT), ('Bz', u.nT),
            ('sigma B', u.nT),
            ('Ion instrument', u.dimensionless_unscaled),
            ('Status', u.dimensionless_unscaled),
            ('Tp_par', u.K), ('Tp_perp', u.K),
            ('carrot', u.dimensionless_unscaled),
            ('r_sun', u.AU), ('clat', u.deg),
            ('clong', u.deg), ('earth_he_angle', u.deg),
            ('n_p', u.cm**-3), ('vp_x', u.km / u.s),
            ('vp_y', u.km / u.s), ('vp_z', u.km / u.s),
            ('vth_p_par', u.km / u.s), ('vth_p_perp', u.km / u.s)]) 

def is_earthlike(self, specs, plan_id, star_ind):
        """Depricated Determine if this planet is Earth-Like or Not, given specs/star id/planet id
        """
        # extract planet and star properties
        Rp_plan = strip_units(specs['Rp'][plan_id])
        a_plan = strip_units(specs['a'][plan_id])
        L_star = specs['L'][star_ind]
        L_plan = L_star / (a_plan**2) # adjust star luminosity by distance^2 in AU
        # its radius (in earth radii) and solar-equivalent luminosity must be
        # between given bounds.  The lower Rp bound is not axis-parallel, but
        # the best axis-parallel bound is 0.90, so that's what we use.
        return (Rp_plan >= 0.90 and Rp_plan <= 1.4) and (L_plan >= 0.3586 and L_plan <= 1.1080) 

def calc_fdmag(self, dMag, smin, smax, subpop=-2):
        """Calculates probability density of dMag by integrating over projected
        separation
        
        Args:
            dMag (float ndarray):
                Planet delta magnitude(s)
            smin (float ndarray):
                Value of minimum projected separation (AU) from instrument
            smax (float ndarray):
                Value of maximum projected separation (AU) from instrument
            subpop (int):
                planet subtype to use for calculation of comp0
                -2 - planet population
                -1 - earthLike population
                (i,j) - kopparapu planet subtypes
        
        Returns:
            float:
                Value of probability density
        
        """
        if subpop == -2:
            f = np.zeros(len(smin))
            for k, dm in enumerate(dMag):
                f[k] = interpolate.InterpolatedUnivariateSpline(self.xnew,self.EVPOCpdf_pop(self.xnew,dm),ext=1).integral(smin[k],smax[k])
        elif subpop == -1:
            f = np.zeros(len(smin))
            for k, dm in enumerate(dMag):
                f[k] = interpolate.InterpolatedUnivariateSpline(self.xnew,self.EVPOCpdf_earthLike(self.xnew,dm),ext=1).integral(smin[k],smax[k])
        else:
            f = np.zeros(len(smin))
            for k, dm in enumerate(dMag):
                f[k] = interpolate.InterpolatedUnivariateSpline(self.xnew,self.EVPOCpdf_hs[subpop[0],subpop[1]](self.xnew,dm),ext=1).integral(smin[k],smax[k])
            
        return f


    #MODIFY THIS TO CREATE A CLASSIFICATION FOR EACH pInd WHICH IS THE CLASSIFICATION NUMBER 

def dcomp_dt(self, intTimes, TL, sInds, fZ, fEZ, WA, mode, C_b=None, C_sp=None):
        """Calculates derivative of completeness with respect to integration time
        
        Args:
            intTimes (astropy Quantity array):
                Integration times
            TL (TargetList module):
                TargetList class object
            sInds (integer ndarray):
                Integer indices of the stars of interest
            fZ (astropy Quantity array):
                Surface brightness of local zodiacal light in units of 1/arcsec2
            fEZ (astropy Quantity array):
                Surface brightness of exo-zodiacal light in units of 1/arcsec2
            WA (astropy Quantity):
                Working angle of the planet of interest in units of arcsec
            mode (dict):
                Selected observing mode
            C_b (astropy Quantity array):
                Background noise electron count rate in units of 1/s (optional)
            C_sp (astropy Quantity array):
                Residual speckle spatial structure (systematic error) in units of 1/s
                (optional)                
                
        Returns:
            astropy Quantity array:
                Derivative of completeness with respect to integration time (units 1/time)
        
        """
        intTimes, sInds, fZ, fEZ, WA, smin, smax, dMag = self.comps_input_reshape(intTimes, TL, sInds, fZ, fEZ, WA, mode, C_b=C_b, C_sp=C_sp)
        
        ddMag = TL.OpticalSystem.ddMag_dt(intTimes, TL, sInds, fZ, fEZ, WA, mode).reshape((len(intTimes),))
        dcomp = self.calc_fdmag(dMag, smin, smax)
        mask = smin>self.PlanetPopulation.rrange[1].to('AU').value
        dcomp[mask] = 0.
        
        return dcomp*ddMag 

def comp_calc(self, smin, smax, dMag, subpop=-2):
        """Calculates completeness for given minimum and maximum separations
        and dMag
        
        Note: this method assumes scaling orbits when scaleOrbits == True has
        already occurred for smin, smax, dMag inputs
        
        Args:
            smin (float ndarray):
                Minimum separation(s) in AU
            smax (float ndarray):
                Maximum separation(s) in AU
            dMag (float ndarray):
                Difference in brightness magnitude
            subpop (int):
                planet subtype to use for calculation of comp0
                -2 - planet population
                -1 - earthLike population
                (i,j) - kopparapu planet subtypes
        
        Returns:
            float ndarray:
                Completeness values
        
        """
        if subpop == -2:
            comp = self.EVPOC_pop(smin, smax, 0., dMag)
        elif subpop == -1:
            comp = self.EVPOC_earthlike(smin, smax, 0., dMag)
        else:
            #for ii,j in itertools.product(np.arange(len(self.Rp_hi)),np.arange(len(self.L_lo))):
            comp = self.EVPOC_hs[subpop[0],subpop[1]](smin, smax, 0., dMag)
        # remove small values
        comp[comp<1e-6] = 0.
        
        return comp 

def comp_per_intTime(self, intTimes, TL, sInds, fZ, fEZ, WA, mode, C_b=None, C_sp=None):
        """Calculates completeness for integration time
        
        Args:
            intTimes (astropy Quantity array):
                Integration times
            TL (TargetList module):
                TargetList class object
            sInds (integer ndarray):
                Integer indices of the stars of interest
            fZ (astropy Quantity array):
                Surface brightness of local zodiacal light in units of 1/arcsec2
            fEZ (astropy Quantity array):
                Surface brightness of exo-zodiacal light in units of 1/arcsec2
            WA (astropy Quantity):
                Working angle of the planet of interest in units of arcsec
            mode (dict):
                Selected observing mode
            C_b (astropy Quantity array):
                Background noise electron count rate in units of 1/s (optional)
            C_sp (astropy Quantity array):
                Residual speckle spatial structure (systematic error) in units of 1/s
                (optional)
                
        Returns:
            flat ndarray:
                Completeness values
        
        """
        intTimes, sInds, fZ, fEZ, WA, smin, smax, dMag = self.comps_input_reshape(intTimes, TL, sInds, fZ, fEZ, WA, mode, C_b=C_b, C_sp=C_sp)
        
        comp = self.comp_calc(smin, smax, dMag)
        mask = smin>self.PlanetPopulation.rrange[1].to('AU').value
        comp[mask] = 0.
        # ensure completeness values are between 0 and 1
        comp = np.clip(comp, 0., 1.)
        
        return comp 

def time_delay_from_location(self, other_location, right_ascension,
                                 declination, t_gps):
        """Return the time delay from the LISA detector to detector for
        a signal with the given sky location. In other words return
        `t1 - t2` where `t1` is the arrival time in this detector and
        `t2` is the arrival time in the other location. Units(AU)
        Parameters
        ----------
        other_location : numpy.ndarray of coordinates in ICRS frame
            A detector instance.
        right_ascension : float
            The right ascension (in rad) of the signal.
        declination : float
            The declination (in rad) of the signal.
        t_gps : float
            The GPS time (in s) of the signal.
        Returns
        -------
        numpy.ndarray
            The arrival time difference between the detectors.
        """
        dx = np.array([self.location[0] - other_location[0],
                       self.location[1] - other_location[1],
                       self.location[2] - other_location[2]])
        cosd = cos(declination)
        e0 = cosd * cos(right_ascension)
        e1 = cosd * -sin(right_ascension)
        e2 = sin(declination)
        ehat = np.array([e0, e1, e2])
        return dx.dot(ehat) / constants.c.value 

def dcomp_dt(self, intTimes, TL, sInds, fZ, fEZ, WA, mode, C_b=None, C_sp=None):
        """Calculates derivative of completeness with respect to integration time
        
        Args:
            intTimes (astropy Quantity array):
                Integration times
            TL (TargetList module):
                TargetList class object
            sInds (integer ndarray):
                Integer indices of the stars of interest
            fZ (astropy Quantity array):
                Surface brightness of local zodiacal light in units of 1/arcsec2
            fEZ (astropy Quantity array):
                Surface brightness of exo-zodiacal light in units of 1/arcsec2
            WA (astropy Quantity):
                Working angle of the planet of interest in units of arcsec
            mode (dict):
                Selected observing mode
            C_b (astropy Quantity array):
                Background noise electron count rate in units of 1/s (optional)
            C_sp (astropy Quantity array):
                Residual speckle spatial structure (systematic error) in units of 1/s
                (optional)                
                
        Returns:
            astropy Quantity array:
                Derivative of completeness with respect to integration time (units 1/time)
        
        """
        intTimes, sInds, fZ, fEZ, WA, smin, smax, dMag = self.comps_input_reshape(intTimes, TL, sInds, fZ, fEZ, WA, mode, C_b=C_b, C_sp=C_sp)
        
        ddMag = TL.OpticalSystem.ddMag_dt(intTimes, TL, sInds, fZ, fEZ, WA, mode).reshape((len(intTimes),))
        dcomp = self.calc_fdmag(dMag, smin, smax)
        mask = smin>self.PlanetPopulation.rrange[1].to('AU').value
        dcomp[mask] = 0.
        
        return dcomp*ddMag 

def comp_per_intTime(self, intTimes, TL, sInds, fZ, fEZ, WA, mode, C_b=None, C_sp=None):
        """Calculates completeness for integration time
        
        Args:
            intTimes (astropy Quantity array):
                Integration times
            TL (TargetList module):
                TargetList class object
            sInds (integer ndarray):
                Integer indices of the stars of interest
            fZ (astropy Quantity array):
                Surface brightness of local zodiacal light in units of 1/arcsec2
            fEZ (astropy Quantity array):
                Surface brightness of exo-zodiacal light in units of 1/arcsec2
            WA (astropy Quantity):
                Working angle of the planet of interest in units of arcsec
            mode (dict):
                Selected observing mode
            C_b (astropy Quantity array):
                Background noise electron count rate in units of 1/s (optional)
            C_sp (astropy Quantity array):
                Residual speckle spatial structure (systematic error) in units of 1/s
                (optional)
                
        Returns:
            flat ndarray:
                Completeness values
        
        """
        intTimes, sInds, fZ, fEZ, WA, smin, smax, dMag = self.comps_input_reshape(intTimes, TL, sInds, fZ, fEZ, WA, mode, C_b=C_b, C_sp=C_sp)
        
        comp = self.comp_calc(smin, smax, dMag)
        mask = smin>self.PlanetPopulation.rrange[1].to('AU').value
        comp[mask] = 0.
        # ensure completeness values are between 0 and 1
        comp = np.clip(comp, 0., 1.)
        
        return comp 

def gen_physical_properties(self, **specs):
        """Generating universe based on Dulz and Plavchan occurrence rate tables.

        """

        PPop = self.PlanetPopulation
        TL = self.TargetList

        # treat eta as the rate parameter of a Poisson distribution
        targetSystems = np.random.poisson(lam=PPop.eta, size=TL.nStars)
        plan2star = []
        for j, n in enumerate(targetSystems):
            plan2star = np.hstack((plan2star, [j] * n))

        a, e, p, Rp = PPop.gen_plan_params(len(plan2star))

        #filter to just Earth candidates
        #modify as needed NB: a has no yet been scaled by \sqrt{L}
        #so effectively everything is at 1 solar luminosity
        inds = (a < 1.67*u.AU) & (a > 0.95*u.AU) & (Rp < 1.4*u.R_earth ) & (Rp >0.8/np.sqrt(a.value)*u.R_earth)  

        self.a = a[inds]
        self.e = e[inds]
        self.p = p[inds]
        self.Rp = Rp[inds]
        plan2star = plan2star[inds]


        self.plan2star = plan2star.astype(int)
        self.sInds = np.unique(self.plan2star)
        self.nPlans = len(self.plan2star)

        if PPop.scaleOrbits:
            self.a *= np.sqrt(TL.L[self.plan2star])


        # sample all of the orbital and physical parameters
        self.I, self.O, self.w = PPop.gen_angles(self.nPlans)
        self.gen_M0()  # initial mean anomaly
        self.Mp = PPop.MfromRp(self.Rp)  # mass