Версия прошивки должна поддерживать функции монитора батареи . + Диаграммы на основе типов записей RT-PWR-BattPack и RT-PWR-BattCell + (См. руководство). Сервер будет хранить их в течении 24 час. .
+Точность измерения 5 мВ по напряжению и 1 Â ° C по температуре. + Аккумулятор Twizy состоит из 42 банок со схемой (layout) 14S3P, распределенных в 7 модулях. + Система поддерживает измерение только 14 3P напряжений и 7 измерений температуры модуля.
+ +График пакета отображает только время с данными, и исключает неактивное время парковки.
+Модуль отправляет обновления на сервер каждые ~60 секунд во время вождения или зарядки + (при наличии сети и соединения с сервером). + Экстремумы и максимальные стандартные отклонения собираются постоянно, текущие значения те, которые отправлены во время обновления.
+При включени переключателя на Twizy и начала зарядки, + передыдущие значения сбрасываются, так что в конце интервала использования они показывают общие экстремальные значения, + которые произошли во время использования.
+Tipp: Проследите линию минимального напряжения, нажмите на нее, чтобы выделить запись, + а затем используйте панель поиска, чтобы переместить маркер
. + +
+ • lower end = минимальное значение
+ • upper end = текущее значение (+текст)
+ • filled box = плохо
+ • bigger box = большее отклонение
+
Это подсветка (candle stick) диаграммы. Максимальное значение это текущее значение (volt_act / temp_act), + Минимальное значение это наименьшее измеренное значение (volt_min / temp_min) + в текущем интервале использования.
+Высота подсветки (candle height) отражает измеренное максимальное отклонение. + Заполнение подсветки (filled candle) используется для отображения «плохого» отклонения, т.е. выше напряжения пробоя или + повышение температуры выше допустимой. +
+Tipp: Для сравнения напряжений выделите банку и увеличте вертикально, + пока вы не можете четко прочитать отклонения других банок от выбранного маркера.
+ ]]>Tipp: to compare voltage details, select a cell and zoom in vertically until you can clearly read the deviations of the other cells from the selection marker.
]]> - + ++ * This allow many double precision numbers to be added together with twice the + * normal precision. Thus the effective precision of the sum is 106 bits or + * about 32 decimal places. + *
+ * The implementation follows J. R. Shewchuk, + * Adaptive Precision + * Floating-Point Arithmetic and Fast Robust Geometric Predicates, + * Discrete & Computational Geometry 18(3) 305–363 (1997). + *
+ * In the documentation of the member functions, sum stands for the value + * currently held in the accumulator. + ***********************************************************************/ +public class Accumulator { + // _s + _t accumulators for the sum. + private double _s, _t; + /** + * Construct from a double. + *
+ * @param y set sum = y. + **********************************************************************/ + public Accumulator(double y) { _s = y; _t = 0; } + /** + * Construct from another Accumulator. + *
+ * @param a set sum = a. + **********************************************************************/ + public Accumulator(Accumulator a) { _s = a._s; _t = a._t; } + /** + * Set the value to a double. + *
+ * @param y set sum = y. + **********************************************************************/ + public void Set(double y) { _s = y; _t = 0; } + /** + * Return the value held in the accumulator. + *
+ * @return sum. + **********************************************************************/ + public double Sum() { return _s; } + /** + * Return the result of adding a number to sum (but don't change + * sum). + *
+ * @param y the number to be added to the sum. + * @return sum + y. + **********************************************************************/ + public double Sum(double y) { + Accumulator a = new Accumulator(this); + a.Add(y); + return a._s; + } + /** + * Add a number to the accumulator. + *
+ * @param y set sum += y. + **********************************************************************/ + public void Add(double y) { + // Here's Shewchuk's solution... + double u; // hold exact sum as [s, t, u] + // Accumulate starting at least significant end + { Pair r = GeoMath.sum(y, _t); y = r.first; u = r.second; } + { Pair r = GeoMath.sum(y, _s); _s = r.first; _t = r.second; } + // Start is _s, _t decreasing and non-adjacent. Sum is now (s + t + u) + // exactly with s, t, u non-adjacent and in decreasing order (except for + // possible zeros). The following code tries to normalize the result. + // Ideally, we want _s = round(s+t+u) and _u = round(s+t+u - _s). The + // following does an approximate job (and maintains the decreasing + // non-adjacent property). Here are two "failures" using 3-bit floats: + // + // Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp + // [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3 + // + // Case 2: _s+_t is not as close to s+t+u as it shold be + // [64, 5] + 4 -> [64, 8, 1] -> [64, 8] = 72 (off by 1) + // should be [80, -7] = 73 (exact) + // + // "Fixing" these problems is probably not worth the expense. The + // representation inevitably leads to small errors in the accumulated + // values. The additional errors illustrated here amount to 1 ulp of the + // less significant word during each addition to the Accumulator and an + // additional possible error of 1 ulp in the reported sum. + // + // Incidentally, the "ideal" representation described above is not + // canonical, because _s = round(_s + _t) may not be true. For example, + // with 3-bit floats: + // + // [128, 16] + 1 -> [160, -16] -- 160 = round(145). + // But [160, 0] - 16 -> [128, 16] -- 128 = round(144). + // + if (_s == 0) // This implies t == 0, + _s = u; // so result is u + else + _t += u; // otherwise just accumulate u to t. + } + /** + * Negate an accumulator. + *
+ * Set sum = −sum.
+ **********************************************************************/
+ public void Negate() { _s = -_s; _t = -_t; }
+}
diff --git a/OpenVehicleApp/src/net/sf/geographiclib/Constants.java b/OpenVehicleApp/src/net/sf/geographiclib/Constants.java
new file mode 100644
index 00000000..04766784
--- /dev/null
+++ b/OpenVehicleApp/src/net/sf/geographiclib/Constants.java
@@ -0,0 +1,26 @@
+/**
+ * Implementation of the net.sf.geographiclib.Constants class
+ *
+ * Copyright (c) Charles Karney (2013)
+ * Define constants specifying the WGS84 ellipsoid.
+ ***********************************************************************/
+public class Constants {
+ /**
+ * The equatorial radius of WGS84 ellipsoid (6378137 m).
+ **********************************************************************/
+ public static final double WGS84_a = 6378137;
+ /**
+ * The flattening of WGS84 ellipsoid (1/298.257223563).
+ **********************************************************************/
+ public static final double WGS84_f = 1/298.257223563;
+
+ private Constants() {}
+}
diff --git a/OpenVehicleApp/src/net/sf/geographiclib/GeoMath.java b/OpenVehicleApp/src/net/sf/geographiclib/GeoMath.java
new file mode 100644
index 00000000..80631eb9
--- /dev/null
+++ b/OpenVehicleApp/src/net/sf/geographiclib/GeoMath.java
@@ -0,0 +1,191 @@
+/**
+ * Implementation of the net.sf.geographiclib.GeoMath class
+ *
+ * Copyright (c) Charles Karney (2013)
+ * Define mathematical functions and constants so that any version of Java
+ * can be used.
+ **********************************************************************/
+public class GeoMath {
+ /**
+ * The number of binary digits in the fraction of a double precision
+ * number (equivalent to C++'s {@code numeric_limits
+ * @param x the argument.
+ * @return x2.
+ **********************************************************************/
+ public static double sq(double x) { return x * x; }
+
+ /**
+ * The hypotenuse function avoiding underflow and overflow. In Java version
+ * 1.5 and later, Math.hypot can be used.
+ *
+ * @param x the first argument.
+ * @param y the second argument.
+ * @return sqrt(x2 + y2).
+ **********************************************************************/
+ public static double hypot(double x, double y) {
+ x = Math.abs(x); y = Math.abs(y);
+ double a = Math.max(x, y), b = Math.min(x, y) / (a != 0 ? a : 1);
+ return a * Math.sqrt(1 + b * b);
+ // For an alternative method see
+ // C. Moler and D. Morrision (1983) http://dx.doi.org/10.1147/rd.276.0577
+ // and A. A. Dubrulle (1983) http://dx.doi.org/10.1147/rd.276.0582
+ }
+
+ /**
+ * log(1 + x) accurate near x = 0. In Java version 1.5 and
+ * later, Math.log1p can be used.
+ *
+ * This is taken from D. Goldberg,
+ * What every computer
+ * scientist should know about floating-point arithmetic (1991),
+ * Theorem 4. See also, N. J. Higham, Accuracy and Stability of Numerical
+ * Algorithms, 2nd Edition (SIAM, 2002), Answer to Problem 1.5, p 528.
+ *
+ * @param x the argument.
+ * @return log(1 + x).
+ **********************************************************************/
+ public static double log1p(double x) {
+ double
+ y = 1 + x,
+ z = y - 1;
+ // Here's the explanation for this magic: y = 1 + z, exactly, and z
+ // approx x, thus log(y)/z (which is nearly constant near z = 0) returns
+ // a good approximation to the true log(1 + x)/x. The multiplication x *
+ // (log(y)/z) introduces little additional error.
+ return z == 0 ? x : x * Math.log(y) / z;
+ }
+
+ /**
+ * The inverse hyperbolic tangent function. This is defined in terms of
+ * GeoMath.log1p(x) in order to maintain accuracy near x = 0.
+ * In addition, the odd parity of the function is enforced.
+ *
+ * @param x the argument.
+ * @return atanh(x).
+ **********************************************************************/
+ public static double atanh(double x) {
+ double y = Math.abs(x); // Enforce odd parity
+ y = Math.log1p(2 * y/(1 - y))/2;
+ return x < 0 ? -y : y;
+ }
+
+ /**
+ * The cube root function. In Java version 1.5 and later, Math.cbrt can be
+ * used.
+ *
+ * @param x the argument.
+ * @return the real cube root of x.
+ **********************************************************************/
+ public static double cbrt(double x) {
+ double y = Math.pow(Math.abs(x), 1/3.0); // Return the real cube root
+ return x < 0 ? -y : y;
+ }
+
+ /**
+ * The error-free sum of two numbers.
+ *
+ * @param u the first number in the sum.
+ * @param v the second number in the sum.
+ * @return Pair(s, t) with s = round(u +
+ * v) and t = u + v - s.
+ *
+ * See D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B.
+ **********************************************************************/
+ public static Pair sum(double u, double v) {
+ double s = u + v;
+ double up = s - v;
+ double vpp = s - up;
+ up -= u;
+ vpp -= v;
+ double t = -(up + vpp);
+ // u + v = s + t
+ // = round(u + v) + t
+ return new Pair(s, t);
+ }
+
+ /**
+ * Normalize an angle (restricted input range).
+ *
+ * @param x the angle in degrees.
+ * @return the angle reduced to the range [−180°, 180°).
+ *
+ * x must lie in [−540°, 540°).
+ **********************************************************************/
+ public static double AngNormalize(double x)
+ { return x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); }
+
+ /**
+ * Normalize an arbitrary angle.
+ *
+ * @param x the angle in degrees.
+ * @return the angle reduced to the range [−180°, 180°).
+ *
+ * The range of x is unrestricted.
+ **********************************************************************/
+ public static double AngNormalize2(double x)
+ { return AngNormalize(x % 360.0); }
+
+ /**
+ * Difference of two angles reduced to [−180°, 180°]
+ *
+ * @param x the first angle in degrees.
+ * @param y the second angle in degrees.
+ * @return y − x, reduced to the range [−180°,
+ * 180°].
+ *
+ * x and y must both lie in [−180°, 180°]. The
+ * result is equivalent to computing the difference exactly, reducing it to
+ * (−180°, 180°] and rounding the result. Note that this
+ * prescription allows −180° to be returned (e.g., if x is
+ * tiny and negative and y = 180°).
+ **********************************************************************/
+ public static double AngDiff(double x, double y) {
+ double d, t;
+ { Pair r = sum(-x, y); d = r.first; t = r.second; }
+ if ((d - 180.0) + t > 0.0) // y - x > 180
+ d -= 360.0; // exact
+ else if ((d + 180.0) + t <= 0.0) // y - x <= -180
+ d += 360.0; // exact
+ return d + t;
+ }
+ /**
+ * Test for finiteness.
+ *
+ * @param x the argument.
+ * @return true if number is finite, false if NaN or infinite.
+ **********************************************************************/
+ public static boolean isfinite(double x) {
+ return Math.abs(x) <= Double.MAX_VALUE;
+ }
+
+ private GeoMath() {}
+}
diff --git a/OpenVehicleApp/src/net/sf/geographiclib/Geodesic.java b/OpenVehicleApp/src/net/sf/geographiclib/Geodesic.java
new file mode 100644
index 00000000..a3ba048a
--- /dev/null
+++ b/OpenVehicleApp/src/net/sf/geographiclib/Geodesic.java
@@ -0,0 +1,1562 @@
+/**
+ * Implementation of the net.sf.geographiclib.Geodesic class
+ *
+ * Copyright (c) Charles Karney (2013)
+ * The shortest path between two points on a ellipsoid at (lat1,
+ * lon1) and (lat2, lon2) is called the geodesic. Its
+ * length is s12 and the geodesic from point 1 to point 2 has azimuths
+ * azi1 and azi2 at the two end points. (The azimuth is the
+ * heading measured clockwise from north. azi2 is the "forward"
+ * azimuth, i.e., the heading that takes you beyond point 2 not back to point
+ * 1.)
+ *
+ * Given lat1, lon1, azi1, and s12, we can
+ * determine lat2, lon2, and azi2. This is the
+ * direct geodesic problem and its solution is given by the function
+ * {@link #Direct Direct}. (If s12 is sufficiently large that the
+ * geodesic wraps more than halfway around the earth, there will be another
+ * geodesic between the points with a smaller s12.)
+ *
+ * Given lat1, lon1, lat2, and lon2, we can
+ * determine azi1, azi2, and s12. This is the
+ * inverse geodesic problem, whose solution is given by {@link #Inverse
+ * Inverse}. Usually, the solution to the inverse problem is unique. In cases
+ * where there are multiple solutions (all with the same s12, of
+ * course), all the solutions can be easily generated once a particular
+ * solution is provided.
+ *
+ * The standard way of specifying the direct problem is the specify the
+ * distance s12 to the second point. However it is sometimes useful
+ * instead to specify the arc length a12 (in degrees) on the auxiliary
+ * sphere. This is a mathematical construct used in solving the geodesic
+ * problems. The solution of the direct problem in this form is provided by
+ * {@link #ArcDirect ArcDirect}. An arc length in excess of 180° indicates
+ * that the geodesic is not a shortest path. In addition, the arc length
+ * between an equatorial crossing and the next extremum of latitude for a
+ * geodesic is 90°.
+ *
+ * This class can also calculate several other quantities related to
+ * geodesics. These are:
+ *
+ * The quantities m12, M12, M21 which all specify the
+ * behavior of nearby geodesics obey addition rules. If points 1, 2, and 3 all
+ * lie on a single geodesic, then the following rules hold:
+ *
+ * The results of the geodesic calculations are bundled up into a {@link
+ * GeodesicData} object which includes the input parameters and all the
+ * computed results, i.e., lat1, lon1, azi1, lat2,
+ * lon2, azi2, s12, a12, m12, M12,
+ * M21, S12.
+ *
+ * The functions {@link #Direct(double, double, double, double, int) Direct},
+ * {@link #ArcDirect(double, double, double, double, int) ArcDirect}, and
+ * {@link #Inverse(double, double, double, double, int) Inverse} include an
+ * optional final argument outmask which allows you specify which
+ * results should be computed and returned. If you omit outmask, then
+ * the "standard" geodesic results are computed (latitudes, longitudes,
+ * azimuths, and distance). outmask is bitor'ed combination of {@link
+ * GeodesicMask} values. For example, if you wish just to compute the distance
+ * between two points you would call, e.g.,
+ *
+ * Additional functionality is provided by the {@link GeodesicLine} class,
+ * which allows a sequence of points along a geodesic to be computed.
+ *
+ * The shortest distance returned by the solution of the inverse problem is
+ * (obviously) uniquely defined. However, in a few special cases there are
+ * multiple azimuths which yield the same shortest distance. Here is a
+ * catalog of those cases:
+ *
+ * The calculations are accurate to better than 15 nm (15 nanometers) for the
+ * WGS84 ellipsoid. See Sec. 9 of
+ * arXiv:1102.1215v1 for
+ * details. The algorithms used by this class are based on series expansions
+ * using the flattening f as a small parameter. These are only accurate
+ * for |f| < 0.02; however reasonably accurate results will be
+ * obtained for |f| < 0.2. Here is a table of the approximate
+ * maximum error (expressed as a distance) for an ellipsoid with the same
+ * major radius as the WGS84 ellipsoid and different values of the
+ * flattening.
+ * The algorithms are described in
+ *
+ * Example of use:
+ *
+ * @param a equatorial radius (meters).
+ * @param f flattening of ellipsoid. Setting f = 0 gives a sphere.
+ * Negative f gives a prolate ellipsoid. If f > 1, set
+ * flattening to 1/f.
+ * @exception GeographicErr if a or (1 − f ) a is
+ * not positive.
+ **********************************************************************/
+ public Geodesic(double a, double f) {
+ _a = a;
+ _f = f <= 1 ? f : 1/f;
+ _f1 = 1 - _f;
+ _e2 = _f * (2 - _f);
+ _ep2 = _e2 / GeoMath.sq(_f1); // e2 / (1 - e2)
+ _n = _f / ( 2 - _f);
+ _b = _a * _f1;
+ _c2 = (GeoMath.sq(_a) + GeoMath.sq(_b) *
+ (_e2 == 0 ? 1 :
+ (_e2 > 0 ? GeoMath.atanh(Math.sqrt(_e2)) :
+ Math.atan(Math.sqrt(-_e2))) /
+ Math.sqrt(Math.abs(_e2))))/2; // authalic radius squared
+ // The sig12 threshold for "really short". Using the auxiliary sphere
+ // solution with dnm computed at (bet1 + bet2) / 2, the relative error in
+ // the azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2.
+ // (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a
+ // given f and sig12, the max error occurs for lines near the pole. If
+ // the old rule for computing dnm = (dn1 + dn2)/2 is used, then the error
+ // increases by a factor of 2.) Setting this equal to epsilon gives
+ // sig12 = etol2. Here 0.1 is a safety factor (error decreased by 100)
+ // and max(0.001, abs(f)) stops etol2 getting too large in the nearly
+ // spherical case.
+ _etol2 = 0.1 * tol2_ /
+ Math.sqrt( Math.max(0.001, Math.abs(_f)) *
+ Math.min(1.0, 1 - _f/2) / 2 );
+ if (!(GeoMath.isfinite(_a) && _a > 0))
+ throw new GeographicErr("Major radius is not positive");
+ if (!(GeoMath.isfinite(_b) && _b > 0))
+ throw new GeographicErr("Minor radius is not positive");
+ _A3x = new double[nA3x_];
+ _C3x = new double[nC3x_];
+ _C4x = new double[nC4x_];
+
+ A3coeff();
+ C3coeff();
+ C4coeff();
+ }
+
+ /**
+ * Solve the direct geodesic problem where the length of the geodesic
+ * is specified in terms of distance.
+ *
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param azi1 azimuth at point 1 (degrees).
+ * @param s12 distance between point 1 and point 2 (meters); it can be
+ * negative.
+ * @return a {@link GeodesicData} object with the following fields:
+ * lat1, lon1, azi1, lat2, lon2,
+ * azi2, s12, a12.
+ *
+ * lat1 should be in the range [−90°, 90°]; lon1
+ * and azi1 should be in the range [−540°, 540°). The
+ * values of lon2 and azi2 returned are in the range
+ * [−180°, 180°).
+ *
+ * If either point is at a pole, the azimuth is defined by keeping the
+ * longitude fixed, writing lat = ±(90° − ε),
+ * and taking the limit ε → 0+. An arc length greater that
+ * 180° signifies a geodesic which is not a shortest path. (For a
+ * prolate ellipsoid, an additional condition is necessary for a shortest
+ * path: the longitudinal extent must not exceed of 180°.)
+ **********************************************************************/
+ public GeodesicData Direct(double lat1, double lon1,
+ double azi1, double s12) {
+ return Direct(lat1, lon1, azi1, false, s12,
+ GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE |
+ GeodesicMask.AZIMUTH);
+ }
+ /**
+ * Solve the direct geodesic problem where the length of the geodesic is
+ * specified in terms of distance and with a subset of the geodesic results
+ * returned.
+ *
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param azi1 azimuth at point 1 (degrees).
+ * @param s12 distance between point 1 and point 2 (meters); it can be
+ * negative.
+ * @param outmask a bitor'ed combination of {@link GeodesicMask} values
+ * specifying which results should be returned.
+ * @return a {@link GeodesicData} object with the fields specified by
+ * outmask computed.
+ *
+ * lat1, lon1, azi1, s12, and a12 are
+ * always included in the returned result.
+ **********************************************************************/
+ public GeodesicData Direct(double lat1, double lon1,
+ double azi1, double s12, int outmask) {
+ return Direct(lat1, lon1, azi1, false, s12, outmask);
+ }
+
+ /**
+ * Solve the direct geodesic problem where the length of the geodesic
+ * is specified in terms of arc length.
+ *
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param azi1 azimuth at point 1 (degrees).
+ * @param a12 arc length between point 1 and point 2 (degrees); it can
+ * be negative.
+ * @return a {@link GeodesicData} object with the following fields:
+ * lat1, lon1, azi1, lat2, lon2,
+ * azi2, s12, a12.
+ *
+ * lat1 should be in the range [−90°, 90°]; lon1
+ * and azi1 should be in the range [−540°, 540°). The
+ * values of lon2 and azi2 returned are in the range
+ * [−180°, 180°).
+ *
+ * If either point is at a pole, the azimuth is defined by keeping the
+ * longitude fixed, writing lat = ±(90° − ε),
+ * and taking the limit ε → 0+. An arc length greater that
+ * 180° signifies a geodesic which is not a shortest path. (For a
+ * prolate ellipsoid, an additional condition is necessary for a shortest
+ * path: the longitudinal extent must not exceed of 180°.)
+ **********************************************************************/
+ public GeodesicData ArcDirect(double lat1, double lon1,
+ double azi1, double a12) {
+ return Direct(lat1, lon1, azi1, true, a12,
+ GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE |
+ GeodesicMask.AZIMUTH | GeodesicMask.DISTANCE);
+ }
+
+ /**
+ * Solve the direct geodesic problem where the length of the geodesic is
+ * specified in terms of arc length and with a subset of the geodesic results
+ * returned.
+ *
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param azi1 azimuth at point 1 (degrees).
+ * @param a12 arc length between point 1 and point 2 (degrees); it can
+ * be negative.
+ * @param outmask a bitor'ed combination of {@link GeodesicMask} values
+ * specifying which results should be returned.
+ * @return a {@link GeodesicData} object with the fields specified by
+ * outmask computed.
+ *
+ * lat1, lon1, azi1, and a12 are always included
+ * in the returned result.
+ **********************************************************************/
+ public GeodesicData ArcDirect(double lat1, double lon1,
+ double azi1, double a12, int outmask) {
+ return Direct(lat1, lon1, azi1, true, a12, outmask);
+ }
+
+ /**
+ * The general direct geodesic problem. {@link #Direct Direct} and
+ * {@link #ArcDirect ArcDirect} are defined in terms of this function.
+ *
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param azi1 azimuth at point 1 (degrees).
+ * @param arcmode boolean flag determining the meaning of the
+ * s12_a12.
+ * @param s12_a12 if arcmode is false, this is the distance between
+ * point 1 and point 2 (meters); otherwise it is the arc length between
+ * point 1 and point 2 (degrees); it can be negative.
+ * @param outmask a bitor'ed combination of {@link GeodesicMask} values
+ * specifying which results should be returned.
+ * @return a {@link GeodesicData} object with the fields specified by
+ * outmask computed.
+ *
+ * The {@link GeodesicMask} values possible for outmask are
+ *
+ * The function value a12 is always computed and returned and this
+ * equals s12_a12 is arcmode is true. If outmask
+ * includes {@link GeodesicMask#DISTANCE} and arcmode is false, then
+ * s12 = s12_a12. It is not necessary to include {@link
+ * GeodesicMask#DISTANCE_IN} in outmask; this is automatically
+ * included is arcmode is false.
+ **********************************************************************/
+ public GeodesicData Direct(double lat1, double lon1, double azi1,
+ boolean arcmode, double s12_a12, int outmask) {
+ return new GeodesicLine(this, lat1, lon1, azi1,
+ // Automatically supply DISTANCE_IN if necessary
+ outmask | (arcmode ? GeodesicMask.NONE :
+ GeodesicMask.DISTANCE_IN))
+ . // Note the dot!
+ Position(arcmode, s12_a12, outmask);
+ }
+
+ /**
+ * Solve the inverse geodesic problem.
+ *
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param lat2 latitude of point 2 (degrees).
+ * @param lon2 longitude of point 2 (degrees).
+ * @return a {@link GeodesicData} object with the following fields:
+ * lat1, lon1, azi1, lat2, lon2,
+ * azi2, s12, a12.
+ *
+ * lat1 and lat2 should be in the range [−90°,
+ * 90°]; lon1 and lon2 should be in the range
+ * [−540°, 540°). The values of azi1 and azi2
+ * returned are in the range [−180°, 180°).
+ *
+ * If either point is at a pole, the azimuth is defined by keeping the
+ * longitude fixed, writing lat = ±(90° − ε),
+ * taking the limit ε → 0+.
+ *
+ * The solution to the inverse problem is found using Newton's method. If
+ * this fails to converge (this is very unlikely in geodetic applications
+ * but does occur for very eccentric ellipsoids), then the bisection method
+ * is used to refine the solution.
+ **********************************************************************/
+ public GeodesicData Inverse(double lat1, double lon1,
+ double lat2, double lon2) {
+ return Inverse(lat1, lon1, lat2, lon2,
+ GeodesicMask.DISTANCE | GeodesicMask.AZIMUTH);
+ }
+
+ /**
+ * Solve the inverse geodesic problem with a subset of the geodesic results
+ * returned.
+ *
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param lat2 latitude of point 2 (degrees).
+ * @param lon2 longitude of point 2 (degrees).
+ * @param outmask a bitor'ed combination of {@link GeodesicMask} values
+ * specifying which results should be returned.
+ * @return a {@link GeodesicData} object with the fields specified by
+ * outmask computed.
+ *
+ * The {@link GeodesicMask} values possible for outmask are
+ *
+ * lat1, lon1, lat2, lon2, and a12 are
+ * always included in the returned result.
+ **********************************************************************/
+ public GeodesicData Inverse(double lat1, double lon1,
+ double lat2, double lon2, int outmask) {
+ outmask &= GeodesicMask.OUT_ALL;
+ GeodesicData r = new GeodesicData();
+ lon1 = GeoMath.AngNormalize(lon1);
+ lon2 = GeoMath.AngNormalize(lon2);
+ // Compute longitude difference (AngDiff does this carefully). Result is
+ // in [-180, 180] but -180 is only for west-going geodesics. 180 is for
+ // east-going and meridional geodesics.
+ double lon12 = GeoMath.AngDiff(lon1, lon2);
+ // If very close to being on the same half-meridian, then make it so.
+ lon12 = AngRound(lon12);
+ // Make longitude difference positive.
+ int lonsign = lon12 >= 0 ? 1 : -1;
+ lon12 *= lonsign;
+ // If really close to the equator, treat as on equator.
+ lat1 = AngRound(lat1);
+ lat2 = AngRound(lat2);
+ // Save input parameters post normalization
+ r.lat1 = lat1; r.lon1 = lon1; r.lat2 = lat2; r.lon2 = lon2;
+ // Swap points so that point with higher (abs) latitude is point 1
+ int swapp = Math.abs(lat1) >= Math.abs(lat2) ? 1 : -1;
+ if (swapp < 0) {
+ lonsign *= -1;
+ { double t = lat1; lat1 = lat2; lat2 = t; }
+ }
+ // Make lat1 <= 0
+ int latsign = lat1 < 0 ? 1 : -1;
+ lat1 *= latsign;
+ lat2 *= latsign;
+ // Now we have
+ //
+ // 0 <= lon12 <= 180
+ // -90 <= lat1 <= 0
+ // lat1 <= lat2 <= -lat1
+ //
+ // longsign, swapp, latsign register the transformation to bring the
+ // coordinates to this canonical form. In all cases, 1 means no change was
+ // made. We make these transformations so that there are few cases to
+ // check, e.g., on verifying quadrants in atan2. In addition, this
+ // enforces some symmetries in the results returned.
+
+ double phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x;
+ s12x = m12x = Double.NaN;
+
+ phi = lat1 * GeoMath.degree;
+ // Ensure cbet1 = +epsilon at poles
+ sbet1 = _f1 * Math.sin(phi);
+ cbet1 = lat1 == -90 ? tiny_ : Math.cos(phi);
+ { Pair p = SinCosNorm(sbet1, cbet1);
+ sbet1 = p.first; cbet1 = p.second; }
+
+ phi = lat2 * GeoMath.degree;
+ // Ensure cbet2 = +epsilon at poles
+ sbet2 = _f1 * Math.sin(phi);
+ cbet2 = Math.abs(lat2) == 90 ? tiny_ : Math.cos(phi);
+ { Pair p = SinCosNorm(sbet2, cbet2);
+ sbet2 = p.first; cbet2 = p.second; }
+
+ // If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
+ // |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is
+ // a better measure. This logic is used in assigning calp2 in Lambda12.
+ // Sometimes these quantities vanish and in that case we force bet2 = +/-
+ // bet1 exactly. An example where is is necessary is the inverse problem
+ // 48.522876735459 0 -48.52287673545898293 179.599720456223079643
+ // which failed with Visual Studio 10 (Release and Debug)
+
+ if (cbet1 < -sbet1) {
+ if (cbet2 == cbet1)
+ sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
+ } else {
+ if (Math.abs(sbet2) == -sbet1)
+ cbet2 = cbet1;
+ }
+
+ double
+ dn1 = Math.sqrt(1 + _ep2 * GeoMath.sq(sbet1)),
+ dn2 = Math.sqrt(1 + _ep2 * GeoMath.sq(sbet2));
+
+ double
+ lam12 = lon12 * GeoMath.degree,
+ slam12 = Math.abs(lon12) == 180 ? 0 : Math.sin(lam12),
+ clam12 = Math.cos(lam12); // lon12 == 90 isn't interesting
+
+ double a12, sig12, calp1, salp1, calp2, salp2;
+ a12 = sig12 = calp1 = salp1 = calp2 = salp2 = Double.NaN;
+ // index zero elements of these arrays are unused
+ double C1a[] = new double[nC1_ + 1];
+ double C2a[] = new double[nC2_ + 1];
+ double C3a[] = new double[nC3_];
+
+ boolean meridian = lat1 == -90 || slam12 == 0;
+
+ if (meridian) {
+
+ // Endpoints are on a single full meridian, so the geodesic might lie on
+ // a meridian.
+
+ calp1 = clam12; salp1 = slam12; // Head to the target longitude
+ calp2 = 1; salp2 = 0; // At the target we're heading north
+
+ double
+ // tan(bet) = tan(sig) * cos(alp)
+ ssig1 = sbet1, csig1 = calp1 * cbet1,
+ ssig2 = sbet2, csig2 = calp2 * cbet2;
+
+ // sig12 = sig2 - sig1
+ sig12 = Math.atan2(Math.max(csig1 * ssig2 - ssig1 * csig2, 0.0),
+ csig1 * csig2 + ssig1 * ssig2);
+ {
+ LengthsV v =
+ Lengths(_n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
+ cbet1, cbet2,
+ (outmask & GeodesicMask.GEODESICSCALE) != 0, C1a, C2a);
+ s12x = v.s12b; m12x = v.m12b;
+ if ((outmask & GeodesicMask.GEODESICSCALE) != 0) {
+ r.M12 = v.M12; r.M21 = v.M21;
+ }
+ }
+ // Add the check for sig12 since zero length geodesics might yield m12 <
+ // 0. Test case was
+ //
+ // echo 20.001 0 20.001 0 | Geod -i
+ //
+ // In fact, we will have sig12 > pi/2 for meridional geodesic which is
+ // not a shortest path.
+ if (sig12 < 1 || m12x >= 0) {
+ m12x *= _b;
+ s12x *= _b;
+ a12 = sig12 / GeoMath.degree;
+ } else
+ // m12 < 0, i.e., prolate and too close to anti-podal
+ meridian = false;
+ }
+
+ double omg12 = Double.NaN;
+ if (!meridian &&
+ sbet1 == 0 && // and sbet2 == 0
+ // Mimic the way Lambda12 works with calp1 = 0
+ (_f <= 0 || lam12 <= Math.PI - _f * Math.PI)) {
+
+ // Geodesic runs along equator
+ calp1 = calp2 = 0; salp1 = salp2 = 1;
+ s12x = _a * lam12;
+ sig12 = omg12 = lam12 / _f1;
+ m12x = _b * Math.sin(sig12);
+ if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
+ r.M12 = r.M21 = Math.cos(sig12);
+ a12 = lon12 / _f1;
+
+ } else if (!meridian) {
+
+ // Now point1 and point2 belong within a hemisphere bounded by a
+ // meridian and geodesic is neither meridional or equatorial.
+
+ // Figure a starting point for Newton's method
+ double dnm;
+ {
+ InverseStartV v =
+ InverseStart(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
+ lam12,
+ C1a, C2a);
+ sig12 = v.sig12;
+ salp1 = v.salp1; calp1 = v.calp1;
+ salp2 = v.salp2; calp2 = v.calp2;
+ dnm = v.dnm;
+ }
+
+ if (sig12 >= 0) {
+ // Short lines (InverseStart sets salp2, calp2, dnm)
+ s12x = sig12 * _b * dnm;
+ m12x = GeoMath.sq(dnm) * _b * Math.sin(sig12 / dnm);
+ if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
+ r.M12 = r.M21 = Math.cos(sig12 / dnm);
+ a12 = sig12 / GeoMath.degree;
+ omg12 = lam12 / (_f1 * dnm);
+ } else {
+
+ // Newton's method. This is a straightforward solution of f(alp1) =
+ // lambda12(alp1) - lam12 = 0 with one wrinkle. f(alp) has exactly one
+ // root in the interval (0, pi) and its derivative is positive at the
+ // root. Thus f(alp) is positive for alp > alp1 and negative for alp <
+ // alp1. During the course of the iteration, a range (alp1a, alp1b) is
+ // maintained which brackets the root and with each evaluation of
+ // f(alp) the range is shrunk, if possible. Newton's method is
+ // restarted whenever the derivative of f is negative (because the new
+ // value of alp1 is then further from the solution) or if the new
+ // estimate of alp1 lies outside (0,pi); in this case, the new starting
+ // guess is taken to be (alp1a + alp1b) / 2.
+ double ssig1, csig1, ssig2, csig2, eps;
+ ssig1 = csig1 = ssig2 = csig2 = eps = Double.NaN;
+ int numit = 0;
+ // Bracketing range
+ double salp1a = tiny_, calp1a = 1, salp1b = tiny_, calp1b = -1;
+ for (boolean tripn = false, tripb = false; numit < maxit2_; ++numit) {
+ // the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
+ // WGS84 and random input: mean = 2.85, sd = 0.60
+ double v, dv;
+ {
+ Lambda12V w =
+ Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
+ numit < maxit1_, C1a, C2a, C3a);
+ v = w.lam12 - lam12;
+ salp2 = w.salp2; calp2 = w.calp2;
+ sig12 = w.sig12;
+ ssig1 = w.ssig1; csig1 = w.csig1;
+ ssig2 = w.ssig2; csig2 = w.csig2;
+ eps = w.eps; omg12 = w.domg12;
+ dv = w.dlam12;
+ }
+ // 2 * tol0 is approximately 1 ulp for a number in [0, pi].
+ // Reversed test to allow escape with NaNs
+ if (tripb || !(Math.abs(v) >= (tripn ? 8 : 2) * tol0_)) break;
+ // Update bracketing values
+ if (v > 0 && (numit > maxit1_ || calp1/salp1 > calp1b/salp1b))
+ { salp1b = salp1; calp1b = calp1; }
+ else if (v < 0 && (numit > maxit1_ || calp1/salp1 < calp1a/salp1a))
+ { salp1a = salp1; calp1a = calp1; }
+ if (numit < maxit1_ && dv > 0) {
+ double
+ dalp1 = -v/dv;
+ double
+ sdalp1 = Math.sin(dalp1), cdalp1 = Math.cos(dalp1),
+ nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
+ if (nsalp1 > 0 && Math.abs(dalp1) < Math.PI) {
+ calp1 = calp1 * cdalp1 - salp1 * sdalp1;
+ salp1 = nsalp1;
+ { Pair p = SinCosNorm(salp1, calp1);
+ salp1 = p.first; calp1 = p.second; }
+ // In some regimes we don't get quadratic convergence because
+ // slope -> 0. So use convergence conditions based on epsilon
+ // instead of sqrt(epsilon).
+ tripn = Math.abs(v) <= 16 * tol0_;
+ continue;
+ }
+ }
+ // Either dv was not postive or updated value was outside legal
+ // range. Use the midpoint of the bracket as the next estimate.
+ // This mechanism is not needed for the WGS84 ellipsoid, but it does
+ // catch problems with more eccentric ellipsoids. Its efficacy is
+ // such for the WGS84 test set with the starting guess set to alp1 =
+ // 90deg:
+ // the WGS84 test set: mean = 5.21, sd = 3.93, max = 24
+ // WGS84 and random input: mean = 4.74, sd = 0.99
+ salp1 = (salp1a + salp1b)/2;
+ calp1 = (calp1a + calp1b)/2;
+ { Pair p = SinCosNorm(salp1, calp1);
+ salp1 = p.first; calp1 = p.second; }
+ tripn = false;
+ tripb = (Math.abs(salp1a - salp1) + (calp1a - calp1) < tolb_ ||
+ Math.abs(salp1 - salp1b) + (calp1 - calp1b) < tolb_);
+ }
+ {
+ LengthsV v =
+ Lengths(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
+ cbet1, cbet2,
+ (outmask & GeodesicMask.GEODESICSCALE) != 0, C1a, C2a);
+ s12x = v.s12b; m12x = v.m12b;
+ if ((outmask & GeodesicMask.GEODESICSCALE) != 0) {
+ r.M12 = v.M12; r.M21 = v.M21;
+ }
+ }
+ m12x *= _b;
+ s12x *= _b;
+ a12 = sig12 / GeoMath.degree;
+ omg12 = lam12 - omg12;
+ }
+ }
+
+ if ((outmask & GeodesicMask.DISTANCE) != 0)
+ r.s12 = 0 + s12x; // Convert -0 to 0
+
+ if ((outmask & GeodesicMask.REDUCEDLENGTH) != 0)
+ r.m12 = 0 + m12x; // Convert -0 to 0
+
+ if ((outmask & GeodesicMask.AREA) != 0) {
+ double
+ // From Lambda12: sin(alp1) * cos(bet1) = sin(alp0)
+ salp0 = salp1 * cbet1,
+ calp0 = GeoMath.hypot(calp1, salp1 * sbet1); // calp0 > 0
+ double alp12;
+ if (calp0 != 0 && salp0 != 0) {
+ double
+ // From Lambda12: tan(bet) = tan(sig) * cos(alp)
+ ssig1 = sbet1, csig1 = calp1 * cbet1,
+ ssig2 = sbet2, csig2 = calp2 * cbet2,
+ k2 = GeoMath.sq(calp0) * _ep2,
+ eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2),
+ // Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0).
+ A4 = GeoMath.sq(_a) * calp0 * salp0 * _e2;
+ { Pair p = SinCosNorm(ssig1, csig1);
+ ssig1 = p.first; csig1 = p.second; }
+ { Pair p = SinCosNorm(ssig2, csig2);
+ ssig2 = p.first; csig2 = p.second; }
+ double C4a[] = new double[nC4_];
+ C4f(eps, C4a);
+ double
+ B41 = SinCosSeries(false, ssig1, csig1, C4a),
+ B42 = SinCosSeries(false, ssig2, csig2, C4a);
+ r.S12 = A4 * (B42 - B41);
+ } else
+ // Avoid problems with indeterminate sig1, sig2 on equator
+ r.S12 = 0;
+
+ if (!meridian &&
+ omg12 < 0.75 * Math.PI && // Long difference too big
+ sbet2 - sbet1 < 1.75) { // Lat difference too big
+ // Use tan(Gamma/2) = tan(omg12/2)
+ // * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2))
+ // with tan(x/2) = sin(x)/(1+cos(x))
+ double
+ somg12 = Math.sin(omg12), domg12 = 1 + Math.cos(omg12),
+ dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
+ alp12 = 2 * Math.atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
+ domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
+ } else {
+ // alp12 = alp2 - alp1, used in atan2 so no need to normalize
+ double
+ salp12 = salp2 * calp1 - calp2 * salp1,
+ calp12 = calp2 * calp1 + salp2 * salp1;
+ // The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
+ // salp12 = -0 and alp12 = -180. However this depends on the sign
+ // being attached to 0 correctly. The following ensures the correct
+ // behavior.
+ if (salp12 == 0 && calp12 < 0) {
+ salp12 = tiny_ * calp1;
+ calp12 = -1;
+ }
+ alp12 = Math.atan2(salp12, calp12);
+ }
+ r.S12 += _c2 * alp12;
+ r.S12 *= swapp * lonsign * latsign;
+ // Convert -0 to 0
+ r.S12 += 0;
+ }
+
+ // Convert calp, salp to azimuth accounting for lonsign, swapp, latsign.
+ if (swapp < 0) {
+ { double t = salp1; salp1 = salp2; salp2 = t; }
+ { double t = calp1; calp1 = calp2; calp2 = t; }
+ if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
+ { double t = r.M12; r.M12 = r.M21; r.M21 = t; }
+ }
+
+ salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
+ salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
+
+ if ((outmask & GeodesicMask.AZIMUTH) != 0) {
+ // minus signs give range [-180, 180). 0- converts -0 to +0.
+ r.azi1 = 0 - Math.atan2(-salp1, calp1) / GeoMath.degree;
+ r.azi2 = 0 - Math.atan2(-salp2, calp2) / GeoMath.degree;
+ }
+ // Returned value in [0, 180]
+ r.a12 = a12;
+ return r;
+ }
+
+ /**
+ * Set up to compute several points on a single geodesic.
+ *
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param azi1 azimuth at point 1 (degrees).
+ * @return a {@link GeodesicLine} object.
+ *
+ * lat1 should be in the range [−90°, 90°]; lon1
+ * and azi1 should be in the range [−540°, 540°). The
+ * full set of capabilities is included.
+ *
+ * If the point is at a pole, the azimuth is defined by keeping the
+ * lon1 fixed, writing lat1 = ±(90 − ε),
+ * taking the limit ε → 0+.
+ **********************************************************************/
+ public GeodesicLine Line(double lat1, double lon1, double azi1) {
+ return Line(lat1, lon1, azi1, GeodesicMask.ALL);
+ }
+ /**
+ * Set up to compute several points on a single geodesic with a subset of the
+ * capabilities included.
+ *
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param azi1 azimuth at point 1 (degrees).
+ * @param caps bitor'ed combination of {@link GeodesicMask} values specifying
+ * the capabilities the {@link GeodesicLine} object should possess, i.e.,
+ * which quantities can be returned in calls to {@link
+ * GeodesicLine#Position GeodesicLine.Position}.
+ * @return a {@link GeodesicLine} object.
+ *
+ * The {@link GeodesicMask} values are
+ *
+ * If the point is at a pole, the azimuth is defined by keeping lon1
+ * fixed, writing lat1 = ±(90 − ε), and taking
+ * the limit ε → 0+.
+ **********************************************************************/
+ public GeodesicLine Line(double lat1, double lon1, double azi1, int caps) {
+ return new GeodesicLine(this, lat1, lon1, azi1, caps);
+ }
+
+ /**
+ * @return a the equatorial radius of the ellipsoid (meters). This is
+ * the value used in the constructor.
+ **********************************************************************/
+ public double MajorRadius() { return _a; }
+
+ /**
+ * @return f the flattening of the ellipsoid. This is the
+ * value used in the constructor.
+ **********************************************************************/
+ public double Flattening() { return _f; }
+
+ /**
+ * @return total area of ellipsoid in meters2. The area of a
+ * polygon encircling a pole can be found by adding EllipsoidArea()/2 to
+ * the sum of S12 for each side of the polygon.
+ **********************************************************************/
+ public double EllipsoidArea() { return 4 * Math.PI * _c2; }
+
+ /**
+ * A global instantiation of Geodesic with the parameters for the WGS84
+ * ellipsoid.
+ **********************************************************************/
+ public static final Geodesic WGS84 =
+ new Geodesic(Constants.WGS84_a, Constants.WGS84_f);
+
+ // This is a reformulation of the geodesic problem. The notation is as
+ // follows:
+ // - at a general point (no suffix or 1 or 2 as suffix)
+ // - phi = latitude
+ // - beta = latitude on auxiliary sphere
+ // - omega = longitude on auxiliary sphere
+ // - lambda = longitude
+ // - alpha = azimuth of great circle
+ // - sigma = arc length along great circle
+ // - s = distance
+ // - tau = scaled distance (= sigma at multiples of pi/2)
+ // - at northwards equator crossing
+ // - beta = phi = 0
+ // - omega = lambda = 0
+ // - alpha = alpha0
+ // - sigma = s = 0
+ // - a 12 suffix means a difference, e.g., s12 = s2 - s1.
+ // - s and c prefixes mean sin and cos
+
+ protected static double SinCosSeries(boolean sinp,
+ double sinx, double cosx,
+ double c[]) {
+ // Evaluate
+ // y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) :
+ // sum(c[i] * cos((2*i+1) * x), i, 0, n-1)
+ // using Clenshaw summation. N.B. c[0] is unused for sin series
+ // Approx operation count = (n + 5) mult and (2 * n + 2) add
+ int k = c.length, n = k - (sinp ? 1 : 0);
+ double
+ ar = 2 * (cosx - sinx) * (cosx + sinx), // 2 * cos(2 * x)
+ y0 = (n & 1) != 0 ? c[--k] : 0, y1 = 0; // accumulators for sum
+ // Now n is even
+ n /= 2;
+ while (n-- != 0) {
+ // Unroll loop x 2, so accumulators return to their original role
+ y1 = ar * y0 - y1 + c[--k];
+ y0 = ar * y1 - y0 + c[--k];
+ }
+ return sinp
+ ? 2 * sinx * cosx * y0 // sin(2 * x) * y0
+ : cosx * (y0 - y1); // cos(x) * (y0 - y1)
+ }
+
+ private class LengthsV {
+ private double s12b, m12b, m0, M12, M21;
+ private LengthsV() {
+ s12b = m12b = m0 = M12 = M21 = Double.NaN;
+ }
+ }
+
+ private LengthsV Lengths(double eps, double sig12,
+ double ssig1, double csig1, double dn1,
+ double ssig2, double csig2, double dn2,
+ double cbet1, double cbet2,
+ boolean scalep,
+ // Scratch areas of the right size
+ double C1a[], double C2a[]) {
+ // Return m12b = (reduced length)/_b; also calculate s12b = distance/_b,
+ // and m0 = coefficient of secular term in expression for reduced length.
+ LengthsV v = new LengthsV(); // To hold s12b, m12b, m0, M12, M21;
+ C1f(eps, C1a);
+ C2f(eps, C2a);
+ double
+ A1m1 = A1m1f(eps),
+ AB1 = (1 + A1m1) * (SinCosSeries(true, ssig2, csig2, C1a) -
+ SinCosSeries(true, ssig1, csig1, C1a)),
+ A2m1 = A2m1f(eps),
+ AB2 = (1 + A2m1) * (SinCosSeries(true, ssig2, csig2, C2a) -
+ SinCosSeries(true, ssig1, csig1, C2a));
+ v.m0 = A1m1 - A2m1;
+ double J12 = v.m0 * sig12 + (AB1 - AB2);
+ // Missing a factor of _b.
+ // Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure accurate
+ // cancellation in the case of coincident points.
+ v.m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
+ csig1 * csig2 * J12;
+ // Missing a factor of _b
+ v.s12b = (1 + A1m1) * sig12 + AB1;
+ if (scalep) {
+ double csig12 = csig1 * csig2 + ssig1 * ssig2;
+ double t = _ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
+ v.M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
+ v.M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
+ }
+ return v;
+ }
+
+ private static double Astroid(double x, double y) {
+ // Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
+ // This solution is adapted from Geocentric::Reverse.
+ double k;
+ double
+ p = GeoMath.sq(x),
+ q = GeoMath.sq(y),
+ r = (p + q - 1) / 6;
+ if ( !(q == 0 && r <= 0) ) {
+ double
+ // Avoid possible division by zero when r = 0 by multiplying equations
+ // for s and t by r^3 and r, resp.
+ S = p * q / 4, // S = r^3 * s
+ r2 = GeoMath.sq(r),
+ r3 = r * r2,
+ // The discrimant of the quadratic equation for T3. This is zero on
+ // the evolute curve p^(1/3)+q^(1/3) = 1
+ disc = S * (S + 2 * r3);
+ double u = r;
+ if (disc >= 0) {
+ double T3 = S + r3;
+ // Pick the sign on the sqrt to maximize abs(T3). This minimizes loss
+ // of precision due to cancellation. The result is unchanged because
+ // of the way the T is used in definition of u.
+ T3 += T3 < 0 ? -Math.sqrt(disc) : Math.sqrt(disc); // T3 = (r * t)^3
+ // N.B. cbrt always returns the double root. cbrt(-8) = -2.
+ double T = GeoMath.cbrt(T3); // T = r * t
+ // T can be zero; but then r2 / T -> 0.
+ u += T + (T != 0 ? r2 / T : 0);
+ } else {
+ // T is complex, but the way u is defined the result is double.
+ double ang = Math.atan2(Math.sqrt(-disc), -(S + r3));
+ // There are three possible cube roots. We choose the root which
+ // avoids cancellation. Note that disc < 0 implies that r < 0.
+ u += 2 * r * Math.cos(ang / 3);
+ }
+ double
+ v = Math.sqrt(GeoMath.sq(u) + q), // guaranteed positive
+ // Avoid loss of accuracy when u < 0.
+ uv = u < 0 ? q / (v - u) : u + v, // u+v, guaranteed positive
+ w = (uv - q) / (2 * v); // positive?
+ // Rearrange expression for k to avoid loss of accuracy due to
+ // subtraction. Division by 0 not possible because uv > 0, w >= 0.
+ k = uv / (Math.sqrt(uv + GeoMath.sq(w)) + w); // guaranteed positive
+ } else { // q == 0 && r <= 0
+ // y = 0 with |x| <= 1. Handle this case directly.
+ // for y small, positive root is k = abs(y)/sqrt(1-x^2)
+ k = 0;
+ }
+ return k;
+ }
+
+ private class InverseStartV {
+ private double sig12, salp1, calp1,
+ // Only updated if return val >= 0
+ salp2, calp2,
+ // Only updated for short lines
+ dnm;
+ private InverseStartV() {
+ sig12 = salp1 = calp1 = salp2 = calp2 = dnm = Double.NaN;
+ }
+ }
+
+ private InverseStartV InverseStart(double sbet1, double cbet1, double dn1,
+ double sbet2, double cbet2, double dn2,
+ double lam12,
+ // Scratch areas of the right size
+ double C1a[], double C2a[]) {
+ // Return a starting point for Newton's method in salp1 and calp1 (function
+ // value is -1). If Newton's method doesn't need to be used, return also
+ // salp2 and calp2 and function value is sig12.
+
+ // To hold sig12, salp1, calp1, salp2, calp2, dnm.
+ InverseStartV w = new InverseStartV();
+ w.sig12 = -1; // Return value
+ double
+ // bet12 = bet2 - bet1 in [0, pi); bet12a = bet2 + bet1 in (-pi, 0]
+ sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
+ cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
+ double sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
+ boolean shortline = cbet12 >= 0 && sbet12 < 0.5 &&
+ cbet2 * lam12 < 0.5;
+ double omg12 = lam12;
+ if (shortline) {
+ double sbetm2 = GeoMath.sq(sbet1 + sbet2);
+ // sin((bet1+bet2)/2)^2
+ // = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2)
+ sbetm2 /= sbetm2 + GeoMath.sq(cbet1 + cbet2);
+ w.dnm = Math.sqrt(1 + _ep2 * sbetm2);
+ omg12 /= _f1 * w.dnm;
+ }
+ double somg12 = Math.sin(omg12), comg12 = Math.cos(omg12);
+
+ w.salp1 = cbet2 * somg12;
+ w.calp1 = comg12 >= 0 ?
+ sbet12 + cbet2 * sbet1 * GeoMath.sq(somg12) / (1 + comg12) :
+ sbet12a - cbet2 * sbet1 * GeoMath.sq(somg12) / (1 - comg12);
+
+ double
+ ssig12 = GeoMath.hypot(w.salp1, w.calp1),
+ csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
+
+ if (shortline && ssig12 < _etol2) {
+ // really short lines
+ w.salp2 = cbet1 * somg12;
+ w.calp2 = sbet12 - cbet1 * sbet2 *
+ (comg12 >= 0 ? GeoMath.sq(somg12) / (1 + comg12) : 1 - comg12);
+ { Pair p = SinCosNorm(w.salp2, w.calp2);
+ w.salp2 = p.first; w.calp2 = p.second; }
+ // Set return value
+ w.sig12 = Math.atan2(ssig12, csig12);
+ } else if (Math.abs(_n) > 0.1 || // Skip astroid calc if too eccentric
+ csig12 >= 0 ||
+ ssig12 >= 6 * Math.abs(_n) * Math.PI * GeoMath.sq(cbet1)) {
+ // Nothing to do, zeroth order spherical approximation is OK
+ } else {
+ // Scale lam12 and bet2 to x, y coordinate system where antipodal point
+ // is at origin and singular point is at y = 0, x = -1.
+ double y, lamscale, betscale;
+ // In C++ volatile declaration needed to fix inverse case
+ // 56.320923501171 0 -56.320923501171 179.664747671772880215
+ // which otherwise fails with g++ 4.4.4 x86 -O3
+ double x;
+ if (_f >= 0) { // In fact f == 0 does not get here
+ // x = dlong, y = dlat
+ {
+ double
+ k2 = GeoMath.sq(sbet1) * _ep2,
+ eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2);
+ lamscale = _f * cbet1 * A3f(eps) * Math.PI;
+ }
+ betscale = lamscale * cbet1;
+
+ x = (lam12 - Math.PI) / lamscale;
+ y = sbet12a / betscale;
+ } else { // _f < 0
+ // x = dlat, y = dlong
+ double
+ cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
+ bet12a = Math.atan2(sbet12a, cbet12a);
+ double m12b, m0;
+ // In the case of lon12 = 180, this repeats a calculation made in
+ // Inverse.
+ LengthsV v =
+ Lengths(_n, Math.PI + bet12a,
+ sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
+ cbet1, cbet2, false, C1a, C2a);
+ m12b = v.m12b; m0 = v.m0;
+
+ x = -1 + m12b / (cbet1 * cbet2 * m0 * Math.PI);
+ betscale = x < -0.01 ? sbet12a / x :
+ -_f * GeoMath.sq(cbet1) * Math.PI;
+ lamscale = betscale / cbet1;
+ y = (lam12 - Math.PI) / lamscale;
+ }
+
+ if (y > -tol1_ && x > -1 - xthresh_) {
+ // strip near cut
+ if (_f >= 0) {
+ w.salp1 = Math.min(1.0, -x);
+ w.calp1 = - Math.sqrt(1 - GeoMath.sq(w.salp1));
+ } else {
+ w.calp1 = Math.max(x > -tol1_ ? 0.0 : -1.0, x);
+ w.salp1 = Math.sqrt(1 - GeoMath.sq(w.calp1));
+ }
+ } else {
+ // Estimate alp1, by solving the astroid problem.
+ //
+ // Could estimate alpha1 = theta + pi/2, directly, i.e.,
+ // calp1 = y/k; salp1 = -x/(1+k); for _f >= 0
+ // calp1 = x/(1+k); salp1 = -y/k; for _f < 0 (need to check)
+ //
+ // However, it's better to estimate omg12 from astroid and use
+ // spherical formula to compute alp1. This reduces the mean number of
+ // Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
+ // (min 0 max 5). The changes in the number of iterations are as
+ // follows:
+ //
+ // change percent
+ // 1 5
+ // 0 78
+ // -1 16
+ // -2 0.6
+ // -3 0.04
+ // -4 0.002
+ //
+ // The histogram of iterations is (m = number of iterations estimating
+ // alp1 directly, n = number of iterations estimating via omg12, total
+ // number of trials = 148605):
+ //
+ // iter m n
+ // 0 148 186
+ // 1 13046 13845
+ // 2 93315 102225
+ // 3 36189 32341
+ // 4 5396 7
+ // 5 455 1
+ // 6 56 0
+ //
+ // Because omg12 is near pi, estimate work with omg12a = pi - omg12
+ double k = Astroid(x, y);
+ double
+ omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
+ somg12 = Math.sin(omg12a); comg12 = -Math.cos(omg12a);
+ // Update spherical estimate of alp1 using omg12 instead of lam12
+ w.salp1 = cbet2 * somg12;
+ w.calp1 = sbet12a - cbet2 * sbet1 * GeoMath.sq(somg12) / (1 - comg12);
+ }
+ }
+ if (w.salp1 > 0) // Sanity check on starting guess
+ { Pair p = SinCosNorm(w.salp1, w.calp1);
+ w.salp1 = p.first; w.calp1 = p.second; }
+ else {
+ w.salp1 = 1; w.calp1 = 0;
+ }
+ return w;
+ }
+
+ private class Lambda12V {
+ private double lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
+ eps, domg12, dlam12;
+ private Lambda12V() {
+ lam12 = salp2 = calp2 = sig12 = ssig1 = csig1 = ssig2 = csig2
+ = eps = domg12 = dlam12 = Double.NaN;
+ }
+ }
+
+ private Lambda12V Lambda12(double sbet1, double cbet1, double dn1,
+ double sbet2, double cbet2, double dn2,
+ double salp1, double calp1,
+ boolean diffp,
+ // Scratch areas of the right size
+ double C1a[], double C2a[], double C3a[]) {
+ // Object to hold lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
+ // eps, domg12, dlam12;
+
+ Lambda12V w = new Lambda12V();
+
+ if (sbet1 == 0 && calp1 == 0)
+ // Break degeneracy of equatorial line. This case has already been
+ // handled.
+ calp1 = -tiny_;
+
+ double
+ // sin(alp1) * cos(bet1) = sin(alp0)
+ salp0 = salp1 * cbet1,
+ calp0 = GeoMath.hypot(calp1, salp1 * sbet1); // calp0 > 0
+
+ double somg1, comg1, somg2, comg2, omg12;
+ // tan(bet1) = tan(sig1) * cos(alp1)
+ // tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
+ w.ssig1 = sbet1; somg1 = salp0 * sbet1;
+ w.csig1 = comg1 = calp1 * cbet1;
+ { Pair p = SinCosNorm(w.ssig1, w.csig1);
+ w.ssig1 = p.first; w.csig1 = p.second; }
+ // SinCosNorm(somg1, comg1); -- don't need to normalize!
+
+ // Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful
+ // about this case, since this can yield singularities in the Newton
+ // iteration.
+ // sin(alp2) * cos(bet2) = sin(alp0)
+ w.salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
+ // calp2 = sqrt(1 - sq(salp2))
+ // = sqrt(sq(calp0) - sq(sbet2)) / cbet2
+ // and subst for calp0 and rearrange to give (choose positive sqrt
+ // to give alp2 in [0, pi/2]).
+ w.calp2 = cbet2 != cbet1 || Math.abs(sbet2) != -sbet1 ?
+ Math.sqrt(GeoMath.sq(calp1 * cbet1) +
+ (cbet1 < -sbet1 ?
+ (cbet2 - cbet1) * (cbet1 + cbet2) :
+ (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
+ Math.abs(calp1);
+ // tan(bet2) = tan(sig2) * cos(alp2)
+ // tan(omg2) = sin(alp0) * tan(sig2).
+ w.ssig2 = sbet2; somg2 = salp0 * sbet2;
+ w.csig2 = comg2 = w.calp2 * cbet2;
+ { Pair p = SinCosNorm(w.ssig2, w.csig2);
+ w.ssig2 = p.first; w.csig2 = p.second; }
+ // SinCosNorm(somg2, comg2); -- don't need to normalize!
+
+ // sig12 = sig2 - sig1, limit to [0, pi]
+ w.sig12 = Math.atan2(Math.max(w.csig1 * w.ssig2 - w.ssig1 * w.csig2, 0.0),
+ w.csig1 * w.csig2 + w.ssig1 * w.ssig2);
+
+ // omg12 = omg2 - omg1, limit to [0, pi]
+ omg12 = Math.atan2(Math.max(comg1 * somg2 - somg1 * comg2, 0.0),
+ comg1 * comg2 + somg1 * somg2);
+ double B312, h0;
+ double k2 = GeoMath.sq(calp0) * _ep2;
+ w.eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2);
+ C3f(w.eps, C3a);
+ B312 = (SinCosSeries(true, w.ssig2, w.csig2, C3a) -
+ SinCosSeries(true, w.ssig1, w.csig1, C3a));
+ h0 = -_f * A3f(w.eps);
+ w.domg12 = salp0 * h0 * (w.sig12 + B312);
+ w.lam12 = omg12 + w.domg12;
+
+ if (diffp) {
+ if (w.calp2 == 0)
+ w.dlam12 = - 2 * _f1 * dn1 / sbet1;
+ else {
+ double dummy;
+ LengthsV v =
+ Lengths(w.eps, w.sig12, w.ssig1, w.csig1, dn1, w.ssig2, w.csig2, dn2,
+ cbet1, cbet2, false, C1a, C2a);
+ w.dlam12 = v.m12b;
+ w.dlam12 *= _f1 / (w.calp2 * cbet2);
+ }
+ }
+
+ return w;
+ }
+
+ protected double A3f(double eps) {
+ // Evaluate sum(_A3x[k] * eps^k, k, 0, nA3x_-1) by Horner's method
+ double v = 0;
+ for (int i = nA3x_; i > 0; )
+ v = eps * v + _A3x[--i];
+ return v;
+ }
+
+ protected void C3f(double eps, double c[]) {
+ // Evaluate C3 coeffs by Horner's method
+ // Elements c[1] thru c[nC3_ - 1] are set
+ for (int j = nC3x_, k = nC3_ - 1; k > 0; ) {
+ double t = 0;
+ for (int i = nC3_ - k; i > 0; --i)
+ t = eps * t + _C3x[--j];
+ c[k--] = t;
+ }
+
+ double mult = 1;
+ for (int k = 1; k < nC3_; ) {
+ mult *= eps;
+ c[k++] *= mult;
+ }
+ }
+
+ protected void C4f(double eps, double c[]) {
+ // Evaluate C4 coeffs by Horner's method
+ // Elements c[0] thru c[nC4_ - 1] are set
+ for (int j = nC4x_, k = nC4_; k > 0; ) {
+ double t = 0;
+ for (int i = nC4_ - k + 1; i > 0; --i)
+ t = eps * t + _C4x[--j];
+ c[--k] = t;
+ }
+
+ double mult = 1;
+ for (int k = 1; k < nC4_; ) {
+ mult *= eps;
+ c[k++] *= mult;
+ }
+ }
+
+ // Generated by Maxima on 2010-09-04 10:26:17-04:00
+
+ // The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
+ protected static double A1m1f(double eps) {
+ double
+ eps2 = GeoMath.sq(eps),
+ t;
+ t = eps2*(eps2*(eps2+4)+64)/256;
+ return (t + eps) / (1 - eps);
+ }
+
+ // The coefficients C1[l] in the Fourier expansion of B1
+ protected static void C1f(double eps, double c[]) {
+ double
+ eps2 = GeoMath.sq(eps),
+ d = eps;
+ c[1] = d*((6-eps2)*eps2-16)/32;
+ d *= eps;
+ c[2] = d*((64-9*eps2)*eps2-128)/2048;
+ d *= eps;
+ c[3] = d*(9*eps2-16)/768;
+ d *= eps;
+ c[4] = d*(3*eps2-5)/512;
+ d *= eps;
+ c[5] = -7*d/1280;
+ d *= eps;
+ c[6] = -7*d/2048;
+ }
+
+ // The coefficients C1p[l] in the Fourier expansion of B1p
+ protected static void C1pf(double eps, double c[]) {
+ double
+ eps2 = GeoMath.sq(eps),
+ d = eps;
+ c[1] = d*(eps2*(205*eps2-432)+768)/1536;
+ d *= eps;
+ c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288;
+ d *= eps;
+ c[3] = d*(116-225*eps2)/384;
+ d *= eps;
+ c[4] = d*(2695-7173*eps2)/7680;
+ d *= eps;
+ c[5] = 3467*d/7680;
+ d *= eps;
+ c[6] = 38081*d/61440;
+ }
+
+ // The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
+ protected static double A2m1f(double eps) {
+ double
+ eps2 = GeoMath.sq(eps),
+ t;
+ t = eps2*(eps2*(25*eps2+36)+64)/256;
+ return t * (1 - eps) - eps;
+ }
+
+ // The coefficients C2[l] in the Fourier expansion of B2
+ protected static void C2f(double eps, double c[]) {
+ double
+ eps2 = GeoMath.sq(eps),
+ d = eps;
+ c[1] = d*(eps2*(eps2+2)+16)/32;
+ d *= eps;
+ c[2] = d*(eps2*(35*eps2+64)+384)/2048;
+ d *= eps;
+ c[3] = d*(15*eps2+80)/768;
+ d *= eps;
+ c[4] = d*(7*eps2+35)/512;
+ d *= eps;
+ c[5] = 63*d/1280;
+ d *= eps;
+ c[6] = 77*d/2048;
+ }
+
+ // The scale factor A3 = mean value of (d/dsigma)I3
+ protected void A3coeff() {
+ _A3x[0] = 1;
+ _A3x[1] = (_n-1)/2;
+ _A3x[2] = (_n*(3*_n-1)-2)/8;
+ _A3x[3] = ((-_n-3)*_n-1)/16;
+ _A3x[4] = (-2*_n-3)/64;
+ _A3x[5] = -3/128.0;
+ }
+
+ // The coefficients C3[l] in the Fourier expansion of B3
+ protected void C3coeff() {
+ _C3x[0] = (1-_n)/4;
+ _C3x[1] = (1-_n*_n)/8;
+ _C3x[2] = ((3-_n)*_n+3)/64;
+ _C3x[3] = (2*_n+5)/128;
+ _C3x[4] = 3/128.0;
+ _C3x[5] = ((_n-3)*_n+2)/32;
+ _C3x[6] = ((-3*_n-2)*_n+3)/64;
+ _C3x[7] = (_n+3)/128;
+ _C3x[8] = 5/256.0;
+ _C3x[9] = (_n*(5*_n-9)+5)/192;
+ _C3x[10] = (9-10*_n)/384;
+ _C3x[11] = 7/512.0;
+ _C3x[12] = (7-14*_n)/512;
+ _C3x[13] = 7/512.0;
+ _C3x[14] = 21/2560.0;
+ }
+
+ // Generated by Maxima on 2012-10-19 08:02:34-04:00
+
+ // The coefficients C4[l] in the Fourier expansion of I4
+ protected void C4coeff() {
+ _C4x[0] = (_n*(_n*(_n*(_n*(100*_n+208)+572)+3432)-12012)+30030)/45045;
+ _C4x[1] = (_n*(_n*(_n*(64*_n+624)-4576)+6864)-3003)/15015;
+ _C4x[2] = (_n*((14144-10656*_n)*_n-4576)-858)/45045;
+ _C4x[3] = ((-224*_n-4784)*_n+1573)/45045;
+ _C4x[4] = (1088*_n+156)/45045;
+ _C4x[5] = 97/15015.0;
+ _C4x[6] = (_n*(_n*((-64*_n-624)*_n+4576)-6864)+3003)/135135;
+ _C4x[7] = (_n*(_n*(5952*_n-11648)+9152)-2574)/135135;
+ _C4x[8] = (_n*(5792*_n+1040)-1287)/135135;
+ _C4x[9] = (468-2944*_n)/135135;
+ _C4x[10] = 1/9009.0;
+ _C4x[11] = (_n*((4160-1440*_n)*_n-4576)+1716)/225225;
+ _C4x[12] = ((4992-8448*_n)*_n-1144)/225225;
+ _C4x[13] = (1856*_n-936)/225225;
+ _C4x[14] = 8/10725.0;
+ _C4x[15] = (_n*(3584*_n-3328)+1144)/315315;
+ _C4x[16] = (1024*_n-208)/105105;
+ _C4x[17] = -136/63063.0;
+ _C4x[18] = (832-2560*_n)/405405;
+ _C4x[19] = -128/135135.0;
+ _C4x[20] = 128/99099.0;
+ }
+}
diff --git a/OpenVehicleApp/src/net/sf/geographiclib/GeodesicData.java b/OpenVehicleApp/src/net/sf/geographiclib/GeodesicData.java
new file mode 100644
index 00000000..1a667a5e
--- /dev/null
+++ b/OpenVehicleApp/src/net/sf/geographiclib/GeodesicData.java
@@ -0,0 +1,78 @@
+/**
+ * Implementation of the net.sf.geographiclib.GeodesicData class
+ *
+ * Copyright (c) Charles Karney (2013)
+ * GeodesicLine facilitates the determination of a series of points on a single
+ * geodesic. The starting point (lat1, lon1) and the azimuth
+ * azi1 are specified in the constructor. {@link #Position Position}
+ * returns the location of point 2 a distance s12 along the geodesic.
+ * Alternatively {@link #ArcPosition ArcPosition} gives the position of point 2
+ * an arc length a12 along the geodesic.
+ *
+ * The calculations are accurate to better than 15 nm (15 nanometers). See
+ * Sec. 9 of
+ * arXiv:1102.1215v1 for
+ * details. The algorithms used by this class are based on series expansions
+ * using the flattening f as a small parameter. These are only accurate
+ * for |f| < 0.02; however reasonably accurate results will be
+ * obtained for |f| < 0.2.
+ *
+ * The algorithms are described in
+ *
+ * Here's an example of using this class
+ *
+ * @param g A {@link Geodesic} object used to compute the necessary
+ * information about the GeodesicLine.
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param azi1 azimuth at point 1 (degrees).
+ *
+ * lat1 should be in the range [−90°, 90°]; lon1
+ * and azi1 should be in the range [−540°, 540°).
+ *
+ * If the point is at a pole, the azimuth is defined by keeping lon1
+ * fixed, writing lat1 = ±(90° − ε), and
+ * taking the limit ε → 0+.
+ **********************************************************************/
+ public GeodesicLine(Geodesic g,
+ double lat1, double lon1, double azi1) {
+ this(g, lat1, lon1, azi1, GeodesicMask.ALL);
+ }
+
+ /**
+ * Constructor for a geodesic line staring at latitude lat1, longitude
+ * lon1, and azimuth azi1 (all in degrees) with a subset of the
+ * capabilities included.
+ *
+ * @param g A {@link Geodesic} object used to compute the necessary
+ * information about the GeodesicLine.
+ * @param lat1 latitude of point 1 (degrees).
+ * @param lon1 longitude of point 1 (degrees).
+ * @param azi1 azimuth at point 1 (degrees).
+ * @param caps bitor'ed combination of {@link GeodesicMask} values
+ * specifying the capabilities the GeodesicLine object should possess,
+ * i.e., which quantities can be returned in calls to {@link #Position
+ * Position}.
+ *
+ * The {@link GeodesicMask} values are
+ *
+ * @param s12 distance between point 1 and point 2 (meters); it can be
+ * negative.
+ * @return a {@link GeodesicData} object with the following fields:
+ * lat1, lon1, azi1, lat2, lon2,
+ * azi2, s12, a12. Some of these results may be
+ * missing if the GeodesicLine did not include the relevant capability.
+ *
+ * The values of lon2 and azi2 returned are in the range
+ * [−180°, 180°).
+ *
+ * The GeodesicLine object must have been constructed with caps
+ * |= {@link GeodesicMask#DISTANCE_IN}; otherwise no parameters are set.
+ **********************************************************************/
+ public GeodesicData Position(double s12) {
+ return Position(false, s12,
+ GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE |
+ GeodesicMask.AZIMUTH);
+ }
+ /**
+ * Compute the position of point 2 which is a distance s12 (meters)
+ * from point 1 and with a subset of the geodesic results returned.
+ *
+ * @param s12 distance between point 1 and point 2 (meters); it can be
+ * negative.
+ * @param outmask a bitor'ed combination of {@link GeodesicMask} values
+ * specifying which results should be returned.
+ * @return a {@link GeodesicData} object including the requested results.
+ *
+ * The GeodesicLine object must have been constructed with caps
+ * |= {@link GeodesicMask#DISTANCE_IN}; otherwise no parameters are set.
+ * Requesting a value which the GeodesicLine object is not capable of
+ * computing is not an error (no parameters will be set).
+ **********************************************************************/
+ public GeodesicData Position(double s12, int outmask) {
+ return Position(false, s12, outmask);
+ }
+
+ /**
+ * Compute the position of point 2 which is an arc length a12
+ * (degrees) from point 1.
+ *
+ * @param a12 arc length between point 1 and point 2 (degrees); it can
+ * be negative.
+ * @return a {@link GeodesicData} object with the following fields:
+ * lat1, lon1, azi1, lat2, lon2,
+ * azi2, s12, a12. Some of these results may be
+ * missing if the GeodesicLine did not include the relevant capability.
+ *
+ * The values of lon2 and azi2 returned are in the range
+ * [−180°, 180°).
+ *
+ * The GeodesicLine object must have been constructed with caps
+ * |= {@link GeodesicMask#DISTANCE_IN}; otherwise no parameters are set.
+ **********************************************************************/
+ public GeodesicData ArcPosition(double a12) {
+ return Position(true, a12,
+ GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE |
+ GeodesicMask.AZIMUTH | GeodesicMask.DISTANCE);
+ }
+ /**
+ * Compute the position of point 2 which is an arc length a12
+ * (degrees) from point 1 and with a subset of the geodesic results returned.
+ *
+ * @param a12 arc length between point 1 and point 2 (degrees); it can
+ * be negative.
+ * @param outmask a bitor'ed combination of {@link GeodesicMask} values
+ * specifying which results should be returned.
+ * @return a {@link GeodesicData} object giving lat1, lon2,
+ * azi2, and a12.
+ *
+ * The GeodesicLine object must have been constructed with caps
+ * |= {@link GeodesicMask#DISTANCE_IN}; otherwise no parameters are set.
+ * Requesting a value which the GeodesicLine object is not capable of
+ * computing is not an error (no parameters will be set).
+ **********************************************************************/
+ public GeodesicData ArcPosition(double a12, int outmask) {
+ return Position(true, a12, outmask);
+ }
+
+ /**
+ * The general position function. {@link #Position(double, int) Position}
+ * and {@link #ArcPosition(double, int) ArcPosition} are defined in terms of
+ * this function.
+ *
+ * @param arcmode boolean flag determining the meaning of the second
+ * parameter; if arcmode is false, then the GeodesicLine object must have
+ * been constructed with caps |= {@link GeodesicMask#DISTANCE_IN}.
+ * @param s12_a12 if arcmode is false, this is the distance between
+ * point 1 and point 2 (meters); otherwise it is the arc length between
+ * point 1 and point 2 (degrees); it can be negative.
+ * @param outmask a bitor'ed combination of {@link GeodesicMask} values
+ * specifying which results should be returned.
+ * @return a {@link GeodesicData} object with the requested results.
+ *
+ * The {@link GeodesicMask} values possible for outmask are
+ *
+ * Requesting a value which the GeodesicLine object is not capable of
+ * computing is not an error; Double.NaN is returned instead.
+ **********************************************************************/
+ public GeodesicData Position(boolean arcmode, double s12_a12,
+ int outmask) {
+ outmask &= _caps & GeodesicMask.OUT_ALL;
+ GeodesicData r = new GeodesicData();
+ if (!( Init() &&
+ (arcmode ||
+ (_caps & GeodesicMask.DISTANCE_IN & GeodesicMask.OUT_ALL) != 0) ))
+ // Uninitialized or impossible distance calculation requested
+ return r;
+ r.lat1 = _lat1; r.lon1 = _lon1; r.azi1 = _azi1;
+
+ // Avoid warning about uninitialized B12.
+ double sig12, ssig12, csig12, B12 = 0, AB1 = 0;
+ if (arcmode) {
+ // Interpret s12_a12 as spherical arc length
+ r.a12 = s12_a12;
+ sig12 = s12_a12 * GeoMath.degree;
+ double s12a = Math.abs(s12_a12);
+ s12a -= 180 * Math.floor(s12a / 180);
+ ssig12 = s12a == 0 ? 0 : Math.sin(sig12);
+ csig12 = s12a == 90 ? 0 : Math.cos(sig12);
+ } else {
+ // Interpret s12_a12 as distance
+ r.s12 = s12_a12;
+ double
+ tau12 = s12_a12 / (_b * (1 + _A1m1)),
+ s = Math.sin(tau12),
+ c = Math.cos(tau12);
+ // tau2 = tau1 + tau12
+ B12 = - Geodesic.SinCosSeries(true,
+ _stau1 * c + _ctau1 * s,
+ _ctau1 * c - _stau1 * s,
+ _C1pa);
+ sig12 = tau12 - (B12 - _B11);
+ r.a12 = sig12 / GeoMath.degree;
+ ssig12 = Math.sin(sig12); csig12 = Math.cos(sig12);
+ if (Math.abs(_f) > 0.01) {
+ // Reverted distance series is inaccurate for |f| > 1/100, so correct
+ // sig12 with 1 Newton iteration. The following table shows the
+ // approximate maximum error for a = WGS_a() and various f relative to
+ // GeodesicExact.
+ // erri = the error in the inverse solution (nm)
+ // errd = the error in the direct solution (series only) (nm)
+ // errda = the error in the direct solution (series + 1 Newton) (nm)
+ //
+ // f erri errd errda
+ // -1/5 12e6 1.2e9 69e6
+ // -1/10 123e3 12e6 765e3
+ // -1/20 1110 108e3 7155
+ // -1/50 18.63 200.9 27.12
+ // -1/100 18.63 23.78 23.37
+ // -1/150 18.63 21.05 20.26
+ // 1/150 22.35 24.73 25.83
+ // 1/100 22.35 25.03 25.31
+ // 1/50 29.80 231.9 30.44
+ // 1/20 5376 146e3 10e3
+ // 1/10 829e3 22e6 1.5e6
+ // 1/5 157e6 3.8e9 280e6
+ double
+ ssig2 = _ssig1 * csig12 + _csig1 * ssig12,
+ csig2 = _csig1 * csig12 - _ssig1 * ssig12;
+ B12 = Geodesic.SinCosSeries(true, ssig2, csig2, _C1a);
+ double serr = (1 + _A1m1) * (sig12 + (B12 - _B11)) - s12_a12 / _b;
+ sig12 = sig12 - serr / Math.sqrt(1 + _k2 * GeoMath.sq(ssig2));
+ ssig12 = Math.sin(sig12); csig12 = Math.cos(sig12);
+ // Update B12 below
+ }
+ }
+
+ double omg12, lam12, lon12;
+ double ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2;
+ // sig2 = sig1 + sig12
+ ssig2 = _ssig1 * csig12 + _csig1 * ssig12;
+ csig2 = _csig1 * csig12 - _ssig1 * ssig12;
+ double dn2 = Math.sqrt(1 + _k2 * GeoMath.sq(ssig2));
+ if ((outmask & (GeodesicMask.DISTANCE | GeodesicMask.REDUCEDLENGTH |
+ GeodesicMask.GEODESICSCALE)) != 0) {
+ if (arcmode || Math.abs(_f) > 0.01)
+ B12 = Geodesic.SinCosSeries(true, ssig2, csig2, _C1a);
+ AB1 = (1 + _A1m1) * (B12 - _B11);
+ }
+ // sin(bet2) = cos(alp0) * sin(sig2)
+ sbet2 = _calp0 * ssig2;
+ // Alt: cbet2 = hypot(csig2, salp0 * ssig2);
+ cbet2 = GeoMath.hypot(_salp0, _calp0 * csig2);
+ if (cbet2 == 0)
+ // I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case
+ cbet2 = csig2 = Geodesic.tiny_;
+ // tan(omg2) = sin(alp0) * tan(sig2)
+ somg2 = _salp0 * ssig2; comg2 = csig2; // No need to normalize
+ // tan(alp0) = cos(sig2)*tan(alp2)
+ salp2 = _salp0; calp2 = _calp0 * csig2; // No need to normalize
+ // omg12 = omg2 - omg1
+ omg12 = Math.atan2(somg2 * _comg1 - comg2 * _somg1,
+ comg2 * _comg1 + somg2 * _somg1);
+
+ if ((outmask & GeodesicMask.DISTANCE) != 0 && arcmode)
+ r.s12 = _b * ((1 + _A1m1) * sig12 + AB1);
+
+ if ((outmask & GeodesicMask.LONGITUDE) != 0) {
+ lam12 = omg12 + _A3c *
+ ( sig12 + (Geodesic.SinCosSeries(true, ssig2, csig2, _C3a)
+ - _B31));
+ lon12 = lam12 / GeoMath.degree;
+ // Use GeoMath.AngNormalize2 because longitude might have wrapped
+ // multiple times.
+ lon12 = GeoMath.AngNormalize2(lon12);
+ r.lon2 = GeoMath.AngNormalize(_lon1 + lon12);
+ }
+
+ if ((outmask & GeodesicMask.LATITUDE) != 0)
+ r.lat2 = Math.atan2(sbet2, _f1 * cbet2) / GeoMath.degree;
+
+ if ((outmask & GeodesicMask.AZIMUTH) != 0)
+ // minus signs give range [-180, 180). 0- converts -0 to +0.
+ r.azi2 = 0 - Math.atan2(-salp2, calp2) / GeoMath.degree;
+
+ if ((outmask &
+ (GeodesicMask.REDUCEDLENGTH | GeodesicMask.GEODESICSCALE)) != 0) {
+ double
+ B22 = Geodesic.SinCosSeries(true, ssig2, csig2, _C2a),
+ AB2 = (1 + _A2m1) * (B22 - _B21),
+ J12 = (_A1m1 - _A2m1) * sig12 + (AB1 - AB2);
+ if ((outmask & GeodesicMask.REDUCEDLENGTH) != 0)
+ // Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure
+ // accurate cancellation in the case of coincident points.
+ r.m12 = _b * ((dn2 * (_csig1 * ssig2) - _dn1 * (_ssig1 * csig2))
+ - _csig1 * csig2 * J12);
+ if ((outmask & GeodesicMask.GEODESICSCALE) != 0) {
+ double t = _k2 * (ssig2 - _ssig1) * (ssig2 + _ssig1) / (_dn1 + dn2);
+ r.M12 = csig12 + (t * ssig2 - csig2 * J12) * _ssig1 / _dn1;
+ r.M21 = csig12 - (t * _ssig1 - _csig1 * J12) * ssig2 / dn2;
+ }
+ }
+
+ if ((outmask & GeodesicMask.AREA) != 0) {
+ double
+ B42 = Geodesic.SinCosSeries(false, ssig2, csig2, _C4a);
+ double salp12, calp12;
+ if (_calp0 == 0 || _salp0 == 0) {
+ // alp12 = alp2 - alp1, used in atan2 so no need to normalized
+ salp12 = salp2 * _calp1 - calp2 * _salp1;
+ calp12 = calp2 * _calp1 + salp2 * _salp1;
+ // The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
+ // salp12 = -0 and alp12 = -180. However this depends on the sign
+ // being attached to 0 correctly. The following ensures the correct
+ // behavior.
+ if (salp12 == 0 && calp12 < 0) {
+ salp12 = Geodesic.tiny_ * _calp1;
+ calp12 = -1;
+ }
+ } else {
+ // tan(alp) = tan(alp0) * sec(sig)
+ // tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
+ // = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
+ // If csig12 > 0, write
+ // csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
+ // else
+ // csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
+ // No need to normalize
+ salp12 = _calp0 * _salp0 *
+ (csig12 <= 0 ? _csig1 * (1 - csig12) + ssig12 * _ssig1 :
+ ssig12 * (_csig1 * ssig12 / (1 + csig12) + _ssig1));
+ calp12 = GeoMath.sq(_salp0) + GeoMath.sq(_calp0) * _csig1 * csig2;
+ }
+ r.S12 = _c2 * Math.atan2(salp12, calp12) + _A4 * (B42 - _B41);
+ }
+
+ return r;
+ }
+
+ /**
+ * @return true if the object has been initialized.
+ **********************************************************************/
+ private boolean Init() { return _caps != 0; }
+
+ /**
+ * @return lat1 the latitude of point 1 (degrees).
+ **********************************************************************/
+ public double Latitude()
+ { return Init() ? _lat1 : Double.NaN; }
+
+ /**
+ * @return lon1 the longitude of point 1 (degrees).
+ **********************************************************************/
+ public double Longitude()
+ { return Init() ? _lon1 : Double.NaN; }
+
+ /**
+ * @return azi1 the azimuth (degrees) of the geodesic line at point 1.
+ **********************************************************************/
+ public double Azimuth()
+ { return Init() ? _azi1 : Double.NaN; }
+
+ /**
+ * @return azi0 the azimuth (degrees) of the geodesic line as it
+ * crosses the equator in a northward direction.
+ **********************************************************************/
+ public double EquatorialAzimuth() {
+ return Init() ?
+ Math.atan2(_salp0, _calp0) / GeoMath.degree : Double.NaN;
+ }
+
+ /**
+ * @return a1 the arc length (degrees) between the northward
+ * equatorial crossing and point 1.
+ **********************************************************************/
+ public double EquatorialArc() {
+ return Init() ?
+ Math.atan2(_ssig1, _csig1) / GeoMath.degree : Double.NaN;
+ }
+
+ /**
+ * @return a the equatorial radius of the ellipsoid (meters). This is
+ * the value inherited from the Geodesic object used in the constructor.
+ **********************************************************************/
+ public double MajorRadius()
+ { return Init() ? _a : Double.NaN; }
+
+ /**
+ * @return f the flattening of the ellipsoid. This is the value
+ * inherited from the Geodesic object used in the constructor.
+ **********************************************************************/
+ public double Flattening()
+ { return Init() ? _f : Double.NaN; }
+
+ /**
+ * @return caps the computational capabilities that this object was
+ * constructed with. LATITUDE and AZIMUTH are always included.
+ **********************************************************************/
+ public int Capabilities() { return _caps; }
+
+ /**
+ * @param testcaps a set of bitor'ed {@link GeodesicMask} values.
+ * @return true if the GeodesicLine object has all these capabilities.
+ **********************************************************************/
+ public boolean Capabilities(int testcaps) {
+ testcaps &= GeodesicMask.OUT_ALL;
+ return (_caps & testcaps) == testcaps;
+ }
+
+}
diff --git a/OpenVehicleApp/src/net/sf/geographiclib/GeodesicMask.java b/OpenVehicleApp/src/net/sf/geographiclib/GeodesicMask.java
new file mode 100644
index 00000000..6b72ea79
--- /dev/null
+++ b/OpenVehicleApp/src/net/sf/geographiclib/GeodesicMask.java
@@ -0,0 +1,77 @@
+/**
+ * Implementation of the net.sf.geographiclib.GeodesicMask class
+ *
+ * Copyright (c) Charles Karney (2013)
+ * These masks do double duty. They specify (via the outmask parameter)
+ * which results to return in the {@link GeodesicData} object returned by the
+ * general routines {@link Geodesic#Direct(double, double, double, double, int)
+ * Geodesic.Direct} and {@link Geodesic#Inverse(double, double, double, double,
+ * int) Geodesic.Inverse} routines. They also signify (via the caps
+ * parameter) to the {@link GeodesicLine#GeodesicLine(Geodesic, double, double,
+ * double, int) GeodesicLine.GeodesicLine} constructor and to {@link
+ * Geodesic#Line(double, double, double, int) Geodesic.Line} what capabilities
+ * should be included in the {@link GeodesicLine} object.
+ **********************************************************************/
+public class GeodesicMask {
+ protected static final int CAP_NONE = 0;
+ protected static final int CAP_C1 = 1<<0;
+ protected static final int CAP_C1p = 1<<1;
+ protected static final int CAP_C2 = 1<<2;
+ protected static final int CAP_C3 = 1<<3;
+ protected static final int CAP_C4 = 1<<4;
+ protected static final int CAP_ALL = 0x1F;
+ protected static final int OUT_ALL = 0x7F80;
+
+ /**
+ * No capabilities, no output.
+ **********************************************************************/
+ public static final int NONE = 0;
+ /**
+ * Calculate latitude lat2. (It's not necessary to include this as a
+ * capability to {@link GeodesicLine} because this is included by default.)
+ **********************************************************************/
+ public static final int LATITUDE = 1<<7 | CAP_NONE;
+ /**
+ * Calculate longitude lon2.
+ **********************************************************************/
+ public static final int LONGITUDE = 1<<8 | CAP_C3;
+ /**
+ * Calculate azimuths azi1 and azi2. (It's not necessary to
+ * include this as a capability to {@link GeodesicLine} because this is
+ * included by default.)
+ **********************************************************************/
+ public static final int AZIMUTH = 1<<9 | CAP_NONE;
+ /**
+ * Calculate distance s12.
+ **********************************************************************/
+ public static final int DISTANCE = 1<<10 | CAP_C1;
+ /**
+ * Allow distance s12 to be used as input in the direct
+ * geodesic problem.
+ **********************************************************************/
+ public static final int DISTANCE_IN = 1<<11 | CAP_C1 | CAP_C1p;
+ /**
+ * Calculate reduced length m12.
+ **********************************************************************/
+ public static final int REDUCEDLENGTH = 1<<12 | CAP_C1 | CAP_C2;
+ /**
+ * Calculate geodesic scales M12 and M21.
+ **********************************************************************/
+ public static final int GEODESICSCALE = 1<<13 | CAP_C1 | CAP_C2;
+ /**
+ * Calculate area S12.
+ **********************************************************************/
+ public static final int AREA = 1<<14 | CAP_C4;
+ /**
+ * All capabilities, calculate everything.
+ **********************************************************************/
+ public static final int ALL = OUT_ALL| CAP_ALL;
+}
diff --git a/OpenVehicleApp/src/net/sf/geographiclib/GeographicErr.java b/OpenVehicleApp/src/net/sf/geographiclib/GeographicErr.java
new file mode 100644
index 00000000..d437a000
--- /dev/null
+++ b/OpenVehicleApp/src/net/sf/geographiclib/GeographicErr.java
@@ -0,0 +1,24 @@
+/**
+ * Implementation of the net.sf.geographiclib.GeographicErr class
+ *
+ * Copyright (c) Charles Karney (2013)
+ * A class to handle exceptions. It's derived from RuntimeException so it
+ * can be caught by the usual catch clauses.
+ **********************************************************************/
+public class GeographicErr extends RuntimeException {
+ /**
+ * Constructor
+ *
+ * @param msg a string message, which is accessible in the catch
+ * clause via getMessage().
+ **********************************************************************/
+ public GeographicErr(String msg) { super(msg); }
+}
diff --git a/OpenVehicleApp/src/net/sf/geographiclib/Pair.java b/OpenVehicleApp/src/net/sf/geographiclib/Pair.java
new file mode 100644
index 00000000..1e1140b4
--- /dev/null
+++ b/OpenVehicleApp/src/net/sf/geographiclib/Pair.java
@@ -0,0 +1,32 @@
+/**
+ * Implementation of the net.sf.geographiclib.Pair class
+ *
+ * Copyright (c) Charles Karney (2013)
+ * This duplicates the C++ class {@code std::pair
+ * @param first the first member of the pair.
+ * @param second the second member of the pair.
+ **********************************************************************/
+ public Pair(double first, double second)
+ { this.first = first; this.second = second; }
+}
diff --git a/OpenVehicleApp/src/net/sf/geographiclib/PolygonArea.java b/OpenVehicleApp/src/net/sf/geographiclib/PolygonArea.java
new file mode 100644
index 00000000..75cb09a4
--- /dev/null
+++ b/OpenVehicleApp/src/net/sf/geographiclib/PolygonArea.java
@@ -0,0 +1,383 @@
+/**
+ * Implementation of the net.sf.geographiclib.PolygonArea class
+ *
+ * Copyright (c) Charles Karney (2013)
+ * This computes the area of a geodesic polygon using the method given
+ * Section 6 of
+ *
+ * This class lets you add vertices one at a time to the polygon. The area
+ * and perimeter are accumulated in two times the standard floating point
+ * precision to guard against the loss of accuracy with many-sided polygons.
+ * At any point you can ask for the perimeter and area so far. There's an
+ * option to treat the points as defining a polyline instead of a polygon; in
+ * that case, only the perimeter is computed.
+ *
+ * Example of use:
+ *
+ * @param earth the Geodesic object to use for geodesic calculations.
+ * @param polyline if true that treat the points as defining a polyline
+ * instead of a polygon.
+ **********************************************************************/
+ public PolygonArea(Geodesic earth, boolean polyline) {
+ _earth = earth;
+ _area0 = _earth.EllipsoidArea();
+ _polyline = polyline;
+ _mask = GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE |
+ GeodesicMask.DISTANCE |
+ (_polyline ? GeodesicMask.NONE : GeodesicMask.AREA);
+ _perimetersum = new Accumulator(0);
+ if (!_polyline)
+ _areasum = new Accumulator(0);
+ Clear();
+ }
+
+ /**
+ * Clear PolygonArea, allowing a new polygon to be started.
+ **********************************************************************/
+ public void Clear() {
+ _num = 0;
+ _crossings = 0;
+ _perimetersum.Set(0);
+ if (!_polyline) _areasum.Set(0);
+ _lat0 = _lon0 = _lat1 = _lon1 = Double.NaN;
+ }
+
+ /**
+ * Add a point to the polygon or polyline.
+ *
+ * @param lat the latitude of the point (degrees).
+ * @param lon the latitude of the point (degrees).
+ *
+ * lat should be in the range [−90°, 90°] and lon
+ * should be in the range [−540°, 540°).
+ **********************************************************************/
+ public void AddPoint(double lat, double lon) {
+ lon = GeoMath.AngNormalize(lon);
+ if (_num == 0) {
+ _lat0 = _lat1 = lat;
+ _lon0 = _lon1 = lon;
+ } else {
+ GeodesicData g = _earth.Inverse(_lat1, _lon1, lat, lon, _mask);
+ _perimetersum.Add(g.s12);
+ if (!_polyline) {
+ _areasum.Add(g.S12);
+ _crossings += transit(_lon1, lon);
+ }
+ _lat1 = lat; _lon1 = lon;
+ }
+ ++_num;
+ }
+
+ /**
+ * Add an edge to the polygon or polyline.
+ *
+ * @param azi azimuth at current point (degrees).
+ * @param s distance from current point to next point (meters).
+ *
+ * azi should be in the range [−540°, 540°). This does
+ * nothing if no points have been added yet. Use PolygonArea.CurrentPoint to
+ * determine the position of the new vertex.
+ **********************************************************************/
+ public void AddEdge(double azi, double s) {
+ if (_num > 0) { // Do nothing if _num is zero
+ GeodesicData g = _earth.Direct(_lat1, _lon1, azi, s, _mask);
+ _perimetersum.Add(g.s12);
+ if (!_polyline) {
+ _areasum.Add(g.S12);
+ _crossings += transit(_lon1, g.lon2);
+ }
+ _lat1 = g.lat2; _lon1 = g.lon2;
+ ++_num;
+ }
+ }
+
+ /**
+ * Return the results so far.
+ *
+ * @return PolygonResult(num, perimeter, area) where
+ * num is the number of vertices, perimeter is the perimeter
+ * of the polygon or the length of the polyline (meters), and area
+ * is the area of the polygon (meters2) or Double.NaN of
+ * polyline is true in the constructor.
+ *
+ * Counter-clockwise traversal counts as a positive area.
+ **********************************************************************/
+ public PolygonResult Compute() { return Compute(false, true); }
+ /**
+ * Return the results so far.
+ *
+ * @param reverse if true then clockwise (instead of counter-clockwise)
+ * traversal counts as a positive area.
+ * @param sign if true then return a signed result for the area if
+ * the polygon is traversed in the "wrong" direction instead of returning
+ * the area for the rest of the earth.
+ * @return PolygonResult(num, perimeter, area) where
+ * num is the number of vertices, perimeter is the perimeter
+ * of the polygon or the length of the polyline (meters), and area
+ * is the area of the polygon (meters2) or Double.NaN of
+ * polyline is true in the constructor.
+ **********************************************************************/
+ public PolygonResult Compute(boolean reverse, boolean sign) {
+ if (_num < 2)
+ return new PolygonResult(_num, 0, _polyline ? Double.NaN : 0);
+ if (_polyline)
+ return new PolygonResult(_num, _perimetersum.Sum(), Double.NaN);
+
+ GeodesicData g = _earth.Inverse(_lat1, _lon1, _lat0, _lon0, _mask);
+ Accumulator tempsum = new Accumulator(_areasum);
+ tempsum.Add(g.S12);
+ int crossings = _crossings + transit(_lon1, _lon0);
+ if ((crossings & 1) != 0)
+ tempsum.Add((tempsum.Sum() < 0 ? 1 : -1) * _area0/2);
+ // area is with the clockwise sense. If !reverse convert to
+ // counter-clockwise convention.
+ if (!reverse)
+ tempsum.Negate();
+ // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
+ if (sign) {
+ if (tempsum.Sum() > _area0/2)
+ tempsum.Add(-_area0);
+ else if (tempsum.Sum() <= -_area0/2)
+ tempsum.Add(+_area0);
+ } else {
+ if (tempsum.Sum() >= _area0)
+ tempsum.Add(-_area0);
+ else if (tempsum.Sum() < 0)
+ tempsum.Add(+_area0);
+ }
+ return new PolygonResult(_num, _perimetersum.Sum(g.s12), 0 + tempsum.Sum());
+ }
+
+ /**
+ * Return the results assuming a tentative final test point is added;
+ * however, the data for the test point is not saved. This lets you report
+ * a running result for the perimeter and area as the user moves the mouse
+ * cursor. Ordinary floating point arithmetic is used to accumulate the
+ * data for the test point; thus the area and perimeter returned are less
+ * accurate than if AddPoint and Compute are used.
+ *
+ * @param lat the latitude of the test point (degrees).
+ * @param lon the longitude of the test point (degrees).
+ * @param reverse if true then clockwise (instead of counter-clockwise)
+ * traversal counts as a positive area.
+ * @param sign if true then return a signed result for the area if
+ * the polygon is traversed in the "wrong" direction instead of returning
+ * the area for the rest of the earth.
+ * @return PolygonResult(num, perimeter, area) where
+ * num is the number of vertices, perimeter is the perimeter
+ * of the polygon or the length of the polyline (meters), and area
+ * is the area of the polygon (meters2) or Double.NaN of
+ * polyline is true in the constructor.
+ *
+ * lat should be in the range [−90°, 90°] and lon
+ * should be in the range [−540°, 540°).
+ **********************************************************************/
+ public PolygonResult TestPoint(double lat, double lon,
+ boolean reverse, boolean sign) {
+ if (_num == 0)
+ return new PolygonResult(1, 0, _polyline ? Double.NaN : 0);
+
+ double perimeter = _perimetersum.Sum();
+ double tempsum = _polyline ? 0 : _areasum.Sum();
+ int crossings = _crossings;
+ int num = _num + 1;
+ for (int i = 0; i < (_polyline ? 1 : 2); ++i) {
+ GeodesicData g =
+ _earth.Inverse(i == 0 ? _lat1 : lat, i == 0 ? _lon1 : lon,
+ i != 0 ? _lat0 : lat, i != 0 ? _lon0 : lon,
+ _mask);
+ perimeter += g.s12;
+ if (!_polyline) {
+ tempsum += g.S12;
+ crossings += transit(i == 0 ? _lon1 : lon,
+ i != 0 ? _lon0 : lon);
+ }
+ }
+
+ if (_polyline)
+ return new PolygonResult(num, perimeter, Double.NaN);
+
+ if ((crossings & 1) != 0)
+ tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
+ // area is with the clockwise sense. If !reverse convert to
+ // counter-clockwise convention.
+ if (!reverse)
+ tempsum *= -1;
+ // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
+ if (sign) {
+ if (tempsum > _area0/2)
+ tempsum -= _area0;
+ else if (tempsum <= -_area0/2)
+ tempsum += _area0;
+ } else {
+ if (tempsum >= _area0)
+ tempsum -= _area0;
+ else if (tempsum < 0)
+ tempsum += _area0;
+ }
+ return new PolygonResult(num, perimeter, 0 + tempsum);
+ }
+
+ /**
+ * Return the results assuming a tentative final test point is added via an
+ * azimuth and distance; however, the data for the test point is not saved.
+ * This lets you report a running result for the perimeter and area as the
+ * user moves the mouse cursor. Ordinary floating point arithmetic is used
+ * to accumulate the data for the test point; thus the area and perimeter
+ * returned are less accurate than if AddPoint and Compute are used.
+ *
+ * @param azi azimuth at current point (degrees).
+ * @param s distance from current point to final test point (meters).
+ * @param reverse if true then clockwise (instead of counter-clockwise)
+ * traversal counts as a positive area.
+ * @param sign if true then return a signed result for the area if
+ * the polygon is traversed in the "wrong" direction instead of returning
+ * the area for the rest of the earth.
+ * @return PolygonResult(num, perimeter, area) where
+ * num is the number of vertices, perimeter is the perimeter
+ * of the polygon or the length of the polyline (meters), and area
+ * is the area of the polygon (meters2) or Double.NaN of
+ * polyline is true in the constructor.
+ *
+ * azi should be in the range [−540°, 540°).
+ **********************************************************************/
+ public PolygonResult TestEdge(double azi, double s,
+ boolean reverse, boolean sign) {
+ if (_num == 0) // we don't have a starting point!
+ return new PolygonResult(0, Double.NaN, Double.NaN);
+
+ int num = _num + 1;
+ double perimeter = _perimetersum.Sum() + s;
+ if (_polyline)
+ return new PolygonResult(num, perimeter, Double.NaN);
+
+ double tempsum = _areasum.Sum();
+ int crossings = _crossings;
+ {
+ double lat, lon, s12, S12, t;
+ GeodesicData g =
+ _earth.Direct(_lat1, _lon1, azi, false, s, _mask);
+ tempsum += g.S12;
+ crossings += transit(_lon1, g.lon2);
+ g = _earth.Inverse(g.lat2, g.lon2, _lat0, _lon0, _mask);
+ perimeter += g.s12;
+ tempsum += g.S12;
+ crossings += transit(g.lon2, _lon0);
+ }
+
+ if ((crossings & 1) != 0)
+ tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
+ // area is with the clockwise sense. If !reverse convert to
+ // counter-clockwise convention.
+ if (!reverse)
+ tempsum *= -1;
+ // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
+ if (sign) {
+ if (tempsum > _area0/2)
+ tempsum -= _area0;
+ else if (tempsum <= -_area0/2)
+ tempsum += _area0;
+ } else {
+ if (tempsum >= _area0)
+ tempsum -= _area0;
+ else if (tempsum < 0)
+ tempsum += _area0;
+ }
+
+ return new PolygonResult(num, perimeter, 0 + tempsum);
+ }
+
+ /**
+ * @return a the equatorial radius of the ellipsoid (meters). This is
+ * the value inherited from the Geodesic object used in the constructor.
+ **********************************************************************/
+
+ public double MajorRadius() { return _earth.MajorRadius(); }
+
+ /**
+ * @return f the flattening of the ellipsoid. This is the value
+ * inherited from the Geodesic object used in the constructor.
+ **********************************************************************/
+ public double Flattening() { return _earth.Flattening(); }
+
+ /**
+ * Report the previous vertex added to the polygon or polyline.
+ *
+ * @return Pair(lat, lon), the current latitude and longitude.
+ *
+ * If no points have been added, then Double.NaN is returned. Otherwise,
+ * lon will be in the range [−180°, 180°).
+ **********************************************************************/
+ public Pair CurrentPoint() { return new Pair(_lat1, _lon1); }
+}
diff --git a/OpenVehicleApp/src/net/sf/geographiclib/PolygonResult.java b/OpenVehicleApp/src/net/sf/geographiclib/PolygonResult.java
new file mode 100644
index 00000000..fdbb7ef1
--- /dev/null
+++ b/OpenVehicleApp/src/net/sf/geographiclib/PolygonResult.java
@@ -0,0 +1,39 @@
+/**
+ * Implementation of the net.sf.geographiclib.PolygonResult class
+ *
+ * Copyright (c) Charles Karney (2013)
+ * @param num the number of vertices in the polygon.
+ * @param perimeter the perimeter of the polygon or the length of the
+ * polyline (meters).
+ * @param area the area of the polygon (meters2).
+ **********************************************************************/
+ public PolygonResult(int num, double perimeter, double area) {
+ this.num = num;
+ this.perimeter = perimeter;
+ this.area = area;
+ }
+}
diff --git a/OpenVehicleApp/src/net/sf/geographiclib/package-info.java b/OpenVehicleApp/src/net/sf/geographiclib/package-info.java
new file mode 100644
index 00000000..6bceb8f6
--- /dev/null
+++ b/OpenVehicleApp/src/net/sf/geographiclib/package-info.java
@@ -0,0 +1,182 @@
+/**
+ *
+ * This is a Java implementation of the geodesic algorithms from GeographicLib. This is a
+ * self-contained library which makes it easy to do geodesic computations
+ * for an ellipsoid of revolution in a Java program. It requires Java
+ * version 1.1 or later.
+ *
+ *
+ * The Java library is part of GeographicLib which available for download at
+ *
+ * as either a compressed tar file (tar.gz) or a zip file. After unpacking
+ * the source, the Java library can be found in GeographicLib-1.31/java. (This
+ * library is completely independent from the rest of GeodegraphicLib.) The
+ * library consists of the files in the src/main/java/net/sf/geographiclib
+ * subdirectory.
+ *
+ *
+ * Also included are 3 small test programs
+ *
+ * Here, for example, is {@code Inverse.java}
+ *
+ *
+ * The important classes are
+ *
+ * The documentation is generated using javadoc when {@code mvn package} is run
+ * (the top of the documentation tree is {@code target/apidocs/index.html}).
+ * This is also available on the web at
+ *
+ * http://geographiclib.sf.net/html/C/index.html.
+ *
+ *
+ *
+ *
+ *
+ *
+ *
+ * {@code
+ * GeodesicData g = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2,
+ * GeodesicMask.DISTANCE); }
+ *
+ *
+ *
+ * |f| error
+ * 0.01 25 nm
+ * 0.02 30 nm
+ * 0.05 10 um
+ * 0.1 1.5 mm
+ * 0.2 300 mm
+ *
+ *
+ *
+ * {@code
+ * // Solve the direct geodesic problem.
+ *
+ * // This program reads in lines with lat1, lon1, azi1, s12 and prints
+ * // out lines with lat2, lon2, azi2 (for the WGS84 ellipsoid).
+ *
+ * import java.util.*;
+ * import net.sf.geographiclib.*;
+ * public class Direct {
+ * public static void main(String[] args) {
+ * try {
+ * Scanner in = new Scanner(System.in);
+ * double lat1, lon1, azi1, s12;
+ * while (true) {
+ * lat1 = in.nextDouble(); lon1 = in.nextDouble();
+ * azi1 = in.nextDouble(); s12 = in.nextDouble();
+ * GeodesicData g = Geodesic.WGS84.Direct(lat1, lon1, azi1, s12);
+ * System.out.println(g.lat2 + " " + g.lon2 + " " + g.azi2);
+ * }
+ * }
+ * catch (Exception e) {}
+ * }
+ * }}
+ **********************************************************************/
+public class Geodesic {
+
+ /**
+ * The order of the expansions used by Geodesic.
+ **********************************************************************/
+ protected static final int GEOGRAPHICLIB_GEODESIC_ORDER = 6;
+
+ protected static final int nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
+ protected static final int nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
+ protected static final int nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER;
+ protected static final int nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
+ protected static final int nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
+ protected static final int nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
+ protected static final int nA3x_ = nA3_;
+ protected static final int nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
+ protected static final int nC3x_ = (nC3_ * (nC3_ - 1)) / 2;
+ protected static final int nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER;
+ protected static final int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
+ private static final int maxit1_ = 20;
+ private static final int maxit2_ = maxit1_ + GeoMath.digits + 10;
+
+ // Underflow guard. We require
+ // tiny_ * epsilon() > 0
+ // tiny_ + epsilon() == epsilon()
+ protected static final double tiny_ = Math.sqrt(GeoMath.min);
+ private static final double tol0_ = GeoMath.epsilon;
+ // Increase multiplier in defn of tol1_ from 100 to 200 to fix inverse case
+ // 52.784459512564 0 -52.784459512563990912 179.634407464943777557
+ // which otherwise failed for Visual Studio 10 (Release and Debug)
+ private static final double tol1_ = 200 * tol0_;
+ private static final double tol2_ = Math.sqrt(tol0_);
+ // Check on bisection interval
+ private static final double tolb_ = tol0_ * tol2_;
+ private static final double xthresh_ = 1000 * tol2_;
+
+ protected static double AngRound(double x) {
+ // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
+ // for reals = 0.7 pm on the earth if x is an angle in degrees. (This
+ // is about 1000 times more resolution than we get with angles around 90
+ // degrees.) We use this to avoid having to deal with near singular
+ // cases when x is non-zero but tiny (e.g., 1.0e-200).
+ final double z = 1/16.0;
+ double y = Math.abs(x);
+ // The compiler mustn't "simplify" z - (z - y) to y
+ y = y < z ? z - (z - y) : y;
+ return x < 0 ? -y : y;
+ }
+
+ protected static Pair SinCosNorm(double sinx, double cosx) {
+ double r = GeoMath.hypot(sinx, cosx);
+ return new Pair(sinx/r, cosx/r);
+ }
+
+ protected double _a, _f, _f1, _e2, _ep2, _b, _c2;
+ private double _n, _etol2;
+ private double _A3x[], _C3x[], _C4x[];
+
+ /**
+ * Constructor for a ellipsoid with
+ *
+ *
+ *
+ *
+ *
+ *
+ *
+ *
+ *
+ * {@code
+ * import net.sf.geographiclib.*;
+ * public class GeodesicLineTest {
+ * public static void main(String[] args) {
+ * // Print waypoints between JFK and SIN
+ * Geodesic geod = Geodesic.WGS84;
+ * double
+ * lat1 = 40.640, lon1 = -73.779, // JFK
+ * lat2 = 1.359, lon2 = 103.989; // SIN
+ * GeodesicData g = geod.Inverse(lat1, lon1, lat2, lon2,
+ * GeodesicMask.DISTANCE | GeodesicMask.AZIMUTH);
+ * GeodesicLine line = new GeodesicLine(geod, lat1, lon1, g.azi1,
+ * GeodesicMask.DISTANCE_IN | GeodesicMask.LONGITUDE);
+ * double
+ * s12 = g.s12,
+ * a12 = g.a12,
+ * ds0 = 500e3; // Nominal distance between points = 500 km
+ * int num = (int)(Math.ceil(s12 / ds0)); // The number of intervals
+ * {
+ * // Use intervals of equal length
+ * double ds = s12 / num;
+ * for (int i = 0; i <= num; ++i) {
+ * g = line.Position(i * ds,
+ * GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE);
+ * System.out.println(i + " " + g.lat2 + " " + g.lon2);
+ * }
+ * }
+ * {
+ * // Slightly faster, use intervals of equal arc length
+ * double da = a12 / num;
+ * for (int i = 0; i <= num; ++i) {
+ * g = line.ArcPosition(i * da,
+ * GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE);
+ * System.out.println(i + " " + g.lat2 + " " + g.lon2);
+ * }
+ * }
+ * }
+ * }}
+ **********************************************************************/
+public class GeodesicLine {
+
+ private static final int nC1_ = Geodesic.nC1_;
+ private static final int nC1p_ = Geodesic.nC1p_;
+ private static final int nC2_ = Geodesic.nC2_;
+ private static final int nC3_ = Geodesic.nC3_;
+ private static final int nC4_ = Geodesic.nC4_;
+
+ private double _lat1, _lon1, _azi1;
+ private double _a, _f, _b, _c2, _f1, _salp0, _calp0, _k2,
+ _salp1, _calp1, _ssig1, _csig1, _dn1, _stau1, _ctau1, _somg1, _comg1,
+ _A1m1, _A2m1, _A3c, _B11, _B21, _B31, _A4, _B41;
+ // index zero elements of _C1a, _C1pa, _C2a, _C3a are unused
+ private double _C1a[], _C1pa[], _C2a[], _C3a[],
+ _C4a[]; // all the elements of _C4a are used
+ private int _caps;
+
+ /**
+ * Constructor for a geodesic line staring at latitude lat1, longitude
+ * lon1, and azimuth azi1 (all in degrees).
+ *
+ *
+ **********************************************************************/
+ public GeodesicLine(Geodesic g,
+ double lat1, double lon1, double azi1,
+ int caps) {
+ _a = g._a;
+ _f = g._f;
+ _b = g._b;
+ _c2 = g._c2;
+ _f1 = g._f1;
+ // Always allow latitude and azimuth
+ _caps = caps | GeodesicMask.LATITUDE | GeodesicMask.AZIMUTH;
+
+ // Guard against underflow in salp0
+ azi1 = Geodesic.AngRound(GeoMath.AngNormalize(azi1));
+ lon1 = GeoMath.AngNormalize(lon1);
+ _lat1 = lat1;
+ _lon1 = lon1;
+ _azi1 = azi1;
+ // alp1 is in [0, pi]
+ double alp1 = azi1 * GeoMath.degree;
+ // Enforce sin(pi) == 0 and cos(pi/2) == 0. Better to face the ensuing
+ // problems directly than to skirt them.
+ _salp1 = azi1 == -180 ? 0 : Math.sin(alp1);
+ _calp1 = Math.abs(azi1) == 90 ? 0 : Math.cos(alp1);
+ double cbet1, sbet1, phi;
+ phi = lat1 * GeoMath.degree;
+ // Ensure cbet1 = +epsilon at poles
+ sbet1 = _f1 * Math.sin(phi);
+ cbet1 = Math.abs(lat1) == 90 ? Geodesic.tiny_ : Math.cos(phi);
+ { Pair p = Geodesic.SinCosNorm(sbet1, cbet1);
+ sbet1 = p.first; cbet1 = p.second; }
+ _dn1 = Math.sqrt(1 + g._ep2 * GeoMath.sq(sbet1));
+
+ // Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
+ _salp0 = _salp1 * cbet1; // alp0 in [0, pi/2 - |bet1|]
+ // Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following
+ // is slightly better (consider the case salp1 = 0).
+ _calp0 = GeoMath.hypot(_calp1, _salp1 * sbet1);
+ // Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
+ // sig = 0 is nearest northward crossing of equator.
+ // With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
+ // With bet1 = pi/2, alp1 = -pi, sig1 = pi/2
+ // With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2
+ // Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
+ // With alp0 in (0, pi/2], quadrants for sig and omg coincide.
+ // No atan2(0,0) ambiguity at poles since cbet1 = +epsilon.
+ // With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
+ _ssig1 = sbet1; _somg1 = _salp0 * sbet1;
+ _csig1 = _comg1 = sbet1 != 0 || _calp1 != 0 ? cbet1 * _calp1 : 1;
+ { Pair p = Geodesic.SinCosNorm(_ssig1, _csig1);
+ _ssig1 = p.first; _csig1 = p.second; } // sig1 in (-pi, pi]
+ // Geodesic.SinCosNorm(_somg1, _comg1); -- don't need to normalize!
+
+ _k2 = GeoMath.sq(_calp0) * g._ep2;
+ double eps = _k2 / (2 * (1 + Math.sqrt(1 + _k2)) + _k2);
+
+ if ((_caps & GeodesicMask.CAP_C1) != 0) {
+ _A1m1 = Geodesic.A1m1f(eps);
+ _C1a = new double[nC1_ + 1];
+ Geodesic.C1f(eps, _C1a);
+ _B11 = Geodesic.SinCosSeries(true, _ssig1, _csig1, _C1a);
+ double s = Math.sin(_B11), c = Math.cos(_B11);
+ // tau1 = sig1 + B11
+ _stau1 = _ssig1 * c + _csig1 * s;
+ _ctau1 = _csig1 * c - _ssig1 * s;
+ // Not necessary because C1pa reverts C1a
+ // _B11 = -SinCosSeries(true, _stau1, _ctau1, _C1pa, nC1p_);
+ }
+
+ if ((_caps & GeodesicMask.CAP_C1p) != 0) {
+ _C1pa = new double[nC1p_ + 1];
+ Geodesic.C1pf(eps, _C1pa);
+ }
+
+ if ((_caps & GeodesicMask.CAP_C2) != 0) {
+ _C2a = new double[nC2_ + 1];
+ _A2m1 = Geodesic.A2m1f(eps);
+ Geodesic.C2f(eps, _C2a);
+ _B21 = Geodesic.SinCosSeries(true, _ssig1, _csig1, _C2a);
+ }
+
+ if ((_caps & GeodesicMask.CAP_C3) != 0) {
+ _C3a = new double[nC3_];
+ g.C3f(eps, _C3a);
+ _A3c = -_f * _salp0 * g.A3f(eps);
+ _B31 = Geodesic.SinCosSeries(true, _ssig1, _csig1, _C3a);
+ }
+
+ if ((_caps & GeodesicMask.CAP_C4) != 0) {
+ _C4a = new double[nC4_];
+ g.C4f(eps, _C4a);
+ // Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
+ _A4 = GeoMath.sq(_a) * _calp0 * _salp0 * g._e2;
+ _B41 = Geodesic.SinCosSeries(false, _ssig1, _csig1, _C4a);
+ }
+ }
+
+ /**
+ * A default constructor. If GeodesicLine.Position is called on the
+ * resulting object, it returns immediately (without doing any calculations).
+ * The object can be set with a call to {@link Geodesic.Line}. Use {@link
+ * Init()} to test whether object is still in this uninitialized state.
+ * (This constructor was useful in C++, e.g., to allow vectors of
+ * GeodesicLine objects. It may not be needed in Java, so make it private.)
+ **********************************************************************/
+ private GeodesicLine() { _caps = 0; }
+
+ /**
+ * Compute the position of point 2 which is a distance s12 (meters)
+ * from point 1.
+ *
+ *
+ *
+ *
+ *
+ * {@code
+ * // Compute the area of a geodesic polygon.
+ *
+ * // This program reads lines with lat, lon for each vertex of a polygon.
+ * // At the end of input, the program prints the number of vertices,
+ * // the perimeter of the polygon and its area (for the WGS84 ellipsoid).
+ *
+ * import java.util.*;
+ * import net.sf.geographiclib.*;
+ *
+ * public class Planimeter {
+ * public static void main(String[] args) {
+ * PolygonArea p = new PolygonArea(Geodesic.WGS84, false);
+ * try {
+ * Scanner in = new Scanner(System.in);
+ * while (true) {
+ * double lat = in.nextDouble(), lon = in.nextDouble();
+ * p.AddPoint(lat, lon);
+ * }
+ * }
+ * catch (Exception e) {}
+ * PolygonResult r = p.Compute();
+ * System.out.println(r.num + " " + r.perimeter + " " + r.area);
+ * }
+ * }}
+ **********************************************************************/
+public class PolygonArea {
+
+ private Geodesic _earth;
+ private double _area0; // Full ellipsoid area
+ private boolean _polyline; // Assume polyline (don't close and skip area)
+ private int _mask;
+ private int _num;
+ private int _crossings;
+ private Accumulator _areasum, _perimetersum;
+ private double _lat0, _lon0, _lat1, _lon1;
+ private static int transit(double lon1, double lon2) {
+ // Return 1 or -1 if crossing prime meridian in east or west direction.
+ // Otherwise return zero.
+ // Compute lon12 the same way as Geodesic.Inverse.
+ lon1 = GeoMath.AngNormalize(lon1);
+ lon2 = GeoMath.AngNormalize(lon2);
+ double lon12 = GeoMath.AngDiff(lon1, lon2);
+ int cross =
+ lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
+ (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
+ return cross;
+ }
+
+ /**
+ * Constructor for PolygonArea.
+ * Geodesic routines from GeographicLib implemented in Java
+ * @author Charles F. F. Karney (charles@karney.com)
+ * @version 1.31
+ *
+ * Abstract
+ * Downloading the source
+ *
+ *
+ * Sample programs
+ *
+ *
+ *
+ * {@code
+ * // Solve the inverse geodesic problem.
+ *
+ * // This program reads in lines with lat1, lon1, lat2, lon2 and prints
+ * // out lines with azi1, azi2, s12 (for the WGS84 ellipsoid).
+ *
+ * import java.util.*;
+ * import net.sf.geographiclib.*;
+ * public class Inverse {
+ * public static void main(String[] args) {
+ * try {
+ * Scanner in = new Scanner(System.in);
+ * double lat1, lon1, lat2, lon2;
+ * while (true) {
+ * lat1 = in.nextDouble(); lon1 = in.nextDouble();
+ * lat2 = in.nextDouble(); lon2 = in.nextDouble();
+ * GeodesicData g = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2);
+ * System.out.println(g.azi1 + " " + g.azi2 + " " + g.s12);
+ * }
+ * }
+ * catch (Exception e) {}
+ * }
+ * }}
+ * Compiling and running a sample program
+ * Three difference ways of compiling and running {@code Inverse.java} are
+ * given. These differ in the degree to which they utilize
+ * maven to manage your Java code and
+ * its dependencies. (Thanks to Skip Breidbach for supplying the maven
+ * support.)
+ * Without using maven
+ * Compile and run as follows
+ * cd inverse/src/main/java
+ * javac -cp .:../../../../src/main/java Inverse.java
+ * echo -30 0 29.5 179.5 | java -cp .:../../../../src/main/java Inverse
+ * Using maven to package GeographicLib
+ * Use maven to create a jar file by
+ * running (in the main java directory)
+ * mvn package
+ * (Your first run of maven may take a long time, because it needs to download
+ * some additional packages to your local repository.) Then compile and run
+ * Inverse.java with
+ * cd inverse/src/main/java
+ * javac -cp .:../../../../target/GeographicLib-1.31.jar Inverse.java
+ * echo -30 0 29.5 179.5 |
+ * java -cp .:../../../../target/GeographicLib-1.31.jar Inverse
+ * Using maven to build and run {@code Inverse.java}
+ * Use maven to install GeographicLib by
+ * running (in the main java directory)
+ * mvn install
+ * (Your first run of maven may take a long time, because it needs to download
+ * some additional packages to your local repository.) Then compile and run
+ * Inverse.java using {@code inverse/pom.xml} with
+ * cd inverse
+ * mvn compile
+ * echo -30 0 29.5 179.5 | mvn -q exec:java
+ * Using the library
+ *
+ *
+ *
+ * import net.sf.geographiclib.*
+ * in your source code.
+ *
+ *
+ * External links
+ *
+ *
+ **********************************************************************/
+package net.sf.geographiclib;