#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of shell archive."
# Contents:  README j2b.bas newcomb.bas precess.bas
# Wrapped by pcl@black on Fri Feb 14 17:33:03 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f README -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"README\"
else
echo shar: Extracting \"README\" \(2719 characters\)
sed "s/^X//" >README <<'END_OF_README'
XThe three programs in this package use high-precision formulae to calculate
Xprecession between equatorial co-ordinates at different epochs.  Since high
Xprecision is necessary, the programs check for the double-precision library
Xbeing loaded and generate a warning if this is not present.  With GW-BASIC,
Xthe double precision library is loaded by specifying the "/D" switch on the
Xcommand line when the BASIC is started.  All three ask for initial
Xco-ordinates (RA and Dec); perform a few sanity checks on the data and then
Xcalculate the transformed co-ordinates and then stop.  Since full source
Xcode is provided, they may be easily modified to return to the input phase
Xfor further calculation.
X
XPRECESS.BAS uses the IAU 1976 precession constants between two epochs,
Xknown and required, specified as Julian dates.
X
XJ2B.BAS calculates transformations between J2000 and B1950 co-ordinates,
Xas these are the two systems most likely to be met with for the next few
Xyears.  The specified position is transformed in two ways: both as if it
Xwere a B1950 position converted to J2000, and as if it were a J2000
Xposition converted to B1950.  Both results are displayed.
X
XNEWCOMB.BAS uses Newcomb's precession constants, which are used by most
Xolder tabulations.  It is provided to assist in comparing observations
Xmade at disparate epochs in the past.  The epochs, in this case, are
Xspecified as years (e.g. 1991.25).
X
XI wish to acknowledge the assistance of Graeme Waddington in the writing
Xof these routines.  Of course, I take full responsibility for any errors
Xthere may be.
X
X
XAll these programs are copyright (c) 1991 by PC Leyland, but permission
Xis given for the free copying and modification of the programs, subject to
Xthe following restrictions:
X
X1) You may charge a fee for the act of transferring the programs and
X   for any media cost; you may not charge any sort of license fee for
X   the programs themselves.
X
X2) You may not remove the copyright notice at the head of each program.
X
X3) If you do modify the programs, you must insert a REM statement
X   specifying the nature of your modification.
X
X4) If you transfer a copy, or copies, of the programs you may not impose
X   conditions on the recipient more restrictive than the present set.
X
X
XRoughly speaking, you can do what you like so long as you don't make a
Xprofit from my work, do give me credit for doing it, don't blame me for your
Xbugs, and don't prevent other people from enjoying whatever benefit they
Xmay get from the work.
X
XNo guarantees are made that these programs are in any way "correct".  I
Xbelieve they're ok, but I would apreciate notification of, and fixes for,
Xany bugs which may surface.
X
XPaul Leyland			3 June 1991
Xpcl@uk.ac.oxford.convex
END_OF_README
if test 2719 -ne `wc -c <README`; then
    echo shar: \"README\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f j2b.bas -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"j2b.bas\"
