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

# 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

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;