Projection: geocentric to geographic3D

Consider a North Sea point with coordinates derived by GPS satellite in the WGS 84 geographical coordinate system with coordinates of:

           latitude    53 deg 48 min 33.82 sec N, 
           longitude 02 deg 07 min 46.38 sec E, 
    and ellipsoidal height 73.0m, 

whose coordinates are required in terms of the ED50 geographical coordinate system which takes the International 1924 ellipsoid. The three parameter datum shift from WGS 84 to ED50 for this North Sea area is given as dX = +84.87m, dY = +96.49m, dZ = +116.95m. 

The WGS 84 geographical coordinates convert to the following geocentric values using the above formulas for X, Y, Z:

   XA = 3771 793.97m
   YA =   140 253.34m
   ZA = 5124 304.35m

Applying the quoted datum shifts to these, we obtain new geocentric values now related to ED50:

   XB = 3771 878.84m
   YB =   140 349.83m
   ZB = 5124 421.30m

These convert to ED50 values on the International 1924 ellipsoid as:
           latitude    53 deg 48 min 36.565 sec N, 
           longitude 02 deg 07 min 51.477 sec E, 
    and ellipsoidal height 28.02 m, 

Note that the derived height is referred to the International 1924 ellipsoidal surface and will need a further correction for the height of the geoid at this point in order to relate it to Mean Sea Level.

Note: These formulas have been transcribed from EPSG Guidance Note #7-2. Users are encouraged to use that document rather than the text which follows as reference because limitations in the transcription will be avoided.

Latitude, P, and Longitude, L, in terms of Geographic Coordinate Reference System A may 
be expressed in terms of a geocentric (earth centred) Cartesian coordinate reference system X, Y, Z 
with the Z axis corresponding with the Polar axis positive northwards, the X axis through 
the intersection of the Greenwich meridian and equator, and the Y axis through the 
intersection of the equator with longitude 90 degrees E. If the prime meridian for geogCRS A is not 
Greewich, longitudes must first be transformed to their Greenwich equivalent. If the earth's 
spheroidal semi major axis is a, semi minor axis  b, and inverse flattening 1/f,  then

   X =   (nu + h) cos P cos L
   Y =   (nu + h) cos P sin L
   Z =  ((1 - e^2) nu + h) sin P

where nu is the prime vertical radius of curvature at latitude P and is equal to 
   nu = a /(1 - e^2*sin^2(P))^0.5,
   P and L are respectively the latitude and longitude (related to Greenwich) of the point 
   h is height above the ellipsoid, (topographic height plus geoidal height), and
   e is the eccentricity of the ellipsoid where e^2 = (a^2 -b^2)/a^2 = 2f -f^2
                                                                                                                                                 
Cartesian coordinates in geocentric coordinate reference system B may be used to derive geographical coordinates in terms of geographic coordinate reference system B by:
           
P = atan[(Z + eta b sin^3 q) / (p – e^2 a cos^3 q)]
L  = atan (Y/X)

where
eta = e^2 / (1 – e^2)
b = a(1 – f)
p = (X^2 + Y^2)^0.5
q = atan[(Z a) / (p b)]

and where L is relative to Greenwich. If the geographic system has a non Greenwich prime meridian, the Greenwich value of the local prime meridian should be applied to longitude.

Then
  	h (p / cos P) – nu

(Note that h is the height above the ellipsoid. This is the height value which is delivered by Transit and GPS satellite observations but is not the topographic height value which is normally used for national mapping and levelling operations. The topographic height is usually the height above mean sea level or an alternative 
level reference for the country. If one starts with a topographic height,  it will be necessary to convert it to an ellipsoid height before using the above transformation formulas. h = N + H, where N is the geoid height above the ellipsoid at the point and is sometimes negative, and H is the height of the point above the geoid. The height above the geoid is often taken to be that above mean sea level, perhaps with a constant correction applied. Geoid heights of points above the nationally used ellipsoid may not be readily available. For the WGS84 ellipsoid the value of N, 
representing the height of the geoid relative to the ellipsoid, can vary between values of -100m in the Sri Lanka area to +80m in the North Atlantic.)