Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# -*- coding: utf-8 -*-
from sgp4.ext import jday
from ephem.angles import Angle
from ephem.planets import earth, moon, mercury, jupiter, EarthLocation
from ephem import J2000
jd_tt = 2450203.5
ggr = EarthLocation('75 W', '45 N', 0.0, temperature=10.0, pressure=1010.0)
ra, dec, dis = ggr(jd_tt).observe(moon).radec(J2000)
print('RA: %s' % Angle(ra).hours())
print('dec: %s' % Angle(dec).degrees())
print('')
print('RA diff: %s' % (Angle(ra).hours() - 11.739849403))
print('dec diff: %s' % (Angle(dec).degrees() - -0.31860323))
print('')
print(dis)
#JD = 2450203.500000 RA = 11.739849403 Dec = -0.31860323 Dis = 0.0026040596
exit(0)
jd = jday(2006, 6, 19, 0, 0, 0)
pos = earth(jd).observe(mercury).radec(J2000) # Astrometric
print(pos.ra)
def HA_to_Dec(self, J2000_coordinate, site):
"""
HA to Dec
"""
assert type(J2000_coordinate) == coords.FK5
assert J2000_coordinate.equinox.value == 'J2000.000'
# Coordinate precessed to Jnow (as an astropy coordinate object)
Jnow_coordinate = J2000_coordinate.precess_to(Time.now())
# Coordinate as a pyephem coordinate (J2000)
RA_string, Dec_string = J2000_coordinate.to_string(
precision=2, sep=':').split(' ')
ephem_coordinate = ephem.FixedBody(
RA_string, Dec_string, epoch=ephem.J2000)
ephem_coordinate = ephem.readdb(
'Polaris,f|M|F7,{},{},2.02,2000'.format(RA_string, Dec_string))
ephem_coordinate.compute(site)
HourAngle = ephem_coordinate.ra - site.sidereal_time()
self.logger.info('Astropy J2000: {}'.format(
J2000_coordinate.to_string(precision=2, sep=':')))
self.logger.info(
'pyephem Jnow: {} {}'.format(ephem_coordinate.ra, ephem_coordinate.dec))
self.logger.info('RA decimal = {:f}'.format(ephem_coordinate.ra))
self.logger.info('LST decimal = {:f}'.format(site.sidereal_time()))
self.logger.info('HA decimal = {:f}'.format(HourAngle))
return HourAngleimport
print "Masking %2.0f%% of sky"% ((1 - msk.sum() / float(len(msk)))*100)
except(AttributeError):
print "Weights not included in file. No mask will be applied."
data.shape = lats.shape
# Generate source locations
if not opts.src is None:
srclist,cutoff,catalogs = a.scripting.parse_srcs(opts.src, opts.cat)
if not opts.cal is None:
cat = a.cal.get_catalog(opts.cal, srclist, cutoff, catalogs)
else:
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()
def __init__(self, ra, dec, name='', epoch=ephem.J2000,
jys=0., index=-1, mfreq=.150,
ionref=(0.,0.), srcshape=(0.,0.,0.), **kwargs):
"""ra = source's right ascension (epoch=J2000)
dec = source's declination (epoch=J2000)
jys = source strength in Janskies at mfreq)
mfreq = frequency (in GHz) where source strength was measured
index = power-law index of source emission vs. freq."""
phs.RadioFixedBody.__init__(self, ra, dec, mfreq=mfreq, name=name,
epoch=epoch, ionref=ionref, srcshape=srcshape)
RadioBody.__init__(self, jys, index)
def __str__(self):
def convert_m(isys, osys, iepoch=e.J2000, oepoch=e.J2000):
"""Return the 3x3 matrix corresponding to a coordinate
system/precssion transformation (see 'convert').
NOTE: To obtain a transformation matrix to dot with isys to get osys,
must reverse isys and osys arguments in convert_m. """
m = np.array([[1,0,0],[0,1,0],[0,0,1]], dtype=np.double)
for i in range(3):
c = convert(m[:,i], isys, osys, iepoch=iepoch, oepoch=oepoch)
m[:,i] = radec2eq(c)
return m
def __init__(self, ra, dec, mfreq=.150, name='', epoch=ephem.J2000,
ionref=(0.,0.), srcshape=(0.,0.,0.), **kwargs):
RadioBody.__init__(self, name, mfreq, ionref, srcshape)
ephem.FixedBody.__init__(self)
self._ra, self._dec = ra, dec
self._epoch = epoch
def __str__(self):
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()
srclist,cutoff,catalogs = a.scripting.parse_srcs(opts.src, opts.cat)
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:
def __init__(self, ra, dec, name='', epoch=ephem.J2000,
jys=0., index=-1, mfreq=.150,
ionref=(0.,0.), srcshape=(0.,0.,0.), **kwargs):
"""ra = source's right ascension (epoch=J2000)
dec = source's declination (epoch=J2000)
jys = source strength in Janskies at mfreq)
mfreq = frequency (in GHz) where source strength was measured
index = power-law index of source emission vs. freq."""
phs.RadioFixedBody.__init__(self, ra, dec, mfreq=mfreq, name=name,
epoch=epoch, ionref=ionref, srcshape=srcshape)
RadioBody.__init__(self, jys, index)
def __str__(self):
DIM = img.shape[0]
RES = 1. / (kwds['d_ra'] * a.img.deg2rad * DIM)
im = a.img.Img(DIM*RES, RES, mf_order=0)
tx,ty,tz = im.get_top(center=(DIM/2,DIM/2))
# Define a weighting for gridding data into the skymap
map_wgts = np.exp(-(tx**2+ty**2) / np.sin(opts.fwidth*a.img.deg2rad)**2)
map_wgts.shape = (map_wgts.size,)
valid = np.logical_not(map_wgts.mask)
map_wgts = map_wgts.compress(valid)
# Get coordinates of image pixels in the epoch of the observation
ex,ey,ez = im.get_eq(ra, dec, center=(DIM/2,DIM/2))
ex = ex.compress(valid); ey = ey.compress(valid); ez = ez.compress(valid)
img = img.flatten(); img = img.compress(valid)
# Precess the pixel coordinates to the (J2000) epoch of the map
m = a.coord.convert_m('eq','eq',
iepoch=kwds['obs_date'], oepoch=ephem.J2000)
ex,ey,ez = np.dot(m, np.array([ex,ey,ez]))
# Put the data into the skymap
skymap.add((ex,ey,ez), map_wgts, img)
skymap.to_fits(opts.map, clobber=True)