TOOLS  FOR  FOUR CARDINAL STARS

How to Compute the Exact Position of a Star

Haluk Akcam - Jan. 01, 2004

Here you will find some of the computation methods used in preparing the research article "Four Cardinal Stars". Since the article itself is a summary for the public interest, the tools presented here for downloading are also in demo version.

Proper Motion

Proper motion is the apparent angular motion of a star on the celestial sphere. It results both from the actual movement of the star in space (peculiar motion), and from the star's motion relative to the solar system. The components of proper motion are m_a and m_d in equatorial coordinate system. For short time intervals such as t<100 years, the difference may be negligible, and can be calculated by applying these approximate formulae:

a_t = a₀ + m_at    and    d_t = d₀ + m_dt    where    t = t - t₀

[Note that m_a is given in the form m_a × cos d in most of the catalogs, and sometimes in mas. Therefore, for the above formula, the given value needs to be divided by cos d factor and then converted into arcsec. See the example below.]

But, if the time interval is covering a long span, these formulae are becoming insufficient, and we need more elaborated rigorous formulae. In that case, we need to define the movement of the star by introducing two additional parameters; the radial velocity and the distance. Radial velocity is the velocity of a star along the line of sight of an observer. A star's velocity in space relative to the sun is the space velocity (V), and has two components; radial velocity (V_R) and tangential velocity (V_T), where V² = V_R² + V_T², which are mostly given in km/s units.

Distance (d) of a star is given usually in parsecs (pc), a unit of length that is the distance at which one AU subtends an angle of one second of arc (i.e. parallax of one second). If the distance is given in light years, we should remember that one parsec is equal to 3.261 563 7941 light years. Some catalogs give only the parallax (p) of the star. Then we need to convert it to distance by using the parallax formula p = 1 / d. (Parallaxes less than about 0.01 arcsecond in older catalogs, or 0.001 arcsecond in the Hipparcos Catalog, are usually poor, and should be used with caution.)

Tangential or transverse velocity (V_T) has two components (proper motions in arcsecs and distance in parsecs):

V_T² = V_TA² + V_TD²    and    V_TA = m_a × d × k₁    and    V_TD = m_d × d × k₁

Conversion constant k₁ is required here to convert all terms into km/s unit.

(conversion constant)    k₁ = 149,597,871.475 / ( 365.25 × 86,400 ) = 4.740 47049

Radial velocity and the two components of tangential velocity help to get the Cartesian velocity coordinates v_x , v_y , and v_z , as follows:

v_x = (V_R cos d₀ cos a₀) - (V_TA sin a₀) - (V_TD sin d₀ cos a₀)

v_y = (V_R cos d₀ sin a₀) + (V_TA cos a₀) - (V_TD sin d₀ sin a₀)

v_z = V_R sin d₀ + V_TD cos d₀

Now we get the three initial Cartesian coordinates for the star at the epoch t₀, using the same coordinate system:

x₀ = d cos d₀ cos a₀    and    y₀ = d cos d₀ sin a₀    and    z₀ = d sin d₀

For practical purposes, we assume that stellar motion is linear. Therefore, after an interval of time t = t - t₀ (in years) , we get the new positions of the star as:

x_t = x₀ + v_x ( t / k₂ )    and    y_t = y₀ + v_y ( t / k₂ )    and    z_t = z₀ + v_z ( t / k₂ )

Here is required the conversion constant k₂ to convert the terms into parsec unit.

(conversion constant)    k₂ = ( 149,597,871.475 cosec 1'' ) / ( 365.25 × 86,400 ) = 977,792.226 75013

Finally, we get the new (epoch t) equatorial coordinate components a_t and d_t by using these equations:

d_xy = (x_t² + y_t²)^½    and    d_t = arctg (z_t / d_xy)    and    a_t = arctg (y_t / x_t)

With these equations, d_t product bears the same sign of z_t, and a_t product designation follows these rules:

x_t > 0 and y_t > 0 results in 1st quadrant; a_t is the same of the product (which is a positive number)

x_t < 0 and y_t > 0 results in 2nd quadrant; a_t is 180° + product (which is a negative number)

x_t < 0 and y_t < 0 results in 3rd quadrant; a_t is 180° + product (which is a positive number)

x_t > 0 and y_t < 0 results in 4th quadrant; a_t is 360° + product (which is a negative number)

Thus we can plot the star on a chart by applying d_t and a_t of the epoch t, but with the fixed reference coordinates of the epoch t₀. To shift the reference coordinates to the same epoch t, we need to compute the reduction amount for the precession of the equinoxes for the given time interval.

Precession of the Equinoxes

The slow change in the direction of the axis of rotation of the Earth due to the application of solar, lunar and planetary forces, is called precession of the equinoxes, with a period of about 25,800 years. The precession results primarily from the gravitational attraction of the Sun and the Moon, and secondarily from the attraction of the planets on the equatorial bulge. This gradual westward shifting of the vernal point consequently changes the stellar coordinates both in equatorial and ecliptic coordinate systems. Therefore, corrections become necessary for any given epoch on the coordinate parameters of stars, to determine their exact positions after the appropriate reference system.

For short time periods, approximate reductions may be sufficient to obtain the coordinates of the star for the new epoch t by applying these formulae (all values are in degrees of arc):

a_t = a₀ + M + N sin a_m tan d_m    and    d_t = d₀ + N cos a_m

a_m = a_t - ½ ( M + N sin a_t tan d_t )    and    d_m = d_t - ½ N cos a_m

a_m = a₀ + ½ ( M + N sin a₀ tan d₀ )    and    d_m = d₀ + ½ N cos a_m

M = 1.281 2323 T + 0.000 3879 T² + 0.000 0101 T³    and    N = 0.556 7530 T - 0.000 1185 T² - 0.000 0116 T³

Where T = ( t - 2000.0 ) / 100 = ( JD - 2,451,545.0 ) / 36,525 and subscript zero refers to epoch J2000.0.

These formulae are derived from L77 model, and good for a few decades, but not sufficient for long term periods. In order to find the exact coordinates, for the reduction of mean equatorial positions from initial epoch t₀ to epoch of date t, and vice versa, the following rigorous formulae must be used:

sin ( a_t - z_A ) cos d_t = sin ( a₀ + z_A ) cos d₀

cos ( a_t - z_A ) cos d_t = cos ( a₀ + z_A ) cos q_A cos d₀ - sin q_A sin d₀

sin d_t = cos ( a₀ + z_A ) sin q_A cos d₀ + cos q_A sin d₀

or

sin ( a₀ + z_A ) cos d₀ = sin ( a_t - z_A ) cos d_t

cos ( a₀ + z_A ) cos d₀ = cos ( a_t - z_A ) cos q_A cos d_t + sin q_A sin d_t

sin d₀ = - cos ( a_t - z_A ) sin q_A cos d_t + cos q_A sin d_t

with respect to the standard epoch t₀ = J2000.0 where

z_A = 2.72767 + 2306.080472 T + 0.3023262 T² + 0.01801752 T³ - 5.708×10^-6 T⁴ - 3.04×10^-7 T⁵ - 1.3×10^-10 T⁶

z_A = -2.72767 + 2306.07607 T + 1.0956768 T² + 0.01826676 T³ - 2.8276×10^-5 T⁴ - 2.486×10^-7 T⁵ - 5×10^-11 T⁶

q_A = 2004.190936 T - 0.426698 T² - 0.04182364 T³ - 7.291×10^-6 T⁴ - 1.127×10^-7 T⁵ + 3.6×10^-10 T⁶ + 9×10^-12 T⁷

and T = ( t - 2000.0 ) / 100 = ( JD - 2,451,545.0 ) / 36,525

Note: These equatorial precession angle expressions (z_A, z_A, and q_A in seconds of arc) are different than L77 (IAU 1976) model of Lieske. They have been first proposed as quantities of precession consistent with the IAU 2000A model, and were provided in the IERS Conventions 2000, based on the MHB corrections to precession rates, which complies with the recommendations of IAU 2000 Resolution B1.6. But in 2003, P. Bretagnon has suggested new terms, which I find more reliable than those proposed after the IAU 2000A model. Therefore what is given here are not those provided by the IERS, but Bretagnon's new expressions.

After getting the arctangent of a_t, we need to determine the actual quadrant of the angle, for it is not easy to estimate it, and also it is a requirement if we use above equations in a program. For this purpose, we follow the procedure:

