'..This program calculates planetary ephemerides.

'..It uses formulae and methods from Peter Duffett-Smith's book:
'.."Practical Astronomy with your Calculator".

'..Single-precision floating point arithmetic is used, so the ephemerides 
'..are very rough (eg: RA and DEC may be out by as much as 30 mins or more).
'..This "roughness" is also partly due to the approximate methods in the 
'..aforementioned book (eg: orbits are assumed to be circular).

'..The calculation code was taken from an AmigaBASIC program I wrote a few
'..years ago. It's terribly monolithic, uses poor variable names, etc.

const latitude=-41.25    '..Launceston's latitude (default)
const RC=57.295779 	 '..Conversion Factor (Deg-Rad)
const OE=23.43  	 '..Obliquity of the Ecliptic
const AU=149597870	 '..Astronomical Unit

'..miscelaneous arrays
DIM Z(2),ZH(2),ZM(2),ZS(2)
DIM X(2),F(2)
DIM IT(2),LT(2),LM(2),LS(2),LT$(2)

'..Planet names
DIM planet$(9)

for i=1 to 9
  read planet$(i)
next

DATA Mercury,Venus,Earth,Mars,Jupiter,Saturn,Uranus,Neptune,Pluto


'...Epoch Data
DIM EP(9,9)

FOR i=1 TO 9  
  FOR J=1 TO 9 
   READ EP(i,J) 
  NEXT J
NEXT i

'..Mercury
DATA .24085,320.66305,77.06645,.205629,.387099
DATA 7.00427,48.03493,6.74,1.918E-6

'..Venus
DATA .61521,310.97453,131.21928,.006785
DATA .723332,3.39438,76.45475,16.92,1.721E-5

'..Earth
DATA 1.00004,99.53431,102.51044,.016720,1,0,0,0,0

'..Mars
DATA 1.88089,249.62919,335.59881,.093382
DATA 1.523691,1.84983,49.36466,9.36,4.539E-6

'..Jupiter
DATA 11.86224,355.21414,13.91992,.048460
DATA 5.202804,1.30450,100.19608,196.74,1.994E-4

'..Saturn
DATA 29.45771,104.17278,92.55833,.055630
DATA 9.538844,2.48933,113.43842,165.6,1.74E-4

'..Uranus
DATA 84.01247,205.78286,170.25472,.04725
DATA 19.181854,.77316,73.87283,65.8,7.768E-5

'..Neptune
DATA 164.79558,249.91462,44.40592,.008586
DATA 30.057960,1.77236,131.50506,62.2,7.597E-5

'..Pluto
DATA 246.378,202.3345,224.258,.246115
DATA 39.29976,17.14451,109.9965,8.2,4.073E-6


'...GST/GMT constant from epoch 1975 to epoch 2000
DIM BV(25)
FOR R=0 TO 25
  READ BV(R)
NEXT R

DATA 17.397221,17.413525,17.363611,17.379644,17.395558,17.411473,17.361678
DATA 17.377592,17.393506,17.409421,17.359625,17.375540,17.391454,17.407368
DATA 17.357573,17.373487,17.389402,17.405316,17.355521,17.371435,17.387349
DATA 17.403264,17.353468,17.369383,17.385297,17.401211


'..trigonometric functions
SUB arcsin(x)
  arcsin = (ATN(X/SQR(-X*X+1)))*RC
END SUB

SUB arccos(x)
  arccos = (-ATN(X/SQR(-X*X+1))+1.5708)*RC
END SUB


{ ** main ** }
string planet_name size 10

'..get the arguments
if (argcount=2 or argcount=3) and ucase$(arg$(2))="TODAY" then
  planet_name = ucase$(arg$(1))
  dt$=date$
  DY = val(mid$(dt$,4,2))
  M = val(mid$(dt$,1,2))
  Y = val(mid$(dt$,7))
  case
    argcount=2 : LAT = latitude
    argcount=3 : LAT = val(arg$(3))
  end case
else
 if argcount=4 then 
  planet_name = ucase$(arg$(1)) 
  DY = val(arg$(2))
  M = val(arg$(3))
  Y = val(arg$(4))
  LAT = latitude
 else
  if argcount=5 then
    planet_name = ucase$(arg$(1)) 
    DY = val(arg$(2))
    M = val(arg$(3))
    Y = val(arg$(4))
    LAT = val(arg$(5))
  else
    print "usage: planet <name> today | <dd> <mm> <yyyy> [<latitude>]"
    print
    print "   eg: planet mars today"
    print "       planet jupiter 2 3 1992"
    print "       planet venus today -35"   
    stop
  end if
 end if
end if

'..known planet?
PL=0
for i=1 to 9
  if planet_name = ucase$(planet$(i)) then PL=i
next

if PL=0 then 
  print "unknown planet!"
  stop
