Precessing positions from B1950 to J2000

[Root ]

Contents

Overview

[Top]

This page provides details of a method for precessing the coordinates of objects from the B1950.0 epoch to the now standard J2000.0 epoch, ignoring the proper motion of the stars. The algorithm used can handle coordinates with declinations up to +90 or -90 degrees (unlike many approximate methods) and the QBASIC program provided will accept input and provide output in the notation used in Burnham's Celestial Handbook, which uses B1950.0 coordinates. I have provided a DOS version of the program compiled with the excellent FirstBasic compiler, for those who need to convert a number of positions.

The celestial sphere can be mapped by using two main systems of coordinates, each referred to a reference plane. The celestial equator is the reference plane for equatorial coordinates (right ascension and declination), and the ecliptic plane is the basis of ecliptic coordinates (usually heliocentric longitude and latitude). The ecliptic and equator are represented as two great circles in the celestial sphere, and they meet at two points on the sphere. One of these points is used as the zero of right ascension. Declination has a natural zero at the celestial equator.

The Earth spins like a top on its own axis as it orbits the Sun. The Moon and the Sun exert tidal forces on the equatorial bulge of the Earth, and there are small effects from the rest of the planets in the solar system. This complex set of forces results in a movement in the celestial equator with time, which causes the zero point for RA to shift back along the ecliptic, and causes the tilt of the equator to change. This combined movement has short and long period components which are ususally broken down into 'corrections' as follows;

The two precession effects are lumped together into the general precession. The motion of the equator means that you must quote a date with the RA and DEC for an object. If you take account of nutation, then your positions are apparent positions. If you leave out nutation, then your positions are referred to the mean equator of the date you have chosen. All observations are referred to the equinox of the date, or epoch, chosen. Three main epochs are in common use;

Star charts, catalogues and observing lists use one of the epochs B1950.0 or J2000.0, and catalogues may list both. Star chart programs can produce coordinates refered to any epoch, but you should make sure you know which epoch the software is working in. The Yale Bright Star catalogue lists B1900.0 and J2000.0 positions. The 'B' and the 'J' refer to the use of the Besselian year and Julian year to count the epochs - the Besselian years are based on a solar year and are not now used. The Julian years are 365.25 days long, and can be related to the Julian day number easily.

The stars are not really fixed on the inside of a large sphere (even though it seems like that on a good dark night). They orbit the galactic centre and drift in their own local groups and clusters. On time scales which we can comprehend, this movement of the nearer stars can be modelled as the proper motion and is given in arcseconds per year (for both RA and DEC). The star with the largest proper motion is 61-cyg, at 4.1 and 3.2 arcseconds per year in RA and DEC. If you ignore the proper motion of this star when converting from B1950 to J2000, then you will have a J2000 RA which is about 0.2 minutes of time wrong, and a DEC which is about 3 arcminutes wrong. Most of the stars in the Bright Star Atlas have proper motions much less than this (about 17 more than 1 arcsec per year in RA out of 9000 odd stars). Deep space objects do not have proper motions that we can measure!

If you have an observing guide or comet ephemeris in B1950.0 coordinates, and you want to plot the positions of the object on a J2000.0 star atlas, you must correct the B1950 coordinates to J2000.0 by calculating the effects of precession. A rough correction, valid for shorter time periods, is often found in astronomical reference books - here is the recipe for B1950.0 to J2000.0;

1.  change B1950.0 RA and DEC to decimal degrees
2.  work out the following formulas

	 RA2000 = RA + 0.640265 + 0.278369 * sin(RA) * tan(DEC)
	DEC2000 = DEC + 0.278369 * cos(RA)
	
3.  convert RA back to hours, and format as Hrs, min, sec as 
    preferred.

The main problem with this recipe is that tan(DEC) will grow very large for stars with a high declination. You can't precess the coordinates of the Pole Star using the recipe above. For objects which stay close to the ecliptic (like planets and the stars in the 'occultation band') the recipe works well. A minor problem is that stars with RA's close to 24h in 1950 coordinates may produce a J2000 result higher than 24h, although that is easily taken care of.