A = sin ( a₀ + z_A ) cos d₀    and    B = cos ( a₀ + z_A ) cos q_A cos d₀ - sin q_A sin d₀

x = ( B × cos z_A ) - ( A × sin z_A )    and    y =  ( A × cos z_A ) + ( B × sin z_A )

where    tan a_t = y / x    and:

x > 0 and y > 0 results in 1st quadrant; a_t is the same of the product (which is a positive number)

x < 0 and y > 0 results in 2nd quadrant; a_t is 180° + product (which is a negative number)

x < 0 and y < 0 results in 3rd quadrant; a_t is 180° + product (which is a positive number)

x > 0 and y < 0 results in 4th quadrant; a_t is 360° + product (which is a negative number)

Since -90° < d_t < +90°, d_t bears the same sign of arcsin(d_t), and there is no need for further explanation. Same procedure is also necessary for a₀ , by applying the second part of above equations, namely by using the z_A parameter:

C = sin ( a_t - z_A ) cos d_t    and    D = cos ( a_t - z_A ) cos q_A cos d_t + sin q_A sin d_t

x = ( D × cos z_A ) + ( C × sin z_A )    and    y = ( C × cos z_A ) - ( D × sin z_A )

where    tan a₀ = y / x    and again:

x > 0 and y > 0 results in 1st quadrant; a₀ is the same of the product (which is a positive number)

x < 0 and y > 0 results in 2nd quadrant; a₀ is 180° + product (which is a negative number)

x < 0 and y < 0 results in 3rd quadrant; a₀ is 180° + product (which is a positive number)

x > 0 and y < 0 results in 4th quadrant; a₀ is 360° + product (which is a negative number)

As above, d₀ bears the same sign of arcsin(d₀).

Apparent Magnitude and Distance

Besides the position of the star, its apparent (as seen from Earth) brightness varies also by the time due to its motion. Magnitude is an arbitrary measure of the brightness of the stars. Apparent magnitude (m) is the brightness of a star as seen from Earth, and absolute magnitude (M) is the true brightness of a star. It is assumed that the apparent magnitude of a star would be equal to the absolute magnitude, if the star were placed at a distance (d) of 10 parsecs from the Earth. (It is the distance at which a star would have a parallax of 0.1 seconds of arc.) The relation between these two values is shown by this formula: M = m + 5 - 5 log₁₀(d) or M = m + 5 + 5 log₁₀(p).

The new distance of the star at the epoch t can be determined by its Cartesian coordinates as it is described above in units of parsecs. Apparent magnitudes are usually given in the catalogs, next to the distances. Thus:

d_t² = x_t² + y_t² + z_t²    and    m_t = m₀ + 5 log₁₀(d_t / d₀)

For instance, a star at 100 pc with V_m = 3.00 would have a new apparent magnitude of  2.77, after approaching 10 pc to our system. But, another star at 20 pc with the same parameters would be then glittering with the new 1.49 V_m.

Ecliptic Coordinates

After getting the equatorial coordinates (a_t, d_t), it is easy to convert them to ecliptic coordinates (l_t, b_t). First, the obliquity of the ecliptic (e) must be computed for the epoch t, and then the following equations must be applied (numerical expressions are after Bretagnon-2003):

e = 84381.4088 - 46.836051 t - 1.667×10^-4 t² + 1.99911×10^-3 t³ - 5.23×10^-7 t⁴ - 2.48×10^-8 t⁵ - 3×10^-11 t⁶

Where e is in seconds of arc and T = ( t - 2000.0 ) / 100 = ( JD - 2,451,545.0 ) / 36,525.

Then, to convert equatorial to ecliptic coordinates:

tan Q = cotan d sin a    and    tan l = ( sin ( Q + e ) / sin Q ) tan a    and    tan b = cotan ( Q + e ) sin l

Or, to convert ecliptic to equatorial coordinates:

tan P = cotan b sin l    and    tan a = ( sin ( P - e ) / sin P ) tan l    and    tan d = cotan ( P - e ) sin a

All terms here are in degrees of arc, and arctangent values need to be placed in appropriate quadrants. For this purpose, we apply the following procedure to find the quadrant of l:

x = cos a    and    y = ( sin ( Q + e ) / sin Q ) sin a    and    tan l = y / x

x > 0 and y > 0 results in 1st quadrant; l is the same of the product (which is a positive number)

x < 0 and y > 0 results in 2nd quadrant; l is 180° + product (which is a negative number)

x < 0 and y < 0 results in 3rd quadrant; l is 180° + product (which is a positive number)

x > 0 and y < 0 results in 4th quadrant; l is 360° + product (which is a negative number)

Where b bears the same sign of tan(b).

Also, we apply the following procedure to find the quadrant of a:

x = cos l    and    y = ( sin ( P - e ) / sin P ) sin l    and    tan a = y / x

x > 0 and y > 0 results in 1st quadrant; a is the same of the product (which is a positive number)

x < 0 and y > 0 results in 2nd quadrant; a is 180° + product (which is a negative number)

x < 0 and y < 0 results in 3rd quadrant; a is 180° + product (which is a positive number)

x > 0 and y < 0 results in 4th quadrant; a is 360° + product (which is a negative number)

Where d bears the same sign of tan(d).

Now, it is time to see the applications with three examples:

Example #1

For the purpose of the article "Four Cardinal Stars", I have compiled appropriate data from 30 recent catalogs for comparison, which you may find here. As an example, Aldebaran is chosen here due to its relatively high radial velocity. The required raw data for the star are:

a = 68.98000195 ±0.85  (arcdegrees, error in 1000 mas and multiplied by cosd) Eq=J2000, Epoch=J1991.25

d = +16.50976164 ±0.53 (arcdegrees, error in 1000 mas) Equator=J2000, Epoch=J1991.25

m_a = +64.7 ±0.5 (in 1000 mas, and multiplied by cosd) Equator=J2000, Epoch=J1991.25

m_d = -187.2 ±1.3  (in 1000 mas)  Equator=J2000, Epoch=J1991.25

V_R = +54.3 ±0.7 km/s      p = 50.09 ±0.95 (in 1000 mas)      m = 0.867 ±0.017

Precise figures of additional constants can be found in any scientific reference source. But, for the convenience of the reader, it is practical to apply pi as p = 3.141 592 653 58979, and cosec(1") = 206,264.806 24790, which is 1/sin(1"). The derived distance of Aldebaran from its trigonometric parallax is in parsecs d = 20.0 ±0.4 pc, or 65.1 ±1.2 ly. As a reminder, it is good to work with error limits, for it gives a more realistic opinion regarding the derived results.

Following these figures, we may now compute the new position of Aldebaran for the spring of the year 10 BC, just because of practical purposes, since the difference between the epoch of catalogs and the date is straight 2000 in number. Further, it must be remembered that there is no calendar year as 0.

The approximate formulae are easy to apply and the results will be:

a_t = 68.98000195 + -2000 × 0.0647 / ( cos (16.50976164) × 3600 ) = 68.94251184 = 68°56'33.04"

d_t = +16.50976164 + -2000 × -0.1872 / 3600  = 16.61376164 = +16°36'49.54"

These are the Eq=J2000 coordinates of Aldebaran for the epoch 10 BC, which are to be converted to the equator of the date. But, first a short reminder: J2000 denotes that 1999 Julian years are passed since the fictitious zero point, and spring of 10 BC denotes that 9.75 Julian years are needed to reach that zero point. Thus, the sum of interval becomes 2008.75 Julian years. 

T = -20.0875    and    M = -25.66209835    and    N = -11.13756821

a_m = 54.56083939    and    d_m = 13.38476810    and    a_t = 41.12119996°    and    d_t = +10.15577456°

Approximation shows that Aldebaran was at 41°07'16.32" R.A. and +10°09'20.79" Declination for the equator and the epoch of 10 BC.

But the rigorous formulae give these results: (Note: Figures are taken from the computer program, which uses the accurate values, so that the reader may see the slight differences.) 

V_TA = 0.0647 × 19.96 × 4.74047049 = 6.12314684

V_TD = -0.1872 × 19.96 × 4.74047049 = -17.71643159

v_x = 18.67405810 - 5.71568380 + 1.80589133 = 14.76426563

v_y = 48.59687955 + 2.19633467 + 4.69960428 = 55.49281850

v_z = 15.43090305 - 16.98600672 = -1.55510367

x₀ = 6.86574784    and    y₀ = 17.86724231    and    z₀ = 5.67336188

x_t = 6.86574784 + 14.76426563 × 2000 / 977792.22675013 = 6.83554866

y_t = 17.86724231 + 55.49281850 × 2000 / 977792.22675013 = 17.75373595

z_t = 5.67336188 - 1.55510367 × 2000 / 977792.22675013 = 5.67654273

d_xy = 19.02419159    and    d_t = 16.61433994 = +16°36'51.62"    and    a_t = 68.94228170 = 68°56'32.21"

When compared with the above approximated results, difference in R.A. for 2000 years is 0.83", but 2.08" in declination, because of the apparent relative trajectory of the star. Therefore, it is always safe to use the rigorous formulae.

Now, it is easy to calculate the new distance and apparent magnitude of Aldebaran, for the year 10 BC. It was slightly closer (0.11 pc) to the Earth with d_t = 19.85 pc, and a little more brighter with m_t = 0.855, some 2000 years ago.

Final stage is to convert these Eq=J2000 coordinates to the equator of the epoch, namely for the year 10 BC.

T = -20.0875    and    z_A = -12.87351578°    and    z_A = -12.78771425°    and    q_A = -11.13699118°

a_t = 41.09935645° = 41°05'57.68"    and    d_t = +10.20822845° = 10°12'29.62"

Here we have a 1.31 arcminutes difference in R.A., and 3.15 arcminutes difference in declination between two computations. Therefore, it is important not to apply approximate formulae for long term calculations.

After securing the exact position of Aldebaran for the epoch and the equator of the date, we may now find the ecliptic coordinates of the star, by applying the above formulae:

T = -20.0875    and    e_t = 23.69609750° = 23°41'45.95"

Q = 74.68042130    and    l_t = 41.82299120° = 41°49'22.77"    and    b_t = -5.60789008° = -05°36'28.40"

Here we see Aldebaran in 12° Taurus, about 6° below the apparent path of the Sun, and the obliquity of the ecliptic being a little steeper.

Example #2

Second example is about the position of Regulus for the date January 23, 5000 AD 12:00 UT. Here, the difference may seem to be 3000 years, 22 days, and 12 hours since the beginning of 2000, but actually it is not. Because, the time parameter of all calculations is expressed in Julian years, which is exactly 365.25 days. On the other hand, the civil or calendar year is 365 days, or 366 days in a leap year (any year divisible by 4 except centenary years not divisible by 400, in the Gregorian calendar, which is the one we use now). Therefore, since January 01, 2000 AD 12:00 UT (=J2000.0), which is expressed as JD 2451545.0, there are 1,095,750.0 days until January 23, 5000 AD 12:00 UT, which is JD 3547295.0. Thus the interval becomes exactly 3000 Julian years ( 3000 × 365.25 = 1095750 ).

In this example, we work with the error limits, and the required raw data for Regulus are:

a = 152.09358075 ±0.71  (arcdegrees, error in 1000 mas and multiplied by cosd) Eq=J2000, Epoch=J1991.25

d = +11.96719513 ±0.49 (arcdegrees, error in 1000 mas) Equator=J2000, Epoch=J1991.25

m_a = -248.7 ±0.4 (in 1000 mas, and multiplied by cosd) Equator=J2000, Epoch=J1991.25

m_d = +5.3 ±0.7  (in 1000 mas)  Equator=J2000, Epoch=J1991.25

V_R = +5.9 ±1.3 km/s      p = 42.09 ±0.79 (in 1000 mas)      m = 1.360 ±0.031

From the parallax (p), we obtain a distance between 23.32 - 24.21 pc, which can be defined as d = 23.8 ±0.4 pc.

Approximation shows that the position of Regulus on January 23, 5000 AD 12:00 UT is supposed to be:

a = 191°05'24.98"    and    d = -04°09'41.01"    or    l = 191°49'35.47"    and    b = +00°28'50.98"

But, the rigorous formulae give these results:

a = 191°05'39.66"    and    d = -04°07'35.16"    or    l = 191°48'53.20"    and    b = +00°30'53.00"

V_m = 1.36    and    d = 23.78 pc    or    d = 77.55 ly

