title: Trilateration with approximated distances using latitude and longitude
body: <p>I can't quite figure this one out.</p>
  <p>I am trying to approximate the location (latitude / longitude) of a beacon based on 3 distance measurements from 3 fixed locations. However the distance readings available may have an error of up to 1 km.</p>
  <p>Similar questions regarding trilateration have been asked <a href="https://stackoverflow.com/questions/2813615/trilateration-using-3-latitude-and-longitude-points-and-3-distances">here (with precise measurements)</a>, <a href="https://gis.stackexchange.com/questions/66/trilateration-using-3-latitude-longitude-points-and-3-distances">here</a>, <a href="https://stackoverflow.com/questions/23957325/trilateration-with-distance-measurement-errors">here (distance measurement errors in Java, but not in lat/lon coordinates and no answers)</a> as well as others. I also managed to dig up <a href="https://inside.mines.edu/%7Ewhereman/papers/Murphy-Hereman-Trilateration-MCS-07-1995.pdf" rel="nofollow noreferrer">this paper</a> dealing with imperfect measurement data, however it for one assumes a cartesian coordinate system and is also rather mathematical than close to a usable implementation.</p>
  <p>So none of above links and answers are really applicable to the following problem:</p>
  <ul>
  <li>All available distance measurements are <strong>approximated</strong> in km (where data most frequently contains readings in-between 1 km and 100 km, in case this matters)</li>
  <li>Measurement errors of up to 1 km are possible.</li>
  <li>3 distance measurements are performed based on 3 fixed (latitude / longitude known) positions.</li>
  <li>target approximation should also be a latitude / longitude combination.</li>
  </ul>
  <p>So far I have adapted <a href="https://gis.stackexchange.com/a/415">this Answer</a> to C#, however I noticed that due to the measurement inaccuracies this algorithm does not work (as the algorithm assumes the 3 distance-circles to perfectly intersect with each other):</p>
  <pre class="lang-csharp prettyprint-override"><code>public static class Trilateration
  {
      public static GeoLocation Compute(DistanceReading point1, DistanceReading point2, DistanceReading point3)
      {
          // not my code :P
          // assuming elevation = 0
          const double earthR = 6371d;
  
          //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
          //   2d. Convert Lat/Long(radians) to ECEF
          double xA = earthR * (Math.Cos(Radians(point1.GeoLocation.Latitude)) * Math.Cos(Radians(point1.GeoLocation.Longitude)));
          double yA = earthR * (Math.Cos(Radians(point1.GeoLocation.Latitude)) * Math.Sin(Radians(point1.GeoLocation.Longitude)));
          double zA = earthR * Math.Sin(Radians(point1.GeoLocation.Latitude));
  
          double xB = earthR * (Math.Cos(Radians(point2.GeoLocation.Latitude)) * Math.Cos(Radians(point2.GeoLocation.Longitude)));
          double yB = earthR * (Math.Cos(Radians(point2.GeoLocation.Latitude)) * Math.Sin(Radians(point2.GeoLocation.Longitude)));
          double zB = earthR * (Math.Sin(Radians(point2.GeoLocation.Latitude)));
  
          double xC = earthR * (Math.Cos(Radians(point3.GeoLocation.Latitude)) * Math.Cos(Radians(point3.GeoLocation.Longitude)));
          double yC = earthR * (Math.Cos(Radians(point3.GeoLocation.Latitude)) * Math.Sin(Radians(point3.GeoLocation.Longitude)));
          double zC = earthR * Math.Sin(Radians(point3.GeoLocation.Latitude));
  
          // a 64 bit Vector3 implementation :)
          Vector3_64 P1 = new(xA, yA, zA);
          Vector3_64 P2 = new(xB, yB, zB);
          Vector3_64 P3 = new(xC, yC, zC);
  
          //from wikipedia
          //transform to get circle 1 at origin
          //ransform to get circle 2d on x axis
          Vector3_64 ex = (P2 - P1).Normalize();
          double i = Vector3_64.Dot(ex, P3 - P1);
          Vector3_64 ey = (P3 - P1 - i * ex).Normalize();
          Vector3_64 ez = Vector3_64.Cross(ex, ey);
          double d = (P2 - P1).Length;
          double j = Vector3_64.Dot(ey, P3 - P1);
  
          //from wikipedia
          //plug and chug using above values
          double x = (Math.Pow(point1.DistanceKm, 2d) - Math.Pow(point2.DistanceKm, 2d) + Math.Pow(d, 2d)) / (2d * d);
          double y = ((Math.Pow(point1.DistanceKm, 2d) - Math.Pow(point3.DistanceKm, 2d) + Math.Pow(i, 2d) + Math.Pow(j, 2d)) / (2d * j)) - ((i / j) * x);
  
          // only one case shown here
          double z = Math.Sqrt(Math.Pow(point1.DistanceKm, 2d) - Math.Pow(x, 2d) - Math.Pow(y, 2d));
  
          //triPt is a vector with ECEF x,y,z of trilateration point
          Vector3_64 triPt = P1 + x * ex + y * ey + z * ez;
  
          //convert back to lat/long from ECEF
          //convert to degrees
          double lat = Degrees(Math.Asin(triPt.Z / earthR));
          double lon = Degrees(Math.Atan2(triPt.Y, triPt.X));
  
          return new GeoLocation(lat, lon);
      }
  
      [MethodImpl(MethodImplOptions.AggressiveInlining)]
      private static double Radians(double degrees) =&gt;
          degrees * Math.Tau / 360d;
  
      [MethodImpl(MethodImplOptions.AggressiveInlining)]
      private static double Degrees(double radians) =&gt;
          radians * 360d / Math.Tau;
  }
  </code></pre>
  <p>Above code most often than not <em>does not</em> work in my case, and instead only returns &quot;Not a number&quot; as it tries to take the square root of a negative number when calculating the final <code>z</code> value (due to measurement inaccuracies).</p>
  <p>In my case measurements may return data like this (visualized with some random online tool):</p>
  <p><a href="https://i.stack.imgur.com/BHXpN.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/BHXpN.png" alt="enter image description here" /></a></p>
  <p>where only 2 or even none of the distance circles intersect:</p>
  <p><a href="https://i.stack.imgur.com/CmEqa.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/CmEqa.png" alt="enter image description here" /></a></p>
  <p>What I am looking for is the an algorithm returning the best possible approximation of the target location based on three distance measurements with a known maximum error of 1 km or further approaches I could take.</p>
  <p>I have also thought of iterating over points on the circles to then determining the minimum average distance to all the points on the other circles but the 3-dimensional sphere geometry of the earth is giving me a headache. Also there's probably a way better and simpler approach to this  which I just can't figure out right now.</p>
  <p>As this is more of an algorithmic problem, rather than any language-specific thing, I appreciate any help in whatever programming language, pseudo code or natural language.</p>