else
echo shar: Extracting \"j2b.bas\" \(6057 characters\)
sed "s/^X//" >j2b.bas <<'END_OF_j2b.bas'
X10  REM Precession between B1950 and J2000
X20  REM Copyright (c) 1991 by PC Leyland.
X30 DEFDBL A-H,O-Z
X40 DEFINT I-N
X50 DIM ROT#(3,3), V#(3), J20#(3), B19#(3)
X60 REM Did we fire up in double precision?
X70 REM The DEFDBL took care of simple variables, but we shall also need
X80 REM the sqr and trig functions.  The code below ensures at least 12
X90 REM digits precision, which should be adequate.  (It's a Newton-Raphson
X100 REM calculation of srq(2), by the way.)
X110 X0 = 1#
X120 FOR I = 1 TO 6 : X1 = 1# / X0 + X0 / 2# : X0 = X1 : NEXT I
X130 X1 = ABS (SQR (2#) - X0)
X140 IF X1 < .000000000001# GOTO 200
X150 PRINT CHR$(7); "Warning -- not double precision"; CHR$ (7)
X160 PRINT "Either re-start BASIC with the /d switch"
X170 PRINT "or bear in mind that your output may not be sufficiently accurate."
X180 PRINT
X190 REM A function to convert degrees, minutes and seconds to radians.
X200 DEF FND2R(D,XM,S) = (((S / 60#) + XM) / 60# + D) / 57.29577951308232# 
X210 REM The same for hours, minutes and seconds to radians
X220 DEF FNH2R(H,XM,S) = FND2R(H,XM,S) * 15#
X230 REM The inverse functions.
X240 DEF FNR2D(R) = R * 57.29577951308232#
X250 DEF FNR2H(R) = FNR2D(R) / 15#
X260 PRINT "Enter R.A (three numbers: hours, minutes seconds)"
X270 INPUT H, XM, S
X280 IF (H < 0) OR (H >= 24) THEN PRINT "***Unusual RA hours***" : REM Sanity Clause
X290 IF (XM < 0) OR  (XM >= 60) THEN PRINT "***Unusual RA minutes***"
X300 IF (S < 0) OR  (S >= 60) THEN PRINT "***Unusual RA seconds***"
X310 PRINT "Enter Dec (three numbers: degrees, minutes seconds)"
X320 INPUT DEG, XMIN, SEC
X330 IF (DEG < -90) OR  (DEG > 90) THEN PRINT "***Unusual Dec degrees***"
X340 ISIGN = SGN (DEG)
X350 IF DEG <> 0 GOTO 390
X360 PRINT "Please confirm sign -- type + or - key"
X370 INPUT S$ : IF S$ <> "+" AND S$ <> "-" GOTO 360
X380 IF S$ = "+" THEN ISIGN = 1 ELSE ISIGN = -1
X390 IF (XMIN < 0) OR  (XMIN >= 60) THEN PRINT "***Unusual Dec minutes***"
X400 IF (SEC < 0) OR  (SEC >= 60) THEN PRINT "***Unusual Dec seconds***"
X410 REM If RA hours is less than zero, negate the minutes and seconds also.
X420 REM This shouldn't be needed, but what the hell.
X430 IF H < 0 THEN XM = - XM : S = -S
X440 REM Likewise of Dec degrees, minutes and seconds.
X450 IF ISIGN < 0 THEN XMIN = - XMIN : SEC = -SEC
X460 REM Geronimo!
X470 ALPHA0 = FNH2R(H, XM, S)
X480 DELTA0 = FND2R(DEG, XMIN, SEC)
X490 REM Rotation Matrix.  As given takes B1950 to J2000; Transpose is inverse.
X500 ROT#(1,1) =  .9999256794956877#
X510 ROT#(1,2) = -.0111814832204662#
X520 ROT#(1,3) = -.0048590038153592#
X530 ROT#(2,1) =  .0111814832391717#
X540 ROT#(2,2) =  .9999374848933135#
X550 ROT#(2,3) = -.0000271625947142#
X560 ROT#(3,1) =  .0048590037723143#
X570 ROT#(3,2) = -.000027170293744#
X580 ROT#(3,3) =  .9999881946023742#
X590 REM Convert co-ordinates to direction cosines, then rotate by precession.
X600 V#(1) = COS(DELTA0) * COS(ALPHA0)
X610 V#(2) = COS(DELTA0) * SIN(ALPHA0)
X620 V#(3) = SIN(DELTA0)
X630 FOR I = 1 TO 3
X640 J20#(I) = 0# : B19#(I) = 0# 
X650 FOR J = 1 TO 3
X660 J20#(I) = J20# (I) + ROT#(I,J) * V#(J)
X670 B19#(I) = B19# (I) + ROT#(J,I) * V#(J)
X680 NEXT J
X690 NEXT I
X700 REM J20# and B19# now holds the direction cosines of the precessed co-ordinates.
X710 REM The new declination is ASIN (xxx#(3)); the Right Ascension is
X720 REM ATAN2 (xxx#(2), xxx#(1)).    All would be simple if only this poxy
X730 REM language had the asin() and  atan2() functions of C or Fortran.
X740 ARG1 = J20#(3) : GOSUB 990 : DELTA20 = RESULT
X750 ARG1 = B19#(3) : GOSUB 990 : DELTA19 = RESULT
X760 ARG1 = J20#(2) : ARG2 = J20#(1) : GOSUB 1130 : ALPHA20 = RESULT
X770 ARG1 = B19#(2) : ARG2 = B19#(1) : GOSUB 1130 : ALPHA19 = RESULT
X780 SIGN20$ = " " : SIGN19$ = " "
X790 IF DELTA20 < 0 THEN SIGN20$ = "-" : DELTA20 = -DELTA20
X800 IF DELTA19 < 0 THEN SIGN19$ = "-" : DELTA19 = -DELTA19
X810 DELTA20 = FNR2D (DELTA20) : ALPHA20 = FNR2H (ALPHA20)
X820 DELTA19 = FNR2D (DELTA19) : ALPHA19 = FNR2H (ALPHA19)
X830 IDEG20 = INT (DELTA20)
X840 IM20 = INT (60# * (DELTA20 - IDEG20))
X850 S20 = 3600# * (DELTA20 - IDEG20  - IM20 / 60#)
X860 IH20 = INT (ALPHA20) : IF IH20 = 24 THEN IH20 = 0 : REM Rounding error possible.
X870 IMIN20 = INT (60# * (ALPHA20 - IH20))
X880 SEC20 = 3600# * (ALPHA20 - IH20  - IMIN20 / 60#)
X890 IDEG19 = INT (DELTA19)
X900 IM19 = INT (60# * (DELTA19 - IDEG19))
X910 S19 = 3600# * (DELTA19 - IDEG19  - IM19 / 60#)
X920 IH19 = INT (ALPHA19) : IF IH19 = 24 THEN IH19 = 0
X930 IMIN19 = INT (60# * (ALPHA19 - IH19))
X940 SEC19 = 3600# * (ALPHA19 - IH19  - IMIN19 / 60#)
X950 PRINT   " h   m   s         o  '  ";CHR$(34);"   (2000)     h   m   s         o  '  ";CHR$(34);"   (1950)"
X960 FORM$ = "##  ##  ##.###   &## ## ##.##          ##  ##  ##.###   &## ## ##.##"
X970 PRINT USING FORM$ ; IH20; IMIN20; SEC20; SIGN20$; IDEG20; IM20; S20; IH19; IMIN19; SEC19; SIGN19$; IDEG19; IM19; S19
X980 END
X990 REM Return arcsin (ARG1) as RESULT in range abs (RESULT) <= pi / 2.
X1000 REM If only this brain-damaged language had proper local variables
X1010 REM and multi-line functions.
X1020 ASSIGN = SGN (ARG1) : ARG1 = ABS (ARG1)
X1030 IF ARG1 <= 1# THEN GOTO 1060
X1040 PRINT "Internal error -- argument to asin > 1"
X1050 STOP
X1060 ASTEMP = SQR (1# - ARG1 * ARG1)
X1070 IF ARG1 > .7071067811865475# THEN GOTO 1100 : REM sqrt (0.5)
X1080 ASTEMP = ATN (ARG1 / ASTEMP)
X1090 GOTO 1110
X1100 ASTEMP = 1.570796326794897# - ATN (ASTEMP / ARG1) : REM pi / 2
X1110 RESULT = ASSIGN * ASTEMP
X1120 RETURN
X1130 REM Return atan2 (ARG1, ARG2) as RESULT in range 0 <= RESULT < 2 * pi
X1140 IF (ARG1 + ARG2) <> ARG1 GOTO 1180
X1150 RESULT = 1.570796326794897# : REM ARG2 zero, return pi/2 or 3pi/2
X1160 IF ARG1 < 0# RESULT = 3 * RESULT
X1170 GOTO 1280
X1180 IF ARG2 > 0# THEN GOTO 1240
X1190 IF ARG1 >= 0# THEN GOTO 1220
X1200 RESULT = 3.141592653589793# + ATN (ARG1 / ARG2) : REM Third quadrant
X1210 GOTO 1280
X1220 RESULT = 3.141592653589793# - ATN (-ARG1 / ARG2) : REM Second quadrant
X1230 GOTO 1280 : REM ARG1 >= 0, ARG2 < 0
X1240 IF ARG1 <= 0# GOTO 1270
X1250 RESULT = ATN (ARG1 / ARG2) : REM First quadrant
X1260 GOTO 1280
X1270 RESULT = 6.283185307179587# - ATN (-ARG1 / ARG2) : REM Fourth quadrant
X1280 RETURN
END_OF_j2b.bas
if test 6057 -ne `wc -c <j2b.bas`; then
    echo shar: \"j2b.bas\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f newcomb.bas -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"newcomb.bas\"
else
echo shar: Extracting \"newcomb.bas\" \(6112 characters\)
sed "s/^X//" >newcomb.bas <<'END_OF_newcomb.bas'
X10 REM Precession with Newcomb's constants.
X20 REM Copyright (c) 1991 by PC Leyland.
X30 DEFDBL A-H,O-Z
X40 DEFINT I-N
X50 DIM ROT#(3,3), V1#(3), V2#(3)
X60 REM Did we fire up in double precision?
X70 REM The DEFDBL took care of simple variables, but we shall also need
X80 REM the sqr and trig functions.  The code below ensures at least 12
X90 REM digits precision, which should be adequate.  (It's a Newton-Raphson
X100 REM calculation of sqr(2), by the way.)
X110 X0 = 1#
X120 FOR I = 1 TO 6 : X1 = 1# / X0 + X0 / 2# : X0 = X1 : NEXT I
X130 X1 = ABS (SQR (2#) - X0)
X140 IF X1 < .000000000001# GOTO 200
X150 PRINT CHR$(7); "Warning -- not double precision"; CHR$ (7)
X160 PRINT "Either re-start BASIC with the /d switch"
X170 PRINT "or bear in mind that your output may not be sufficiently accurate."
X180 PRINT
X190 REM A function to convert degrees, minutes and seconds to radians.
X200 DEF FND2R(D,XM,S) = (((S / 60#) + XM) / 60# + D) / 57.29577951308232# 
X210 REM The same for hours, minutes and seconds to radians
X220 DEF FNH2R(H,XM,S) = FND2R(H,XM,S) * 15#
X230 REM The inverse functions.
X240 DEF FNR2D(R) = R * 57.29577951308232#
X250 DEF FNR2H(R) = FNR2D(R) / 15#
X260 PRINT "Enter year of known co-ordinates";
X270 INPUT E0
X280 IF (E0 < 1750#) OR (E0 > 2050#) THEN PRINT "***Unusual epoch***" : REM Sanity clause.
X290 PRINT "Enter required epoch";
X300 INPUT E1
X310 PRINT "Enter known R.A (three numbers: hours, minutes seconds)"
X320 INPUT H, XM, S
X330 IF (H < 0) OR (H >= 24) THEN PRINT "***Unusual RA hours***"
X340 IF (XM < 0) OR  (XM >= 60) THEN PRINT "***Unusual RA minutes***"
X350 IF (S < 0) OR  (S >= 60) THEN PRINT "***Unusual RA seconds***"
X360 PRINT "Enter known Dec (three numbers: degrees, minutes seconds)"
X370 INPUT DEG, XMIN, SEC
X380 IF (DEG < -90) OR  (DEG >= 90) THEN PRINT "***Unusual Dec degrees***"
X390 ISIGN = SGN (DEG)
X400 IF DEG <> 0 GOTO 440
X410 PRINT "Please confirm sign -- type + or - key"
X420 INPUT S$ : IF S$ <> "+" AND S$ <> "-" GOTO 410
X430 IF S$ = "+" THEN ISIGN = 1 ELSE ISIGN = -1
X440 IF (XMIN < 0) OR  (XMIN >= 60) THEN PRINT "***Unusual Dec minutes***"
X450 IF (SEC < 0) OR  (SEC >= 60) THEN PRINT "***Unusual Dec seconds***"
X460 REM If RA hours is less than zero, negate the minutes and seconds also.
X470 REM This shouldn't be needed, but what the hell.
X480 IF H < 0 THEN XM = - XM : S = -S
X490 REM Likewise of Dec degrees, minutes and seconds.
X500 IF ISIGN < 0 THEN XMIN = - XMIN : SEC = -SEC
X510 REM Geronimo!
X520 T1 = (E0 - 1900#) / 100# : T2 = (E1 - 1900#) / 100#
X530 REM Convert to radians.
X540 ALPHA0 = FNH2R(H, XM, S)
X550 DELTA0 = FND2R(DEG, XMIN, SEC)
X560 REM Precession constants are in centuries ...
X570 DT = T2 - T1
X580 C1 = 2304.25# : C2 = 1.396#
X590 C3 = .302# : C4 = .018#
X600 C5 = .791#
X610 C6 = 2004.682# : C7 = -.853#
X620 C8 = -.426# : C9 = -.042#
X630 ZETA = ((C1 + C2 * T1) + (C3 + C4 * DT) * DT ) * DT
X640 Z = ZETA + C5 * DT * DT
X650 THETA = ((C6 + C7 * T1) + (C8 + C9 * DT) * DT) * DT
X660 REM ... and are in arcsec, so convert to radians.
X670 ZETA = FND2R (0#, 0#, ZETA)
X680 Z = FND2R (0#, 0#, Z)
X690 THETA = FND2R (0#, 0#, THETA)
X700 REM Rotation matrix formulation for simplicity.
X710 CS1 = COS (Z) : CS2 = COS (THETA) : CS3 = COS (ZETA)
X720 SN1 = SIN (Z) : SN2 = SIN (THETA) : SN3 = SIN (ZETA)
X730 ROT#(1,1) = -SN1 * SN3 + CS1 * CS2 * CS3
X740 ROT#(1,2) = -SN1 * CS3 - CS1 * CS2 * SN3
X750 ROT#(1,3) = -CS1 * SN2
X760 ROT#(2,1) = +CS1 * SN3 + SN1 * CS2 * CS3
X770 ROT#(2,2) = +CS1 * CS3 - SN1 * CS2 * SN3
X780 ROT#(2,3) = -SN1 * SN2
X790 ROT#(3,1) = +SN2 * CS3
X800 ROT#(3,2) = -SN2 * SN3
X810 ROT#(3,3) = +CS2
X820 REM Convert co-ordinates to direction cosines, then rotate by precession.
X830 V1#(1) = COS(DELTA0) * COS(ALPHA0)
X840 V1#(2) = COS(DELTA0) * SIN(ALPHA0)
X850 V1#(3) = SIN(DELTA0)
X860 FOR I = 1 TO 3
X870 V2#(I) = 0#
X880 FOR J = 1 TO 3
X890 V2#(I) = V2# (I) + ROT#(I,J) * V1#(J)
X900 NEXT J
X910 NEXT I
X920 REM V2# now holds the direction cosines of the precessed co-ordinates.
X930 REM The new declination is ASIN(V2#(3)); the Right Ascension is
X940 REM ATAN2 (V2#(2), V2#(1)).  All would be simple if only this poxy
X950 REM language had the asin() and  atan2() functions of C or Fortran.
X960 ARG1 = V2(3) : GOSUB 1110 : DELTA1 = RESULT
X970 ARG1 = V2(2) : ARG2 = V2(1) : GOSUB 1250 : ALPHA1 = RESULT
X980 SIGN1$ = " "
X990 IF DELTA1 < 0 THEN SIGN1$ = "-" : DELTA1 = -DELTA1
X1000 DELTA1 = FNR2D (DELTA1) : ALPHA1 = FNR2H (ALPHA1)
X1010 IDEG1 = INT (DELTA1)
X1020 IM1 = INT (60# * (DELTA1 - IDEG1))
X1030 S1 = 3600# * (DELTA1 - IDEG1  - IM1 / 60#)
X1040 IH1 = INT (ALPHA1) : IF IH1 = 24 THEN IH1 = 0 : REM Rounding error possible.
X1050 IMIN1 = INT (60# * (ALPHA1 - IH1))
X1060 SEC1 = 3600# * (ALPHA1 - IH1  - IMIN1 / 60#)
X1070 PRINT   " h   m   s         o  '  ";CHR$(34);"
X1080 FORM$ = "##  ##  ##.###   &## ## ##.##"
X1090 PRINT USING FORM$ ; IH1; IMIN1; SEC1; SIGN1$; IDEG1; IM1; S1
X1100 END
X1110 REM Return arcsin (ARG1) as RESULT in range abs (RESULT) <= pi / 2.
X1120 REM If only this brain-damaged language had proper local variables
X1130 REM and multi-line functions.
X1140 ASSIGN = SGN (ARG1) : ARG1 = ABS (ARG1)
X1150 IF ARG1 <= 1# THEN GOTO 1180
X1160 PRINT "Internal error -- argument to asin > 1"
X1170 STOP
X1180 ASTEMP = SQR (1# - ARG1 * ARG1)
X1190 IF ARG1 > .7071067811865475# THEN GOTO 1220 : REM sqrt (0.5)
X1200 ASTEMP = ATN (ARG1 / ASTEMP)
X1210 GOTO 1230
X1220 ASTEMP = 1.570796326794897# - ATN (ASTEMP / ARG1) : REM pi / 2
X1230 RESULT = ASSIGN * ASTEMP
X1240 RETURN
X1250 REM Return atan2 (ARG1, ARG2) as RESULT in range 0 <= RESULT < 2 * pi
X1260 IF (ARG1 + ARG2) <> ARG1 GOTO 1300
X1270 RESULT = 1.570796326794897# : REM ARG2 zero, return pi/2 or 3pi/2
X1280 IF ARG1 < 0# RESULT = 3 * RESULT
X1290 GOTO 1400
X1300 IF ARG2 > 0# THEN GOTO 1360
X1310 IF ARG1 >= 0# THEN GOTO 1340
X1320 RESULT = 3.141592653589793# + ATN (ARG1 / ARG2) : REM Third quadrant
X1330 GOTO 1400
X1340 RESULT = 3.141592653589793# - ATN (-ARG1 / ARG2) : REM Second quadrant
X1350 GOTO 1400 : REM ARG1 >= 0, ARG2 < 0
X1360 IF ARG1 <= 0# GOTO 1390
X1370 RESULT = ATN (ARG1 / ARG2) : REM First quadrant
X1380 GOTO 1400
X1390 RESULT = 6.283185307179587# - ATN (-ARG1 / ARG2) : REM Fourth quadrant
X1400 RETURN
END_OF_newcomb.bas
if test 6112 -ne `wc -c <newcomb.bas`; then
    echo shar: \"newcomb.bas\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f precess.bas -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"precess.bas\"
else
echo shar: Extracting \"precess.bas\" \(6351 characters\)
sed "s/^X//" >precess.bas <<'END_OF_precess.bas'
X10 REM Precession with IAU 1976 constants.
X20 REM Copyright (c) 1991 by PC Leyland.
X30 DEFDBL A-H,O-Z
X40 DEFINT I-N
X50 DIM ROT#(3,3), V1#(3), V2#(3)
X60 REM Did we fire up in double precision?
X70 REM The DEFDBL took care of simple variables, but we shall also need
X80 REM the sqr and trig functions.  The code below ensures at least 12
X90 REM digits precision, which should be adequate.  (It's a Newton-Raphson
X100 REM calculation of sqr(2), by the way.)
X110 X0 = 1#
X120 FOR I = 1 TO 6 : X1 = 1# / X0 + X0 / 2# : X0 = X1 : NEXT I
X130 X1 = ABS (SQR (2#) - X0)
X140 IF X1 < .000000000001# GOTO 200
X150 PRINT CHR$(7); "Warning -- not double precision"; CHR$ (7)
X160 PRINT "Either re-start BASIC with the /d switch"
X170 PRINT "or bear in mind that your output may not be sufficiently accurate."
X180 PRINT
X190 REM A function to convert degrees, minutes and seconds to radians.
X200 DEF FND2R(D,XM,S) = (((S / 60#) + XM) / 60# + D) / 57.29577951308232# 
X210 REM The same for hours, minutes and seconds to radians
X220 DEF FNH2R(H,XM,S) = FND2R(H,XM,S) * 15#
X230 REM The inverse functions.
X240 DEF FNR2D(R) = R * 57.29577951308232#
X250 DEF FNR2H(R) = FNR2D(R) / 15#
X260 PRINT "Enter Julian date of known co-ordinates";
X270 INPUT E0
X280 IF (E0 < 2341972#) OR (E0 > 2469808#) THEN PRINT "***Unusual epoch***" : REM Sanity clause, 1700 to 2050AD.
X290 PRINT "Enter required epoch";
X300 INPUT E1
X310 PRINT "Enter known R.A (three numbers: hours, minutes seconds)"
X320 INPUT H, XM, S
X330 IF (H < 0) OR (H >= 24) THEN PRINT "***Unusual RA hours***"
X340 IF (XM < 0) OR  (XM >= 60) THEN PRINT "***Unusual RA minutes***"
X350 IF (S < 0) OR  (S >= 60) THEN PRINT "***Unusual RA seconds***"
X360 PRINT "Enter known Dec (three numbers: degrees, minutes seconds)"
X370 INPUT DEG, XMIN, SEC
X380 IF (DEG < -90) OR  (DEG >= 90) THEN PRINT "***Unusual Dec degrees***"
X390 ISIGN = SGN (DEG)
X400 IF DEG <> 0 GOTO 440
X410 PRINT "Please confirm sign -- type + or - key"
X420 INPUT S$ : IF S$ <> "+" AND S$ <> "-" GOTO 410
X430 IF S$ = "+" THEN ISIGN = 1 ELSE ISIGN = -1
X440 IF (XMIN < 0) OR  (XMIN >= 60) THEN PRINT "***Unusual Dec minutes***"
X450 IF (SEC < 0) OR  (SEC >= 60) THEN PRINT "***Unusual Dec seconds***"
X460 REM If RA hours is less than zero, negate the minutes and seconds also.
X470 REM This shouldn't be needed, but what the hell.
X480 IF H < 0 THEN XM = - XM : S = -S
X490 REM Likewise of Dec degrees, minutes and seconds.
X500 IF ISIGN < 0 THEN XMIN = - XMIN : SEC = -SEC
X510 REM Geronimo!
X520 T1 = (E0 - 2451545#) / 36525# : T2 = (E1 - 2451545#) / 36525#
X530 REM Convert to radians.
X540 ALPHA0 = FNH2R(H, XM, S)
X550 DELTA0 = FND2R(DEG, XMIN, SEC)
X560 REM Precession constants are in centuries ...
X570 DT = T2 - T1
X580 C1 = 2306.2181# : C2 = 1.39656# : C3 = -.000139#
X590 C4 = .30188# : C5 = -.000345# : C6 = .017998#
X600 C7 = .7928# : C8 = .000411# : C9 = .000205#
X610 C10 = 2004.3109# : C11 = -.8533# : C12 = -.000217#
X620 C13 = -.42665# : C14 = -.000217# : C15 = -.041833#
X630 ZETA = ((C1 + (C2 + C3 * T1) * T1) + ((C4 + C5 * T1) + C6 * DT) * DT ) * DT
X640 Z = ZETA + ((C7 + C8 * T1) + C9 * DT) * DT * DT
X650 THETA = ((C10 + (C11 + C12 * T1) * T1) + ((C13 + C14 * T1) + C15 * DT) * DT) * DT
X660 REM ... and are in arcsec, so convert to radians.
X670 ZETA = FND2R (0#, 0#, ZETA)
X680 Z = FND2R (0#, 0#, Z)
X690 THETA = FND2R (0#, 0#, THETA)
X700 REM Rotation matrix formulation for simplicity.
X710 CS1 = COS (Z) : CS2 = COS (THETA) : CS3 = COS (ZETA)
X720 SN1 = SIN (Z) : SN2 = SIN (THETA) : SN3 = SIN (ZETA)
X730 ROT#(1,1) = -SN1 * SN3 + CS1 * CS2 * CS3
X740 ROT#(1,2) = -SN1 * CS3 - CS1 * CS2 * SN3
X750 ROT#(1,3) = -CS1 * SN2
X760 ROT#(2,1) = +CS1 * SN3 + SN1 * CS2 * CS3
X770 ROT#(2,2) = +CS1 * CS3 - SN1 * CS2 * SN3
X780 ROT#(2,3) = -SN1 * SN2
X790 ROT#(3,1) = +SN2 * CS3
X800 ROT#(3,2) = -SN2 * SN3
X810 ROT#(3,3) = +CS2
X820 REM Convert co-ordinates to direction cosines, then rotate by precession.
X830 V1#(1) = COS(DELTA0) * COS(ALPHA0)
X840 V1#(2) = COS(DELTA0) * SIN(ALPHA0)
X850 V1#(3) = SIN(DELTA0)
X860 FOR I = 1 TO 3
X870 V2#(I) = 0#
X880 FOR J = 1 TO 3
X890 V2#(I) = V2# (I) + ROT#(I,J) * V1#(J)
X900 NEXT J
X910 NEXT I
X920 REM V2# now holds the direction cosines of the precessed co-ordinates.
X930 REM The new declination is ASIN(V2#(3)); the Right Ascension is
X940 REM ATAN2 (V2#(2), V2#(1)).  All would be simple if only this poxy
X950 REM language had the asin() and  atan2() functions of C or Fortran.
X960 ARG1 = V2(3) : GOSUB 1110 : DELTA1 = RESULT
X970 ARG1 = V2(2) : ARG2 = V2(1) : GOSUB 1250 : ALPHA1 = RESULT
X980 SIGN1$ = " "
X990 IF DELTA1 < 0 THEN SIGN1$ = "-" : DELTA1 = -DELTA1
X1000 DELTA1 = FNR2D (DELTA1) : ALPHA1 = FNR2H (ALPHA1)
X1010 IDEG1 = INT (DELTA1)
X1020 IM1 = INT (60# * (DELTA1 - IDEG1))
X1030 S1 = 3600# * (DELTA1 - IDEG1  - IM1 / 60#)
X1040 IH1 = INT (ALPHA1) : IF IH1 = 24 THEN IH1 = 0 : REM Rounding error possible.
X1050 IMIN1 = INT (60# * (ALPHA1 - IH1))
X1060 SEC1 = 3600# * (ALPHA1 - IH1  - IMIN1 / 60#)
X1070 PRINT   " h   m   s         o  '  ";CHR$(34);"
X1080 FORM$ = "##  ##  ##.###   &## ## ##.##"
X1090 PRINT USING FORM$ ; IH1; IMIN1; SEC1; SIGN1$; IDEG1; IM1; S1
X1100 END
X1110 REM Return arcsin (ARG1) as RESULT in range abs (RESULT) <= pi / 2.
X1120 REM If only this brain-damaged language had proper local variables
X1130 REM and multi-line functions.
X1140 ASSIGN = SGN (ARG1) : ARG1 = ABS (ARG1)
X1150 IF ARG1 <= 1# THEN GOTO 1180
X1160 PRINT "Internal error -- argument to asin > 1"
X1170 STOP
X1180 ASTEMP = SQR (1# - ARG1 * ARG1)
X1190 IF ARG1 > .7071067811865475# THEN GOTO 1220 : REM sqrt (0.5)
X1200 ASTEMP = ATN (ARG1 / ASTEMP)
X1210 GOTO 1230
X1220 ASTEMP = 1.570796326794897# - ATN (ASTEMP / ARG1) : REM pi / 2
X1230 RESULT = ASSIGN * ASTEMP
X1240 RETURN
X1250 REM Return atan2 (ARG1, ARG2) as RESULT in range 0 <= RESULT < 2 * pi
X1260 IF (ARG1 + ARG2) <> ARG1 GOTO 1300
X1270 RESULT = 1.570796326794897# : REM ARG2 zero, return pi/2 or 3pi/2
X1280 IF ARG1 < 0# RESULT = 3 * RESULT
X1290 GOTO 1400
X1300 IF ARG2 > 0# THEN GOTO 1360
X1310 IF ARG1 >= 0# THEN GOTO 1340
X1320 RESULT = 3.141592653589793# + ATN (ARG1 / ARG2) : REM Third quadrant
X1330 GOTO 1400
X1340 RESULT = 3.141592653589793# - ATN (-ARG1 / ARG2) : REM Second quadrant
X1350 GOTO 1400 : REM ARG1 >= 0, ARG2 < 0
X1360 IF ARG1 <= 0# GOTO 1390
X1370 RESULT = ATN (ARG1 / ARG2) : REM First quadrant
X1380 GOTO 1400
X1390 RESULT = 6.283185307179587# - ATN (-ARG1 / ARG2) : REM Fourth quadrant
X1400 RETURN
END_OF_precess.bas
if test 6351 -ne `wc -c <precess.bas`; then
    echo shar: \"precess.bas\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of shell archive.
exit 0
