Земная долгота/широта под спутником (декартовы координаты) в определенную эпоху

С того момента, как вы находитесь, вам необходимо:

1. Рассчитайте координаты XYZ вашей наземной станции
2. Найдите линию в пространстве между спутником и наземной станцией
3. Найдите координаты XYZ точки на пересечении этой линии с эллипсоидом, определяемым всеми точками с высотами 300 км.
4. Найдите широту и долготу этой точки.

Для (1) вы можете использовать эту функцию 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);

Для (2) и (3) эта функция поможет

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;

А для (4) вам просто нужна функция, которая вычисляет координаты широты и долготы из позиций XYZ, и это сделает это:

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;

И тогда вы готовы, и у вас есть наземное положение точки, где сигнал GPS, достигающий базовой станции, пересекает высоту 300 км.

Я сам написал эти функции и тщательно их протестировал, поэтому знаю, что они работают, но имейте в виду, что последняя функция, несмотря на то, что она имеет очень хорошую точность, по причине, по которой я мало исследовал, похоже, дает большие ошибки на широтах. +-45°.

Последнее, что нужно отметить. Если вы хотите большей точности и не рассчитываете координаты XYZ спутников GPS, вы можете использовать файлы эфемерид sp3. Их легко скачать отсюда . Существует один файл в день, и в каждом файле координаты XYZ каждого спутника заносятся в таблицу каждые 15 минут. Для получения позиций в любое время сплайн-интерполяция дает очень хорошее решение, почти такое же, как рекомендуемая полиномиальная интерполяция Невилла.