This answer is divided into multiple sections:
Analysis and Reduction of the Problem, showing how to find the desired point with "canned" routines.
Illustration: a Working Prototype, giving working code.
Example, showing examples of the solutions.
Pitfalls, discussing potential problems and how to cope with them.
ArcGIS implementation, comments about creating a custom ArcGIS tool and where to obtain the needed routines.
Analysis and Reduction of the Problem
Let's begin by observing that in the (perfectly round) spherical model there will always be a solution--in fact, exactly two solutions. Given base points A, B, and C, each pair determines its "perpendicular bisector," which is the set of points equidistant from the two given points. This bisector is a geodesic (great circle). Spherical geometry is elliptic: any two geodesics intersect (in two unique points). Thus, the intersection points of AB's bisector and BC's bisector are--by definition--equidistant from A, B, and C, thereby solving the problem. (See the first figure below.)
Things look more complicated on an ellipsoid, but because it is a small perturbation of the sphere, we can expect similar behavior. (The analysis of this would take us too far afield.) The complicated formulas used (internally within a GIS) to compute accurate distances on an ellipsoid are not a conceptual complication, though: the problem is basically the same. To see how simple the problem really is, let's state it somewhat abstractly. In this statement, "d(U,V)" refers to the true, fully accurate distance between points U and V.
Given three points A, B, C (as lat-lon pairs) on an ellipsoid, find a point X for which (1) d(X,A) = d(X,B) = d(X,C) and (2) this common distance is as small as possible.
These three distances all depend on the unknown X. Thus the differences in distances u(X) = d(X,A) - d(X,B) and v(X) = d(X,B) - d(X,C) are Real-valued functions of X. Again, somewhat abstractly, we may assemble these differences into an ordered pair. We will also use (lat, lon) as coordinates for X, allowing us to consider it as an ordered pair, too, say X = (phi, lambda). In this setup, the function
F(phi, lambda) = (u(X), v(X))
is a function from a portion of a two-dimensional space taking values in two-dimensional space and our problem reduces to
Find all possible (phi, lambda) for which F(phi, lambda) = (0,0).
Here is where the abstraction pays off: plenty of great software exists to solve this (purely numerical multidimensional root-finding) problem. The way it works is that you write a routine to compute F, then you pass it along to the software along with any information about restrictions on its input (phi must lie between -90 and 90 degrees and lambda must lie between -180 and 180 degrees). It cranks away for a fraction of a second and returns (typically) just one value of (phi, lambda), if it can find one.
There are details to handle, because there's an art to this: there are various solution methods to choose from, depending on how F "behaves"; it helps to "steer" the software by giving it a reasonable starting point for its search (this is one way we can obtain the nearest solution, rather than any other one); and you usually need to specify how accurate you would like the solution to be (so it knows when to stop the search). (For more on what GIS analysts need to know about such details, which come up a lot in GIS problems, please visit Recommend topics to be included in a Computer Science for Geospatial Technologies course and look in the "Miscellany" section near the end.)
Illustration: a Working Prototype
The analysis shows we need to program two things: a crude initial estimate of the solution and the calculation of F itself.
The initial estimate can be made by a "spherical average" of the three base points. This is obtained by representing them in geocentric Cartesian (x,y,z) coordinates, averaging those coordinates, and projecting that average back to the sphere and re-expressing it in latitude and longitude. The size of the sphere is immaterial and the calculations are thereby made straightforward: because this is just a starting point, we don't need ellipsoidal calculations.
For this working prototype I used Mathematica 8.
sphericalMean[points_] := Module[{sToC, cToS, cMean},
sToC[{f_, l_}] := {Cos[f] Cos[l], Cos[f] Sin[l], Sin[f]};
cToS[{x_, y_, z_}] := {ArcTan[x, y], ArcTan[Norm[{x, y}], z]};
cMean = Mean[sToC /@ (points Degree)];
If[Norm[Most@cMean] < 10^(-8), Mean[points], cToS[cMean]] / Degree
]
(The final If condition tests whether the average might fail clearly to indicate a longitude; if so, it falls back to a straight arithmetic mean of the latitudes and longitudes of its input--maybe not a great choice, but at least a valid one. For those using this code for implementation guidance, note that the arguments of Mathematica's ArcTan are reversed compared to most other implementations: its first argument is the x-coordinate, its second is the y-coordinate, and it returns the angle made by the vector (x,y).)
As far as the second part goes, because Mathematica--like ArcGIS and almost all other GISes--contains code to compute accurate distances on the ellipsoid, there's almost nothing to write. We just call the root-finding routine:
tri[a_, b_, c_] := Block[{d = sphericalMean[{a, b, c}], sol, f, q},
sol = FindRoot[{GeoDistance[{Mod[f, 180, -90], Mod[q, 360, -180]}, a] ==
GeoDistance[{Mod[f, 180, -90], Mod[q, 360, -180]}, b] ==
GeoDistance[{Mod[f, 180, -90], Mod[q, 360, -180]}, c]},
{{f, d[[1]]}, {q, d[[2]]}},
MaxIterations -> 1000, AccuracyGoal -> Infinity, PrecisionGoal -> 8];
{Mod[f, 180, -90], Mod[q, 360, -180]} /. sol
];
The most noteworthy aspect of this implementation is how it dodges the need to constrain the latitude (f) and longitude (q) by always computing them modulo 180 and 360 degrees, respectively. This avoids having to constrain the problem (which often creates complications). The control parameters MaxIterations etc. are tweaked to make this code provide the greatest possible accuracy it can.
To see it in action, let's apply it to the three base points given in a related question:
sol = tri @@ (bases = {{-6.28530175, 106.9004975375}, {-6.28955287, 106.89573839}, {-6.28388865789474, 106.908087643421}})
{-6.29692, 106.907}
The computed distances between this solution and the three points are
{1450.23206979, 1450.23206979, 1450.23206978}
(these are meters). They agree through the eleventh significant digit (which is too precise, actually, since distances are rarely accurate to better than a millimeter or so). Here's a picture of these three points (black), their three mutual bisectors, and the solution (red):
Example
To test this implementation and get a better understanding of how the problem behaves, here is a contour plot of the root mean square discrepancy in distances for three widely spaced base points. (The RMS discrepancy is obtained by computing all three differences d(X,A)-d(X,B), d(X,B)-d(X,C), and d(X,C)-d(X,A), averaging their squares, and taking the square root. It equals zero when X solves the problem and otherwise increases as X moves away from a solution, and so measures how "close" we are to being a solution at any location.)
The base points (60,-120), (10,-40), and (45,10) are shown in red in this Plate Carree projection; the solution (49.2644488, -49.9052992)--which required 0.03 seconds to compute--is in yellow. Its RMS discrepancy is less than three nanometers, despite all relevant distances being thousands of kilometers. The dark areas show small values of the RMS and the light areas show high values.
This map clearly shows another solution lies near (-49.2018206, 130.0297177) (computed to an RMS of two nanometers by setting the initial search value diametrically opposite the first solution.)
Pitfalls
Numerical instability
When the base points are nearly collinear and close together, all solutions will be nearly half a world away and extremely difficult to pin down accurately. The reason is that small changes in a location across the globe--moving it towards or away from the base points--will induce only incredibly tiny changes in the differences of distances. There's just not enough accuracy and precision built into the usual computation of geodetic distances to pin down the result.
For example, starting with the base points at (45.001, 0), (45, 0), and (44.999,0), which are separated along the Prime Meridian by only 111 meters between each pair, tri obtains the solution (11.8213, 77.745). The distances from it to the base points are 8,127,964.998 77; 8,127,964.998 41; and 8,127,964.998 65 meters, respectively. They agree to the nearest millimeter! I'm not sure how accurate this result may be, but would not be in the least surprised if other implementations returned locations far from this one, showing almost as good equality of the three distances.
Computation time
These calculations, because they involve considerable searching using complicated distance computations, are not fast, usually requiring a noticeable fraction of a second. Real-time applications need to be aware of this.
ArcGIS implementation
Python is the preferred scripting environment for ArcGIS (beginning with version 9). The scipy.optimize package has a multivariate rootfinder root which should do what FindRoot does in the Mathematica code. Of course ArcGIS itself offers accurate ellipsoidal distance calculations. The rest, then, is all implementation details: decide how the base point coordinates will be obtained (from a layer? typed in by the user? from a text file? from the mouse?) and how the output will be presented (as coordinates displayed on screen? as a graphic point? as a new point object in a layer?), write that interface, port the Mathematica code shown here (straightforward), and you will be all set.
