Programa en basic de calculo de efemerides planetarias. (poca aproximacion)

'  EFE.BAS
cls


dim b$(32)
c%=360
pi#=4*atn(1)

cd#=360/365.2422 ' grados recorridos en un dia medio
ra#=180/pi#      ' factor para transformar radianes en grados
rr#=360/pi#      ' el doble de ra#


for f=1 to 8
  read a$(f),t(f),w(f),p(f),e(f),s(f),i(f),no(f)
next f




principio:
'input "Fecha aaaa,mm,dd", y,m,d
y=1984:m=7:d=12
'y=1980:m=1:d=1


' calculo del periodo de tiempo en dias desde el 1 de enero de 1980
if m<3 then k=m else k=0
j=int(365.25*(y-k))+int(30.4375*(m+1+12*k))+d-int(y/100)+_
  int(int((y/100)/4))+1720996.5


ed%=j-2444238.5+1     ' periodo de dias desde el dia con datos iniciales
print ed%
'                   POSICION DE LA TIE
  'ed%=1655
  we#=98.83354  ' longitud del periodo
  ' 102.596403=   longitud del perihelio
  ee#=0.016718  ' excentricidad de la tie


  ' cd#*ed%   =   grados que recorre en los dias del periodo considerado
  '               en su orbita, referido al
  ' 1.00004   =   periodo orbital de la tie en a¤os sidereos
  z#=cd#*ed%/1.00004   :gosub rango:ne#=z#


  'le sumamos la longitud del punto de partida y le restamos la del perihelio.
  ' y tenemos la distancia desde el punto a calcular hasta el perihelio
  z#=ne#+we#-102.596403 :gosub rango:me#=z#


  z#=ne#+rr#*ee#*sin(me#/ra#)+we#:gosub rango:le#=z#
  ve#=le#-102.596403                   ' anomalia verdadera
  re#=(1-ee#^2)/(1+ee#*cos(ve#/ra#))   ' radio vector




