Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def toEphem(self):
if str(self.epoch).lower() == str(Epoch.J2000).lower():
epoch = ephem.J2000
elif str(self.epoch).lower() == str(Epoch.B1950).lower():
epoch = ephem.B1950
else:
epoch = ephem.now()
return ephem.Equatorial(self.ra.R, self.dec.R, epoch=epoch)
def to_fits(ftag,i,src,cnt,history=''):
filename = fname(ftag,cnt)
print 'Saving data to', filename
while len(i.shape) < 4: i.shape = i.shape + (1,)
cen = ephem.Equatorial(src.ra, src.dec, epoch=aa.epoch)
# We precess the coordinates of the center of the image here to
# J2000, just to have a well-defined epoch for them. For image coords to
# be accurately reconstructed, precession needs to be applied per pixel
# and not just per phase-center because ra/dec axes aren't necessarily
# aligned between epochs. When reading these images, to be 100% accurate,
# one should precess the ra/dec coordinates back to the date of the
# observation, infer the coordinates of all the pixels, and then
# precess the coordinates for each pixel independently.
cen = ephem.Equatorial(cen, epoch=ephem.J2000)
a.img.to_fits(filename, i, clobber=True,
object=src.src_name, obs_date=str(aa.date),
ra=cen.ra*a.img.rad2deg, dec=cen.dec*a.img.rad2deg, epoch=2000.,
d_ra=L[-1,-1]*a.img.rad2deg, d_dec=M[1,1]*a.img.rad2deg,
freq=np.average(aa[0].beam.afreqs),history=history)
Compute topocentric libration angles and the position angle of the moon's rotational axis.
:return: -
"""
# Astrometric libration angles and the selenographic co-longitude of the sun are provided
# by PyEphem.
self.astrometric_lib_lat = self.moon.libration_lat
self.astrometric_lib_long = self.moon.libration_long
self.colong = self.moon.colong
# Compute topocentric libration values. For algorithmic details, refer to the user's guide.
j2000 = ephem.Date(datetime.datetime(2000, 1, 1, 12, 0))
t = (ephem.Date(self.location_time.date) - j2000) / 36525.
epsilon = radians(23.439281 - 0.013002575 * t)
ma = ephem.Equatorial(self.ra, self.de, epoch=ephem.Date(self.location_time.date))
me = ephem.Ecliptic(ma)
cap_i = radians(1.54266)
omega = radians(125.044555 - 1934.13626194 * t)
lm = radians(218.31664563 + 481267.88119575 * t - 0.00146638 * t ** 2)
i = acos(cos(cap_i) * cos(epsilon) + sin(cap_i) * sin(epsilon) * cos(omega))
omega_prime = atan2(-sin(cap_i) * sin(omega) / sin(i), (
cos(cap_i) * sin(epsilon) - sin(cap_i) * cos(epsilon) * cos(omega)) / sin(i))
self.pos_rot_north = atan2(-sin(i) * cos(omega_prime - self.ra),
cos(self.de) * cos(i) - sin(self.de) * sin(i) * sin(
omega_prime - self.ra))
self.topocentric_lib_lat = asin(
-sin(cap_i) * cos(me.lat) * sin(me.lon - omega) - cos(cap_i) * sin(me.lat))
self.topocentric_lib_long = (atan2(
cos(cap_i) * cos(me.lat) * sin(me.lon - omega) - sin(cap_i) * sin(me.lat),
cos(me.lat) * cos(me.lon - omega)) - lm + omega) % (2. * pi)
def to_fits(ftag,i,src,cnt,history=''):
filename = fname(ftag,cnt)
print 'Saving data to', filename
while len(i.shape) < 4: i.shape = i.shape + (1,)
cen = ephem.Equatorial(src.ra, src.dec, epoch=aa.epoch)
# We precess the coordinates of the center of the image here to
# J2000, just to have a well-defined epoch for them. For image coords to
# be accurately reconstructed, precession needs to be applied per pixel
# and not just per phase-center because ra/dec axes aren't necessarily
# aligned between epochs. When reading these images, to be 100% accurate,
# one should precess the ra/dec coordinates back to the date of the
# observation, infer the coordinates of all the pixels, and then
# precess the coordinates for each pixel independently.
cen = ephem.Equatorial(cen, epoch=ephem.J2000)
a.img.to_fits(filename, i, clobber=True,
object=src.src_name, obs_date=str(aa.date),
ra=cen.ra*a.img.rad2deg, dec=cen.dec*a.img.rad2deg, epoch=2000.,
d_ra=L[-1,-1]*a.img.rad2deg, d_dec=M[1,1]*a.img.rad2deg,
freq=np.average(aa[0].beam.afreqs),history=history)
if not opts.cal is None:
cat = a.cal.get_catalog(opts.cal, srclist, cutoff, catalogs)
aa = a.cal.get_aa(opts.cal, .1, opts.freq, 1)
aa.set_jultime(opts.juldate)
cat.compute(aa)
else:
cat = a.src.get_catalog(srclist, cutoff, catalogs)
m = a.coord.convert_m('eq', opts.osys, oepoch=opts.oepoch)
ths,phis = h.px2crd(np.arange(h.npix()), ncrd=2)
ras,decs = phis, np.pi/2 - ths
afreq = np.array([opts.freq])
for srcname in cat:
src = cat[srcname]
eq = e.Equatorial(src._ra, src._dec, epoch=e.J2000)
eq = e.Equatorial(eq, epoch=opts.oepoch)
ra,dec = eq.get()
# Account for size of source
a1,a2,th = src.srcshape
dras,ddecs = ras - ra, decs - dec
print '--------------------------------------------------'
src.update_jys(afreq)
strength = src.get_jys()[0] * opts.sscale
print 'Adding', srcname, 'with strength %f Jy' % strength,
print 'and index', src.index
print 'Source shape: a1=%f, a2=%f, th=%f' % (a1, a2, th)
da1 = dras*np.cos(th) - ddecs*np.sin(th)
da2 = dras*np.sin(th) + ddecs*np.cos(th)
delta = (da1/a1)**2 + (da2/a2)**2
px = np.where(delta <= 1)[0]
if len(px) == 0:
print 'Treating as point source.'
def precess(self, epoch=Epoch.NOW):
if str(epoch).lower() == str(Epoch.J2000).lower():
epoch = ephem.J2000
elif str(epoch).lower() == str(Epoch.B1950).lower():
epoch = ephem.B1950
elif str(epoch).lower() == str(Epoch.NOW).lower():
epoch = ephem.now()
j2000 = self.toEphem()
now = ephem.Equatorial(j2000, epoch=epoch)
return Position.fromRaDec(
Coord.fromR(now.ra), Coord.fromR(now.dec), epoch=Epoch.NOW)
cat = a.src.get_catalog(srclist, cutoff, catalogs)
o = ephem.Observer()
if opts.juldate is None:
o.date = ephem.J2000
o.epoch = o.date
try: del(cat['Sun'])
except(KeyError): pass
else:
o.date = a.phs.juldate2ephem(opts.juldate)
o.epoch = o.date
for s in cat.values():
try: a.phs.RadioFixedBody.compute(s, o)
except(TypeError): a.phs.RadioSpecial.compute(s, o)
#cat.compute(o)
# lat/lon coordinates of sources
scrds = [ephem.Equatorial(s.ra,s.dec,epoch=o.epoch) for s in cat.values()]
afreqs = np.array([.150])
cat.update_jys(afreqs)
sflxs = cat.get_jys().squeeze()
snams = cat.keys()
if opts.osys == 'ga':
scrds = [ephem.Galactic(s, epoch=opts.oepoch) for s in scrds]
elif opts.osys == 'ec':
scrds = [ephem.Ecliptic(s, epoch=opts.oepoch) for s in scrds]
slats = np.array([float(s.get()[1]) for s in scrds]) * a.img.rad2deg
slons = np.array([float(s.get()[0]) for s in scrds]) * a.img.rad2deg
if opts.osys == 'eq': slons = 360 - slons
slons = np.where(slons < -180, slons + 360, slons)
slons = np.where(slons >= 180, slons - 360, slons)
# Generate map grid/outline
map.drawmapboundary()
o.add_option('--fwidth', dest='fwidth', type='float', default=10,
help='Width of gaussian, in degrees, used to downweight data away from pointing center for each facet. Default 10.')
opts, args = o.parse_args(sys.argv[1:])
# Open skymap
if os.path.exists(opts.map): skymap = a.map.Map(fromfits=opts.map)
else: skymap = a.map.Map(nside=opts.nside)
skymap.set_interpol(opts.interpolate)
prev_dra = None
for i, filename in enumerate(args):
img, kwds = a.img.from_fits(filename)
img = img.squeeze()
# Read ra/dec of image center, which are stored in J2000
assert(kwds['epoch'] == 2000)
s = ephem.Equatorial(kwds['ra']*a.img.deg2rad, kwds['dec']*a.img.deg2rad,
epoch=ephem.J2000)
# To precess the entire image to J2000, we actually need to precess the
# coords of the center back to the epoch of the observation (obs_date),
# get the pixel coordinates in that epoch, and then precess the coords of
# each pixel. This is because extrapolating pixel coords from the J2000
# center assumes the J2000 ra/dec axes, which may be tilted relative to
# the ra/dec axes of the epoch of the image.
s = ephem.Equatorial(s, epoch=kwds['obs_date'])
ra, dec = s.get()
print '-----------------------------------------------------------'
print 'Reading file %s (%d / %d)' % (filename, i + 1, len(args))
print kwds
print 'Pointing (ra, dec):', ra, dec
print 'Image Power:', np.abs(img).sum()
if prev_dra != kwds['d_ra']:
prev_dra = kwds['d_ra']
"""
An astronomy library for precessing coordinates between epochs and converting
between topocentric (z = up, x = east), ecliptic (heliocentric), equatorial
(celestial), and galactic coordinate systems. Vectors are 3 dimensional
with magnitude 1, representing a point on the unit sphere. Includes generic
3 vector rotation code.
"""
import numpy as np, ephem as e
sys_dict = {
'eq': e.Equatorial,
'ec': e.Ecliptic,
'ga': e.Galactic
}
def convert(crd, isys, osys, iepoch=e.J2000, oepoch=e.J2000):
"""Convert 'crd' from coordinate system isys to osys, including
epoch precession. Valid coordinate systems are 'ec' (Ecliptic), 'eq'
(Equatorial), and 'ga' (Galactic). Epochs may be date strings, or
numerical ephemeris times."""
if len(crd) == 3: crd = eq2radec(crd)
c1 = sys_dict[isys[:2].lower()](crd[0], crd[1], epoch=iepoch)
return sys_dict[osys[:2].lower()](c1, epoch=oepoch).get()
def convert_m(isys, osys, iepoch=e.J2000, oepoch=e.J2000):
"""Return the 3x3 matrix corresponding to a coordinate
system/precssion transformation (see 'convert').