Coordenadas geocentricas del sol. Programa en Basic

'coordenadas del sol  (basado en las formulas de Jean Meeus "Astronomical algorithms")
screen 9
cls
pi=4*atn(1)

ra=2*pi/360

'dj=2448543.5 ' dia juliano de ejemplo
dj=2448843.5

t=(dj-2451545)/36525

'Longitud geometrica media del sol
Lms=280.46646+36000.76983*t+0.0003032*t^2
print "Lms = ";Lms
call cuadrante (Lms)
print "Lms = ";Lms

'Anomalia media del sol
Ams=357.52911+35999.05029*t-0.0001537*t^2
call cuadrante (Ams)
call cuadrante (2*Ams)
call cuadrante (3*Ams)
'Ecuacion del centro del sol
Ecs=(1.914602-0.004817*t-0.000014*t^2)*sin(Ams*ra)+_
     (0.019993-0.000101*t)*sin(2*Ams*ra)+_
     (0.000289)*sin(3*Ams*ra)
call cuadrante(Ecs)
'Longitud verdadera del sol
Lvs=Lms+Ecs
call cuadrante(Lvs)
'correccion
u=125.04-1934.136*t
'Longitd aparente del sol
Las=Lvs-0.00569-0.00478*sin(u*ra)
call cuadrante(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$(cint(mi))
print "lo.sol: ";Las#;lo$;"  =  ";
print lo1$


end


sub s
while inkey$="":wend
end sub

sub cuadrante(k#)
  do
  if k#<0   then k#=k#+360 else exit loop
  loop
  do
  if k#>360 then k#=k#-360 else exit loop
  loop
end sub