I'm working with some code that does a mathematical conversion from Earth Centered Earth Fixed Coordinates to Latitude/Longitude/Altitude Coordinates. One iteration of this function is written in Python- another in C. The C method came after the Python method was written and is supposed to be a direct translation. However, I'm finding that the Python method, although it's using floating point values for it's calculations, gets a correct answer, while the C method, even when using doubles, seems to only be rounding after only 5 decimal places and therefore gets an incorrect answer when attempting to convert ECEF to LLA coordinates. Here is the Python method for the conversion, along with an example coordinate:
import math
import numpy as np
def ecef2lla(eceflist):
x = float(eceflist[0])
y = float(eceflist[1])
z = float(eceflist[2])
a = 6378137 # radius
e = 8.1819190842622e-2 # eccentricity
e_string = "E value:" + " " + str(e)
print(e_string)
asq = math.pow(a,2)
asq_string = "asq value:" + " " + str(asq)
print(asq_string)
esq=math.pow(e,2)
esq_string = "esq value:" + " " + str(esq)
print(esq_string)
b = math.sqrt( asq * (1-esq) )
b_string = "b value:" + " " + str(b)
print(b_string)
bsq = math.pow(b,2)
ep = math.sqrt( (asq - bsq)/bsq)
ep_string = "ep value:" + " " + str(ep)
print(ep_string)
p = math.sqrt( math.pow(x,2) + math.pow(y,2) )
p_string = "p value:" + " " + str(p)
print(p_string)
th = math.atan2(a*z, b*p)
th_string = "th value:" + " " + str(th)
print(th_string)
lon = math.atan2(y,x)
lon_string = "lon value:" + " " + str(lon)
print(lon_string)
lat = math.atan2( (z + math.pow(ep,2)*b*math.pow(math.sin(th),3) ), (p - esq*a*math.pow(math.cos(th),3)))
lat_string = "lat value:" + " " + str(lat)
print(lat_string)
N = a/( math.sqrt(1-esq*math.pow(math.sin(lat),2)) )
N_divisor = math.sqrt(1-esq*math.pow(math.sin(lat),2))
ND_string = "N divisor value is" + " " + str(N_divisor)
print(ND_string)
N_string = "N value:" + " " + str(N)
print(N_string)
alt = p / math.cos(lat) - N
alt_string = "alt value:" + " " + str(alt)
print(alt_string)
lon = lon % (2*math.pi)
ret = [lat*180/math.pi, lon*180/math.pi, alt]
return ret
if __name__ == '__main__':
lla_coords = ecef2lla([2155172.63, 2966340.65, 5201390])
print(lla_coords)
The following method is the translated C method, using double data types, with thinking being this would ensure maximum precision:
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <unistd.h>
void ecef_to_lla(double *coordinate){
double ecef_latitude = coordinate[0];
double ecef_longitude = coordinate[1];
double ecef_altitude = coordinate[2];
double a = 6378137;
//double e = 8.1819190842622e-2;
double e = 0.081819190842622;
printf("E value: %.16f\n",e);
double asq = pow(a,2);
printf("asq string: %.16f\n",asq);
double esq = pow(e,2);
printf("esq string: %.16f\n",esq);
double b = sqrt(asq * (1-esq));
printf("b string: %.16f\n",b);
double bsq = pow(b,2);
double ep = sqrt((asq-bsq)/bsq);
printf("ep string: %.16f\n",ep);
double p = sqrt(pow(ecef_latitude,2) + pow(ecef_longitude,2));
printf("p string: %.16f\n",p);
double th = atan2(a*ecef_altitude,b*p);
printf("th string: %.16f\n",th);
double lla_longitude = atan2(ecef_longitude,ecef_latitude);
double lla_latitude = atan2((ecef_altitude + pow(ep,2)*b*pow(sin(th),2)),p-esq*a*pow(cos(th),3));
printf("Lat is %.16f\n",lla_latitude);
double n = a/(sqrt(1-esq*pow(sin(lla_latitude),2)));
double n_divisor = sqrt(1-esq*pow(sin(lla_latitude),2));
printf("n_divisor is %.16f\n",n_divisor);
printf("n is %.16f\n",n);
double lla_altitude = p / cos(lla_latitude) - n;
printf("alt is %.16f\n",lla_altitude);
lla_longitude = fmod(lla_longitude ,(2 * M_PI));
lla_latitude = fmod(lla_latitude ,(2 * M_PI));
lla_latitude = lla_latitude * 180/M_PI;
lla_longitude = lla_longitude * 180/M_PI;
coordinate[0] = lla_latitude;
coordinate[1] = lla_longitude;
coordinate[2] = lla_altitude;
}
int main(void){
double coordinates[3] = {2155172.63, 2966340.65, 5201390};
ecef_to_lla(coordinates);
}
The LLA coordinates returned by the Python method are:
54.99999538240099, 54.00000006037918, 8.26665045041591
while those returned by the C method are:
55.0268382213345930,54.0000000603791790,4279.4874338442459702
So the altitude calculationseems to be the issue
There's a tiny transcription error. Change the line
double lla_latitude = atan2((ecef_altitude + pow(ep,2)*b*pow(sin(th),2)),p-esq*a*pow(cos(th),3));
to
double lla_latitude = atan2((ecef_altitude + pow(ep,2)*b*pow(sin(th),3)),p-esq*a*pow(cos(th),3));
When I make that change, I get 8.2666504495 for the altitude.