Posición del Sol en el Sistema de Referencia Local

Un Sistema de Referencia está formado por un punto O denominado origen y una base formada por tres vectores unitarios i, j y k mutuamente perpendiculares. Las direcciones de dichos vectores son los ejes X, Y y Z, respectivamente.

El vector posición r del punto P, se expresa r=xi+yj+zk, o bien, en términos de los ángulos θ y φ

x=r·sinθ·cosφ
y=r·sinθ·sinφ
z=r·cosθ

r=r·sinθ·cosφ i+r·sinθ·sinφ j+r·cosθ k

Sistema de Referencia en el centro de la Tierra

Situamos el observador en el centro de la Tierra que observa como el Sol se mueve en una órbita inclinada ε=23.45° respecto del plano ecuatorial.

En el Sistema de Referencia de la eclíptica la posición del Sol se define con respecto a los ejes X', Y' Z'. La posición angular del Sol está determinada por el ángulo φ que varía de 0° en el equinocio de primavera (en Marzo) hasta completar 360° un año después. El eje X=X' tiene la dirección de la línea de los equinocios

El vector unitario S que apunta desde el centro de la Tierra hacia el Sol se expresa en este SR (Sistema de Referencia)

S=cosφ i'+sinφ j'

Donde i' y j' son vectores unitarios en las direcciones X' y Y' respectivamente.

Sistema de Referencia girado un ángulo ε

Vamos a expresar el vector S en el SR X,Y,Z, donde el plano XY coincide con el plano del Ecuador. El SR Ecuatorial está relacionado con el SR de la Eclíptica por una rotación de ángulo ε alrededor del eje X=X'. Debido a que el eje de la Tierra forma un ángulo de ε=23.45° con la perpendicular al plano de la eclíptica, la relación entre los vectores unitarios puede verse en la figura..

i'=i
j'=cos ε j+sin ε k
k'=-sin ε j+cos ε k

El vector S se expresa en el SR Ecuatorial X,Y,Z del siguiente modo

S=cosφ i+sinφ·cos ε j+sinφ·sin ε k

En el SR Ecuatorial X,Y,Z la posición del Sol se puede escribir utilizando dos ángulos: el ángulo de declinación δ y el ángulo RA, tal como vemos en la figura.

S=cosδ·cosRA i+cosδ·sinRA j+sinδ k

Es también útil considerar un nuevo Sistema de Referencia X" Y" Z" en el cual X" e Y" están en el plano ecuatorial X Y y el eje Z" coincide con Z, el eje de la Tierra. En este SR, el vector S que apunta al Sol siempre está en el plano X" Z"

Sistema de Referencia Local

Consideramos ahora un punto situado sobre la superficie de la Tierra cuya latitud λ y cuyo ángulo horario ω. El Sistema de Referencia Local está formado por un origen situado en el punto y una base definida por tes vectores unitarios E, N, Z mutuamente perpendiculares. Las direcciones de dichos vectores son, la dirección Oeste-Este (tangente al paralelo en el punto), la dirección Sur-Norte (tangente al meridiano en el punto), y la dirección radial.

Expresamos los vectores unitarios E, N, Z, en términos de los vectores unitarios i", j", k".

Como vemos en la figura, E=-sinω i"+cosω j"

Al principio de esta página ya vimos cómo se expresaba el vector r en coordenadas esféricas

Z=cosλ·cosω i"+cosλ·sinω j"+sinλ k"

Calculamos el vector N a partir de Z y E, fijándonos que el producto vectorial Z×E=N, del mismo modo que se cumple que i×j=k

N ^ = Z ^ × E ^ =| i ^ " j ^ " k ^ " cosλcosω cosλsinω sinλ sinω cosω 0 |= cosωsinλ· i ^ "sinωcosλ· j ^ "+cosλ· k ^ "

El vector S está en el plano X" Z", y forma un ángulo δ (declinación) con el eje X"

S=cosδ·i"+sinδ· k"

Altitud y Acimut

Expresamos el vector S en el SR Local, para ello, calculamos las componentes del vector S a lo largo de las direcciones de los vectores unitarios Z, N y E mediante el producto escalar (Recordar que la componente X" del vector S es el producto escalar Sx=S·i")

N ^ =cosωsinλ· i ^ "sinωcosλ· j ^ "+cosλ· k ^ " E ^ =sinω· i ^ "+cosω· j ^ " Z ^ =cosωcosλ· i ^ "+sinωcosλ· j ^ "+sinλ· k ^ " S =cosδ· i ^ "+sinδ· k ^ " S N = S · N ^ =cosωsinλcosδ+cosλsinδ S E = S · E ^ =sinωcosδ S Z = S · Z ^ =cosωcosλcosδ+sinλsinδ S =sinωcosδ· E ^ +(cosλsinδcosωsinλcosδ)· N ^ + (cosωcosλcosδ+sinλsinδ)· Z ^

