function [lat,lon,hae]=toGeodetic(x,y, z) 
    c_wgs84_a = 6378137.0;
    c_wgs84_e2 = 0.00669437999013;

    p = sqrt(x * x + y * y);

    lon = atan2(y, x);
    lat = atan2(z / p, 0.01);
    n = computeRn(lat);
    hae = p / cos(lat) - n;
    old_hae = -1e-9;
    num = z / p;

    while (abs(hae - old_hae) > 1e-4) 
        old_hae = hae;
        den = 1 - c_wgs84_e2 * n / (n + hae);
        lat = atan2(num, den);
        n = computeRn(lat);
        hae = p / cos(lat) - n;
    end

    lat = rad2deg(lat);
    lon = rad2deg(lon);