end if

if PL=3 then 
  print "can't give ephemerides for earth!"
  stop
end if

'..legal day?
if DY<1 or DY>31 then
  print "day must be from 1 to 31."
  stop
end if

'..legal month?
if M<1 or M>12 then 
  print "month must be from 1 to 12."
  stop
end if

'..legal year?
if Y<1975 then
  print "year must be 1975 or later."
  stop
end if

'..legal latitude?
if LAT < -70 or LAT > 65 then 
  print "latitude must be from -70 to +65 degrees."
  stop
end if

'...Calculate days since Jan 1 1975
L=0
D=0

IF M<=2 THEN 
  D=E-1 
ELSE 
  D=D+1
  D=D*30.6
  D=INT(D)
END IF

'...Leap year?
Y$=STR$(Y)
YN$=RIGHT$(Y$,2)
YN=VAL(YN$)
IF YN/4=INT(YN/4) THEN L=1 '...YES

IF L=1 AND M<=2 THEN D=D*62 ELSE IF L=0 AND M<=2 THEN D=D*63
IF L=1 AND M>2 THEN D=E-62 ELSE IF L=0 AND M>2 THEN D=E-63

IF M<=2 THEN D=D/2:D=INT(D)

D=D+DY '...add days elapsed in month

'...Days since 1975
TD=0
N=75
days_in_year=0

WHILE N<YN
  IF N/4<>INT(N/4) THEN days_in_year=365 ELSE days_in_year=366
  TD=TD+days_in_year
  N=N+1
WEND

TD=TD+D 	'...Total Days elapsed since epoch 1975
DS=D  		'...days elapsed in current year

'print "** calculated days"

'...Planet Calc's
NP=(360/365.25)*(TD/EP(PL,1))

while NP>360
  NP=NP-360
wend

L1=(360/3.1415927)*EP(PL,4)
L2=NP+EP(PL,2)-EP(PL,3)
L3=SIN(L2/RC)
L4=L1*L3
L=L4+NP+EP(PL,2)

IF L>360 THEN L=L-360
IF L<0 THEN L=L+360

VP=L-EP(PL,3)

R1=EP(PL,5)*(1-EP(PL,4)^2)
R2=EP(PL,4):R3=COS(VP/RC)
R4=1+(R2*R3)
R=R1/R4

'print "** done planet calcs"

'...Earth Calc's
NE=(360/365.25)*(TD/EP(3,1))

while NE>360
  NE=NE-360
wend

LE1=(360/3.1415927)*EP(3,4)
LE2=NE+EP(3,2)-EP(3,3)
LE3=SIN(LE2/RC)
LE4=LE1*LE3
le=LE4+NE+EP(3,2)

IF le>360 THEN le=le-360
IF le<0 THEN le=le+360

VE=le-EP(3,3)

RE1=1-EP(3,4)^2
RE2=EP(3,4)
RE3=COS(VE/RC)
RE4=1+(RE2*RE3)
re=RE1/RE4

'print "** done earth calcs"

'...Ecliptic Coordinate Calc's
PS1=SIN((L-EP(PL,7))/RC)*SIN(EP(PL,6)/RC)

PS=arcsin(PS1)

LL1=TAN((L-EP(PL,7))/RC)*COS(EP(PL,6)/RC)
LL=ATN(LL1)*RC+EP(PL,7)

'...L & LL should be very nearly the same
while LL<L-2
  LL=LL+180
wend

while LL>L+2 
  LL=LL-180
wend

RR=R*COS(PS/RC)

IF PL>3 THEN 
  '..Outer planet calcs
  LD=(LL-le)/RC
  A1=(re*SIN(LD))/(RR-re*COS(LD))
  A=(ATN(A1)*RC)+LL
  IF A>360 THEN A=A-360
  IF A<0 THEN A=A+360
ELSE
  '...Inner planet calcs
  LD=(le-LL)/RC
  A1=ATN((RR*SIN(LD))/(re-RR*COS(LD)))*RC
  A=180+le+A1
  IF A>360 THEN A=A-360
  IF A<0 THEN A=A+360
END IF

'print "** done ecliptic coordinate calcs"

AD=(A-LL)/RC
LRD=(LL-le)/RC
B=ATN((RR*TAN(PS/RC)*SIN(AD))/(re*SIN(LRD)))*RC

'...Equatorial Coordinate Conversion

'...Declination...
DE1=SIN(B/RC)*COS(OE/RC)
DE2=COS(B/RC)*SIN(OE/RC)*SIN(A/RC)
DE=arcsin(DE1+DE2)

'...Right Ascension...
RA1=TAN(A/RC)*COS(OE/RC)
RA2=(TAN(B/RC)*SIN(OE/RC))/COS(A/RC)
RA=ATN(RA1-RA2)*RC

'...Find Quadrants of A & A' & make equivalent...
QUADRANT:

X(1)=INT(A)
X(2)=INT(RA)

FOR i=1 TO 2
  Xdegs=X(i)
  IF Xdegs<=90 THEN F(i)=1
  IF Xdegs>90 AND Xdegs<=180 THEN F(i)=2
  IF Xdegs>180 AND Xdegs<=270 THEN F(i)=3
  IF Xdegs>270 THEN F(i)=4
NEXT i

IF F(2)>F(1) THEN RA=RA-180:X(2)=RA:GOTO QUADRANT
IF F(2)<F(1) THEN RA=RA+180:X(2)=RA:GOTO QUADRANT

RA=RA/15
IF RA>24 THEN RA=RA-24
RA=ABS(RA)

'...Convert Decimal Deg's/Hrs to Deg's/Hrs, Mins, Secs
Z(1)=RA
Z(2)=DE

FOR i=1 TO 2
  Zdegs=Z(i)
  ZH(i)=INT(Zdegs)
  Zmins=(Zdegs-ZH(i))*60
  ZM(i)=INT(Zmins)
  Zsecs=(Zmins-ZM(i))*60
  ZS(i)=INT(Zsecs)
NEXT

'print "** done RA/DEC conversion"

'...Display RA & Dec
title$ = "Ephemeral data for "+planet$(PL)+" on"+str$(DY)+" -"+str$(M)+" -"
title$ = title$+str$(Y)

dashed_line$=""
for i%=1 to len(title$)
  dashed_line$ = dashed_line$ + "-"
next

print
print dashed_line$
print title$
print dashed_line$

PRINT "    Right Ascension:";TAB(22);
PRINT str$(ZH(1));"h";
PRINT str$(ZM(1));"m";
PRINT str$(ZS(1));"s"

PRINT "        Declination:";TAB(22);
PRINT str$(ZH(2));"d";
PRINT str$(ZM(2));"m";
PRINT str$(ZS(2));"s"

'...Distance...
P1=re^2+R^2-2*re*R*COS((L-le)/RC)
P=SQR(P1)
km$=str$(clng((AU*P)/1e6))

PRINT "Distance from Earth:";TAB(22);left$(str$(P),5);" AU";
PRINT " (~";right$(km$,len(km$)-1);" million km)"

'...Angular Size...
TH=EP(PL,8)/P
TH=TH

PRINT "       Angular Size:";TAB(22);left$(str$(TH),5);" arc seconds"

'...Phase...
D=A-L
phase=.5*(1+COS(D/RC))
phase=phase*100
IF phase>100 THEN phase=phase/10 '...adjust for high values of phase

PRINT "              Phase:";TAB(22);left$(str$(phase),6);"%"

'...Apparent Magnitude

'...decrease number of decimal places
M1=R*P
M2=EP(PL,9)*SQR(phase)
M=M1/M2
M=5*LOG(M)/LOG(10)-26.7

M=M+4	'..this is a kludge in order to give a sane answer!

M$=STR$(M)
M$=LEFT$(M$,5)

PRINT " Apparent Magnitude:";TAB(22);left$(str$(M),5)

'...Rise & Set Times
AV=.065709:CV=1.002743:DV=.997257
IF Y>1975 AND Y<=2000 THEN BV!=BV(Y-1975) ELSE BV!=17.4

T1=-TAN(LAT/RC)*TAN(DE/RC)
T2=arccos(T1)
T=T2*(1/15)

'...Find Greenwich Sidereal Time (GST)
GR=24+RA-T
IF GR>24 THEN GR=GR-24
IF GR<0 THEN GR=GR+24
GS=RA+T
IF GS>24 THEN GS=GS-24
IF GS<0 THEN GS=GS+24

'...Convert GST to GMT...
DSS=DS*AV

T0=DSS-BV!
GR=GR-T0
IF GR<0 THEN GR=GR+24
IF GR>24 THEN GR=GR-24
GR=GR*DV
GS=GS-T0
IF GS<0 THEN GS=GS+24
IF GS>24 THEN GS=GS-24
GS=GS*DV

'...Convert to H M S...
LT(1)=GR:LT(2)=GS

FOR i=1 TO 2
  IT(i)=INT(LT(i))
  FL=LT(i)-IT(i)
  L1=FL*60
  LM(i)=INT(L1) '...mins...
  FL=L1-LM(i)
  LS(i)=INT(FL*60) '...secs...
  LT$(i)=STR$(IT(i))+"h"+STR$(LM(i))+"m"+STR$(LS(i))+"s"
NEXT i

GR$=LT$(1)
GS$=LT$(2)

PRINT "               Rise:";TAB(22);GR$
PRINT "                Set:";TAB(22);GS$

END