Ground longitude/latitude under a satellite (cartesian coordinates) at a specfic epoch

Question

The script I'm wanting to develop uses the cartesian coordinates (XYZ) from a satellite, and in conjunction with the range, elevation and azimuth from a location, I then take a satellite’s orbital information and get the ground longitude/latitude under that satellite at a given time.

One step further from this: imagine the signal from a satellite piercing the atmosphere at exactly 300km above sea level. At this particular point when altitude is 300km, I need to calculate the ground longitude/latitude.

In the pyephem module there appears to be already a method (ephem.readtle) that can achieve this, but for TLE (two line element) data only. I'd like to use a satellite's cartesian coordinates to develop this. Is there such a method already out there? Or perhaps somebody with experience in this domain can point me in the right direction.

A similar question already exists referring to ECEF from Azimuth, Elevation, Range and Observer Lat,Lon,Alt, but it's not the same problem.

Here's what I have developed already:

• satellite cartesian coordinates, XYZ
• azimuth, elevation and range of satellite from ground station
• ground station coordinates in lat, long, height above sea level

Here's what I need: ground longitude/latitude under a satellite at a specific epoch, and in particular where the piercing point in the atmosphere (the point which the signal from the satellite pierces the atmosphere) is 300km altitude.

---

Ok, so I didn't resolve yet the issue with how to resolve for ground tracks for altitudes of 300km, but I believe the method I've wrote converting XYZ to ellipsoidal is complete:

 def cartesian_to_ellipsoidal(self):
    x = 4433469.9438
    y = 362672.7267
    z = 4556211.6409
    r = np.sqrt(x**2 + y**2 + z**2)

    # WGS-84 PARAMETERS, semimajor and semiminor axis
    a = 6378137.0
    b = 6356752.314 

    # Eccentricity
    e_squared = (a**2 - b**2) / a**2

    # Auxiliary quantities
    p = np.sqrt(x**2 + y**2)

    # Latitude (phi) & Longitude (lam)
    phi = np.rad2deg(np.arctan(z / ((1- e_squared) * p)))
    lam = np.rad2deg(np.arctan(y/x))

    # Radius of curvature in prime vertical               
    N = a / np.sqrt(1 - e_squared * (np.sin(np.deg2rad(phi)))**2)

    # Altitude
    h = (p / np.cos(np.deg2rad(phi))) - N

    return lam, phi, h

The XYZ coordinates are taken as a sample from this Matlab worked example Matlab worked example. The results are to an accuracy I'm satisfied with. What I don't understand and perhaps a dumb question, but when a ground track is calculated along with the altitude, why is the altitude not zero? One would expect since a ground track is mapped to the surface of the earth, the altitude should be zero.

One last point: the desired ground track was for where the line betweent the satellite and ground station pierces the amosphere 300km above the earth. Considering distances of 20000km (26000 km radius), then would an adjustment to my code - just to compensate for the desired 300km altitude scenario - make much a difference to outcome? If not, ground track data may just suffice.

Answer 1

From the point you are, you need to:

1. Calculate the XYZ coordinates of your ground station
2. Find the line in space between the satellite and the ground station
3. Find the XYZ coordinates of the point at the intersection of that line with the ellipsoid defined by all the points with elevations of 300 km
4. Find the latitude and longitude of that point.

For (1) you can use this matlab function:

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);

For (2) and (3) this function will do the trick

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;

And for (4) you just need a function that calculate lat lon coordinates from XYZ positions, and this one would do it:

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;

And then you are ready and you have the ground position of the point were the GPS signal reaching the base station crosses the 300 km altitude.

I wrote these functions myself and I've tested them extensively so I know they work, but be aware that the last function, despite it has very good accuracy, for a reason I haven't investigated much, seem to produce large errors at latitudes of +-45°.

Last thing to note. If you want more accuracy and avoid calculating the XYZ coordinates of the GPS satellites, you can use the sp3 ephemeris files. They are easy to download from here. There is one file per day, and each file have the XYZ positions of each satellite tabulated every 15 minutes. To get positions any time spline interpolation gives a very good solution, pretty much the same as the recommended Neville's polynomial interpolation.