First, we detect an error of 14.68" in R.A., and 125.85" in declination between two computations.

Now, we need to find the effect of error margins on the obtained results. We have two variations of six parameters to compute the limits both for R.A. and declination: Namely; a₀, d₀, m_a, m_d, V_R, and p. It means that we have to repeat the same procedure 2⁶ = 64 times, which is a time consuming operation, but we may create a simple computer program for it.

As a result, we obtain the following limits:

a is valid between these ranges:    191°05'38.23"    and    191°05'41.10"    ( 191°05'39.67" ±01.43" )

d is valid between these ranges:    -04°07'37.33"    and    -04°07'33.00"    ( -04°07'35.16" ±02.16" )

Further, we apply the same procedure for apparent magnitude and distance. But, here we have seven parameters (including m₀) with two variants, and the application should be therefore repeated 2⁷ = 128 times:

V_m is valid between these ranges:    1.393    and    1.330    ( 1.362 ±0.031 )

d is valid between these ranges:    23.335    and    24.235    ( 23.785 ±0.450 )    parsecs.

Example #3

In this final example, we work with error limits to find out the angular distance between two stars at a given time, yet by using long terms. The date is April 02, 126843 BC 07:30 UT, which is corresponding to ( JD-44,607,891.1875 ), or ( J-122129.75 ). Yet, numerical quantities about precession and the obliquity of the ecliptic are not precise enough even to estimate the positions for that remote date, due to lack of efficient theories. Therefore, we can calculate only the angular distances, but without any positioning for the equator of the date. Because of this incapability, all computations here are based on the J2000.0 equator reference system. Another problem is about the long-term validity of derived figures out of relative short time observations. It must be remembered that catalog data such as proper motions or radial velocities are the results of the last hundred years or so, which are mostly based on primitive ground observations. But, as an example for the application of error margins, we will try this remote date. Our stars are Aldebaran and Regulus.

Approximation shows:

a = 66°39'12.05"    and    d = +22°57'50.59"    for Aldebaran

a = 160°51'31.58"    and    d = +11°47'40.60"    for Regulus.

Rigorous formulae show:

a = 65°14'26.86"    and    d = +26°13'24.91"    for Aldebaran

a = 161°03'49.48"    and    d = +11°38'18.28"    for Regulus

As it is seen here, for long term computations, there is a good possibility to encounter tremendous errors when approximation is applied. Now let's see the error limits for rigorous formulae results:

Aldebaran:    a = 65.23806399 ±0.10391352    and    d = 26.22954131 ±0.24455669

where    V_m = -0.018 ±0.057    and    d = 13.29 ±0.50 pc

Regulus:    a = 161.06557094 ±0.08366046    and    d = 11.63812194 ±0.02865638

where    V_m = 1.316 ±0.048    and    d = 23.29 ±0.62 pc

Which means that we find Aldebaran in a spherical rectangle of about 11.2' × 29.3' size, and Regulus in 9.8' × 3.4'. That is, uncertainty about Aldebaran's position is ten times more than of Regulus, mostly due to relative high vertical motion and radial velocity of Aldebaran. Regarding to magnitude and distance; Aldebaran seems more brighter and closer than today, but Regulus shows no big difference. 

The angular distance between two stars is:

u = arccos [ sin ( d₁ + Dd₁ ) × sin ( d₂ + Dd₂ ) + cos ( d₁ + Dd₁ ) × cos ( d₂ + Dd₂ ) × cos ( a₁ + Da₁ -a₂ - Da₂ ) ]

Or in a more simplified but approximated way:

u = arccos [ sin d₁ × sin d₂ + cos d₁ × cos d₂ × cos ( a₁ -a₂ ) ]

Thus, the angular distance between two stars is found almost exactly 90°, without error limits. But, we need to apply the same formula four times to find the possible limits, which ends up with this result:

u = 90°00'11.45" ±12'21.67"

The reverse method is to find the temporal limits for the expected angular distance 90°, which is:

t_J = -122215 ± 2519

Finally, we understand that due to the error margins, the expected angular distance between two stars is not exactly in 126843 BC, but can be found some time in an interval of five thousand years.

Copyright © 2004-2008 Haluk Akcam. All rights reserved.

[back to list of tools] [main page] [table of contents]