'*****************************************************
'*              Position of planets                  *
'* ------------------------------------------------- *
'* This program calculates the ephemeris of a planet *
'* in our solar system to adjust an equatorial tele- *
'* scope.                                            *
'* ------------------------------------------------- *
'* Ref.: "Mathematiques par l'informatique indivi-   *
'*        duelle - Programmes en BASIC, MASSON,      *
'*        Paris, 1982, p 100 - 105" [BIBLI 06].      *
'* ------------------------------------------------- *
'* Inputs:                                           *
'*   Date: day,month,year                            *
'*   Hour UT: hour                                   *
'*   Planet number: 1 to 8                           *
'*   (Mercury: 1 Venus : 2 Mars   : 4 Jupiter: 5     *
'*    Saturn : 6 Uranus: 7 Neptune: 8)               *
'* Outputs:                                          *
'*   Ascent in hours (0 to 24)                       *
'*   Declination in deg. (-90 to 90 North or South)  *
'*                                                   *
'* SAMPLE RUN:                                       *
'* (Find ascent and declination of planet Mars on    *
'*  March 10th, 1982 at 6h UT)                       *
'*                                                   *
'* Date (D,M,Y): 10,3,1982                           *
'* Hour UT: 6                                        *
'* Mercury: 1 Venus : 2 Mars   : 4 Jupiter: 5        *
'* Saturn : 6 Uranus: 7 Neptune: 8                   *
'* Planet number: 4                                  *
'*                                                   *
'* Ascent     : 13 H  8 MN                           *
'* Declination:  3 ø 45 MN S                         *
'*                                                   *
'* This program allows building the following table: *
'*                                                   *
'* Date: March 10th, 1982 at 6H UT.                  *
'*                                                   *
'*       Planet     Ascent     Declination           *
'*       Mercury    21 H 51    14 ø 45 mn S          *
'*       Venus      20 H 26    14 ø 57 mn S          *
'*       Mars       13 H 08     3 ø 45 mn S          *
'*       Jupiter    14 H 32    13 ø 30 mn S          *
'*       Saturn     13 H 22     5 ø 42 mn S          *
'*                                                   *
'*             English Version By J-P Moreau, Paris. *
'*                        (www.jpmoreau.fr)          *
'*****************************************************
defdbl a-h,o-z
defint i-n
'
dim A(9,8),B(12)
'Init planetary constants (fill table A by columns)
for j=1 to 8
  for i=1 to 9
    read A(i,j)
  next i
next j
data 4.01166,.071425454,1.32493,.000000742289,.823045
data .000000566185,.205615,.122225,.387099
data 3.60861,.027963119,2.271616,.00000065572,1.32291
data .000000436681,.006816,.0592301,.723332
data 1.72727,.0172028,1.76688,.000000818559,0
data 0,.016751,0,1
data 2.17756,.0091467658,5.83378,.000000879297,.851616
data .000000371232,.093309,.0322939,1.5236
data 4.68279,.00145099,.2289,.000000857,1.73578
data .000000482933,.048376,.0228418,5.202799
data 4.8567,.00058484,1.5974,.000000412,1.96856
data .000000417308,.054311,.0435026,9.552098
data 4.3224,.000205424,2.9523,.000000762,1.2825
data .000000238237,.047319,.013482,19.21694
data 1.5223,.000105061,.7637,.000000393,2.28102
data .00000052517,.008262,.0310536,30.12912
'calendar constants
for i=1 to 12
  read B(i)
next i
data 0,31,59,25,90,25,120,25,151,25,181,25,212
data 25,243,25,273,25,304,25,334,25
cls
print
'input hour, date (DD/MM/YY) and planet number
input " Date (D,M,Y): ", j,m,aa
input " Hour UT: ", h
print " Mercury: 1 Venus : 2 Mars   : 4 Jupiter: 5"
print " Saturn : 6 Uranus: 7 Neptune: 8"
input " Planet number: ", n
'calculate time t
t=365.25*(aa-1901)+B(m)+j
t=INT(t)+h/24
'calculate earth coordinates
j=3 : gosub 1000 'planet #3 coordinates
g=x : h=y 'save earth coordinates
'calculate coordinates of planet #n
j=n : gosub 1000
'calculate geocentric equatorial coordinates
x=x-g : y=y-h
t=y*.917484#-z*.397772#
z=y*.397772#+z*.917484# : y=t
'calculate ascent and declination
PI=4#*ATN(1#)
xx=x:yy=y:gosub 2000  'aa=ATN2(y,x)
aa=zz
yy=z:xx=SQR(x*x+y*y):gosub 2000 'd=ATN2(z,SQR(x*x+y*y))
d=zz
'conversion
aa=aa*12#/PI
if aa<0 then aa=24#+aa
h=INT(aa) : m=INT(60*(aa-h))
d=d*180#/PI
A$="N" : if d<0# then A$="S"
d=ABS(d) : b=INT(d) : c=INT(60#*(d-b))
'print results
print
print " Ascent     : ";
print using "##"; h; : print " H ";
print using "##"; m; : print " MN"
print " Declination: ";
print using "##"; b; : print " ø ";
print using "##"; c; : print " MN "; A$
END

1000 'planet coordinates
'calculate planetary constants
  XL=A(1,j)+A(2,j)*t
  O=A(3,j)+A(4,j)*t
  P=A(5,j)+A(6,j)*t
  E=A(7,j)
  XI=A(8,j)
  XA=A(9,j)
'solve Kepler's equation
  xm=XL-O : u=xm
  for j=1 to 10
    u=xm+E*sin(u)
  next j
'formula (3) of reference book
 r=XA*(cos(u)-E)
 v=XA*SQR(1-E*E)*sin(u)
 O=O-P : xm=sin(O) : O=cos(O)
 q=sin(P) : P=cos(P)
 xj=sin(XI) : XI=cos(XI)
 x=(P*O-XI*q*xm)*r+(-P*xm-XI*q*O)*v
 y=(q*O+XI*P*xm)*r+(-q*xm+XI*P*O)*v
 z=xj*xm*r+xj*O*v
return

2000 'calculate zz=ATN2(yy,xx)
if ABS(xx) < 1e-12 then
  if yy>0 then
    zz=PI/2#
  elseif ABS(yy) < 1e-12 then
    zz=-PI/2#
  else
    zz=0#
  end if
else
  zz=ATN(yy/xx)
  if xx<0# then zz=zz+PI
end if
return

'end of file planets.bas