Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
"orbit": [
26553.4 * u.km,
0.741 * u.one,
63.4 * u.deg,
0.0 * u.deg,
-10.12921 * u.deg,
0.0 * u.rad,
],
"period": 28 * u.day,
}
moon_leo = {
"body": Moon,
"tof": 60 * u.day,
"raan": -2.18 * 1e-4 * u.deg,
"argp": 15.0 * 1e-3 * u.deg,
"inc": 6.0 * 1e-4 * u.deg,
"orbit": [
6678.126 * u.km,
0.01 * u.one,
28.5 * u.deg,
0.0 * u.deg,
0.0 * u.deg,
0.0 * u.rad,
],
"period": 28 * u.day,
}
moon_geo = {
"body": Moon,
"tof": 60 * u.day,
"raan": 6.0 * u.deg,
if x_range > 300.:
xw_min = coord_wraps[coord_index] - 360
xw_max = coord_wraps[coord_index] - np.spacing(360.)
elif xw_min < 0.:
xw_min = max(-180., xw_min - 0.1 * x_range)
xw_max = min(+180., xw_max + 0.1 * x_range)
else:
xw_min = max(0., xw_min - 0.1 * x_range)
xw_max = min(360., xw_max + 0.1 * x_range)
elif coord_type == 'latitude':
xw_min = max(-90., xw_min - 0.1 * x_range)
xw_max = min(+90., xw_max + 0.1 * x_range)
if coord_type in LONLAT:
xw_min *= u.deg.to(unit)
xw_max *= u.deg.to(unit)
ranges.append((xw_min, xw_max))
return ranges
else:
ApertureClass = photutils.SkyRectangularAperture
else:
raise ValueError('invalid aperture_type')
# If on-sky, converts aperture parameters to quantities
if coord_type == 'sky':
if aperture_type == 'circular':
n_params = 1
else:
n_params = 3
assert len(aperture_params) == n_params, 'invalid number of parameters'
units = [astropy.units.arcsec, astropy.units.arcsec, astropy.units.deg]
for ii in range(n_params):
if not isinstance(aperture_params[ii], astropy.units.Quantity):
aperture_params[ii] *= units[ii]
aperture = ApertureClass(coords, *aperture_params)
# Overrides the aperture class so that it inherits from MarvinAperture and
# can gain the methods we defined there. Sets the parent to self.
aperture.__class__ = type('MarvinAperture', (ApertureClass, MarvinAperture), {})
aperture.parent = self
return aperture
else:
data = np.loadtxt(self.multiFile)
else:
data = [self.pos]
columns = ['Gaia_ID', 'TIC_ID', 'RA', 'Dec', 'Gaia_sep', 'TIC_sep', 'Gmag', 'Tmag', 'pmra', 'pmdec', 'parallax']
t = Table(np.zeros(11), names=columns)
t['RA'].unit, t['Dec'].unit, t['Gaia_sep'].unit, t['TIC_sep'].unit = u.arcsec, u.arcsec, u.arcsec, u.arcsec
t['pmra'].unit, t['pmdec'].unit = u.mas/u.year, u.mas/u.year
for i in range(len(data)):
pos = data[i]
gaia = self.crossmatch_by_position(0.005, 'Mast.GaiaDR2.Crossmatch', pos)[0]
tess = self.crossmatch_by_position(0.5, 'Mast.Tic.Crossmatch', pos)[0]
pos[0], pos[1] = pos[0]*u.deg, pos[1]*u.deg
gaiaPos = [gaia['MatchRA'], gaia['MatchDEC']]
sepGaia = self.crossmatch_distance(pos, gaiaPos)
tessPos = [tess['MatchRA'], tess['MatchDEC']]
sepTess = self.crossmatch_distance(pos, tessPos)
t.add_row([gaia['MatchID'], tess['MatchID'], pos[0], pos[1], sepGaia, sepTess, gaia['phot_g_mean_mag'],
tess['Tmag'], gaia['pmra'], gaia['pmdec'], gaia['parallax']])
t.remove_row(0)
return t
def _check_beam_areas(self, threshold, mean_beam, mask=None):
"""
Check that the beam areas are the same to within some threshold
"""
if mask is not None:
assert len(mask) == len(self.beams)
mask = np.array(mask, dtype='bool')
else:
mask = np.ones(len(self.beams), dtype='bool')
qtys = dict(sr=u.Quantity(self.beams, u.sr),
major=u.Quantity([bm.major for bm in self.beams], u.deg),
minor=u.Quantity([bm.minor for bm in self.beams], u.deg),
# position angles are not really comparable
#pa=u.Quantity([bm.pa for bm in self.beams], u.deg),
)
errormessage = ""
for (qtyname, qty) in (qtys.items()):
minv = qty[mask].min().value
maxv = qty[mask].max().value
mn = getattr(mean_beam, qtyname).value
maxdiff = np.max(np.abs((maxv-mn, minv-mn)))/mn
if isinstance(threshold, dict):
th = threshold[qtyname]
else:
th = threshold
cos_b = cos_c * cos_a + sin_c * sin_a * cos_B
# sin_b = np.sqrt(1 - cos_b**2)
# sine rule and cosine rule for A (using both lets arctan2 pick quadrant).
# multiplying both sin_A and cos_A by x=sin_b * sin_c prevents /0 errors
# at poles. Correct for the x=0 multiplication a few lines down.
# sin_A/sin_a == sin_B/sin_b # Sine rule
xsin_A = sin_a * sin_B * sin_c
# cos_a == cos_b * cos_c + sin_b * sin_c * cos_A # cosine rule
xcos_A = cos_a - cos_b * cos_c
A = Angle(np.arctan2(xsin_A, xcos_A), u.radian)
# Treat the poles as if they are infinitesimally far from pole but at given lon
small_sin_c = sin_c < 1e-12
if small_sin_c.any():
# For south pole (cos_c = -1), A = posang; for North pole, A=180 deg - posang
A_pole = (90*u.deg + cos_c*(90*u.deg-Angle(posang, u.radian))).to(u.rad)
if A.shape:
# broadcast to ensure the shape is like that of A, which is also
# affected by the (possible) shapes of lat, posang, and distance.
small_sin_c = np.broadcast_to(small_sin_c, A.shape)
A[small_sin_c] = A_pole[small_sin_c]
else:
A = A_pole
outlon = (Angle(lon, u.radian) + A).wrap_at(360.0*u.deg).to(u.deg)
outlat = Angle(np.arcsin(cos_b), u.radian).to(u.deg)
return outlon, outlat
release.
"""
if not distlimit.isscalar:
raise ValueError('distlimit must be a scalar in search_around_3d')
if coords1.isscalar or coords2.isscalar:
raise ValueError('One of the inputs to search_around_3d is a scalar. '
'search_around_3d is intended for use with array '
'coordinates, not scalars. Instead, use '
'``coord1.separation_3d(coord2) < distlimit`` to find '
'the coordinates near a scalar coordinate.')
if len(coords1) == 0 or len(coords2) == 0:
# Empty array input: return empty match
return (np.array([], dtype=int), np.array([], dtype=int),
Angle([], u.deg),
u.Quantity([], coords1.distance.unit))
kdt2 = _get_cartesian_kdtree(coords2, storekdtree)
cunit = coords2.cartesian.x.unit
# we convert coord1 to match coord2's frame. We do it this way
# so that if the conversion does happen, the KD tree of coord2 at least gets
# saved. (by convention, coord2 is the "catalog" if that makes sense)
coords1 = coords1.transform_to(coords2)
kdt1 = _get_cartesian_kdtree(coords1, storekdtree, forceunit=cunit)
# this is the *cartesian* 3D distance that corresponds to the given angle
d = distlimit.to_value(cunit)
idxs1 = []
def gal2cel_angle(glon,glat,angle,offset=1e-7):
from astropy.coordinates import SkyCoord
import astropy.units as u
origin = SkyCoord(glon,glat,unit=u.deg,frame='galactic')
return estimate_angle(angle,origin,'fk5',offset)
xp_out, yp_out = np.meshgrid(x, y, indexing='xy', sparse=False, copy=False)
world_out = wcs_out.pixel_to_world(xp_out, yp_out)
# Convert the input world coordinates to the frame of the output world
# coordinates.
world_in = world_in.transform_to(world_out.frame)
# Finally, compute the pixel positions in the *output* image of the pixels
# from the *input* image.
xp_inout, yp_inout = wcs_out.world_to_pixel(world_in)
world_in_unitsph = world_in.represent_as('unitspherical')
xw_in, yw_in = world_in_unitsph.lon.to_value(u.deg), world_in_unitsph.lat.to_value(u.deg)
world_out_unitsph = world_out.represent_as('unitspherical')
xw_out, yw_out = world_out_unitsph.lon.to_value(u.deg), world_out_unitsph.lat.to_value(u.deg)
# Put together the parameters common both to the serial and parallel implementations. The aca
# function is needed to enforce that the array will be contiguous when passed to the low-level
# raw C function, otherwise Cython might complain.
aca = np.ascontiguousarray
common_func_par = [0, ny_in, nx_out, ny_out, aca(xp_inout), aca(yp_inout),
aca(xw_in), aca(yw_in), aca(xw_out), aca(yw_out), aca(array),
shape_out]
if nproc == 1:
array_new, weights = _reproject_slice([0, nx_in] + common_func_par)
west_depression = Angle(np.sqrt(2. * ( abs(self.height_above_west)
/ (thorconsts.EQUAT_RAD))) * u.rad)
# print "west_depression",west_depression
if self.height_above_west > 0. :
self.setalt = Angle(-0.833,unit=u.deg) - west_depression
# zd = 90 deg 50 arcmin
elif self.height_above_west <= 0. :
self.setalt = Angle(-0.833,unit=u.deg) + west_depression
east_depression = Angle(np.sqrt(2. * (abs(self.height_above_east)
/ (thorconsts.EQUAT_RAD))) * u.rad)
# print "east_depression",east_depression
if self.height_above_east > 0. :
self.risealt = Angle(-0.833 * u.deg) - east_depression
# zd = 90 deg 50 arcmin
elif self.height_above_east <= 0. :
self.risealt = Angle(-0.833 * u.deg) + east_depression
# print "self.risealt = ",self.risealt
# compute and store the speed of diurnal rotation at this site
# for later use in barycentric velocity correction.
axisdist = np.sqrt(self.location.x ** 2 + self.location.y ** 2)
axisdist = axisdist.to(u.km)
# one sidereal day is about 86164.09 seconds.
self.diurnalspeed = (2. * np.pi * axisdist) / (86164.09 * u.s)