A better method for use with programmable calculators or any form of computer program or spreadsheet is based on the idea of rotations around rectangular axes, so that no overflow can occur. This second method is developed in the next section.

As a practical application, I have written some format conversion routines which can use the notation found in Burnham's Celestial Handbook, a much loved observation guide for variable and double stars and deep space objects. Burnham's uses B1950.0 coordinates and does not include comprehensive star maps. If you want to plot the objects mentioned on a J2000.0 star atlas, then you may find the QBASIC program here useful.

Algorithm

[Top]

The method here is found in Duffett-Smith, and is also explained in Meeus and (in much more detail) in the Explanatory Supplement.

As a concrete example, take the variable star R-lyr for which Burnham's gives the B1950 coordinates as 18538n4353. Translating from this compact notation, this is RA 18h 53.8 minutes and DEC 43 deg 53 minutes North. Wil Tirion's The Cambridge Star Atlas gives the J2000.0 coordinates of this star as 18 h 55.3 min and +43 57, so we can check the calculation to moderate accuracy. Other variable stars in the Burnham's list for Lyra are not shown on Tirion's map, so we might need to plot them from our converted coordinates. The basic steps in the algorithm are;

  1. change the B1950 RA and DEC to decimal degrees or radians as needed
  2. work out the rectangular coordinates in a set of axes chosen to relate to the celestial pole and the zero of RA
  3. rotate the coordinates through three angles designed to bring the 1950 pole to the 2000 pole, and the 1950 RA zero point to the 2000 zero point. These angles depend on the time period between B1950 and J2000.
  4. work out the new RA and DEC
  5. change RA back to hours, and convert to the preferred format
The formulas for each step are outlined below, and the values are worked out for our example of the variable star R-lyr. I have used an ordinary scientific calculator to work out the values.

1. Change to decimal degrees

This is very simple;
     ra = (hrs + min/60) * 15   
	 
    dec = degrees + minutes/60, 
          dec is negative if 's' is in the position as given in
          Burnham's and positive if 'n'.
For our example;
    Burnham's gives 18538n4353 as the 1950 coordinates, so

     ra = (18 + 53.8/60) * 15 = 18.896667 *15 = 283.45 degrees
    dec = 43 + 53/60   = 43.883333 degrees

2. Convert to rectangular coordinates

The system of rectangular coordinates used has;

These are the formulas for converting to rectangular coordinates;

    X = cos(ra) * cos(dec)
    Y = sin(ra) * cos(dec)
    Z = sin(dec)

You can derive these formulas using nothing more than simple trigonometry applied to two right angled triangles. An explanation is available.

Applying these formulas to our example we have;

    X = cos(283.45) * cos(43.883333) =  0.1676447
    Y = sin(283.45) * cos(43.883333) = -0.7009849
    Z = sin(43.883333)               =  0.6931922

3. Combined rotation through three angles

To correct for the effects of precession, we need to rotate the coordinate axes. The star has remained in the same place, except for any proper motion. The rotations are applied as follows;

Rotations can be modelled using matrix algebra, each rotation is represented by a 3 by 3 matrix, and the three matrices are multiplied to calculate a final 'precession matrix'. Duffett-Smith has a good explanation of this method. Meeus approaches the whole task using formulas involving the coordinates. The point about precession is that the 'precession matrix' depends only on the difference in time between the two epochs.

Below are the entries for the precession matrix for conversion from B1950.0 to J2000.0;


NOTE: 1999 July 15. Eric Flesch has pointed out that the 
middle column bottom row value should be -0.000027 not
0.000027. This will change subsequent results. Correction
needs to be made to the QBASIC code as well.

              --                                 --
              |  0.999926   -0.011179   -0.004859  |
              |  0.011179    0.999938   -0.000027  |
              |  0.004859    0.000027    0.999988  | 
              --                                 --

You can multiply the column vector representing the B1950 components by this precession matrix to generate the column vector of J2000 components. As neither QBASIC or my TI-80 calculator have matrix operations, below are the formulas for each component written out longhand;

x' = 0.999926*x -0.011179*y -0.004859*z  
y' = 0.011179*x +0.999938*y -0.000027*z 
z' = 0.004859*x +0.000027*y +0.999988*z  

