cls
'calculamos el valor de pi y el coeficiente ra que pasa grados
'a radianes
pi=4*atn(1):ra=2*pi/360
'damos el valor de la primera fecha de la serie que
'queremos calcular. a es el año, m es el mes y d es el dia.
'Si d es un numero fraccionario indica la fraccion decimal de ese dia
'a partir de las 00.00 h, es decir la media noche. Un dia con fraccion
'de 0.5 dara dias julianos enteros pues estos se miden desde las
'12.00 del medio dia.
print "dj = Dia juliano "
print "t = Tiempo T en centurias julianas desde J2000"
print "Lms = Longitud geometrica media del sol "
print "Ams = Anomalia media del sol"
print "Ecs = Ecuacion del centro del sol "
print "Lvs = Longitud verdadera del sol"
print "Las = Longitud aparente del sol"
print "eot = excentricidad de la orbita de la tierra"
print ""
print " dj t Lms Ams Ecs Lvs Las eot Avs rvs oe arc ds Las Las"
print "____________________________________________________________________"
'calculamos el dia juliano como se explico en otro post.
'a=2000:m=1:d=1.5
a=1992:m=10:d=13
an=a:me=m
if m<3 then
m=m+12
a=a-1
end if
A=int(a/100)
B=2-A+int(A/4)
dj = int(365.25 * (a + 4716)) + int(30.6001 * (m+1)) + d + B - 1524.5
'si queremos hacer la serie a partir de un dia juliano determinado
'hay que poner ese dato aqui.
'dj=2448976.5
'dj=2448908.5
'................................
'Hacemos un bucle a partir del dia inicial, en este caso para 20 dias mas.
for n=1 to 20
'damos los valores sucesivos
dj=dj+n-1
print dj;" ";
'..............
' t es el numero de periodos de 36525 dias desde la epoca J2000,
' es decir desde el 1.5 de enero de ese año.
t=(dj-2451545)/36525
print using ("###.#######",t);" ";
'..............
'despues calculamos la longitud media del sol con referencia
'al equinocio de cada fecha
Lms=280.46646+36000.76983*t+0.0003032*t^2
'print using ("###.#####",Lms);" ";
call cuadrante Lms
print using ("###.#####",Lms);" ";
'..............
'a continuaion se calcula la anomalia media del sol, que
'es igual a la anomalia media de la tierra.
Ams=357.52911+35999.05029*t-0.0001537*t^ 2
call cuadrante Ams
print using ("###.#####",Ams);" ";
Ams2=2*Ams
call cuadrante (Ams2)
Ams3=3*Ams
call cuadrante (Ams3)
'...................
'el siguiente paso es calcular la ecuacion del centro del sol
Ecs=(1.914602-0.004817*t-0.000014*t^2)*sin(Ams*ra)+_
(0.019993-0.000101*t)*sin(Ams2*ra)+_
(0.000289)*sin(Ams3*ra)
call cuadrante Ecs
print using ("###.#####",Ecs);" ";
'..................
'Sumando a la ecuacion anterior la longitud media calculada
'antes tenemos la longitud verdadera del sol
Lvs=Lms+Ecs
call cuadrante Lvs
print using ("###.#####",Lvs);" ";
'...................
'hacemos una correccion por la nutacion y la aberracion:
u=125.04-1934.136*t
call cuadrante u
'print "Correccion u = ";u
'y ya tenemos la longitud aparente del sol, referida al
'verdadero equinocio de la fecha considerada:
Las=Lvs-0.00569-0.00478*sin(u*ra)
call cuadrante Las
print using ("###.#####",Las);" ";
lo=int(Las)
si=int(lo/30)+1
gr=lo-(si-1)*30
if gr=30 then si=si+1:gr=0
if si>12 then si=si-12
mi=(Las-lo)*60
m=int(mi)
s=int((mi-m)*60)
lo$=str$(lo)+"g :"+str$(m)+"m :"+str$(s)+"s"
lo1$=str$(gr)+"("+str$(si)+")"+str$(int(mi))
'---------------
'A continuacion calculamos el radio vector, es decir la distancia
'entre el centro del sol y el centro de la tierra.
'(en unidades astronomicas)
'Para este calculo tenemos primero que conocer la excentricidad de la
' tierra y la anomalia verdadera del sol:
'La excentricidad de la orbita de la tierra es:
eot=0.016708634 - 0.000042037*t - 0.0000001267*t^2
print using ("##.########",eot);" ";
'................
'para saber la anomalia verdadera del sol basta con sumar
'a la anomalia media su ecuacion del centro:
Avs=Ecs+Ams
call cuadrante Avs
print using ("###.#####",Ams);" ";
'-------------------
' y el radio vector del sol sera:
rvs=1.000001018*(1-eot^2)/(1+eot*cos(Avs*ra))
'print 1.000001018*(1-eot^2)
print using ("##.########",rvs);" ";
'............
'Por ultimo, para calcular la ascension reta y la declinacion, debemos
'primero calcular la oblicuidad de la ecliptica:
oe=(23+26/60+21.448/3600)-46.8150/3600*t-0.00059/3600*t^2+0.001813/3600*t^3
print using ("##.########",oe);" ";
'..............
'y por tanto la ascension recta del sol sera:
oeu=oe+0.00256*cos(u*ra)
call cuadrante oeu
arc=atn(cos(ra*oeu)*sin(Las*ra)/cos(Las*ra))/ra+180
call cuadrante arc
print using ("###.#####",arc);" ";
'................
'y la declinacion
ds=asn(sin(ra*oeu)*sin(Las*ra))/ra
print using ("###.#####",ds);" ";
'..............
print using ("###",lo);"º ";using ("##",m);"m ";using ("##",s);"s ";" = ";lo1$
next n
end
'-----------------
sub s
while inkey$="":wend
end sub
sub cuadrante byref k
if k<0 then
do: k=k+360: loop until k>0
end if
if k>360 then
do: k=k-360: loop until k<360
end if
end sub
'Antonio Martinez
No hay comentarios:
Publicar un comentario