I am trying to find the destination point, given a starting point lat/long, bearing & distance. The calculator from this website below gives me the desired results.
http://www.movable-type.co.uk/scripts/latlong.html
When I try to implement the same through code, I don't get the right results.
Below is my code -
private GLatLng pointRadialDistance(double lat1, double lon1,
double radianBearing, double radialDistance)
{
double rEarth = 6371.01;
lat1 = DegreeToRadian(lat1);
lon1 = DegreeToRadian(lon1);
radianBearing = DegreeToRadian(radianBearing);
radialDistance = radialDistance / rEarth;
double lat = Math.Asin(Math.Sin(lat1) * Math.Cos(radialDistance) + Math.Cos(lat1)
* Math.Sin(radialDistance) * Math.Cos(radianBearing));
double lon;
if (Math.Cos(lat) == 0)
{ // Endpoint a pole
lon = lon1;
}
else
{
lon = ((lon1 - Math.Asin(Math.Sin(radianBearing) * Math.Sin(radialDistance) / Math.Cos(lat))
+ Math.PI) % (2 * Math.PI)) - Math.PI;
}
lat = RadianToDegree(lat);
lon = RadianToDegree(lon);
GLatLng newLatLng = new GLatLng(lat, lon);
return newLatLng;
}
public double Bearing(double lat1, double long1, double lat2, double long2)
{
//Convert input values to radians
lat1 = DegreeToRadian(lat1);
long1 = DegreeToRadian(long1);
lat2 = DegreeToRadian(lat2);
long2 = DegreeToRadian(long2);
double deltaLong = long2 - long1;
double y = Math.Sin(deltaLong) * Math.Cos(lat2);
double x = Math.Cos(lat1) * Math.Sin(lat2) -
Math.Sin(lat1) * Math.Cos(lat2) * Math.Cos(deltaLong);
double bearing = Math.Atan2(y, x);
return bearing;
}
public double DegreeToRadian(double angle)
{
return Math.PI * angle / 180.0;
}
public double RadianToDegree(double angle)
{
return 180.0 * angle / Math.PI;
}
From the main program, I call the sub procedures as follows -
double bearing = Bearing(-41.294444, 174.814444, -40.90521, 175.6604);
GLatLng endLatLng = pointRadialDistance(-41.294444, 174.814444, bearing, 80);
I get below results -
Bearing=1.02749621782165
endLatLng=-40.5751022737927,174.797458881699
The answer I expect is -40.939722,175.646389
(from website link above).
Can anyone suggest what mistake I am making in the code here?
Here's some code that achieves what you want to do.
public static GeoLocation FindPointAtDistanceFrom(GeoLocation startPoint, double initialBearingRadians, double distanceKilometres)
{
const double radiusEarthKilometres = 6371.01;
var distRatio = distanceKilometres / radiusEarthKilometres;
var distRatioSine = Math.Sin(distRatio);
var distRatioCosine = Math.Cos(distRatio);
var startLatRad = DegreesToRadians(startPoint.Latitude);
var startLonRad = DegreesToRadians(startPoint.Longitude);
var startLatCos = Math.Cos(startLatRad);
var startLatSin = Math.Sin(startLatRad);
var endLatRads = Math.Asin((startLatSin * distRatioCosine) + (startLatCos * distRatioSine * Math.Cos(initialBearingRadians)));
var endLonRads = startLonRad
+ Math.Atan2(
Math.Sin(initialBearingRadians) * distRatioSine * startLatCos,
distRatioCosine - startLatSin * Math.Sin(endLatRads));
return new GeoLocation
{
Latitude = RadiansToDegrees(endLatRads),
Longitude = RadiansToDegrees(endLonRads)
};
}
public struct GeoLocation
{
public double Latitude { get; set; }
public double Longitude { get; set; }
}
public static double DegreesToRadians(double degrees)
{
const double degToRadFactor = Math.PI / 180;
return degrees * degToRadFactor;
}
public static double RadiansToDegrees(double radians)
{
const double radToDegFactor = 180 / Math.PI;
return radians * radToDegFactor;
}
very simple solution in geometry library (V3), if you do not have a problema with using the google maps api V3 (depending on the application - realtime asset tracking, for example - the free license is not applicable OR you might not want to refactor from V2 to V3).
1st: declare an extra library ALONG with your current declaration:
<script type="text/javascript" src="http://maps.google.com/maps/api/js?libraries=geometry&sensor=false"></script>
2nd: establish starting point, heading and distance
var nyc = new google.maps.LatLng(40.715, -74.002);
var distance = 5576673;
var heading = 51.2145;
3rd: go there
var endPoint = google.maps.geometry.spherical.computeOffset(nyc, distance, heading);
var london = new google.maps.Marker({
position: endPoint,
map: map
});
done, you are now in London town. for more info regarding computeDistance, computeHeading and computeArea:
http://www.svennerberg.com/2011/04/calculating-distances-and-areas-in-google-maps-api-3/
http://code.google.com/intl/en/apis/maps/documentation/javascript/geometry.html
Below is the implementation of the JavaScript code at http://www.movable-type.co.uk/scripts/latlong.html which I wrote for myself and use in my own projects. You can implement it to your project if you will.
Note: Coordinate is a class with X (longitude), Y (latitude), Z (altitude) properties. ToDegree() and ToRadian() are extensions for Double type. Finally, GetTarget() is an extension for a Coordinate instance.
/// <summary>Calculates the destination coordinate by given angle and distance.</summary>
/// <param name="origin">Origin.</param>
/// <param name="bearing">Azimuth.</param>
/// <param name="distance">Distance (km).</param>
/// <returns>Coordinate.</returns>
public static Coordinate GetTarget(
this Coordinate origin, double bearing, double distance, double altitude = 0)
{
var d = distance / 6371;
var rlat = origin.Y.ToRadian();
var rlon = origin.X.ToRadian();
var rbearing = bearing.ToRadian();
var lat2 = rlat + (d * Math.Cos(rbearing));
var dlat = lat2 - rlat;
var dphi = Math.Log((Math.Tan((lat2 / 2) + (Math.PI / 4))) / (Math.Tan((rlat / 2) + (Math.PI / 4))));
var q =
Math.Abs(dlat) > 0.0000000001
? dlat / dphi
: Math.Cos(rlat);
var dlon = (d * Math.Sin(rbearing)) / q;
if (Math.Abs(lat2) > Math.PI / 2)
{
lat2 = lat2 > 0 ? Math.PI : Math.PI - lat2;
}
var lon2 = (rlon + dlon + Math.PI) % (2 * Math.PI) - Math.PI;
return new Coordinate
{
X = lon2.ToDegree(),
Y = lat2.ToDegree(),
Z = origin.Z
};
}