How to use the ephem.J2000 function in ephem

# -*- 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)