How to use the skyfield.api.Topos function in skyfield

def calc_comet(comet_df, obstime, earthcoords, numdays=0, alt_table=False):
    # Generating a position.
    cometvec = sun + mpc.comet_orbit(comet_df, ts, GM_SUN)

    cometpos = earth.at(obstime).observe(cometvec)
    ra, dec, distance = cometpos.radec()
    print("RA", ra, "   DEC", dec, "   Distance", distance)

    if earthcoords:
        if len(earthcoords) > 2:
            elev = earthcoords[2]
        else:
            elev = 0
        obstopos = Topos(latitude_degrees=earthcoords[0],
                         longitude_degrees=earthcoords[1],
                         elevation_m=elev)
        print("\nObserver at",
               obstopos.latitude, "N  ", obstopos.longitude, "E  ",
              "Elevation", obstopos.elevation.m, "m")
        obsvec = earth + obstopos

        alt, az, distance = \
            obsvec.at(obstime).observe(cometvec).apparent().altaz()
        print("Altitude", alt, "     Azumuth", az, distance)

        alm_twilights = almanac.dark_twilight_day(eph, obstopos)
        dawn, dusk = find_twilights(obstime, alm_twilights)
        print(WHICH_TWILIGHT, ": Dawn", svt2str(dawn), "Dusk", svt2str(dusk))

        if numdays:

# alternatively under a commercial licence.
#
# The terms of the AGPL v3 license can be found in the main directory of this
# repository.

# Generate a list of ephemerides using Skyfield.
# The output of this script is used in the test in src/frame.c

from skyfield.api import load
from skyfield.api import Topos

planets = load('de421.bsp')
ts = load.timescale()
t = ts.utc(2018, 11, 28, 0, 0, 0)
earth = planets['earth'].at(t)
atlanta_p = (planets['earth'] + Topos('33.7490 N', '84.3880 W'))
atlanta = atlanta_p.at(t)

print('TAI = %r' % (t.tai - 2400000.5))
print('TT  = %r' % (t.tt  - 2400000.5))
print('UT1 = %r' % (t.ut1 - 2400000.5))

def dump_pv(p):
    print(' {{%.12f, %.12f, %.12f},\n'
          '  {%.12f, %.12f, %.12f}},' %
          (p.position.au[0], p.position.au[1], p.position.au[2],
          p.velocity.au_per_d[0], p.velocity.au_per_d[1], p.velocity.au_per_d[2]))

def dump_altazd(p):
    print(' {%.12f, %.12f, %.12f},' %
          (p[0].degrees, p[1].degrees, p[2].au))

for i in range(0,3):
            # (obs_id[i], object_number[i], ra_obs[i], decl_obs[i], station_number[i], lat[i], lon[i], elev[i], obs_name[i], obs_time_sql[i]) = observation_list[i]
            obs_id.append( observation_list[i]['obs_id'] )
            object_number.append( observation_list[i]['object_number'] )
            ra_obs.append( observation_list[i]['ra'] ) 
            decl_obs.append( observation_list[i]['declination'] ) 
            ra_obs_rad.append( radians(observation_list[i]['ra'] )) 
            decl_obs_rad.append( radians(observation_list[i]['declination'] )) 
            station_number.append( observation_list[i]['station_number'] )
            gslat.append( radians(observation_list[i]['latitude']) )
            gslon.append( radians(observation_list[i]['longitude']) )
            gsalt.append( observation_list[i]['elevation_m'] )
            obs_name.append( observation_list[i]['obs_name'] )
            obs_time_sql.append( observation_list[i]['obs_time'] )
            iod_string.append( observation_list[i]['iod_string'] )
            observer_location.append( Topos(
                latitude_degrees  = observation_list[i]['latitude'], 
                longitude_degrees = observation_list[i]['longitude'], 
                elevation_m       = observation_list[i]['elevation_m']) )

        try:
            epoch_time = obs_time_sql[0].strftime('%Y-%m-%d %H:%M:%S')

            TLE = db.selectTLEEpochBeforeDate(obs_time_sql[0], object_number[0])
            sat_name   = TLE[1]
            tle_line_1 = TLE[2]
            tle_line_2 = TLE[3]
            print("{}\n{}\n{}".format(sat_name,tle_line_1,tle_line_2))
        except:
            log.warning("NO TLE found for object '{}' with EPOCH before {}".format(object_number[0], obs_time_sql[0]))
            continue

def main():
    ts = load.timescale()
    t = ts.tt(2018, 1, 22, 9, 9, 20)
    trial_angles = 10, 20, 30, 40  # the error peaks around 20 degrees
    trial_elevations = 0, 6000, AU_M
    print(__doc__)
    for n in range(1, 5):
        print('=== {} iterations ==='.format(n))
        print('')
        for elevation_m in trial_elevations:
            for degrees in trial_angles:
                top = Topos(latitude_degrees=degrees, longitude_degrees=123,
                            elevation_m=elevation_m)
                xyz_au = top.at(t).position.au
                xyz_au = einsum('ij...,j...->i...', t.M, xyz_au)
                lat, lon, elev = reverse_terra(xyz_au, t.gast, n)
                lat = lat / DEG2RAD
                error_mas = 60.0 * 60.0 * 1000.0 * abs(degrees - lat)
                print('latitude {} degrees, elevation {} m'
                      ' -> error of {:.2f} mas'
                      .format(degrees, elevation_m, error_mas))
        print('')
    print("""\
Given that iterations=3 pushes the maximum error from tens of mas ("mas"

def main():

	#planets = load('de421.bsp')
	#earth = planets['earth']

	Site = CosparSite("data/sites.txt")
	duvall = Topos('47.77166666666667 N', '121.90416666666667 W')



	ts = load.timescale()

	#t = ts.now()

	#stations_url = 'http://celestrak.com/NORAD/elements/stations.txt'
	satellites_file = '../tle/bulk.tle'
	satellites = load.tle(satellites_file)

	#satellite = satellites['ARKYD 6A']
	satellite = satellites['TDRS 11']


	#astrometric = observer_location.at(t).observe(mars)