Coordinate conversions using TPM and PyTPM

Please read the TPM manual to learn the details of using TPM for performing coordinate conversions.

For a detailed list of functions and constants defined in PyTPM, see the section on Functions and constants in PyTPM. For detailed information on the data structures defined in PyTPM, see the section on Data Structures.

See Comparison of PyTPM with SLALIB for a comparison of results obtained using PyTPM, with those obtained using SLALIB.

In this document several examples of using PyTPM module for performing coordinate conversions will be presented.

In the following examples coordinates of M100 is converted between different systems. TPM can convert both positions and velocities. For M100 only the position are significant, except for the FK5-FK4 conversion. The section Converting positions and velocities has examples of converting positions and velocities, with data from HIPPARCOS catalog. Also, the examples here use the convenience function convert.convertv6(), instead of directly using the underlying TPM functions.

As mentioned before, facilities in TPM that are wrapped by PyTPM are provided in the module pytpm.tpm. The pytpm.convert module contains a few convenience functions for performing coordinate conversions. We can load both of these modules in the following manner.

>>> from pytpm import tpm, convert

Function convert.convertv6

The convert.convertv6() function aims to present an easy to use interface to TPM that allows easy coordinate conversions. It takes a large set of input, sets up TPM, performs the conversions, and returns the final coordinates. The source code of this function also serves as an example of using TPM.

This function takes the following arguments. Not all of these need to be specified for a given conversion.

Parameter Description
v6 one V6C vector or a list of V6C vectors.
s1 start state
s2 end state
epoch epoch of the coordinates as a Julian date
equinox equinox of the coordinates as Julian date
utc time of “observation” as a Julian date; exact meaning depends on the type of conversion; defaults to the epoch J2000.0
delta_ut UT1 - UTC in seconds.
delta_at TAI - UTC in seconds.
lon geodetic longitude in degrees
lat geodetic latitude in degrees
alt altitude in meters
xpole ploar motion in radians
ypole ploar motion in radians
T temperature in kelvin
P pressure in milli-bars
H relative humidity (0-1)
wavelength wavelength of observation in microns

A tpm.V6C object is a 6-D vector that stores positions and velocities in Cartesian coordinates.

The function convert.cat2v6() can be used to create a V6C object from catalog data. This function will accept scalar or list of coordinates, and will return a list of V6C vectors.

Function convert.v62cat() can convert V6C objects to catalog coordinates. The coordinates are returned as a dictionary. This function will also take a list of V6C objects and will return a list of dictionaries.

Coordinate systems in TPM are specified using integers or integer constants. These are referred to as states.The following lists the states defined in TPM. Note that each of these states are just integers: TPM_S02 == 2, TPM_S13 == 13 and so on.

State Description
TPM_S00 Null
TPM_S01 Heliocentric mean FK4 system, any equinox
TPM_S02 Heliocentric mean FK5 system, any equinox
TPM_S03 IAU 1980 Ecliptic system
TPM_S04 IAU 1958 Galactic system
TPM_S05 Heliocentric mean FK4 system, B1950 equinox
TPM_S06 Heliocentric mean FK5 system, J2000 equinox
TPM_S07 Geocentric mean FK5 system, J2000 equinox
TPM_S08 TPM_S07 + light deflection
TPM_S09 TPM_S08 + Aberration
TPM_S10 TPM_S09 + precession
TPM_S11 Geocentric apparent FK5, current equinox
TPM_S12 Topocentric mean FK5, J2000 equinox
TPM_S13 TPM_S12 + light definition
TPM_S14 TPM_S13 + aberration
TPM_S15 TPM_S14 + precession
TPM_S16 Topocentric apparent FK5, current equinox
TPM_S17 Topocentric apparent FK5, current equnix
TPM_S18 Topocentric apparent (Hour Angle, Declination)
TPM_S19 Topocentric observed (Azimuth, Elevation)
TPM_S20 Topocentric observed (Hour Angle, Declination)
TPM_S21 Topocentric observed WHAM (longitude, latitude)

Some of these states have additional special names.

Name State
TARGET_FK4 TPM_S01
TARGET_FK5 TPM_S02
TARGET_ECL TPM_S03
TARGET_GAL TPM_S04
TARGET_APP_HADEC TPM_S17
TARGET_OBS_HADEC TPM_S20
TARGET_APP_AZEL TPM_S18
TARGET_OBS_AZEL TPM_S19
TARGET_OBS_WHAM TPM_S21

