¿Conversión de coordenadas ECEF en LLA geodésico?

Question

Estoy intentando convertir las coordenadas ECEF en LLA geodésicas.

Estas son las coordenadas:

X = 1199816.63380
Y = -4812201.6785
Z = 4016140.62425

He comprobado las coordenadas utilizando este y devuelve valores correctos:

Lat: 39.111677542
Long: -76
Alt: 12068.6198

Pero usando mi código, obtengo estos valores:

Lat: 179.79495170382071
Long: -75.99999999999999
Alt: NaN

Y este es mi código:

const RADIUS_EQUATOR = 6378137.0;
const RADIUS_POLES = 6356752.3;
const e = Math.sqrt((RADIUS_EQUATOR*RADIUS_EQUATOR / RADIUS_POLES*RADIUS_POLES) / RADIUS_EQUATOR*RADIUS_EQUATOR);
const eprime = Math.sqrt((RADIUS_EQUATOR*RADIUS_EQUATOR - RADIUS_POLES*RADIUS_POLES) / RADIUS_POLES*RADIUS_POLES);

function radianToDegree(r) {
  return r * (180 / Math.PI);
}

function ECRtoLLA(ECR) {
  let X = ECR.x;
  let Y = ECR.y;
  let Z = ECR.z;

  let a = RADIUS_EQUATOR;
  let b = RADIUS_POLES;

  //Auxiliary values first
  let p = Math.sqrt(X*X + Y*Y);
  let theta = Math.atan((Z*a) / (p*b));

  let sintheta = Math.sin(theta);
  let costheta = Math.cos(theta);

  let num = Z + eprime * eprime * b * sintheta * sintheta * sintheta;
  let denom = p - e * e * a * costheta * costheta * costheta;

  //Now calculate LLA
  let latitude  = Math.atan2(num, denom);
  let longitude = Math.atan2(Y, X);
  let N = getN(latitude, a);
  let altitude  = (p / Math.cos(latitude)) - N;

  if (X < 0 && Y < 0) {
    longitude = longitude - Math.PI;
  }

  if (X < 0 && Y > 0) {
    longitude = longitude + Math.PI;
  }

  latitude = radianToDegree(latitude);
  longitude = radianToDegree(longitude);

  return {
    latitude: latitude,
    longitude: longitude,
    altitude: altitude
  };
}

console.info(ECRtoLLA({
  x: 1199816.63380,
  y: -4812201.6785,
  z: 4016140.62425
}));

¿Alguna indicación o ayuda?

Answer 1

Parece que he calculado mal e y eprime .

Después de arreglarlos, obtuve los resultados correctos. Lo siento, esta es la corrección

const RADIUS_EQUATOR = 6378137.0;
const a = RADIUS_EQUATOR;
const asqr = a*a;

const FLATTENING_DENOM = 298.257223563
const FLATTENING = 1 / FLATTENING_DENOM;
const RADIUS_POLES = RADIUS_EQUATOR*(1-FLATTENING);
const b = RADIUS_POLES;
const bsqr = b*b;

const e = Math.sqrt((asqr - bsqr) / asqr);
const eprime = Math.sqrt((asqr - bsqr) / bsqr);