print "                     EFEMERIDES "
print
ca$=  "     p          Ascension Recta           Declinacion"
d1$=  "\           \  ### \\  ## \\  ## \\   ### \\  ## \\  ## \\  ## ## ##"
ca1$= "     p             Latitud      Longitud"
d2$=  "\           \    ### \\  ## \\   ## ## ##"
print ca1$
for n=1 to 8
  z#=cd#*ed%/t(n)                  :gosub rango:np#=z#
  z#=np#+w(n)-p(n)                 :gosub rango:mp#=z#
  z#=np#+rr#*e(n)*sin(mp#/ra#)+w(n):gosub posit:lp#=z#
  vp#=lp#-p(n)
  rp#=s(n)*(1-e(n)^2)/(1+e(n)*cos(vp#/ra#))


  pps=sin((lp#-no(n))/ra#)*sin(i(n)/ra#)
  ppc=sqr(1-pps^2)
  psi=ra#*atn(pps/ppc)
  yx=lp#-no(n)
  y=sin(yx/ra#)*cos(i(n)/ra#)
  x#=cos(yx/ra#)
  ct=ra#*atn(y/x#)
  q=ct: gosub cuadrante: ct=q
  l1=ct+no(n)
  r1=rp#*cos(psi/ra#)


                               ' longitud ecliptica
  if n>2 then goto 430         ' p. interiores mer y ven
    ll=le#-l1
    a=ra#*atn((r1*sin(ll/ra#))/(re#-r1*cos(ll/ra#)))
    z#=180+le#+a:gosub posit:lam=z#
    goto 470
    430                       ' p. exteriores jup - plu
    ll=l1-le#
    z#=l1+ra#*atn((re#*sin(ll/ra#)/(r1-re#*cos(ll/ra#)))):gosub posit:lam=z#


    470                       ' latitud ecliptica
    bet=ra#*atn(r1*tan(psi/ra#)*sin((lam-l1)/ra#)/(re#*sin((l1-le#)/ra#)))


'  q=lam : gosub cuadrante: lam=q
  ' formato en signos, grados y minutos
  si=int(lam/30)+1
  grr=lam-(si-1)*30
  gr=int(grr)
  mi=int((grr-gr)*60)
  ' declinacion en formato g.m.s.
  sig=1:if bet<0 then sig=-1
  dc1=abs(bet)
  dh =int(dc1)
  dm=int(60*(dc1-dh))
  dh=sig*dh
' print ca1$


  print using d2$;a$(n);dh;"§";dm;"m";si;gr;mi


  goto fin


  ' coordenadas ecliptico-ecuatoriales
  la=lam/ra#
  ec=23.441884/ra#
  be=bet/ra#
  y=sin(la)*cos(ec)-tan(be)*sin(ec)
  x#=cos(la)
  rh#=ra#*atn((sin(la)*cos(ec)-tan(be)*sin(ec))/(cos(la)))


  ' va al cuadrante
  q=rh# : gosub cuadrante: rh#=q
  ' formato en signos, grados y minutos
  si=int(rh#/30)+1
  grr=rh#-(si-1)*30
  gr=int(grr)
  mi=(grr-gr)*60
  ' ascension recta en formato h.m.s.
  rx=rh#/15
  ry=int(rx)
  mx=60*(rx-ry)
  my=int(mx)
  sx=cint(60*(mx-my))
  ' declinacion en formato g.m.s.
  pps#=sin(be)*cos(ec)+cos(be)*sin(ec)*sin(la)
  ppc#=sqr(1-pps#^2)
  dc=ra#*atn(pps#/ppc#)
  sig=1:if dc<0 then sig=-1
  dc1=abs(dc)
  dh=int(dc1)
  dm=60*(dc1-dh)
  dmi=int(dm)
  ds=cint(60*(dm-dmi))
  dh=sig*dh
'  print ca$
' print using d1$;a$(n);ry;"h";my;"m";sx;"s";dh;"§";dmi;"m";ds;"s";si;gr;mi
fin:
next n


' calculo de longitud y latitud
for n=1 to 8


next


' dibuja los planetas en un atlas estelar
print "t para terminar , s para otra fecha"
a$=input$(1)
cls
if a$="t" then end else goto principio


end








' rango de 0-360 grados
rango:
  z#=360*(z#/360-int(z#/360))
return


' numeros positivos
posit:
  if z#>360 then z#=z#-360
  if z#<0 then z#=z#+360
return


cuadrante:
  if x#<0 then q=q+180:return
  if y>0 then return
  q=q+360
return


sub s
  print "para seguir pulsa una tecla"
  while inkey$="" :wend
end sub


'                 (1)     (2)          (3)        (4)       (5)       (6)       (7)
'                 t       we#           p         ee#        s         i        n
'              periodo long. peri. l.perihelio  excentr. semie ma. inclinac.  nodo asc.
data "Mercu",  0.24085,231.2973 ,   77.1442128,.2056306,0.3870986,7.0043579 ,48.0941733
data "Venus",  0.61521,355.73352,  131.2895792,.0067826,0.7233316,3.394435  ,76.4997524
data "Marte",  1.88089,126.30783,  335.6908166,.0933865,1.5236883,1.8498011 ,49.4032001
data "Jupit", 11.86224,146.966365,  14.0095493,.0484658,5.202561 ,1.3041819 ,100.2520175
data "Satur", 29.45771,165.322242,  92.6653974,.0556155,9.554747 ,2.4893741 ,113.4888341
data "Urano", 84.01247,228.070855, 172.7363288,.0463232,19.21814 ,0.7729895 ,73.8768642
data "Neptu",164.79558,260.3578998, 47.8672148,.0090021,30.10957 , 1.7716017,131.5606494
data "Pluto",250.9    ,209.439    ,222.972    ,.25387  ,39.78459 ,17.137    ,109.941
'     tierra   1.00004  98.83354   102.596403  .016718   1          0          0


' (1) periodo orbital en a¤os sidereos
' (2) longitud el 1 de enero de 1980
' (3) longitud del perihelio
' (4) excentricidad de la orbita
' (5) valor medio del semieje mayor en relacion con el valor de la tierra
' (6) declinacion de la orbita a la ecliptica
' (7) nodo ascendente