FK5 equinox and epoch J2000.0, to FK4 equinox and epoch B1950.0

First obtain the FK5 equinox J2000.0 and epoch J2000.0 RA and Dec coordinates in radians.

>>> ra_j2000 = tpm.HMS(hh=12, mm=22, ss=54.899).to_radians()
>>> dec_j2000 = tpm.DMS(dd=15, mm=49, ss=20.57).to_radians()

Create a tpm.V6C vector for the object. Note that convert.cat2v6() will accept a list of coordinates as well.

>>> v6 = convert.cat2v6(ra_j2000, dec_j2000)

Now convert to FK4 equinox B1950.0, but remaining at epoch J2000.0.

In the following 6 stands for FK5 equinox J2000.0 coordinates and 5 stands for FK4 equinox B1950.0. The epoch and equinox are specified using epoch and equinox keywords. They can be interpreted in different ways depending on the exact conversion requested. They are also used only for certain calculations. In the present case, they are applicable to the input coordinates. Read the TPM manual for more information.

>>> v6_fk4 = convert.convertv6(v6, s1=6, s2=5, epoch=tpm.J2000,
   ...: equinox=tpm.J2000)

Convert V6C to catalog data and print results. Function convert.v62cat() will also accept a list of V6C vectors.

>>> d = convert.v62cat(v6_fk4, C=tpm.CB)
>>> print tpm.HMS(r=d['alpha'])
 12H 20M 22.935S
>>> print tpm.DMS(r=d['delta'])
+16D 05' 58.024"

The parameter C is the number of days in a century. The velocities in tpm.V6C vectors are stored AU/day. These Cartesian velocities must be converted to catalog proper motions. The latter are expressed in arc-secs/century.

In the Besselian system a century has approximately 36524.22 days, where as in the Julian system a century has exactly 36525.0 days. The former is used in FK4 and the latter is used in FK5. The default value is set to 36525.0. These two values are provided as the constants tpm.CB and tpm.CJ, respectively; see Constants.

Note that the results above do not agree with the FK4 values given by SIMBAD. This is because the results are for the epoch J2000.0, i.e., FK4 equinox B1950.0 epoch J2000.0. Even though the object doesn’t have proper motion, the FK4 system is rotating with respect to FK5. This results in a fictitious proper motion in the FK4 system. We must apply proper motion from epoch J2000.0 to epoch B1950.0 to get the final result, i.e., FK4 equinox B1950.0 epoch B1950.0.

>>> v6_fk4_ep1950 = convert.proper_motion(v6_fk4, tpm.B1950, tpm.J2000)

Finally convert V6C to catalog data and print results. The final result is in FK4 equinox and epoch B1950.0. The final results agree with the values given by SIMBAD.

>>> d = convert.v62cat(v6_fk4_ep1950, C=tpm.CB)
>>> print tpm.HMS(r=d['alpha'])
 12H 20M 22.943S
>>> print tpm.DMS(r=d['delta'])
+16D 05' 58.241"

FK5 equinox and epoch J2000.0, to IAU 1958 Galactic System

The IAU 1958 galactic system is represented using state 4. The result below is for the epoch J2000.0, since the input coordinates are at epoch J2000.0. The epoch of the galactic coordinates given by SIMBAD is J2000.0. So the result obtained below is what we need, i.e., we don’t need to apply any proper motion corrections.

>>> ra_j2000 = tpm.HMS(hh=12, mm=22, ss=54.899).to_radians()
>>> dec_j2000 = tpm.DMS(dd=15, mm=49, ss=20.57).to_radians()
>>> v6 = convert.cat2v6(ra_j2000, dec_j2000)

>>> v6_gal = convert.convertv6(v6, s1=6, s2=4, epoch=tpm.J2000,
   ...: equinox=tpm.J2000)

>>> d = convert.v62cat(v6_gal)
>>> print tpm.r2d(d['alpha'])
271.136139562
>>> print tpm.r2d(d['delta'])
76.8988689751

IAU 1958 Galactic, to FK5 equinox and epoch J2000.0

Here we set the starting state to galactic i.e., 4 and the end state to FK5 equinox. Since the input coordinates are at epoch J2000.0, the final results will also be at epoch J2000.0, i.e., FK5 equinox and epoch J2000.0.

>>> gal_lon = tpm.d2r(271.1361)
>>> gal_lat = tpm.d2r(76.8989)
>>> v6 = convert.cat2v6(gal_lon, gal_lat)