Applying this formula to our example;

The B1950 coordinates are;

x =  0.1676447
y = -0.7009849
z =  0.6931922

Remember to do the multiplications before adding the terms, and that positive
and negative numbers multiply to give negative numbers! 

x' = (0.999926 * 0.1676447) - (0.011179 * -0.7009849) - (0.004859 * 0.6931922)  
y' = (0.011179 * 0.1676447) + (0.999938 * -0.7009849) - (0.000027 * 0.6931922) 
z' = (0.004859 * 0.1676447) + (0.000027 * -0.7009849) + (0.999988 * 0.6931922) 

The result of each product is shown below

x' = 0.1676323 + 0.0078363 - 0.0033682 =  0.1721004 
y' = 0.0018741 - 0.7009414 - 0.0000187 = -0.6990860 
z' = 0.0008146 - 0.0000189 + 0.6931839 =  0.6939796

and we have the J2000 components;

x' =  0.1721004 
y' = -0.6990860 
z' =  0.6939796

The new values are close to the old values, precession is a small effect, so we can be reasonably confident that the calculations are correct. An error with a sign usually makes a large impact in the result!

4. Convert back to RA and DEC referred to J2000

You convert the rectangular components back to angular measure using the following formulas;

    r = arctan(y/x)
    
    if x < 0 then ra2000 = r + 180
    if y < 0 and x > 0 then ra2000 = r + 360
	
    dec2000 = arcsin(z)
Applying these formulas to our example data;
    x' =  0.1721004 
    y' = -0.6990860 
    z' =  0.6939796

    r = atan(-0.6990860/0.1721004)
      = atan(-4.0620824)
      = -76.169982 degrees
  
      y' < 0 and x > 0, so add 360;
 
    ra2000 = 283.83002 degs
           = 18.922001 hrs
	   
   dec2000 = arcsin(0.6939796)
           = 43.94596 degs  
We can write this in various forms, including the notation used in Burnham's and in degrees, minutes and seconds form.

5. Change to desired notation

Notations are important, and software is notorious for giving the wrong notation when you need results quickly! To convert the J2000.0 (or any) RA and DEC to degrees minutes and seconds form, use the following recipe;

    RA = 18.922001 hrs
	 
    subtract the whole number of hours;
	
    18.922001 - 18 = 0.922001 hrs
	
    multiply the fractional part by 60;
	
    0.922001 * 60 = 55.32006
	
    the number of minutes is 55. Subtract this from what remains;
	
    55.32006 - 55 = 0.32006
	
    multiply by 60 again;
	
    0.32006 * 60 = 19.20 seconds
	
    The RA is thus, 18;55;19.2 
	
    The same method on the degrees gives;
	
    43.94596 degs = 43;56;45	

RA is meaured in hours of time, each of which is equal to 15 degrees. For this reason, RA is usually given to slightly more accuracy than DEC.

Many star atlases and observing guides list the positions to 0.1 minute of time for RA and 1 arcminute for DEC. Burnham's uses a very compact notation for this level of accuracy as explained earlier. To convert our example results into Burnham's notation, use the following recipe;

    RA = 18.922001 hrs
    DEC = 43.94596 degs
	
    subtract the hours from the RA, and multiply remainder by 600
    to get the tenths of minutes;
	
    18.922001 - 18 = 0.922001
	
    0.922001 * 600 = 553.2006
	
    round number to the nearest whole unit, and join the hrs
    figure to the minutes figure;
	
    18553
	
    Now subtract the whole degrees from the DEC figure and multiply
    the remainder by 60, and round the number to the nearest whole figure;
	
    43.94596 - 43 = 0.94596
	
    0.94596 * 60 = 56.7576 = 57
	
    join the degrees and arcminute figures together;
	
    4357
	
    Now join both figures together, but with the letter 'n' in between
    for positive DEC, and 's' for negative dec;
	
    18553n4357

It is now easy to see that the J2000.0 position for R-lyr coincides with the position given in the The Cambridge Star Atlas.

QBASIC code for algorithm

[Top]

At the end of this section, you will find a full QBASIC program which will prompt the user for a position in B1950.0 coordinates using the notation in Burnham's Celestial Handbook. The actual calculation is done using the lines below;

