Python函数和C函数在精度上的差异

Question

我正在使用一些代码，从以地球为中心的地球固定坐标到纬度/经度/高度坐标进行数学转换。这个函数的一次迭代是用Python编写的，另一次是用C编写的。C方法是在Python方法编写之后出现的，应该是直接翻译。然而，我发现Python方法虽然使用浮点值进行计算，但得到了正确的答案，而C方法，即使使用双精度，似乎只在小数点后5位取整，因此在尝试将ECEF转换为LLA坐标时得到了错误的答案。以下是转换的Python方法，以及一个坐标示例：

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)

以下方法是转换的C方法，使用双数据类型，考虑到这将确保最大精度：

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

Python方法返回的LLA坐标为：

54.99999538240099, 54.00000006037918, 8.26665045041591

而C方法返回的是：

55.0268382213345930,54.0000000603791790,4279.4874338442459702

所以海拔高度的计算似乎是个问题

Answer 1

有一个小小的转录错误。换线

double lla_latitude = atan2((ecef_altitude + pow(ep,2)*b*pow(sin(th),2)),p-esq*a*pow(cos(th),3));

到

double lla_latitude = atan2((ecef_altitude + pow(ep,2)*b*pow(sin(th),3)),p-esq*a*pow(cos(th),3));

当我做出改变时，我得到的高度是8.2666504495