>>> v6_fk5 = convert.convertv6(v6, s1=4, s2=6, epoch=tpm.J2000)

>>> d = convert.v62cat(v6_fk5)
>>> print tpm.HMS(r=d['alpha'])
 12H 22M 54.900S
>>> print tpm.DMS(r=d['delta'])
+15D 49' 20.683"

The results are consistent with the accuracy of the input galactic coordinates. When converting from Galactic to FK5, we pass through the FK4 system. The lack of proper motion information in the Galactic system affects the accuracy here. Usually, in catalogs we are interested in converting FK5 to Galactic, and we do not notice this.

FK5 equinox and epoch J2000.0, to IAU 1980 Ecliptic system

The ecliptic system is indicated using the state 3. Here the epoch of the output ecliptic coordinates will be J2000.0, as the input coordinates are at epoch J2000.0.

>>> ra_j2000 = tpm.HMS(hh=12, mm=22, ss=54.899).to_radians()
>>> dec_j2000 = tpm.DMS(dd=15, mm=49, ss=20.57).to_radians()
>>> v6 = convert.cat2v6(ra_j2000, dec_j2000)

>>> v6_ecl = convert.convertv6(v6, s1=6, s2=3, epoch=tpm.J2000,
   ...: equinox=tpm.J2000)

>>> d = convert.v62cat(v6_ecl)
>>> print tpm.r2d(d['alpha'])
178.78256462
>>> print tpm.r2d(d['delta'])
16.7597002513

The results agree with the results from the SLALIB (pyslalib) routine sla_eqecl.

See Comparison of PyTPM with SLALIB for a more detailed comparison between PyTPM and SLALIB.

IAU 1980 Ecliptic system, to FK5 equinox and epoch J2000.0

The starting state is set to 3 for ecliptic and the end state is set to 6 for FK5 equinox J2000.0.

>>> ecl_lon = tpm.d2r(178.78256462)
>>> ecl_lat = tpm.d2r(16.7597002513)
>>> v6 = convert.cat2v6(ecl_lon, ecl_lat)

>>> v6_fk5 = convert.convertv6(v6, s1=3, s2=6, epoch=tpm.J2000)

>>> d = convert.v62cat(v6_fk5)
>>> print tpm.HMS(r=d['alpha'])
 12H 22M 54.898S
>>> print tpm.DMS(r=d['delta'])
+15D 49' 20.570"

FK5 equinox and epoch J2000.0, to Geocentric apparent

Geocentric apparent RA & Dec. for midnight of 2010/1/1 can be calculated as shown below. The state identification number for geocentric apparent position is 11, as shown in the table above.

Obtain UTC and TDB time for the time of observation.

>>> utc = tpm.gcal2j(2010, 1, 1) - 0.5  # midnight
>>> tdb = tpm.utc2tdb(utc)

Obtain coordinates and V6C vector.

>>> ra_j2000 = tpm.HMS(hh=12, mm=22, ss=54.899).to_radians()
>>> dec_j2000 = tpm.DMS(dd=15, mm=49, ss=20.57).to_radians()
>>> v6 = convert.cat2v6(ra_j2000, dec_j2000)

Apply proper motion from epoch J2000.0 to epoch of observation. In this example, this is not needed since proper motion is zero. But we do this for completeness. The result is FK5 J2000 current epoch.

>>> v6 = convert.proper_motion(v6, tt, tpm.J2000)

Convert coordinates from FK5 equinox J2000, current epoch to FK5 equinox and epoch of date.

>>> v6_gc = convert.convertv6(v6, s1=6, s2=11, utc=utc)
>>> d = convert.v62cat(v6_gc)
>>> print tpm.r2d(d['alpha'])
185.860038856
>>> print tpm.r2d(d['delta'])
15.7631353482

The result from SLALIB (pyslalib) for the equivalent conversion, using the sla_map function is given below.

>>> utc = slalib.sla_caldj(2010, 1, 1)[0]  # midnight
>>> tt = slalib.sla_dtt(utc) / 86400.0 + utc

>>> r, d = slalib.sla_map(ra_j2000, dec_j2000, 0, 0, 0, 0.0, 2000.0,
   ...: tt)

>>> tpm.r2d(r)
185.86002229414245
>>> tpm.r2d(d)
15.763142468669891

The difference is about 0.06 arc-sec in RA and about 0.03 arc-sec in Dec.:

>>> (tpm.r2d(r) - 185.860038856) * 3600.0
-0.059622687126648088
>>> (tpm.r2d(d) - 15.7631353482) * 3600.0
0.025633691604554087

See Comparison of PyTPM with SLALIB for a more detailed comparison between PyTPM and SLALIB.

FK5 equinox and epoch J2000.0, to topocentric observed

Topocentric observed azimuth and elevation (and zenith distance) for an observer at the default location (KPNO) is calculated for 2010/1/1 mid-day. The final state i.e., apparent topocentric Az & El, is 19.

For midnight 2010/1/1 this object is below the horizon and hence the refraction calculations are not reliable. So we use mid-day for the following example.

>>> utc = tpm.gcal2j(2010, 1, 1)  # mid-day
>>> tt = tpm.utc2tdb(utc)

>>> ra_j2000 = tpm.HMS(hh=12, mm=22, ss=54.899).to_radians()
>>> dec_j2000 = tpm.DMS(dd=15, mm=49, ss=20.57).to_radians()
>>> v6 = convert.cat2v6(ra_j2000, dec_j2000)

>>> v6 = convert.proper_motion(v6, tt, tpm.J2000)

>>> v6_app = convert.convertv6(v6, s1=6, s2=19, utc=utc)

>>> d = convert.v62cat(v6_app)
>>> print tpm.r2d(d['alpha']), 90 - tpm.r2d(d['delta'])
133.49820871 22.0162437585

To calculate the observed hour angle and declination the v6_app vector obtained above can be used as input. We don’t need to go back to the FK5 equinox and epoch J2000.0 values. The input state is now 19 and the output, i.e., topocentric observed HA & Dec, is 20.

>>> v6_hadec = convert.convertv6(v6_app, s1=19, s2=20, utc=utc)

>>> d = convert.v62cat(v6_hadec)
>>> print tpm.r2d(d['alpha'])
343.586827647
>>> print tpm.r2d(d['delta'])
15.7683070508

To calculate the observed RA we need to find the LAST, since TPM only provides apparent RA. The observed RA can be found by subtracting hour angle from LAST. This is one situation where we need to access the underlying TPM machinery provided in pytpm.tpm.

>>> tstate = tpm.TSTATE()
>>> tpm.tpm_data(tstate, tpm.TPM_INIT)
>>> tstate.utc = utc
>>> tstate.delta_ut = tpm.delta_UT(utc)
>>> tstate.delta_at = tpm.delta_AT(utc)
>>> tstate.lon = tpm.d2r(-111.598333)
>>> tstate.lat = tpm.d2r(31.956389)
>>> tpm.tpm_data(tstate, tpm.TPM_ALL)
>>> last = tpm.r2d(tpm.r2r(tstate.last))
>>> last - tpm.r2d(d['alpha']) + 360.0
185.85569737491355

The same calculation with SLALIB, using sla_aop produces results that agree with PyTPM.

>>> dut = tpm.delta_UT(tpm.gcal2j(2010, 1, 1))  # DUT for mid-day.
>>> utc = slalib.sla_caldj(2010, 1, 1)[0] + 0.5  # mid-day.
>>> tt = slalib.sla_dtt(utc) / 86400.0 + utc

>>> r, d = slalib.sla_map(ra_j2000, dec_j2000, 0, 0, 0, 0.0, 2000.0,
   ...: tt)

>>> lon = tpm.d2r(-111.598333)
>>> lat = tpm.d2r(31.956389)

>>> az, zd, ha, dec, ra = slalib.sla_aop(r, d, utc, dut, lon, lat,
   ...: 2093.093, 0, 0, 273.15, 1013.25, 0, 0.550, 0.0065)

>>> tpm.r2d(tpm.r2r(az)), tpm.r2d(tpm.r2r(zd))
133.498195532 22.0162383595

The hour angle, declination and right ascension are:

>>> print tpm.r2d(tpm.r2r(ha))
343.586827289
>>> print tpm.r2d(tpm.r2r(dec))
15.7683143606
>>> print tpm.r2d(tpm.r2r(ra))
185.855680678

See Comparison of PyTPM with SLALIB for a more detailed comparison between PyTPM and SLALIB.

Converting positions and velocities

Converting position and velocities involve the exact same procedures as in the above sections. The only difference is that we provide velocities in addition to positions and the output will have both the converted positions as well as the converted velocities.

