PyTPM 0.7.1

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 V6C vector for the object. Note that pytpm.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 and epoch J2000.0 coordinates and 5 stands for FK4 equinox and epoch B1950.0. The epoch and equinox are specified using epoch and equinox keywords. But they can be interpreted in different ways depending on the exact conversion requested. In this case, they are applicable to the input coordinates.

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

Convert V6C to catalog data and print results. Function pytpm.convert.v62cat will also accept a list of V6C objects.

>>> 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 AU/day must be converted into “/century. In the Besselian system, a century has approximately 36524.22 days, where as in the Julian system a century has 36525.0 days. The former is used in FK4 and the latter is used in FK5. The default value is set to 36525.0.

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

>>> 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 to IAU 1958 Galactic System

The IAU 1958 galactic system is represented using state 4. The result below is for the 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.

FK5 equinox and epoch J2000 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.

>>> 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 form the SLALIB (pyslalib) routine sla_eqecl.

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 and epoch 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 to Geocentric apparent

Geocentric apparent RA & Dec. for midnight of 2010/1/1 is calculated as shown below. The state identification number for geocentric apparent position is 11.

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

FK5 equinox and epoch J2000 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. Please consult the TPM manual and the PyTPM documentation for more information.

>>> 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

Consult the appropriate section of the PyTPM documentation for a detailed comparison between PyTPM and SLALIB.

Converting positions and velocities

Converting positions and velocities follow exactly the same procedure as the examples shown above. The convert.cat2v6 function will take proper motions, radial velocity and parallax in addition to position. The returned dictionary will have appropriate fields for final proper motions, radial velocity and parallax.

See the file doc/examples/conversions.py for a full example. The file is also included with the HTML documentation and with the source distribution.

For example if tab is a table that contains full 6-D coordinates with keys ra, dec, pma, pmd, px and rv, then a full V6C vector can be constructed as:

>>> v6 = convert.cat2v6(tab['ra'], tab['dec'], tab['pma'],
   ...: tab['pmd'], tab['px'], rv, tpm.CJ)

See docstring of the convert.convertv6 function for the required units for each of these.

To convert this from, say FK5 to Ecliptic, at the same epoch, we can use:

>>> v6o = convert.convertv6(v6, s1=6, s2=3)
>>> cat = convert.v62cat(v6o)

The variable cat will contain a dictionary, or a list of dictionaries, with the relevant catalog quantities. See the docstring of this convert.v62cat for units of output quantities.