El vector S se puede expresar como vemos en la figura en términos de dos ángulos αs (altitud del Sol) y γs (acimut del Sol)

S=cosαs·sinγsE+cosαs·cosγsN+sinαsZ

Igualando las componentes, obtenemos tres relaciones:

cosαs·sinγs=-sinω·cosδ
cosαs·cosγs=cosλ·sinδ-sinλ·cosω·cosδ
sinαs=cosλ·cosω·cosδ+sinλ·sinδ

De las cuales, solamente precisamos dos:

  1. Obtenemos αs la altitud del Sol

  2. sinαs=cosλ·cosω·cosδ+sinλ·sinδ

  3. Obtenemos γs el acimut del Sol medido desde el Norte (positivo en el sentido de las agujas del reloj, hacia el Este)

  4. cos γ s = cosλsinδcosωsinλcosδ cos α s

Para calcular γsse precisa calcular el arco coseno que devuelve un ángulo entre 0-180°. Pero el acimut del Sol puede tomar un valor en el intervalo 0-360°. Mediante la siguiente regla obtenemos el ángulo adecuado:

Calculamos γc que está en el intervalo 0≤γc≤180° mediante

γ c =arccos( cosλsinδcosωsinλcosδ cos α s )

Como vemos en la figura:

Altura máxima

A mediodía, el ángulo horario ω=0, el Sol está exactamente en el meridiano local, apuntando hacia el Sur, por lo que el acimut γs=180° y el Sol tiene la máxima altitud αsm.

sinαsm=cosλ·cosδ+sinλ·sinδ

por lo que αsm=90-(λ-δ)=90-λ+δ

Amanecer y Anochecer (Orto y Ocaso)

Cuando la altitud αs=0 se dice que amanece o se pone el Sol.

cosω= sinλsinδ cosλcosδ =tanλ·tanδ Δh= 1 15 arccos( tanλ·tanδ )

El arco coseno devuelve el ángulo en grados y 15 grados equivale a una hora. El ángulo horario ω=0 al mediodía. El ángulo horario al amanecer es negativo ω<0 y al anochecer es positivo ω>0. Luego la hora del Orto hsr=-Δh y la hora del Ocaso es hssh. La duración del día es 2·Δh

Vamos a determinar la hora local del Orto y Ocaso en Bilbao el día 21 de junio de 2013 (N=172 es el número de día). Bilbao tiene una latitud de λ=43.3° y en el día del solsticio de verano δ=23.45°, por lo que Δh=7.61h

La hora del Orto es -7.61h antes y la del Ocaso +7.61h después medidas desde el mediodía. La duración del día es de 15h 13min

En verano la hora de Bilbao es UTC+2, por lo que el mediodía se produce a las 14h hora local, a la que hemos de añadir la corrección debida a la ecuación del tiempo Δt y a la longitud de -2.94° hacia el Oeste (cuatro minutos por grado ó una hora cada 15 grados). Fijarse que amanece y anochece antes en Gerona (extremo Este) que en La Coruña (extremo Oeste), por tanto hemos de sumar los minutos correspondientes a la longitud.

La ecuación del tiempo para N=172 y h=12, es Δt=-1.39 minutos

La corrección es eqTime-4·longitude=-1.39-4·(-2.94)=10.37 minutos

El resultado es que

Las tablas de la la web del Ministerio de Fomento, nos proporcionan los siguientes datos para el día 21 de junio: 6h 32min y 21h 55min, respectivamente.

En invierno la hora de Bilbao es UTC+1 por lo que el mediodía se produce a las 13h local, a la que hemos de añadir la corrección debida a la ecuación del tiempo Δt y a la longitud de -2.94° hacia el Oeste. El 21 de diciembre el ángulo de declinación δ=-23.45° y se denomina solsticio de invierno por lo que Δh=4.39h.

La hora del Orto es -4.39h antes y la del Ocaso +4.39h después medidas desde el mediodía. La duración del día es de 8h 47min

La ecuación del tiempo para N=355 y h=12, es Δt=2.79 minutos

La corrección es eqTime-4·longitude=2.79-4·(-2.94)=14.55 minutos

El resultado es que

Las tablas de la la web del Ministerio de Fomento, nos proporcionan los siguientes datos para el día 21 de diciembre: 8h 41min y 17h 39min, respectivamente.