All these files mentioned below are available in the examples/ directory of the documentation source code (_downloads in the HTML distribution).

The examples shown here are available in the examples/conversions.py file.

The file examples/hip_full.txt contains data from the HIPPARCOS catalog. This file was generated using the examples/hip_full.py script. We will convert this data between different coordinate systems.

The get_hipdata function below is one way of reading in the data in hip_full.txt. With table reading packages such as asciitable, asciidata and atpy, or the numpy.loadtxt function, this function is not needed. The calculations are much easier to perform with Numpy. PyTPM does not need Numpy or table reading software. Hence I want to have examples that don’t use these, even though the calculations are harder.

import csv
import math
from pytpm import tpm, convert

def get_hipdata():
    """Return data in hip_full.txt.

    The data was created with hip_full.py file. Assumes that the file
    hip_full.txt is in the current directory.

    A dictionary is returned. All positions are in radians. Proper
    motions are in arc-sec per Julian century. Parallax is in
    arc-sec. The keys are:

      ra_icrs: ICRS RA
      dec_icrs: ICRS Dec
      raj2: RA J2000
      decj2: Dec J2000
      rab1: RA B1950
      decb1: Dec B1950
      glon: Galactic longitude (epoch J2000)
      glat: Galactic latitude (epoch J2000)
      elon2: Ecliptic longitude (epoch J2000)
      elat2: Ecliptic latitude (epoch J2000)
      pma: proper motion in RA (without cos(dec) factor)
      pmd: proper motion in Dec
      px: parallax.

    """
    f = open("hip_full.txt", "r")
    s = csv.reader(f, quoting=csv.QUOTE_NONNUMERIC,
                   delimiter=" ", skipinitialspace=True)
    d = dict(ra_icrs=[], dec_icrs=[], px=[], pma=[], pmd=[],
             raj2=[], decj2=[], rab1=[], decb1=[], glon=[],
             glat=[], elon2=[], elat2=[])

    for i in s:
        d["ra_icrs"].append(math.radians(i[0]))
        d["dec_icrs"].append(math.radians(i[1]))
        d["raj2"].append(math.radians(i[5]))
        d["decj2"].append(math.radians(i[6]))
        d["rab1"].append(math.radians(i[7]))
        d["decb1"].append(math.radians(i[8]))
        d["glon"].append(math.radians(i[9]))
        d["glat"].append(math.radians(i[10]))
        d["elon2"].append(math.radians(i[11]))
        d["elat2"].append(math.radians(i[12]))

        # milli-arc-sec/jul yr to arc-sec per Jul. cent.
        # And take out cos in RA proper motion.
        d["pma"].append(i[3] / math.cos(d["decj2"][-1]) / 1000.0 * 100.0)
        d["pmd"].append(i[4] / 1000.0 * 100.0)

        # milli-arsec to arc-sec
        d["px"].append(i[2] / 1000.0)

    f.close()
    return d

First read in the data and create tpm.V6C vectors. Since we don’t have radial velocities, we will set it to 0.0.

hip_tab = get_hipdata()
# Dummy radial velocities.
rv = [0.0 for i in range(len(hip_tab['px']))]
# Create V6C vectors.
v6 = convert.cat2v6(hip_tab['raj2'], hip_tab['decj2'], hip_tab['pma'],
                hip_tab['pmd'], hip_tab['px'], rv, tpm.CJ)

Subsequent steps are identical to those for converting positions. For example for the FK5 to FK4 conversions we can do:

# Convert from FK5 equinox and epoch J2000.0 to FK4 equinox B1950, but
# at the given epoch i.e., J2000.0.
v6o = convert.convertv6(v6, s1=6, s2=5, epoch=tpm.J2000)

# Apply proper motion from J2000.0 to B1950.0. Objects with zero
# velocity in FK5 will have a fictitious proper motion in FK4.
v6o = convert.proper_motion(v6o, tpm.B1950, tpm.J2000)

# Convert V6C vectors into a list of dictionaries, each of which
# contain the 6-D Fk4 B1950.0 coordinates.
cat = convert.v62cat(v6o, tpm.CB)

The variable cat will be a list of dictionaries, each containing the FK4 coordinates of the corresponding entry in the hip_full.txt file.

We can print out the first of these and compare the result with the FK4 data in hip_full.txt.

>>> print tpm.HMS(r=hip_tab['rab1'][0])  # HIPPAROCS FK4 from Vizier.
 23H 57M 26.475S
