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

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.

• When you say cartesian, do you mean Earth-Centered Inertial coordinates? And you want to find the lat/lon directly under the satellite? Just drawing a line from the center of the earth up to the satellite, find out where the line intersects the surface (and 300km above the surface)? If XYZ come from the TLE's, what is the ground station for? – uhoh Aug 25 '16 at 13:19
• Check out questions about Skyfield here and here and also the website. It's easy to use! – uhoh Aug 25 '16 at 13:22
• @uhoh, yes ECEF. I'd like to find the ground track of the satellite in the lon/lat frame, but rather than the satellite itself, the point where the line from the satellite to the receive intersects 300km above the Earth. XYZ don't come from the TLE's but my own python script by way of satellite ephemerides. – pymat Aug 25 '16 at 13:40
• Oh, I understand now - where the line from satellite to ground station pierces 300km atlitude. So the XYZ of the satellite is in Earth Centered Earth Fixed (ECEF=rotating) coordinates, or Earth Centered Inertial (ECI=not-rotating)? – uhoh Aug 25 '16 at 15:36
• Can you try to clean this up into a more well (narrowly) defined and shorter question that can have a clear answer? There may be too much going on here to fit the stackexchange Q&A format. At least to me this seems more like "thinking out loud", – uhoh Aug 26 '16 at 15:14

## 1 Answer

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.

• When you say "the ellipsoid defined by all the points with elevations of 300 km" how do you "grow" the WGS84 ellipsoid by 300 km in altitude? What I'm curious is if each point on WGS84 moves on a local ellipsoid normal, or radially away from the Earth's center, or along a local gravity gradient ("up")? Also, if you'd like to add an answer to What is the Fischer 1960 Mercury Ellipsoid, and why is it called that? with some historical perspective, that would be great! I love your answers, and her story is quite interesting! – uhoh May 13 '18 at 15:33
• @uhoh As you can see in the code, it will just add the 300 km to the semi-major and semi-minor axes that define the reference ellipsoid. And I think the resulting ellipsoid is equivalent to the surface defined by all point at 300 km of elipsoidal elevation relative to the first ellipsoid. Elipsoidal elevations in turn, are measured perpendicular to the ellipsoid surface (i.e. min. distance), NOT radially from the ellipsoid center NOR considering gravity gradients in any way. Interesting question you pointed! But I'll refrain to answer this time just due to a lack of spare time to do it. Cheers – Camilo Rada May 13 '18 at 16:43
• Ah! I finally realized what it is that keeps puzzling me about this. "...with the ellipsoid defined by..." The resulting surface is not an ellipsoid. For small changes, it would probably look like one, but it can't be called one. Quickest way to check would be to take the reference ellipsoid, add 300 km to major and minor radii, calculate the shape of a proper ellipsoid with those radii, then compare the difference. Probably only a few tens of meters. Still, it's better not to call it an ellipsoid if it's not one. – uhoh May 19 '18 at 6:56