title: Finding Earth Coordinates (latitude,longitude), Distance (meters) and Bearing (angle)
body: <p>I need to deal with Earth coordinates in various ways. There is no function in C/C++ which does this straight away.<br>
  Referred below questions:</p>
  
  <ol>
  <li><a href="https://stackoverflow.com/questions/7222382/get-lat-long-given-current-point-distance-and-bearing">Python: Get lat/long given current point, distance and
  bearing</a></li>
  <li><a href="https://stackoverflow.com/questions/30002372/given-point-of-latitude-longitude-distance-and-bearing-how-to-get-the-new-la">C: Given point of (latitude,longitude), distance and bearing, How
  to get the new latitude and
  longitude</a></li>
  </ol>
  
  <p>From 1st one and <a href="http://www.movable-type.co.uk/scripts/latlong.html" rel="nofollow noreferrer">movable type scripts website</a>, I found that below are the formulas:  </p>
  
  <p><strong>Find Bearing(angle) between 2 coordinates</strong></p>
  
  <pre><code>x = cos(lat1Rad)*sin(lat2Rad) - sin(lat1Rad)*cos(lat2Rad)*cos(lon2Rad-lon1Rad);
  y = sin(lon2Rad-lon1Rad) * cos(lat2Rad);
  bearing = atan2(y, x);  // In radians; 
  // Convert to degrees and for -ve add 360
  </code></pre>
  
  <p><strong>Find Distance(meters) between 2 coordinates</strong></p>
  
  <pre><code>PI = 3.14159265358979323846, earthDiameterMeters = 2*6371*1000;
  x = sin((lat2Rad-lat1Rad) / 2);
  y = sin((lon2Rad-lon1Rad) / 2);
  meters = earthDiameterMeters * asin(sqrt(x*x + y*y*cos(lat1Rad)*cos(lat2Rad)));
  </code></pre>
  
  <p><strong>Find Coordinate from Coordinate+Distance+Angle</strong></p>
  
  <pre><code>meters *= 2 / earthDiameterMeters;
  lat2Rad = asin(sin(lat1Rad)*cos(meters) + cos(lat1Rad)*sin(meters)*cos(bearing));
  lon2Rad = lon1Rad + atan2(sin(bearing)*sin(meters)*cos(lat1Rad), 
                            cos(meters) - sin(lat1Rad)*sin(lat2Rad));
  </code></pre>
  
  <p>Below pseudo code should verify above 3 equations mutually:  </p>
  
  <pre><code>struct Coordinate { double lat, lon; } c1, c2;  
  auto degree = FindBearing(c1, c2);
  auto meters = FindDistance(c1, c2);
  auto cX = FindCoordiante(c1, degree, meters);
  </code></pre>
  
  <p>Now actually the answer comes <em>almost</em> near but not the correct. i.e. cX is Not equal to c2!<br>
  There is always a difference of <code>0.0005</code> difference in longitude value.
  e.g.</p>
  
  <pre><code>c1 = (12.968460,77.641308)  
  c2 = (12.967862,77.653130)  
  angle = 92.97         ^^^
  distance = 1282.74  
  cX = (12.967862,77.653613)
                        ^^^
  </code></pre>
  
  <p>I don't have much knowledge of mathematics' <a href="https://en.wikipedia.org/wiki/Haversine_formula" rel="nofollow noreferrer">Havesine Forumla</a>. But what I know is that from <a href="https://www.fcc.gov/encyclopedia/find-terminal-coordinates-given-bearing-and-distance" rel="nofollow noreferrer">the website fcc.gov</a>, the answer always comes correct.</p>
  
  <p><strong>What am I doing wrong?</strong>  </p>
  
  <h2><a href="http://ideone.com/I4sT0f" rel="nofollow noreferrer">Code</a> only for reference</h2>
  
  <p>Though the syntax is in C++, but all the math functions are from C and is easily portable in C as well (hence tagged for both)</p>
  
  <pre><code>#include&lt;iostream&gt;
  #include&lt;iomanip&gt;
  #include&lt;cmath&gt;
  
  
  // Source: http://www.movable-type.co.uk/scripts/latlong.html
  
  static const double PI = 3.14159265358979323846, earthDiameterMeters = 6371.0 * 2 * 1000;
  
  double degreeToRadian (const double degree) { return (degree * PI / 180); };
  double radianToDegree (const double radian) { return (radian * 180 / PI); };
  
  double CoordinatesToAngle (double latitude1,
                             const double longitude1,
                             double latitude2,
                             const double longitude2)
  {
    const auto longitudeDifference = degreeToRadian(longitude2 - longitude1);
    latitude1 = degreeToRadian(latitude1);
    latitude2 = degreeToRadian(latitude2);
  
    using namespace std;
    const auto x = (cos(latitude1) * sin(latitude2)) -
                   (sin(latitude1) * cos(latitude2) * cos(longitudeDifference));
    const auto y = sin(longitudeDifference) * cos(latitude2);
  
    const auto degree = radianToDegree(atan2(y, x));
    return (degree &gt;= 0)? degree : (degree + 360);
  }
  
  double CoordinatesToMeters (double latitude1,
                              double longitude1,
                              double latitude2,
                              double longitude2)
  {
    latitude1 = degreeToRadian(latitude1);
    longitude1 = degreeToRadian(longitude1);
    latitude2 = degreeToRadian(latitude2);
    longitude2 = degreeToRadian(longitude2);
  
    using namespace std;
    auto x = sin((latitude2 - latitude1) / 2), y = sin((longitude2 - longitude1) / 2);
  #if 1
    return earthDiameterMeters * asin(sqrt((x * x) + (cos(latitude1) * cos(latitude2) * y * y)));
  #else
    auto value = (x * x) + (cos(latitude1) * cos(latitude2) * y * y);
    return earthDiameterMeters * atan2(sqrt(value), sqrt(1 - value));
  #endif
  }
  
  std::pair&lt;double,double&gt; CoordinateToCoordinate (double latitude,
                                                   double longitude,
                                                   double angle,
                                                   double meters)
  {
    latitude = degreeToRadian(latitude);
    longitude = degreeToRadian(longitude);
    angle = degreeToRadian(angle);
    meters *= 2 / earthDiameterMeters;
  
    using namespace std;
    pair&lt;double,double&gt; coordinate;
  
    coordinate.first = radianToDegree(asin((sin(latitude) * cos(meters))
                               + (cos(latitude) * sin(meters) * cos(angle))));
    coordinate.second = radianToDegree(longitude
                      + atan2((sin(angle) * sin(meters) * cos(latitude)),
                      cos(meters) - (sin(latitude) * sin(coordinate.first))));
  
    return coordinate;
  }
  
  int main ()
  {
    using namespace std;
    const auto latitude1 = 12.968460, longitude1 = 77.641308,
               latitude2 = 12.967862, longitude2 = 77.653130;
  
    cout &lt;&lt; std::setprecision(10);
    cout &lt;&lt; "(" &lt;&lt; latitude1 &lt;&lt; "," &lt;&lt; longitude1 &lt;&lt; ") --- "
            "(" &lt;&lt; latitude2 &lt;&lt; "," &lt;&lt; longitude2 &lt;&lt; ")\n";
  
    auto angle = CoordinatesToAngle(latitude1, longitude1, latitude2, longitude2);
    cout &lt;&lt; "Angle =  " &lt;&lt; angle &lt;&lt; endl;
  
    auto meters = CoordinatesToMeters(latitude1, longitude1, latitude2, longitude2);
    cout &lt;&lt; "Meters = " &lt;&lt; meters &lt;&lt; endl;
  
    auto coordinate = CoordinateToCoordinate(latitude1, longitude1, angle, meters);
    cout &lt;&lt; "Destination = (" &lt;&lt; coordinate.first &lt;&lt; "," &lt;&lt; coordinate.second &lt;&lt; ")\n";
  }
  </code></pre>