>>> print tpm.HMS(r=c['alpha'])  # PyTPM output.
 23H 57M 26.474S
>>> print tpm.DMS(r=hip_tab['decb1'][0])
+00D 48' 38.264"
>>> print tpm.DMS(r=c['delta'])
+00D 48' 38.257"

Similarly FK5 to ecliptic conversion can be performed in the following manner.

# FK5 equinox and epoch J2000.0, to IAU 1980 ecliptic J2000.0
v6o = convert.convertv6(v6, s1=6, s2=3)
# Convert V6C vectors into a list of dictionaries, each of which
# contain the 6-D Fk4 B1950.0 coordinates.
cat = convert.v62cat(v6o, tpm.CJ)

We can compare the result with that in the HIPPARCOS catalog returned by Vizier.

>>> print tpm.HMS(r=hip_tab['elon2'][0])  # Ecliptic data from Vizier cat.
 00H 01M 44.172S
>>> print tpm.HMS(r=c['alpha'])  # PyTPM conversion of FK5.
 00H 01M 44.172S
>>> print tpm.DMS(r=hip_tab['elat2'][0])
+00D 59' 55.604"
>>> print tpm.DMS(r=c['delta'])
+00D 59' 55.604"

All other conversions can be performed in the same manner.

PyTPM and TPM

All the facilities of TPM can be used form within Python, using functions and classes.

The following are C and Python code that perform equivalent conversions. The code takes the FK5 J2000.0 equinox and epoch coordinates of M100 through all coordinate systems defined in TPM.

These, and the source code for the convert module, should allow translation of TPM C code to PyTPM Python code.

These code fragments can be found in examples/conversion_example.c and examples/conversion_example.py, respectively.

Note that the conversions themselves are not accurate, since we do not apply any proper motions. They are meant to compare C and Python code.

C code

#include "tpm/astro.h"
#include <stdio.h>
#include <math.h>

/* Take a coordinate through all states. */
/* Coordinates for M100 from SIMBAD. */

int main(){
  double ra = (12+22/60.0+54.899/3600.0) * (2*M_PI/24.0);
  double de = (15+49/60.0+20.57/3600.0) * (2*M_PI/360.0);
  double ra1, ra1_d, de1, de1_d;
  double ep = J2000;
  double eq = J2000;
  V6 v6;
  V6 pvec[N_TPM_STATES];
  TPM_TSTATE tstate;
  int s1 = TPM_S06; /* Heliocentric mean J2000 FK5 ~~ ICRS */
  int s2 = TPM_S00; /* Assign required states. */

  for(int i=TPM_S00; i < N_TPM_STATES; i ++){
    tpm_data(&tstate, TPM_INIT);
    tstate.utc = J2000;
    tstate.lon = d2r(-111.598333);
    tstate.lat = d2r(31.956389);
    tstate.alt = 2093.093;
    tstate.delta_ut = delta_UT(tstate.utc);
    tpm_data(&tstate, TPM_ALL);

    v6 = v6init(SPHERICAL);
    v6SetR(v6, 1e9);
    v6SetAlpha(v6, ra);
    v6SetDelta(v6, de);

    pvec[s1] = v6s2c(v6);
    s2 = i;
    tpm(pvec, s1, s2, ep, eq, &tstate);
    v6 = v6c2s(pvec[s2]);

    ra1 = v6GetAlpha(v6);
    de1 = v6GetDelta(v6);
    ra1_d = r2d(ra1);
    if (ra1_d < 0.0) ra1_d += 360.0;
    de1_d = r2d(de1);
    if (de1_d < 0.0) de1_d += 360.0;

    printf("%02d-%02d %-17s %s %s %8.4f %8.4f\n", s1, s2,
      tpm_state(s2), fmt_alpha(ra1), fmt_delta(de1), ra1_d, de1_d);
  }
  return 0;
}

PyTPM code

# Take coordinates of M100 through all states.
from pytpm import tpm

ra = tpm.h2r(12+22/60.0+54.899/3600.0)
de = tpm.d2r(15+49/60.0+20.57/3600.0)
ep = tpm.J2000
eq = tpm.J2000
s1 = tpm.TPM_S06
s2 = tpm.TPM_S00
tstate = tpm.TSTATE()
pvec = tpm.PVEC()

