algorithm - Trilateration using 3 latitude/longitude points, and 3 distances?

I want to find out an unknown target location (latitude and longitude co-ordinates). There are 3 known points (latitude and longitude co-ordinate pairs) and for each point a distance in kilometers to the target location. How can I calculate the co-ordinates of the target location?

For example, say I have the following data points

37.418436,-121.963477   0.265710701754km
37.417243,-121.961889   0.234592423446km
37.418692,-121.960194   0.0548954278262km

What I'd like is what is the math for a function that takes that as input and returns 37.417959,-121.961954 as output.

I understand how to calculate the distance between two points, from http://www.movable-type.co.uk/scripts/latlong.html I understand the general principle that with three circles like these you get exactly one point of overlap. What I'm hazy on is the math needed to calculate that point with this input.

Answer

After some looking around at Wikipedia and the same question/answer at StackOverflow, I figured I would take a stab at it, and try to fill in the gaps.

First off, Not sure where you got the output, but it appears to be wrong. I plotted the points in ArcMap, buffered them to the distances specified, ran intersect to on the buffers, and then captured the vertex of intersection to get the solutions. Your proposed output is the point in green. I calculated the value in the callout box, which is about 3 meters of what ArcMap gave for solution derived from the intersect.

The math on the wikipedia page isn't too bad, just need to covert your geodetic coordinates to the cartesian ECEF, which can be found here. the a/x +h terms can be replaced by the authalic sphere radius, if you aren't using an ellipsoid.

Probably easiest just give you some well(?) documented code, so here it is in python

import math
import numpy

#assuming elevation = 0
earthR = 6371

LatA = 37.418436
LonA = -121.963477
DistA = 0.265710701754
LatB = 37.417243
LonB = -121.961889
DistB = 0.234592423446
LatC = 37.418692
LonC = -121.960194
DistC = 0.0548954278262


#using authalic sphere
#if using an ellipsoid this step is slightly different
#Convert geodetic Lat/Long to ECEF xyz
#   1. Convert Lat/Long to radians
#   2. Convert Lat/Long(radians) to ECEF
xA = earthR *(math.cos(math.radians(LatA)) * math.cos(math.radians(LonA)))
yA = earthR *(math.cos(math.radians(LatA)) * math.sin(math.radians(LonA)))
zA = earthR *(math.sin(math.radians(LatA)))

xB = earthR *(math.cos(math.radians(LatB)) * math.cos(math.radians(LonB)))

yB = earthR *(math.cos(math.radians(LatB)) * math.sin(math.radians(LonB)))
zB = earthR *(math.sin(math.radians(LatB)))

xC = earthR *(math.cos(math.radians(LatC)) * math.cos(math.radians(LonC)))
yC = earthR *(math.cos(math.radians(LatC)) * math.sin(math.radians(LonC)))
zC = earthR *(math.sin(math.radians(LatC)))

P1 = numpy.array([xA, yA, zA])
P2 = numpy.array([xB, yB, zB])
P3 = numpy.array([xC, yC, zC])


#from wikipedia
#transform to get circle 1 at origin
#transform to get circle 2 on x axis
ex = (P2 - P1)/(numpy.linalg.norm(P2 - P1))
i = numpy.dot(ex, P3 - P1)
ey = (P3 - P1 - i*ex)/(numpy.linalg.norm(P3 - P1 - i*ex))
ez = numpy.cross(ex,ey)
d = numpy.linalg.norm(P2 - P1)
j = numpy.dot(ey, P3 - P1)


#from wikipedia
#plug and chug using above values
x = (pow(DistA,2) - pow(DistB,2) + pow(d,2))/(2*d)
y = ((pow(DistA,2) - pow(DistC,2) + pow(i,2) + pow(j,2))/(2*j)) - ((i/j)*x)

# only one case shown here
z = numpy.sqrt(pow(DistA,2) - pow(x,2) - pow(y,2))

#triPt is an array with ECEF x,y,z of trilateration point

triPt = P1 + x*ex + y*ey + z*ez

#convert back to lat/long from ECEF
#convert to degrees
lat = math.degrees(math.asin(triPt[2] / earthR))
lon = math.degrees(math.atan2(triPt[1],triPt[0]))

print lat, lon