'   Corrected precession coefficients.
'   form direction cosines (i.e. rectangular coordinates assuming
'   radius vector of 1 unit)
x = COS(ra) * COS(dec)
y = SIN(ra) * COS(dec)
z = SIN(dec)
'   apply the precession matrix
x2 = .999925708# * x - .0111789372# * y - .0048590035# * z
y2 = .0111789372# * x + .9999375134# * y - .0000271626# * z
z2 = .0048590036# * x - .0000271579# * y + .9999881946# * z
'   convert the new direction cosines to RA, DEC
ra2 = FNatn2(y2, x2)
dec2 = FNasin(z2)

This should look similar to the formulas presented in the last section. Some programmable calculators (and spreadsheets such as Microsoft Excel and Lotus 1-2-3) have built-in matrix multiplication routines, so a more compact code may be used.

The rest of the code is concerned with defining some useful functions, and processing input and output strings. I have used the multiline version of the DEF FN statement to define the functions. By doing this the code will compile using the FirstBasic shareware compiler available from PowerBasic.

The easiest way to transfer the listing into QBASIC would be to;

'*********************************************************************
CLS
'    Make all variables double precision (about 16 significant figs)
DEFDBL A-Z
'    Define some useful constants in radians
pi = 4 * ATN(1)
tpi = 2 * pi
degs = 180 / pi
rads = pi / 180
'	Returns arctan in correct quadrant
DEF FNatn2 (y, x)
    a = ATN(y / x)
    IF x < 0 THEN a = a + pi
    IF (y < 0) AND (x > 0) THEN a = a + tpi
    FNatn2 = a
END DEF
'	Returns arcsine, which QBASIC does not have
DEF FNasin (x)
    c = SQR(1 - x * x)
    FNasin = ATN(x / c)
END DEF
'	returns a string which represents the RA as
'	hhmmm, where the last figure is 1/10 minute of time 
DEF FNra$ (r)
a = INT(r)
b = 600 * (r - a)
r$ = STR$(INT(a * 1000 + b + .5))
r$ = RIGHT$(r$, LEN(r$) - 1)
WHILE LEN(r$) < 5
	r$ = "0" + r$
WEND
FNra$ = r$
END DEF
'	returns a string which represents the DEC as ddmm
DEF FNdec$ (d)
IF d < 0 THEN
	ds$ = "s"
ELSE
	ds$ = "n"
END IF
c = ABS(d)
a = INT(c)
b = 60 * (c - a)
d$ = STR$(INT(a * 100 + b + .5))
d$ = RIGHT$(d$, LEN(d$) - 1)
WHILE LEN(d$) < 4
	d$ = "0" + d$
WEND
FNdec$ = ds$ + d$
END DEF
'   Start of main program
PRINT
PRINT "Precess B1950 to J2000"
PRINT "----------------------"
PRINT
PRINT "type 'q' to quit"
PRINT
10 PRINT   'start of loop
'   Get input string from user
INPUT "B1950 : ", inpos$
'   some simple error checking
IF inpos$ = "q" THEN 20
IF LEN(inpos$) < 10 THEN
	GOTO 10
END IF
'   get the numerical values of RA and DEC
ra = VAL(LEFT$(inpos$, 5)) / 1000
dec = VAL(RIGHT$(inpos$, 4)) / 100
'   get sign of the DEC
SELECT CASE MID$(inpos$, 6, 1)
CASE "s"
  decsgn = -1
CASE "n"
  decsgn = 1
CASE ELSE   'detects wrong input format
	GOTO 10