for i in range(tpm.N_TPM_STATES):
    tpm.tpm_data(tstate, tpm.TPM_INIT)
    tstate.utc = tpm.J2000
    tstate.lon = tpm.d2r(-111.598333)
    tstate.lat = tpm.d2r(31.956389)
    tstate.alt = 2093.093
    tstate.delta_ut = tpm.delta_UT(tstate.utc)
    tpm.tpm_data(tstate, tpm.TPM_ALL)

    v6 = tpm.V6S()
    v6.r = 1e9
    v6.alpha = ra
    v6.delta = de


    pvec[s1] = v6.s2c()
    s2 = i
    tpm.tpm(pvec, s1, s2, ep, eq, tstate)
    v6 = pvec[s2].c2s()

    ra1 = v6.alpha
    de1 = v6.delta
    ra1_d = tpm.r2d(ra1)
    if ra1_d < 0.0 : ra1_d += 360.0
    de1_d = tpm.r2d(de1)
    if de1_d < 0.0 : de1_d += 360.0

    s = "{0:02d}-{1:02d} {2:<17s} {3:s} {4:s} {5:8.4f} {6:8.4f}"
    print s.format(s1, s2, tpm.tpm_state(s2),
                   tpm.fmt_alpha(ra1), tpm.fmt_delta(de1), ra1_d,
                   de1_d)

We create a state structure, tpm.TSTATE, and initialize it by calling tpm.tpm_data() with TPM_INIT. Then we assign values to independent parameters of the state data structure. We then calculate all dependent state properties by calling tpm_data() and passing TPM_ALL. We then create an array of tpm.V6C vectors, pvec (tpm.PVEC), create a V6C vector for our object, and assign it to the desired location in the array, based on the starting state. We then call tpm.tpm() with the state structure and the array of V6C vectors, along with the starting and ending state numbers. Finally we retrieve the appropriate V6C vector from the array, which will give us the final coordinates.

The result from running the above code is given below:

06-00 null               12H 22M 54.898S +15D 49' 20.570" 185.7287  15.8224
06-01 Helio. mean FK4    12H 22M 54.824S +15D 49' 20.447" 185.7284  15.8223
06-02 Helio. mean FK5    12H 22M 54.898S +15D 49' 20.570" 185.7287  15.8224
06-03 IAU 1980 Ecliptic  11H 55M 07.815S +16D 45' 34.920" 178.7826  16.7597
06-04 IAU 1958 Galactic  18H 04M 32.673S +76D 53' 55.928" 271.1361  76.8989
06-05 Helio. mean FK4    12H 20M 22.935S +16D 05' 58.024" 185.0956  16.0995
06-06 Helio. mean FK5    12H 22M 54.898S +15D 49' 20.570" 185.7287  15.8224
06-07 Geoc. mean FK5     12H 22M 54.899S +15D 49' 20.569" 185.7287  15.8224
06-08 S07 + Light Defl.  12H 22M 54.898S +15D 49' 20.571" 185.7287  15.8224
06-09 S08 + Aberration   12H 22M 54.995S +15D 49' 13.474" 185.7291  15.8204
06-10 S09 + Precession   12H 22M 54.995S +15D 49' 13.474" 185.7291  15.8204
06-11 Geoc. app. FK5     12H 22M 54.045S +15D 49' 19.561" 185.7252  15.8221
06-12 Topo. mean FK5     12H 22M 54.899S +15D 49' 20.569" 185.7287  15.8224
06-13 S12 + Light Defl.  12H 22M 54.898S +15D 49' 20.571" 185.7287  15.8224
06-14 S13 + Aberration   12H 22M 55.013S +15D 49' 13.452" 185.7292  15.8204
06-15 S14 + Precession   12H 22M 55.013S +15D 49' 13.452" 185.7292  15.8204
06-16 Topo. app. FK5     12H 22M 54.063S +15D 49' 19.539" 185.7253  15.8221
06-17 Topo. app. HA/Dec  22H 52M 35.524S +15D 49' 19.539" 343.1480  15.8221
06-18 Topo. app. Az/El   08H 50M 11.837S +67D 45' 09.683" 132.5493  67.7527
06-19 Topo. obs. Az/El   08H 50M 11.837S +67D 45' 34.371" 132.5493  67.7595
06-20 Topo. obs. HA/Dec  22H 52M 36.636S +15D 49' 38.307" 343.1527  15.8273
06-21 Topo. obs. WHAM    22H 52M 56.457S -14D 49' 46.993" 343.2352 345.1703