Cálculo de Coordenadas Geodésicas, Distancias y Tiempos Astronómicos en MATLAB

Cálculo de Coordenadas Geodésicas, Distancias y Tiempos Astronómicos

El siguiente código en MATLAB implementa un conjunto de cálculos de geodesia y astronomía de posición: distancias entre estaciones, transformaciones de coordenadas entre marcos de referencia y cálculos de tiempos astronómicos.

1. Longitud de arco paralelo entre estaciones

Se convierte la latitud de la primera estación a radianes y se calcula el radio de curvatura en el primer vertical N1. La diferencia de longitudes se obtiene como lonT = lon1 - lon2.

lat1 = deg2rad(dms2degrees(latY));
lon1 y2 ... lonT = lon1 - lon2;
N1 = a / sqrt(1 - e2*sin(lat1)^2);

Distancia del paralelo: pdist1 = N1*cos(lat1).*lonT;

DifLongP = abs(pdist1 - pdist2);

2. Longitud de arco de meridiano entre estaciones

Se usa la fórmula de longitud de arco de meridiano Sm para cada estación y se calcula la diferencia: DifLongM = abs(Sm1 - Sm2);

3. Distancia espacial entre estaciones

Se obtienen las coordenadas cartesianas geocéntricas de cada estación:

x1 = (N1 + h1)*cos(lat1)*cos(lon1);
y1 = (N1 + h1)*cos(lat1)*sin(lon1);
z1 = (N1*(1 - e2) + h1)*sin(lat1);
p1 = [x1 y1 z1];
Dg = norm(p2 - p1);

4. Coordenadas locales ENU

Las componentes Este (e), Norte (n) y Up (u) se obtienen mediante la matriz de rotación local:

incx = x2 - x1;
e = -sin(lon1)*incx + cos(lon1)*incy;
n = sin(lat1)*cos(lon1)*incx - sin(lat1)*sin(lon1)*incy + cos(lat1)*incz;
u = cos(lat1)*cos(lon1)*incx + cos(lat1)*sin(lon1)*incy + sin(lat1)*incz;

5. Transformación de marco de referencia ITRF93 → ITRF2020

Se aplica una transformación de Helmert de 7 parámetros con variación temporal, desde la época de 2005.35 hasta 2025.49:

t0 = 2015.0;
t1 = 2025.49;
t2 = 2005.35;
it1 = t2 - t1;     % Diferencia de años para mover la estación
x2005 = x1 + VxY*it1;
p1 = [x2005; y2005; z2005];
it0 = t2 - t0;

% Traslaciones (m)
Tx = (-65.8/1000) + (-2.8/1000)*it0;
% Ty, Tz

% Factor de escala (ppb)
D = (4.47/10^9) + (0.12/10^9)*it0;

% Rotaciones (rad)
Rx = ((-3.36/1000/3600)*(pi/180)) + ((-0.11/1000/3600)*(pi/180))*it0;
% Ry, Rz

% Matriz de transformación
M = [D, -Rz, Ry; Rz, D, -Rx; -Ry, Rx, D];
T = [Tx; Ty; Tz];
p2 = p1 + T + M*p1;

6. Coordenadas geodésicas (método iterativo)

Se calcula la longitud directamente y la latitud y altura mediante un proceso iterativo que converge cuando las correcciones son menores que la tolerancia establecida:

x = p2(1);
lonF = atan(y ./ x);
pb = sqrt(x^2 + y^2);
latF = atan(z ./ (1 - e2));
corl = 1; corh = 1; h = 0; count = 0;
while (abs(corl) > 1e-11 && abs(corh) > 1e-4)
    N = a ./ sqrt(1 - e2*sin(latF)^2);
    hn = (z ./ sin(latF)) - N.*(1 - e2);
    latn = atan((z.*(N + hn)) ./ (pb.*(N.*(1 - e2) + hn)));
    corl = latF - latn;
    corh = h - hn;
    count = count + 1;
    latF = latn;
    h = hn;
end

% Conversión rad → DMS
latDeg = degrees2dms(latF*(180/pi));
lonDeg = degrees2dms(lonF*(180/pi));

7. Ascensión recta, declinación y tiempo sidéreo

Datos de partida: ascensión recta α = 11h 12m 34s (junio, CIV, oficial, Este), declinación δ igual a la latitud de la estación 2.

dec = lat2;
alfa_h = 11 + 12/60 + 34/3600;
alfa_rad = deg2rad(alfa_h * 15);

% Ángulo horario al orto (ESTE = negativo)
cosH = tan(dec) / tan(lat1);
H_rad = -acos(cosH);
H_horas = rad2deg(H_rad) / 15;   % negativo, +24 h

% Tiempo sidéreo local
theta_local = alfa_h + H_horas;
theta_G = theta_local - (rad2deg(lon1)/15);
if theta_local < 0, theta_local = theta_local + 24; end
thetaG = theta_local - (rad2deg(lon1)/15);
if thetaG >= 24, thetaG = thetaG - 24; end
if thetaG < 0, thetaG = thetaG + 24; end

8. Tiempo civil y hora oficial

theta0_h = 17 + 49/60 + 15.81/3600;
UT1 = (thetaG - theta0_h) / K;
if UT1 < 0, UT1 = UT1 + 24/K; end
ut_h = floor(UT1);
ut_m = floor((UT1 - ut_h)*60);
ut_s = (UT1 - ut_h - ut_m/60)*3600;
fprintf(' %dh %dm %.1fs\n', ut_h, ut_m, ut_s);

% Hora oficial (UTC+2)
hora_oficial = UT1 + 2;
h = floor(hora_oficial);
m = floor((hora_oficial - h)*60);
s = (hora_oficial - h - m/60)*3600;

9. Diferencia de tiempo sidéreo en el ocaso

Se calcula el ángulo horario de puesta de sol en ambas estaciones y la diferencia entre ellas:

cosH_yebe = -tan(dec)*tan(lat1);
H_yebe = acos(cosH_yebe);
H_yebe_h = rad2deg(H_yebe) / 15;
theta_yebe = alfa_h + H_yebe_h;

cosH_vill = -tan(dec) * tan(lat2);
H_vill = acos(cosH_vill);
H_vill_h = rad2deg(H_vill) / 15;
theta_vill = alfa_h + H_vill_h;

dif_H = abs(H_yebe_h - H_vill_h);
h_dif = floor(dif_H);
m_dif = floor((dif_H - h_dif)*60);
s_dif = (dif_H - h_dif - m_dif/60)*3600;

Nota: el código original contenía variables sin definir (por ejemplo, K, VxY, Ty, Tz, Ry, Rz, y2005, z2005), sintaxis incompleta (como lon1y2..lonT, incy, incz) y fórmulas fragmentadas (como Sm), que deben completarse antes de ejecutar el programa.