Monday, August 17, 2009

Aprendiendo a manejar la información

Quiero ver como perturban la luna y el sol la gravedad en la superficie de la tierra, y que efecto tiene eso sobre los terremotos...
así que debo ver como manejar la información astronómica y geográfica para obtener cuanto afectan esos cuerpos la aceleración de gravedad de donde ocurren los sismos y que diferencia hay con otros puntos.

Cuidado

No conozco mucho sobre como manejar la información astronómica y geográfica...
busqué en Internet los conceptos necesarios para lo que quiero calcular y anoté aquí lo que entendí. Podría ser bastante diferente a lo que se usa normalmente en esos campos y además contener algunos errores.

Valores a calcular

En principio parece bastante sencillo...
Calcularía la aceleración de gravedad debida a ambos cuerpos celestes sobre un punto de la tierra usando la fórmula:

(1) a = G * M / d^2

Donde:
a = Aceleración
G = constante gravitacional
M = masa del cuerpo celeste
d = distancia del punto sobre la tierra al cuerpo celeste

Me interesa ver como varía el valor obtenido sobre la superficie de la tierra (no el valor como tal) y como la variación debe ser pequeña comparada con el valor (el diámetro de la tierra es pequeño comparado con las distancias a ambos cuerpos), debería restarle un cierto valor promedio para poder apreciar mejor esa diferencia. Creo que el valor en el centro de la tierra sería una buena referencia.
Es importante notar que esa aceleración es un vector y debe tenerse en cuenta su dirección, que llamaré û. Va en la dirección del punto en que se mide hacia el cuerpo que la causa, normalmente se toma el vector en dirección contraria pero se usa un signo menos.

Como la fórmula se usará para varios cuerpos (C) en diferentes puntos (P) los identificaré con algunos parámetros entre corchetes:

(2) ā[CP] = - û[CP] G M[C] / d[CP]^2

Donde:
ā[CP] = aceleración debida al cuerpo C en el punto P
û[CP] = vector unitario en la dirección del centro del cuerpo C al punto P
M[C] = masa del cuerpo C
d[CP] = distancia entre el centro del cuerpo C y el punto P

(ā y û son vectores)

Los puntos de interés son los siguientes:
O: Centro de la tierra
T: Punto sobre la superficie de la tierra
L: Luna
S: Sol

Me interesa entonces una aceleración total sobre el punto T, calculada como:
(3) āT = ā[LT] - ā[LO] + ā[ST] - ā[SO]

La dirección de este vector no resultaría muy útil, a menos que sea fácil de ubicar en un mapa.
Así que debería descomponerlo en componentes adecuados, por ejemplo Este (E), Norte (E) y Vertical (V).
Para hacerlo, multiplicaría āT por vectores unitarios en esa 3 direcciones:
(4) ā[E] = āT . û[E]
(5) ā[N] = āT . û[N]
(6) ā[V] = āT . û[V]

Donde:
ā[E] = componente de la aceleración en la dirección Este (signo negativo indicaría Oeste)
ā[N] = componente de la aceleración en la dirección Norte (signo negativo indicaría Sur)
ā[V] = componente vertical de la aceleración (positivo indica saliendo de la tierra, negativo indica entrando a la misma)
û[E] = vector unitario en la dirección Este
û[N] = vector unitario en la dirección Norte
û[V] = vector unitario en dirección vertical (radial, saliendo de la tierra)

El problema pareciera estar resuelto, pero falta algo más... el sistema de coordenadas a usar.

Sistema de referencia geodésico

Los puntos sobre la tierra normalmente se representan en el Sistema de referencia geodésico:
http://es.wikipedia.org/wiki/Sistema_de_referencia_geod%C3%A9sico

Está mejor explicado en la versión en Inglés:
http://en.wikipedia.org/wiki/Geodetic_system

La forma mas común de ubicar los puntos es por longitud, latitud y altitud, representadas por las letras "λ", "φ" y "h" respectivamente.
En inglés también se denominan coordenadas ENU (East North Up)

Es mas fácil para trabajar con vectores usar coordenadas cartesianas.
En inglés llaman a un sistema de coordenadas cartesianas fijo al centro de la tierra ECEF (Earth-centred Earth-fixed).
Su centro es el centro de la tierra (por lo que O sería X=0, Y=0, Z=0), el eje X pasa por el cruce del ecuador con el meridiano de Greenwich, Y pasa por el ecuador a 90° Este, y Z por el polo Norte.

Efemérides

Los cuerpos celestes se ubican respecto al sistema de Efemérides:
http://es.wikipedia.org/wiki/Efemerides#Efem.C3.A9rides_cient.C3.ADfica

Las principales variables son la "ascensión recta" (similar a la longitud) y la "declinación" (similar a la latitud)
La mayor diferencia es que este sistema mantiene una posición fija respecto a las estrellas mientras el geodésico está fijo respecto a la tierra.
La "ascensión recta" suele medirse en horas en vez de grados.

No encontré por ningún lado una explicación directa de como convertir entre ambos sistemas, pero un ejemplo de conversión de "ascensión recta" y "declinación" a "Altura" y "Acimut" me ayudó a entenderlo...
Lo encontré en:
http://www.stargazing.net/kepler/altaz.html

El sistema de efemérides se mide respecto a una época, actualmente se usa el J2000.0,
en que el tiempo se mide en días respecto a las 12:00 GMT del 1 de Enero de 2000.