Creamos un script para dibujar el acimut y la altitud del Sol a cada hora, entre el orto y el ocaso. En el script se realizan las siguientes tareas:

  1. Introduce mediante el comando input

  2. Calcula la declinación δ mediante las fórmula (tómese h=12)

  3. δ=0.006918-0.399912·cos(x)+0.070257·sin(x)-0.006758·cos(2x)+0.000907·sin(2x)-0.002697·cos(3x)
    +0.001480·sin(3x).

    x = 2 π 365 ( N 1 + h 12 24 )

  4. Define la función anónima que calcula la altitud del Sol αs mediante la fórmula

  5. αs=arcsin(cosλ·cosω·cosδ+sinλ·sinδ)

  6. Define la función anónima que calcula el azimut del Sol 0≤γc≤180° mediante la fórmula

  7. γ c = arccos ( cos λ sin δ cos ω sin λ cos δ cos α s )

  8. Calcula el ángulo horario del Orto o del Ocaso tomando el medidodía como cero.

  9. Δh=arccos(-tanλ·tan δ)

    Para convertir el ángulo en horas se tiene en cuenta que 15° equivalen a 1 hora. Se multiplica el ángulo en radianes por el factor de conversión 12/π.

  10. Se crea un vector w de ángulos correspondientes a las horas que van desde el valor entero de -12·Δh/π (Orto) al valor entero de 12·Δh/π (Ocaso)

  11. Para Bilbao que tiene una latitud de λ=43.3°, el día 21 de junio, N=172, la salida del Sol se produce a las -7.61.

    Creamos un vector horas con el operador : cuyos elementos son [-7 -6 -5....0 ...5 6 7]. y convertimos el vector horas en un vector w de ángulos horarios dividiendo por el factor 12/π.

    Se calcula la altitud del Sol αs para cada uno de los valores de dicho vector

    Se calcula el acimut γs teniendo en cuenta que:

    En la gráfica hemos representado el ángulo suplementario de γs, 180-γs, de modo que -180≤γs≤180°

  12. Proporciona los datos de

lambda=input('Latitud: ');
lambda=lambda*pi/180;
N=input('Número de día: ');
h=12;
%declinación
x=2*pi*(N-1-(h-12)/24)/365;
delta=0.006918-0.399912*cos(x)+0.070257*sin(x)-0.006758*cos(2*x)
+0.000907*sin(2*x)-0.002697*cos(3*x)+0.001480*sin(3*x);

%ángulo en radianes del orto u ocaso
hr=acos(-tan(lambda)*tan(delta)); 
%altidud
alfa_s=@(w) asin(cos(lambda)*cos(w)*cos(delta)+sin(lambda)*sin(delta));
%acimut
gamma_c=@(w) acos((cos(lambda)*sin(delta)-cos(w)*sin(lambda)*cos(delta))
./cos(alfa_s(w)) );
hold on
hMax=floor(hr*12/pi); %hora del orto u ocaso, sin decimales

%desde el orto al mediodía
h=-hMax:0;
w=h*pi/12; %convierte hora en radianes
a_s=alfa_s(w)*180/pi;
g_c=180-gamma_c(w)*180/pi;
plot(g_c,a_s,'-ro','markersize',4, 'markerfacecolor','y');

%desde el mediodía al ocaso
h=0:hMax;
w=h*pi/12;
a_s=alfa_s(w)*180/pi;
g_c=180-(2*pi-gamma_c(w))*180/pi;
plot(g_c,a_s,'-ro','markersize',4, 'markerfacecolor','y');
axis([-150 150 0 90]);

%muestra datos 
str=sprintf('duración día: %1.3fh',2*hr*12/pi);
text(-150,87,str)
str=sprintf('altitud máxima: %1.3f°',90-(lambda-delta)*180/pi);
text(-150,82,str)
str=sprintf('declinación: %1.3f°',delta*180/pi);
text(-150,78,str)
grid on
hold off

Corremos el script en la ventana de comandos

>> orto_ocaso
Latitud: 43.3
Número de día: 172

A partir del script orto-ocaso_1 creamos el script orto-ocaso que es una representación tridimensional de la trayectoria del Sol empleando la función plot3. Conocidos la altitud del Sol αs y el acimut γs para cada una de las horas. Calculamos el vector unitario de componentes (x,y,z) que señala la orientación del Sol en la esfera celeste.

S=cosαs·sinγsE+cosαs·cosγsN+sinαsZ

x=cosαs·sinγs
y=cosαs·cosγs
z=sinαs

Para ver un dibujo tridimensional desde distintos puntos de vista, se actúa sobre el botón en la barra del menú titulado Rotate 3D. Se va cambiando el valor de los parámetros del comando view tal como puede apreciarse en la parte inferior izquierda de la ventana.

lambda=input('Latitud: ');
lambda=lambda*pi/180;
N=input('Número de día: ');
h=12;

