EPSG:1076

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.

LMP requires that the input local horizontal coordinates are referenced to true north, and that depths (D) are referenced to the reference surface of the vertical CRS. Depths are required as input but are not changed by the algorithm. They are used to compute the depth correction factor and therefore a depth of zero must be close to Mean Sea Level (strictly speaking to the reference ellipsoid).

Required inputs are:

Latitude of the Well Reference Point (WRP) in the target 2D geographic CRS;
Longitude of the Well Reference Point in the target 2D geographic CRS;
dcf flag: 0 = no depth correction is applied; 2 = a variable depth correction factor along the path is applied;
n[]: Array of length num_points containing local northing coordinates, oriented to true north;
e[]: Array of length num_points containing local easting coordinates, oriented to true east;
D[]: Array of length num_points containing depths below the Vertical Reference Surface at survey stations down the wellbore (referred to as TVD in the SPE paper).

The output of LMP is ellipsoidal coordinates of the wellbore trajectory in a geographic 2D CRS. If required, projected coordinates of the wellbore may be calculated from these.

The following equations assume that the local coordinates (n,e) and depth D are in metres.  If these use a different unit of measure, then as a first step they should be converted to metres.

When dcf flag = 2, output arrays LAT, LON (in radians) are computed as:

	for i=1:num_points
		ρ = a (1 –  e^2) / (1 –  e^2 sin^2(LAT[i-1]))^3/2
		ν = a / (1 –  e^2 sin^2(LAT[i-1]))^1/2
		LAT[i] = LAT[i-1] + (n[i] – n[i-1]) / (ρ – D[i]) 
		LON[i] = LON[i-1] + (e[i] – e[i-1]) / (ν – D[i]) / cos(LAT[i-1])

where for i=1, for the first station: LAT[0] and LON[0] are the geographic coordinates of the WRP in radians; n[0] = 0, e[0] = 0.

and
	a is the semi-major axis of the reference ellipsoid;
	e is the eccentricity of the reference ellipsoid (not to be confused with local easting);
	ρ is the radius of curvature in the meridian;
	ν is the radius of curvature in the prime vertical.

When dcf flag = 0, the equations above are modified to:
		LAT[i] = LAT[i-1] + (n[i] – n[i-1]) / ρ 
		LON[i] = LON[i-1] + (e[i] – e[i-1]) / {ν cos(LAT[i-1])}

Reverse case: 
To compute (n,e) from (LAT,LON,D):
	For i=1:num_points
		ρ = a (1 –  e^2) / (1 –  e^2 sin^2(LAT[i-1]))^3/2
		ν = a / (1 –  e^2 sin^2(LAT[i-1]))^1/2
		n[i] = n[i-1] + (LAT[i] – LAT[i-1]) * (ρ – D[i])
		e[i] = e[i-1] + (LON[i] – LON[i-1]) * (ν – D[i]) * cos(LAT[i-1])

where for the first station i=0 variables are initialised as in the forward case. Note that LAT[i-1], the latitude of the previous point, is used in the reverse formulas to compute the radii of curvature in order to be internally consistent in round trip calculations.

When dcf flag = 0, the equations above are modified to:
	n[i] = n[i-1] + (LAT[i] – LAT[i-1]) * ρ 
	e[i] = e[i-1] + (LON[i] – LON[i-1]) * ν * cos(LAT[i-1])