Celestial Programming : Latitude/Longitude to Geocentric XYZ

This takes geodedic coordinates (e.g. latitude and longitude from a GPS receiver) and converts it into rectangular coordinates (xyz). This is useful for computing topocentric coordinates. Before using it rotations need to be applied for Earth Rotation Angle (or GMST), and possibly precession, nutation, and polar motion depending on the accuracy required.

x = 3928199.4164297534m
y = 3097042.887854688m
z = 3917430.301268257m

Explanatory Supplement to the Astronomical Almanac 7.131

r=[xyz]=[ρcosθcosλρcosθsinλρsinθ]=[(aC+h)cosθcosλ(aC+h)cosθsinλ(aS+h)sinθ] r= \begin{bmatrix} x\\ y\\ z \end{bmatrix} = \begin{bmatrix} \rho\cos\theta'\cos\lambda\\ \rho\cos\theta'\sin\lambda\\ \rho\sin\theta' \end{bmatrix} = \begin{bmatrix} (aC+h)\cos\theta\cos\lambda \\ (aC+h)\cos\theta\sin\lambda \\ (aS+h)\sin\theta \end{bmatrix} C=(cos2θ+(1f)2sin2θ)12S=(1f)2Ca=6378136.6mf=1/298.25642 \begin{align*} C&=(\cos^2\theta + (1-f)^2\sin^2\theta)^{-\frac{1}{2}} \\ S&=(1-f)^2C \\ a&=6378136.6m\\ f&=1/298.25642\\ \end{align*} θ\theta' Geocentric latitude
θ\theta Geodedic latitude
λ\lambda Longitude
ρ\rho Geocentric distance
aa Earth's equitorial radius
ff Earth's flatening ratio
Components of rr will be in the same units as aa, meters in this case.

//Greg Miller (gmiller@gregmiller.net) 2022
//Released as public domain
//www.celestialprogramming.com

//Convert Geodedic Lat Lon to geocentric XYZ position vector
//All angles are input as radians
function convertGeodedicLatLonToITRFXYZ(lat,lon,height){
    //Algorithm from Explanatory Supplement to the Astronomical Almanac 3rd ed. P294
    const a=6378136.6;
    const f=1/298.25642;

    const C=Math.sqrt(((Math.cos(lat)*Math.cos(lat)) + (1.0-f)*(1.0-f) * (Math.sin(lat)*Math.sin(lat))));

    const S=(1-f)*(1-f)*C;
    
    const h=height;

    let r=new Array();
    r[0]=(a*C+h) * Math.cos(lat) * Math.cos(lon);
    r[1]=(a*C+h) * Math.cos(lat) * Math.sin(lon);
    r[2]=(a*S+h) * Math.sin(lat);
    
    return r;
}