How to use the ephem.Equatorial function in ephem

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').