Longitud/latitud terrestre bajo un satélite (coordenadas cartesianas) en una época específica

Desde el punto en que te encuentras, necesitas:

1. Calcula las coordenadas XYZ de tu estación terrestre
2. Encuentre la línea en el espacio entre el satélite y la estación terrestre
3. Encuentre las coordenadas XYZ del punto en la intersección de esa línea con el elipsoide definido por todos los puntos con elevaciones de 300 km
4. Encuentra la latitud y la longitud de ese punto.

Para (1) puede usar esta función de matlab:

function [X,Y,Z] = ll2xyz(lat,lon,alt)
    b=6356752.3141;
    a=6378137.0;
    lat=lat*pi/180;
    %transformation between geografic and geocentric latitude
    gclat=atan(b^2*tan(lat)/a^2);
    lon=lon*pi/180;
    R=sqrt(1./(tan(gclat).^2/b^2 + 1/a^2));
    Z=R.*tan(gclat);
    Z=Z+alt.*sin(lat);
    R=R+alt.*cos(lat);
    X=R.*cos(lon);
    Y=R.*sin(lon);

Para (2) y (3) esta función hará el truco

function [X Y Z] = rectaxelip(x1,x2,y1,y2,z1,z2,alt)
%rectaxelip calculates the point XYZ were a straight line that pass by points P1
%and P2, intersects a ellipsoid alt meters bigger than WGS84 
%ellispsoid (on each axis)

%The director cosines of the line are 
    d=sqrt((x1-x2).^2.*(y1-y2).^2.*(z1-z2).^2);
    L=(x1-x2)./d;
    M=(y1-y2)./d;
    N=(z1-z2)./d;
    %its parametric  for is
    %x = x1 + L * t
    %y = y1 + M * t
    %z = z1 + N * t
    %Now we look for the intersection
    %the ellipsoid would have axes bp and ap

    b=6356752.3141; 
    a=6378137.0;
    bp=b+alt; 
    ap=a+alt;

    %The equation of such ellipsod would be
    %  x^2/ap^2+y^2/ap^2+z^2/bp^2=1
    %subtituting and reorganizng you get a quadratic equation like
    %   A x^2 + B x + C = 0
    %with
    A=(L./ap).^2+(M./ap).^2+(N./bp).^2;
    B=2*((x1.*L./(ap.^2))+(y1.*M./(ap.^2))+(z1.*N./(bp.^2)));
    C=((x1./ap).^2+(y1./ap).^2+(z1./bp).^2)-1;
    ta=(-B+sqrt(B.^2 - (4*A.*C)))./(2*A);
    tb=(-B-sqrt(B.^2 - (4*A.*C)))./(2*A);

    %then
    x1a = x1 + L .* ta;
    y1a = y1 + M .* ta;
    z1a = z1 + N .* ta;
    x1b = x1 + L .* tb;
    y1b = y1 + M .* tb;
    z1b = z1 + N .* tb;
    %Now the distance from each point to the satellite is
    da=sqrt((x1a-x2).^2+(y1a-y2).^2+(z1a-z2).^2);
    db=sqrt((x1b-x2).^2+(y1b-y2).^2+(z1b-z2).^2);
    ESta=da<db;
    t=ta.*ESta+tb.*(~ESta);

    %then
    X = x1 + L .* t;
    Y = y1 + M .* t;
    Z = z1 + N .* t;

Y para (4) solo necesita una función que calcule las coordenadas de latitud y longitud desde las posiciones XYZ, y esta lo haría:

function [lat,lon,alt] = xyz2ll(x,y,z1)
%For some reason the error spikes at +-45
a=6356752.3141; 
b=6378137.0;

signos=sign(z1);
z1=-abs(z1);

lon=atan2(y,x)*180/pi; %Longitude
d1=sqrt(x.^2+y.^2);
z2=-sqrt(1./((1/(a^2))+(1./(b^2*(z1./d1).^2))));
d2=z2./(z1./d1);
p2=-a*d2./(b^2*sqrt(1-(d2.^2/(b^2)))); %slope of the tangent to the ellipse in the point
for i=1:5
    %points in the palne
    dp=(z1-z2-(d1./p2)+(p2.*d2))./(p2-(1./p2));
    zp=(p2.*dp)+z2-(p2.*d2);
    %points in the ellipse
    z2=-sqrt(1./((1/(a^2))+(1./(b^2*(zp./dp).^2))));
    d2=z2./(zp./dp);
    p2=-a*d2./(b^2*sqrt(1-(d2.^2/(b^2))));
end
lat=-90-(atan(p2)*180/pi);   %geographic latitude
lat=abs(lat).*signos;
ecuador=find(abs(z1)<0.001);
lat(ecuador)=0;
alt=sqrt(d1.^2+z1.^2)-sqrt(d2.^2+z2.^2);
alt(ecuador)=d1(ecuador)-a;

Y luego está listo y tiene la posición en el suelo del punto donde la señal GPS que llega a la estación base cruza la altitud de 300 km.

Escribí estas funciones yo mismo y las probé extensamente, así que sé que funcionan, pero tenga en cuenta que la última función, a pesar de que tiene muy buena precisión, por una razón que no he investigado mucho, parece producir grandes errores en latitudes. de +-45°.

Lo último a tener en cuenta. Si quieres más precisión y evitar calcular las coordenadas XYZ de los satélites GPS, puedes utilizar los archivos de efemérides sp3. Son fáciles de descargar desde aquí . Hay un archivo por día, y cada archivo tiene las posiciones XYZ de cada satélite tabuladas cada 15 minutos. Para obtener posiciones en cualquier momento, la interpolación spline brinda una muy buena solución, más o menos lo mismo que la interpolación polinomial de Neville recomendada.