Por cierto, para ver código fuente en el navegador, con las líneas largas ocupando una sola línea, es muy recomendable desactivar los estilos (suele ser en el menú Ver/Estilo).
EDTinvocar_todo.txt
#source("c:/octave/EDT/EDTinvocar_todo.txt");
source("c:/octave/EDT/EDTfunciones.txt");
source("c:/octave/EDT/EDTcalculo_principal.txt");
source("c:/octave/EDT/EDTgrafica_EDT_EEX_EIE.txt");pause(5);
source("c:/octave/EDT/EDTgrafica_lambda_alfa_gr.txt");pause(5);
source("c:/octave/EDT/EDTgrafica_M_nu_gr.txt");pause(5);
source("c:/octave/EDT/EDTgrafica_EDT_DiaSolar.txt");pause(5);subplot(1,1,1);
source("c:/octave/EDT/EDTgrafica_analema.txt");pause(2);
EDTfunciones.txt
#source("c:/octave/ecuTiempo.txt")
format long
global incliECLI=23.43928108*pi/180 #23°26'21".4119 (sin(incliECLI)=0.397776995) inclinación de la eclíptica
global anyotropico=365.242190402 # en días solares medios
global anyosidereo = 365.256363004 # en días solares medios, en la época J2000
global excentricidadTierra = 0.016729
global factorEDTsegrad = (24*60*60)/(2*pi) # escala vertical dif. angular (rad) <-> tiempo (s)
global factorEDTminrad = (24*60)/(2*pi) # escala vertical dif. angular (rad) <-> tiempo (min)
global omega_T = 2*pi / anyotropico # fecuencia angular del Sol Ecuatorial Medio, sistema equinoccial, rad/dia
global omega_S = 2*pi / anyosidereo # fecuencia angular de la anomalía media M, sistema orbital, rad/dia
global omega_P = omega_S*(1-(anyotropico/anyosidereo)); #frecuencia angular de precesión de los equinoccios, rad/dia
#{
La anomalía Excentrica E se calcula en base a la anomalía Media M,
y a la excentricidad e de la órbita, mediante la ecuación de Kepler
M = E - e sen(E)
#}
function ret = anomaliaExcentrica (anomaliaMedia, excentricidad)
fsol = @(x) (x - excentricidad *sin(x) - anomaliaMedia);
ret = fsolve(fsol,anomaliaMedia);
endfunction
#{
La anomalía verdadera, nu, se relaciona con la anomalia Excentrica E,
y la excentricidad e de la órbita, mediante
tan(nu/2)=sqrt( tan(E/2) * (1+e)/(1-e) )
Aqui se obtiene en función de la anomalía media M.
nu-M es periodica de periodo 2*pi.
Si M0=M-n*2*pi tal que 0<=M0<2*pi entonces
nu = nu0-M0+M = nu0 + n*2*pi
#}
function retv = anomaliaVerdadera (M, exc)
retv = zeros(1,length(M));
for i = [1:length(M)]
n = sign(M(i))*abs(floor(M(i)/(2*pi)));
M0=M(i)-n*2*pi;
E = anomaliaExcentrica(M0, exc);
if (E == pi) ret0 = pi;
elseif (E < pi) ret0 = 2*atan(sqrt((1+exc)/(1-exc))*tan(E/2));
elseif (E > pi) E = 2*pi-E; ret0 = 2*(pi - atan(sqrt((1+exc)/(1-exc))*tan(E/2)));
endif
retv(i) = ret0 + n*2*pi;
endfor
endfunction
#{
Anomalía media correspondiente a una anomalía verdadera dada. Función inversa de anomaliaVerdadera().
#}
function retv = MdeAnomaliaVerdadera (nu, exc)
retv = zeros(1,length(nu));
for i = [1:length(nu)]
n = sign(nu(i))*abs(floor(nu(i)/(2*pi)));
nu0=nu(i)-n*2*pi;
if (nu0 == pi) E = pi;
elseif (nu0 < pi) E = 2*atan(sqrt((1-exc)/(1+exc))*tan(nu0/2));
elseif (nu0 > pi) E = 2*(pi - atan(sqrt((1-exc)/(1+exc))*tan((2*pi-nu0)/2)));
endif
M = E - exc * sin(E);
retv(i) = M + n*2*pi;
endfor
endfunction
#{
Devuelve 'y' (alfa) Ascension Recta de un punto en la eclíptica con Longitud Eclíptica 'x' (lambda)
y una inclinación 'incli' de la ecliptica.
Por tanto x medida desde el punto de Aries
La diferencia x-y es una función periodica de periodo pi, impar respecto a 0, par respecto a pi/4.
Si x = X+n*pi, con 0<=X<=pi
x-y(x) = X-y(X) --> y(x)=y(X)+n*pi , n = sign(x)*abs(floor(x/pi))
#}
function retv = ARdeLE(xv, incli)
retv = zeros(1,length(xv));
for i = [1:length(xv)]
x = xv(i);
n = sign(x)*abs(floor(x/pi));
x0=x-n*pi; # x0 es X
ret = acos(cos(x0)/sqrt((cos(x0))^2+(sin(x0)*cos(incli))^2)) + n*pi;
retv(i)=ret;
endfor
endfunction
#{
DeClinación (delta) de un punto con Longitud ecliptica 'x' (lambda), para inclinación 'incli' de la eclíptica.
sin(delta)=sin(lambda)*sin(incli)
#}
function retv = DCdeLE(xv, incli)
retv = zeros(1,length(xv));
for i = [1:length(xv)]
x = xv(i);
ret = asin(sin(x)*sin(incli));
retv(i)=ret;
endfor
endfunction
EDTcalculo_principal.txt
#invocada desde EDTinvocarTodo.txt
#{
d: tiempo medido en dias solares medios, desde el 1 de enero de 2009 (d=0 a las 0h de ese día)
sufijo _e para valores en el equinoccio de marzo y sufijo _p para valores en el perihelio
nu anomalía verdadera del Sol, medida en sistema orbital, origen en el perihelio
L Longitud eclíptica del Sol, medida en sistema equinoccial, origen en el punto de Aries
AR_Sol ascension recta del Sol, medida en sistema equinoccial, origen en el punto de Aries
M anomalía media, medida en sistema orbtial, origen en perihelio. Unida al discurrir uniforme del tiempo, y a ECLI.
lamb longitud eclíptica de ECLI (igual en valor a la ascensión recta de ECUA), sistema equinoccial, origen en el punto de Aries.
AR_lamb ascensión recta de ECLI, sistema equinoccial, origen en el punto de Aries
#}
d_ini=-10; # incluiremos unos pocos dias del año anterior
d_fin=370; # y otros pocos del siguiente
dias=[d_ini:d_fin];
# En 2009 el equinoccio de primavera se producirá a las 11h 44m UT del 20 de marzo
d_e=31+28+19+((11+(44/60))/24);
# En http://aom.giss.nasa.gov/srorbpar.html obtenemos para la distancia angular desde el perihelio al nodo ascendente, en grados
AOM2009_gr = 283.05
# El angulo medido (en sistema orbital) desde el perihelio hasta el punto de Aries, en radianes,
nu_e = (360 - AOM2009_gr)*pi/180;
# A partir de este angulo determinamos cuando se pasa por el perihelio, d_p.
# Este angulo coincide con nu_e, la anomalia verdadera en el equinoccio de marzo, cuando el Sol pasa por el punto de Aries
# La anomalia media en el equinoccio es
M_e = MdeAnomaliaVerdadera(nu_e, excentricidadTierra);
# y M_e=(d_e-d_p)*omega_S de donde
d_p = d_e - (M_e/omega_S);
# anomalia media
M = omega_S*(dias-d_p);
nu_e = anomaliaVerdadera(M_e, excentricidadTierra);
nu = anomaliaVerdadera(M, excentricidadTierra);
# La diferencia angular en d= d_e entre ECLI y Sol en el S.R. orbital es igual a(M_e - nu_e),
# y coincide con la diferencia angular vista desde el S.R. equinoccial (lamb_e - L_e), siendo L_e=0
# Por tanto lamb_e es negativo, e igual a
lamb_e=M_e-nu_e
lamb_p = lamb_e + omega_T*(d_p - d_e);
# es importante recordar que lamb se refiere a ECLI visto en el S. R. Equinoccial, moviendose al ritmo omega_T.
lamb = lamb_e + omega_T * (dias - d_e);
# Precisamente en el S.R. Ecuatorial la AR (ascensión recta) de ECUA es igual a lamb de ECLI
AR_ECUA = lamb;
# La longitud ecliptica del Sol, L, seria (nu - nu_e) de no haber precesión. Pero como el equinoccio se mueve
# en sentido contrario al Sol, hay un pequeñisimo incremento de L adicional de modo que
L = nu - nu_e + omega_P * (dias - d_e);
L_e = 0;
L_p = 0 - nu_e + omega_P * (d_p - d_e);
#ascensión recta del Sol
AR_Sol = ARdeLE(L, incliECLI);
AR_Sol_e = ARdeLE(L_e, incliECLI);
AR_Sol_p = ARdeLE(L_p, incliECLI);
#ECUACION DEL TIEMPO EDT
valorEDT = (AR_ECUA - AR_Sol) * factorEDTminrad;
# ascensión recta de ECLI
AR_ECLI = ARdeLE(lamb, incliECLI);
#Ecuación de la inclinación de la eclíptica
valorEIE = -(AR_ECLI - AR_ECUA) * factorEDTminrad;
#Ecuación del centro o de la excentricidad de la órbita
valorEEX = -(AR_Sol - AR_ECLI) * factorEDTminrad;
EDTgrafica_EDT_EEX_EIE.txt
#invocada desde EDTinvocar_todo.txt
plot(dias, valorEDT,"-k;EDT;","linewidth",3, dias,valorEIE,"-b;EIE;",dias,valorEEX,"-r;EEX;",d_e,0,"+b",d_p,0,"+r");
xlabel("días desde el inicio del año 2009");
ylabel("minutos");
title("Gráfica 1 - ecuación del tiempo EDT");
grid("minor");
set (gca (), "xlim", [d_ini-5, d_fin+5]);
set(gca(),"position",[0.15,0.15,0.75,0.75])
legend("show","location","northwest")
set(text(36,1,"equinoccio"),"color","blue")
set(text(0,2,"perihelio"),"color","red")
fid = fopen ("c:/Octave/EDT/tabla_EDT_EEX_EIE.txt", "w");
fdisp(fid,"indice,d,fecha,EDT,EEX,EIE");
for i = [1:length(dias)]
fprintf(fid,"%d,%d,%s,%f,%f,%f\n",i,dias(i),datestr(datenum(2009,1,1)+dias(i),29),valorEDT(i),valorEEX(i),valorEIE(i));
endfor
fclose(fid);
EDTgrafica_lambda_alfa_gr.txt
#invocada desde EDTinvocar_todo.txt
lambd=[0:360];
alfa = ARdeLE(lambd*pi/180, incliECLI)*180/pi;
plot(lambd, lambd-alfa,"-k",0,0,"+k");
xlabel('{\fontsize{12} \lambda (grados)}');
ylabel('{\fontsize{12} \lambda - \alpha (grados)}');
title("Gráfica 2");
text(60,2.75,'{\fontsize{12}oblicuidad de la eclíptica = 23,439º}');
text(0,-0.25,'{\fontsize{12}punto de Aries}');
grid("minor");
set (gca (), "xlim", [-5, 365]);
set(gca(),"position",[0.15,0.15,0.75,0.75])
#print -dpng "c:/octave/EDT/grafica_lambda_alfa_gr.png";
fid = fopen ("c:/Octave/EDT/tabla_lambda_alfa_gr.txt", "w");
fdisp(fid,"lambda_gr,alfa_rad,lambda_rad,alfa_gr");
for i = [1:length(lambd)]
fprintf(fid,"%d,%f,%f,%f\n",lambd(i),alfa(i),lambd(i)*pi/180,alfa(i)*pi/180);
endfor
fclose(fid);
EDTgrafica_M_nu_gr.txt
#invocada desde EDTinvocar_todo.txt
M=[0:360];
nu = anomaliaVerdadera (M*pi/180, excentricidadTierra)*180/pi;
plot(M, M-nu,"-k",0,0,"+k");
xlabel('{\fontsize{12} M (grados)}');
ylabel('{\fontsize{12} M - \nu (grados)}');
title("Gráfica 3");
text(50,1.5,'{\fontsize{12}excentricidad = 0,016729}');
text(0,0.2,'{\fontsize{12}perihelio}');
grid("minor");
set (gca (), "xlim", [-5, 365]);
set(gca(),"position",[0.15,0.15,0.75,0.75])
#print -dpng "c:/octave/EDT/grafica_M_nu_gr.png";
fid = fopen ("c:/Octave/EDT/tabla_M_nu_gr.txt", "w");
fdisp(fid,"M_gr,nu_rad,M_rad,nu_gr");
for i = [1:length(M)]
fprintf(fid,"%d,%f,%f,%f\n",M(i),nu(i),M(i)*pi/180,nu(i)*pi/180);
endfor
fclose(fid);
EDTgrafica_EDT_DiaSolar.txt
#invocada desde EDTinvocar_todo.txt
difer = zeros(1,length(dias));
for i = [2:length(difer)]
difer(i)=valorEDT(i)-valorEDT(i-1);
endfor
difer(1)=difer(366);
difer=-difer*60;#diferencia dia a dia en segundos CAMBIADA DE SIGNO
diasEnMes=[31,28,31,30,31,30,31,31,30,31,30,31];
nombreMes=["Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dic","Ene"];
diaIniMes=[0,cumsum(diasEnMes)];
subplot(2,1,1);
plot(dias, valorEDT,diaIniMes,valorEDT(-d_ini+1+diaIniMes),'+');
set (gca (), "xlim", [d_ini-5, d_fin+5]);
ylabel("EDT (minutos)");
#set(text(50,18,"ecuación del tiempo"),"color","black");
for i=[0:12] set(text(diaIniMes(i+1)-4,valorEDT(-d_ini+1+diaIniMes(i+1))+5,nombreMes(3*i+1:3*i+3)),"color","black") ;endfor
grid("on");
subplot(2,1,2);
plot(dias,difer,diaIniMes,difer(-d_ini+1+diaIniMes),'+');
set (gca (), "xlim", [d_ini-5, d_fin+5]);
xlabel("días desde el inicio del año 2009");
ylabel(strcat('\Delta D_v (segundos)'));
set(text(52,26,"exceso del día solar verdadero sobre el día solar medio"),"color","black")
set(text(52,23,"valores positivos: días más largos, negativos: días más cortos"),"color","black")
for i=[0:12] set(text(diaIniMes(i+1)-4,difer(-d_ini+1+diaIniMes(i+1))-5,nombreMes(3*i+1:3*i+3)),"color","black") ;endfor
grid("on");
fid = fopen ("c:/Octave/EDT/tabla_EDT_DiaSolar.txt", "w");
fdisp(fid,"indice,d,fecha,EDT_min,difer_s");
for i = [1:length(dias)]
fprintf(fid,"%d,%d,%s,%f,%f\n",i,dias(i),datestr(datenum(2009,1,1)+dias(i),29),valorEDT(i),difer(i));
endfor
fclose(fid);
EDTgrafica_analema.txt
#invocada desde EDTinvocar_todo.txt
#declinacion del Sol en grados
DC_Sol = DCdeLE(L, incliECLI)*180/pi;
DC_Sol_e = DCdeLE(L_e, incliECLI)*180/pi;
DC_Sol_p = DCdeLE(L_p, incliECLI)*180/pi;
# ángulo horario: diferencia angular entre la ascension recta del Sol y de ECUA, Horizontal, en grados
xH=(lamb-AR_Sol)*180/pi;
xH_e=(lamb_e-AR_Sol_e)*180/pi;
xH_p=(lamb_p-AR_Sol_p)*180/pi;
# plot(dias,DC_Sol);pause(2);#print -dpng "c:/octave/EDT/grafica_dias_DC_gr.png";
diasEnMes=[31,28,31,30,31,30,31,31,30,31,30,31];
nombreMes=["Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dic","Ene"];
diaIniMes=[0,cumsum(diasEnMes)];
#plot(xH,DC_Sol,'-k',"linewidth",1,xH_p,DC_Sol_p,'+r',xH_e,DC_Sol_e,'+b')
latN=0 #latitud Norte del punto de observacion
plot(xH,DC_Sol+latN,'-k',"linewidth",1,xH_p,DC_Sol_p+latN,'+r',xH_e,DC_Sol_e+latN,'+b',xH(-d_ini+1+diaIniMes),DC_Sol(-d_ini+1+diaIniMes),'+g')
offx=[-1,-3,-3,-3,+1,+1,-3,-3,+1,+1,+1,+1];
offy=[-1.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5];
for i=[0:11] set(text(xH(-d_ini+1+diaIniMes(i+1))+offx(i+1),DC_Sol(-d_ini+1+diaIniMes(i+1))+offy(i+1),nombreMes(3*i+1:3*i+3)),"color","black") ;endfor
xlabel("Ángulo entre Sol y meridiano (grados)");
ylabel("Declinación del Sol (grados)");
title("Analema en el Ecuador");
grid("on");
#grid("minor");
set (gca (), "xlim", [-15, 15]);
set (gca (), "ylim", [-30+latN, 30+latN]);
set (gca (), "dataaspectratio", [1, 2, 1]);
set(gca(),"position",[0.15,0.15,0.75,0.75])
set(text(-9,0.5+latN,"equinoccio"),"color","blue")
set(text(-7,-23+latN,"perihelio"),"color","red")
set(text(10,0.5+latN,"OESTE"),"color","black")
set(text(-14,0.5+latN,"ESTE"),"color","black")
set(text(-2,27+latN,"NORTE"),"color","black")
set(text(-1,-27+latN,"SUR"),"color","black")
#print -dpng "c:/octave/EDT/grafica_analema.png";
pause(5);
#-----------------------------------------------------------------------------------
# grafica de los analemas a diferentes horas vistos desde una latitud no ecuatorial
source("c:/octave/EDT/EDTfunciones_cambio_coordenadas.txt")
lat=40; #latitud de la observacion
for i=[0:23]
[z,A]= zAhoriz(DC_Sol*pi/180,(xH+15*i)*pi/180,lat*pi/180);z=z*180/pi;A=A*180/pi;
zz{i+1}=90-z;AA{i+1}=A;#zz distancia al horizonte en vez de al cenit, en grados
endfor
plot(AA{1},zz{1},"k",AA{2},zz{2},AA{3},zz{3},AA{4},zz{4},AA{5},zz{5},AA{6},zz{6},AA{7},zz{7},AA{8},zz{8},AA{9},zz{9},AA{10},zz{10},AA{11},zz{11},AA{12},zz{12},AA{13},zz{13},AA{14},zz{14},AA{15},zz{15},AA{16},zz{16},AA{17},zz{17},AA{18},zz{18},AA{19},zz{19},AA{20},zz{20},AA{21},zz{21},AA{22},zz{22},AA{23},zz{23},AA{24},zz{24});
title("analema cada hora a 40º latitud Norte");
xlabel("Azimut desde el Sur (+ Oeste, - Este)(grados)");
ylabel("Altura sobre el horizonte (grados)");
grid("on");
set (gca (), "xlim", [-180, 180]);
set (gca (), "ylim", [-90, 90]);
#print -dpng "c:/octave/EDT/grafica_analemas40N.png";
pause(5);
#-----------------------------------------------------------------------------------
lat=-40; #latitud de la observacion
for i=[0:23]
[z,A]= zAhoriz(DC_Sol*pi/180,(xH+15*i)*pi/180,lat*pi/180);z=z*180/pi;A=A*180/pi;
zz{i+1}=90-z;AA{i+1}=A;#zz distancia al horizonte en vez de al cenit, en grados
endfor
plot(AA{1},zz{1},"k",AA{2},zz{2},AA{3},zz{3},AA{4},zz{4},AA{5},zz{5},AA{6},zz{6},AA{7},zz{7},AA{8},zz{8},AA{9},zz{9},AA{10},zz{10},AA{11},zz{11},AA{12},zz{12},AA{13},zz{13},AA{14},zz{14},AA{15},zz{15},AA{16},zz{16},AA{17},zz{17},AA{18},zz{18},AA{19},zz{19},AA{20},zz{20},AA{21},zz{21},AA{22},zz{22},AA{23},zz{23},AA{24},zz{24});
title("analema cada hora a 40º latitud Sur");
xlabel("Azimut desde el Norte (+ Oeste, - Este)(grados)");
ylabel("Altura sobre el horizonte (grados)");
grid("on");
set (gca (), "xlim", [-180, 180]);
set (gca (), "ylim", [-90, 90]);
#print -dpng "c:/octave/EDT/grafica_analemas40S.png";
pause(5);
#-----------------------------------------------------------------------------------
lat=0; #latitud de la observacion
for i=[0:23]
[z,A]= zAhoriz(DC_Sol*pi/180,(xH+15*i)*pi/180,lat*pi/180);z=z*180/pi;A=A*180/pi;
zz{i+1}=90-z;AA{i+1}=A;#zz distancia al horizonte en vez de al cenit, en grados
endfor
plot(AA{1},zz{1},"k",AA{2},zz{2},AA{3},zz{3},AA{4},zz{4},AA{5},zz{5},AA{6},zz{6},AA{7},zz{7},AA{8},zz{8},AA{9},zz{9},AA{10},zz{10},AA{11},zz{11},AA{12},zz{12},AA{13},zz{13},AA{14},zz{14},AA{15},zz{15},AA{16},zz{16},AA{17},zz{17},AA{18},zz{18},AA{19},zz{19},AA{20},zz{20},AA{21},zz{21},AA{22},zz{22},AA{23},zz{23},AA{24},zz{24});
title("analema cada hora en el ecuador");
xlabel("Azimut desde el Sur (+ Oeste, - Este)(grados)");
ylabel("Altura sobre el horizonte (grados)");
grid("on");
set (gca (), "xlim", [-180, 180]);
set (gca (), "ylim", [-90, 90]);
#print -dpng "c:/octave/EDT/grafica_analemasEcuador.png";
fid = fopen ("c:/Octave/EDT/tabla_analema.txt", "w");
fdisp(fid,"indice,d,fecha,difAR_gr,DC_gr,AR_Sol_gr,AR_ECUA_gr");
for i = [1:length(dias)]
fprintf(fid,"%d,%d,%s,%f,%f,%f,%f\n",i,dias(i),datestr(datenum(2009,1,1)+dias(i),29),xH(i),DC_Sol(i),AR_Sol(i)*180/pi,AR_ECUA(i)*180/pi);
endfor
fclose(fid);
EDTfunciones_cambio_coordenadas.txt
#{
paso de coordenadas ecuatoriales declinación delta, y ángulo horario t,
a coordenadas horizontales distancia cenital z y Azimut A, para una latitud lat
delta y t pueden ser vectores, pero de la misma dimensión. lat se considera escalar.
#}
function [z, A] = zAhoriz(delta, t, lat)
z = zeros(1,length(delta))';
A = zeros(1,length(delta))';
slat=lat;#valor original con signo
lat=abs(lat);#valor absoluto
for i = [1:length(z)]
if (slat<0) t(i)=t(i)+pi; #para latitudes negativas hay desfase pi en el angulo horario
endif
if (t(i)<0 || t(i)>2*pi)
t(i)=mod(t(i),2*pi); # reducimos SIEMPRE el ángulo horario t al intervalo [0,2*pi]
endif
z(i)=acos(sin(lat)*sin(delta(i))+cos(lat)*cos(delta(i))*cos(t(i)));
if (abs(z(i)) < 0.0001) sA=0; # para z=0 el valor de A es indefinido, elegimos A=0.
else sA=cos(delta(i))*sin(t(i))/sin(z(i)); # sA es sin(A)
cA=(-cos(lat)*sin(delta(i))+sin(lat)*cos(delta(i))*cos(t(i)))/sin(z(i)); #cA es cos(A)
endif
# para obtener A(i) en [0,2*pi]
# cuadrante coseno seno para angulo inverso usar
# 1 + + acos o asin
# 2 - + acos o pi-asin
# 3 - - 2*pi-acos o pi-asin
# 4 + - 2*pi-acos o 2*pi+asin
# En octave asin toma valores en [-pi/2,pi/2], y acos toma valores en [0,pi], para el dominio [-1,1].
if (sA>0) A(i)=acos(cA);
elseif (sA<0) A(i)=2*pi-acos(cA);
else A(i)=0;
endif
#Pero nos interesa A(i) en [-pi,pi]
if (A(i)>pi) A(i)=A(i)-2*pi;
endif
if (slat<0) z(i)=pi-z(i);
# ahora se deshace el desfase pi en latitudes negativas.
# Queremos seguir teniendo las 12h del mediodia a 0º de acimut, que ahora se cuenta desde el Norte, no desde el Sur
# y asi la grafica de analemas para latitudes sur queda bien
if (A(i)<=0) A(i)=A(i)+pi;
else A(i)=A(i)-pi;
endif
endif
endfor
endfunction