%declinación
x=2*pi*(N-1-(h-12)/24)/365;
delta=0.006918-0.399912*cos(x)+0.070257*sin(x)-0.006758*cos(2*x)
+0.000907*sin(2*x)-0.002697*cos(3*x)+0.001480*sin(3*x); hr=acos(-tan(lambda)*tan(delta)); %angulo horario del orto u ocaso alfa_s=@(w) asin(cos(lambda)*cos(w)*cos(delta)+sin(lambda)*sin(delta)); gamma_c=@(w) acos((cos(lambda)*sin(delta)-cos(w)*sin(lambda)*cos(delta)) ./cos(alfa_s(w))); hold on hMax=floor(hr*12/pi); %hora del orto u ocaso, sin decimales %desde el orto al mediodía h=-hMax:0; w=h*pi/12; %convierte hora en radianes a_s=alfa_s(w); g_c=gamma_c(w); x=cos(a_s).*sin(g_c); y=cos(a_s).*cos(g_c); z=sin(a_s); plot3(x,y,z,'-ro','markersize',6, 'markerfacecolor','y'); %desde el mediodía al ocaso h=0:hMax; w=h*pi/12; a_s=alfa_s(w); g_c=2*pi-gamma_c(w); x=cos(a_s).*sin(g_c); y=cos(a_s).*cos(g_c); z=sin(a_s); plot3(x,y,z,'-ro','markersize',6, 'markerfacecolor','y'); plot3([0 0],[-0.4 0.6],[0 0],'k') xlabel('x') ylabel('y') zlabel('z') title('Posición del Sol') hold off view(-126,32)

Analema

Si dibujamos la altitud del Sol en función del acimut en un determinado momento y en una localidad dada, obtendremos una figura similar a un ocho.

Como hemos visto, la ecuación del tiempo nos da la diferencia entre el movimiento del Sol (medio) y el actual, en el cual se tiene en cuenta que la órbita de la Tierra es elíptica y el eje de la Tierra está inclinado respecto a la perpendicular a la eclíptica.

A las 12 del mediodía el Sol alcanza su altura máxima αsm=90-λ+δ, cuando su acimut γs es 180º. A lo largo de los 365 días del año, la hora local del mediodía no es constante debido a la corrección que introduce la ecuación del tiempo. (No tenemos en cuenta los adelantos horarios para ahorrar energía, que cambian la hora local en una hora durante medio año).

Si situamos nuestra cámara fotográfica en la misma posición y tomamos una fotografía del Sol a la misma hora local, superponemos las fotografías, veremos la figura de un ocho. Convertimos el tiempo Δt en minutos que nos suministra la ecuación del tiempo en desviación angular en grados teniendo en cuenta que 60 minutos, o una hora son 15 grados.

Δt=229.18·(0.00075+0.001868·cos(x)-0.032077·sin(x)-0.014615·cos(2x)-0.040849·sin(2x))

Donde x se define en función del número de día N y la hora h

x = 2 π 365 ( N 1 + h 12 24 )

El script realiza las siguientes tareas:

  1. Introduce mediante el comando input la latitud λ del lugar
  2. Define una función anónima que calcula el tiempo Δt en minutos, cuando se le pasa el parámetro x relacionado con el número de día N a lo largo del año.
  3. Define una función anónima que calcula la declinación δ, cuando se le pasa el parámetro x relacionado con el número de día N a lo largo del año.
  4. Determina la desviación angular en grados a los que equivale el tiempo Δt
  5. Determima la máxima altitud del Sol al mediodía αsm=90-λ+δ,
  6. Representa mediante el comando plot, la desviación angular en el eje horizontal y la altitud máxima en el vertical

En la figura, se muestra la analema para la latitud λ=43.3° de Bilbao

lambda=input('Latitud: ');
lambda=lambda*pi/180;
h=12;

%ecuación del tiempo en minutos
eq_time=@(x) 229.18*(0.00075+0.001868*cos(x)-0.032077*sin(x)
-0.014615*cos(2*x)-0.040849*sin(2*x));
%declinación
delta=@(x) 0.006918-0.399912*cos(x)+0.070257*sin(x)-0.006758*cos(2*x)
+0.000907*sin(2*x)-0.002697*cos(3*x)+0.001480*sin(3*x);

%convierte minutos a grados
N=1:365;
x=2*pi*(N-1-(h-12)/24)/365;
alfa=90-(lambda-delta(x))*180/pi;
gamma=eq_time(x)*15/60; %60 minutos, 1 hora son 15 grados
plot(gamma,alfa,'-ro','markersize',4, 'markerfacecolor','y');

ylim([0 90])
xlabel('ángulo')
ylabel('Altitud')
title('Analema')

Referencia

Sproul A. B. Derivation of solar geometric relationship using vector analysis. Renewable Energy. 32 (2007) 1187-1205.