2011/09/23

Calculo aproximado de efemerides planetarias. Programa en Basic



'  EFEMERIDES PLANETARIAS
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 multiplicando
                 '                         grados en radianes dividiendo
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

' 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 TIERRA
  'ed%=1655
  we#=98.83354  ' longitud del periodo
  ' 102.596403=   longitud del perihelio
  ee#=0.016718  ' excentricidad de la tierra

  ' cd#*ed%   =   grados que recorre en los dias del periodo considerado
  '               en su orbita, referido al
  ' 1.00004   =   periodo orbital de la tierra 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 PLANETARIAS"
print
print "  Planeta       Ascension Recta           Declinacion"
d1$=  "\           \  ### \\  ## \\  ## \\   ### \\  ## \\  ## \\  ## ## ##"

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#)

  if n>2 then goto 430         ' Planetas interiores Mercurio y venus
    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 ' planetas exteriores   jupiter - pluton
  ll=l1-le#
  ' longitud ecliptica
  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#)))

  ' 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 en horas, m. y s.
  rx=rh#/15
  ry=int(rx)
  mx=60*(rx-ry)
  my=int(mx)
  sx=cint(60*(mx-my))
  ' declinacion en formato en grados, m. y 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 using d1$;a$(n);ry;"h";my;"m";sx;"s";dh;"§";dmi;"m";ds;"s";si;gr;mi
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

No hay comentarios: