offset method
Vincenty inverse calculation
More on Wikipedia
Implementation
@override
GeoLatLng offset(final GeoLatLng from, final double distanceInMeter,
final double bearing) {
final double equatorialRadius = EQUATOR_RADIUS;
final double polarRadius = POLAR_RADIUS;
final double flattening = FLATTENING; // WGS-84 ellipsoid params
final double latitude = from.latitudeInRad;
final double longitude = from.longitudeInRad;
final double alpha1 = degToRadian(bearing);
final double sinAlpha1 = math.sin(alpha1);
final double cosAlpha1 = math.cos(alpha1);
final double tanU1 = (1 - flattening) * math.tan(latitude);
final double cosU1 = 1 / math.sqrt((1 + tanU1 * tanU1));
final double sinU1 = tanU1 * cosU1;
final double sigma1 = math.atan2(tanU1, cosAlpha1);
final double sinAlpha = cosU1 * sinAlpha1;
final double cosSqAlpha = 1 - sinAlpha * sinAlpha;
final double dfUSq = cosSqAlpha *
(equatorialRadius * equatorialRadius - polarRadius * polarRadius) /
(polarRadius * polarRadius);
final double a = 1 +
dfUSq / 16384 * (4096 + dfUSq * (-768 + dfUSq * (320 - 175 * dfUSq)));
final double b =
dfUSq / 1024 * (256 + dfUSq * (-128 + dfUSq * (74 - 47 * dfUSq)));
double sigma = distanceInMeter / (polarRadius * a);
double sigmaP = 2 * PI;
double sinSigma = 0.0;
double cosSigma = 0.0;
double cos2SigmaM = 0.0;
double deltaSigma;
int maxIterations = 200;
do {
cos2SigmaM = math.cos(2 * sigma1 + sigma);
sinSigma = math.sin(sigma);
cosSigma = math.cos(sigma);
deltaSigma = b *
sinSigma *
(cos2SigmaM +
b /
4 *
(cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
b /
6 *
cos2SigmaM *
(-3 + 4 * sinSigma * sinSigma) *
(-3 + 4 * cos2SigmaM * cos2SigmaM)));
sigmaP = sigma;
sigma = distanceInMeter / (polarRadius * a) + deltaSigma;
} while ((sigma - sigmaP).abs() > 1e-12 && --maxIterations > 0);
if (maxIterations == 0) {
throw new StateError("offset calculation faild to converge!");
}
final double tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;
final double lat2 = math.atan2(
sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
(1 - flattening) * math.sqrt(sinAlpha * sinAlpha + tmp * tmp));
final double lambda = math.atan2(
sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);
final double c =
flattening / 16 * cosSqAlpha * (4 + flattening * (4 - 3 * cosSqAlpha));
final double l = lambda -
(1 - c) *
flattening *
sinAlpha *
(sigma +
c *
sinSigma *
(cos2SigmaM +
c * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
double lon2 = longitude + l;
// print("LA ${radianToDeg(lat2)}, LO ${radianToDeg(lon2)}");
if (lon2 > PI) {
lon2 = lon2 - 2 * PI;
}
if (lon2 < -1 * PI) {
lon2 = lon2 + 2 * PI;
}
return new GeoLatLng(radianToDeg(lat2), radianToDeg(lon2));
}