END SELECT
'   extract decimal hrs/degs values from input format
rahrs = INT(ra)
ramin = (ra - rahrs) * 1000
decdeg = INT(dec)
decmin = (dec - decdeg) * 100
ra = rahrs + ramin / 600
dec = decsgn * (decdeg + decmin / 60)
'   convert to radians
ra = ra * rads * 15
dec = dec * rads
'   form direction cosines (i.e. rectangular coordinates assuming
'   radius vector of 1 unit)
x = COS(ra) * COS(dec)
y = SIN(ra) * COS(dec)
z = SIN(dec)
'   apply the precession matrix
x2 = .999925708# * x - .0111789372# * y - .0048590035# * z
y2 = .0111789372# * x + .9999375134# * y - .0000271626# * z
z2 = .0048590036# * x - .0000271579# * y + .9999881946# * z
'   convert the new direction cosines to RA, DEC
ra2 = FNatn2(y2, x2)
dec2 = FNasin(z2)
'   call functions to generate the Burnham notation and print
PRINT "J2000 : "; FNra$(ra2 * degs / 15) + FNdec$(dec2 * degs)
'   loop round for next input
GOTO 10
PRINT
20 PRINT : PRINT "bye"
END
'********************************************************************

The output from a typical session is shown below, I added the comments later! The position 01487n8902 means 1h 48.7 min RA and +89 degs 2 min DEC, and is the 1950 position of the Pole Star!

B1950 : 01487n8902        'pole star
J2000 : 02318n8916

B1950 : 02470n5541        'eta perseus
J2000 : 02507n5553

B1950 : 04330n1625        'alpha taurus
J2000 : 04359n1631

B1950 : 10057n1212        'alpha leo
J2000 : 10084n1157

B1950 : 15573s2229        'delta scorpius
J2000 : 16003s2238

B1950 : 01359s5730        'alpha eridanus
J2000 : 01377s5715

The 1950 positions are apparent positions, including nutation, taken from Chris Mariott's Skymap v2.1 star chart program. I set the program to the date 1949 December 31st 2209 UT, which corresponds to the instant of B1950. The positions produced by this program coincide with the J2000.0 positions of the stars very closely. There is the odd difference here and there, almost certainly due to the nutation correction.

References

[Top]

Duffett-Smith, Peter
Practical Astronomy with your calculator
Cambridge University Press
3rd edition 1988
ISBN 0-521-35699-7

Straight forward methods, easy to understand, helped me through my University astronomy module! Cheap and recommended if you are working with a programmable calculator or fancy writing some simple Basic programs.

Meeus, Jean
Astronomical Algorithms
Willmann-Bell
1st English edition, 1991
ISBN 0-943396-35-2

Comprehensive cookbook of astronomical calculation. Presents formulas and algorithms - you have to make your own program code, although early chapters do explain some common pitfalls and techniques. Not much in the way of derivation or linking of methods to dynamics. Each copy has a serial number, and you can buy a disc supplement with program code in Pascal, QuickBasic or C.

Seidelmann, P. Kenneth (ed)
Explanatory Supplement to the Astronomical Almanac
University Science Books
1992, completely revised
ISBN 0-935702-68-7

Definitive reference on all aspects of the ephemeris and associated calculations. No hostages taken, no example calculations and modern vectoral notation used throughout. Very thorough treatment of high precision precession calculations in section 3.21.

Burnham, Robert Jr
Burnham's Celestial Handbook: An observer's guide to the universe beyond the solar system
Dover Publications Inc
1977, revised and expanded edition, three volumes
ISBN 0-486-23567-X (Vol 1, paper)
ISBN 0-486-23568-8 (Vol 2, paper)
ISBN 0-486-23673-0 (Vol 3, paper)

The original hitchiker's guide to the Galaxy: detailed descriptions of objects observable in amateur telescopes, and references to mythical stories and historical observations.

Hoffleit E.D., Warren Jr. W.H.
The Bright Star Catalog
Yale University, National Space Science Data Center
1991, 5th Revised Ed.

Machine readable version is available for download from the NASA Astronomical Data Center. Occupies about 2 Mb of disc space in text form (800K download), and provides basic data on the 'naked eye stars'.

Tirion, Wil
The Cambridge Star Atlas
Cambridge University Press
2nd edition 1996
ISBN 0-521-56098-5

Nice clear maps based on the Bright Star Catalog. Covers the sky in 20 maps, each taking one complete page. Opposite page lists details of notable objects including their J2000 coordinates. Cheaper than Norton's Star Atas in the UK, and the maps don't get spread over the page fold.

[Root ]


Last Modified 15th July 1999
Keith Burnett

kburnett@btinternet.com