La "Altura" y el "Acimut" como tales no me interesan...
pero como la altura es el ángulo entre la linea del observador al cuerpo celeste y el horizonte, el punto que de latitud y longitud que me interesa es aquel en que la altura sea de 90º

En el ejemplo calculan la altura por la formula:

sin(ALT) = sin(DEC)*sin(LAT)+cos(DEC)*cos(LAT)*cos(HA)

Donde:
ALT = Altura
DEC = Declinación
LAT = Latitud
HA = Ángulo horario (ingles: Hour angle)

Para que la altura sea de 90º, sin(ALT) debe ser 1, lo que se logra si DEC=LAT y HA=0
Esto me da la latitud... falta la longitud...
El Ángulo horario se calcula como:

HA = LST - RA

Donde:
HA = Ángulo horario (ingles: Hour angle)
LST = Tiempo local sideral (ingles: Local Sideral Time)
RA = Ascensión Recta (inglés: Right ascension)

Para calcular el Tiempo local sideral se usa la siguiente fórmula:
LST = 100.46º + 0.985647 * d + long + 15 * UT

Donde:
LST = Tiempo local sideral (en grados)
d = cantidad de días desde J2000.0
long = longitud (en grados)
UT = Hora Universal (en horas) (inglés: Universal Time)
El resultado se lleva al rango de 0º a 360º
Como me interesa HA=0, LST = RA, de donde:

(7) long = RA - (100.46º + 0.985647 * d + 15 * UT)

Podríamos obtener UT de d, dado que d se mide en días desde las 12 del mediodía, habría que sumarle 0.5, tomar la parte fraccionaria y multiplicarlo por 24,
como el resultado final se llevará al rango de 0º a 360º, el obtener la parte fraccionaria es redundante... bastaría:

long = RA - (100.46º + 0.985647 * d + 15 * 24 * (d + 0.5))
long = RA - (280.46º + 360.985647 * d)

Como el DateTime usado en Windows se mide también en días, basta restarle a la fecha que nos interese el valor de la fecha "01/01/2000 12:00"

Cálculos

Para estos cálculos se requiere la información de los cuerpos celestes (Sol y Luna) y de un punto sobre la Tierra en un momento t.

Datos:
t = fecha, medida en días a partir de la fecha "01/01/2000 12:00"

Los datos de cada cuerpo C se pueden encontrar en Internet la posición de cada cuerpos de interés C en un momento t como efemérides.
Por ejemplo, se puede obtener para fechas recientes en:
http://aa.usno.navy.mil/data/docs/geocentric.php

Datos de C:
RA[C] = Ascensión Recta del cuerpo C en el instante t.
DEC[C] = Declinación de C en ese instante.
R[C] = Distancia de C al centro de la tierra.

Para llevarlo a coordenadas geodésicas usaría:

(8) Long[C] = RA[C] - (280.46º + 360.985647 * t)
(9) LAT[C] = DEC[C]

Donde:
Long[C] = Longitud del punto de la tierra donde C se observa verticalmente arriba (a 90º del horizonte)
LAT[C] = Latitud del mismo punto.

Para el punto P, sobre la tierra, generalmente se puede conseguir su longitud y latitud, que deberían ser llevadas a una forma cartesiana.
Por ejemplo, En:
http://earthquake.usgs.gov/eqcenter/recenteqsww/Quakes/quakes_big.php
se pueden encontrar los epicentros de terremotos recientes.

Datos de P:
Long[P] = Longitud del punto P.
LAT[P] = Latitud del punto P.

Para estos cálculos asumiré que puedo representar la tierra aproximadamente como esférica, por lo que dada la Long y LAT de un punto P obtendría:

(10) X[P] = R * cos(Long[P]) * cos(LAT[P])
(11) Y[P] = R * sin(Long[P]) * cos(LAT[P])
(12) Z[P] = R * sin(LAT[P])

Donde:
R = Radio de la tierra;

Los vectores unitarios en las direcciones que interesan serían:
(13) û[V] = (cos(Long[P]) * cos(LAT[P]), sin(Long[P]) * cos(LAT[P]), sin(LAT[P]))
(14) û[E] = (-sin(Long[P]), cos(Long[P]), 0)
(15) û[N] = (-cos(Long[P]) * sin(LAT[P]), -sin(Long[P]) * sin(LAT[P]), cos(LAT[P]))

Ya tengo las ecuaciones necesarias para pasar de los datos al resultado...

La aceleración āT la obtengo de ecuación (3) donde los diferentes términos se calculan usando la (2).
Los valores requeridos por (2) son datos, o se obtienen de las ecuaciones (9) a (12).
Para calcular sus componentes uso (4), (5) y (6), que requieren los valores de û calculados con (13), (14) y (15)

Resultados

Con todas estas ecuaciones puedo obtener la perturbación en la aceleración de gravedad debida al Sol y la Luna sobre un punto P en que ocurrió un terremoto...
pero el número como tal no indica mucho...
debo poder compararlo con algo, por ejemplo, valores en puntos cercanos o en tiempos diferentes.
Podría intentar obtener algún gráfico de como varía esta aceleración con la posición o el tiempo,
pero la aceleración es un vector (3 variables escalares) que depende de otras 2 coordenadas espaciales (latitud y longitud, ya que la calculo sobre la superficie de la tierra)
y del tiempo.