001    /*
002     * Zmanim Java API
003     * Copyright (C) 2004-2008 Eliyahu Hershfeld
004     *
005     * This program is free software; you can redistribute it and/or modify it under the terms of the
006     * GNU General Public License as published by the Free Software Foundation; either version 2 of the
007     * License, or (at your option) any later version.
008     *
009     * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
010     * even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
011     * General Public License for more details.
012     *
013     * You should have received a copy of the GNU General Public License along with this program; if
014     * not, write to the Free Software Foundation, Inc. 59 Temple Place - Suite 330, Boston, MA
015     * 02111-1307, USA or connect to: http://www.fsf.org/copyleft/gpl.html
016     */
017    package net.sourceforge.zmanim.util;
018    
019    /**
020     * A class for various location calculations
021     * Most of the code in this class is ported from <a href="http://www.movable-type.co.uk/">Chris Veness'</a>
022     * <a href="http://www.fsf.org/licensing/licenses/lgpl.html">LGPL</a> Javascript Implementation
023     *
024     * @author &copy; Eliyahu Hershfeld 2009
025     * @version 0.1
026     */
027    public class GeoLocationUtils {
028            private static int DISTANCE = 0;
029            private static int INITIAL_BEARING = 1;
030            private static int FINAL_BEARING = 2;
031    
032            /**
033             * Calculate the initial <a
034             * href="http://en.wikipedia.org/wiki/Great_circle">geodesic</a> bearing
035             * between this Object and a second Object passed to this method using <a
036             * href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a>
037             * inverse formula See T Vincenty, "<a
038             * href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse
039             * Solutions of Geodesics on the Ellipsoid with application of nested
040             * equations</a>", Survey Review, vol XXII no 176, 1975.
041             *
042             * @param location
043             *            the destination location
044             */
045            public static double getGeodesicInitialBearing(GeoLocation location, GeoLocation destination) {
046                    return vincentyFormula(location, destination, INITIAL_BEARING);
047            }
048    
049            /**
050             * Calculate the final <a
051             * href="http://en.wikipedia.org/wiki/Great_circle">geodesic</a> bearing
052             * between this Object and a second Object passed to this method using <a
053             * href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a>
054             * inverse formula See T Vincenty, "<a
055             * href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse
056             * Solutions of Geodesics on the Ellipsoid with application of nested
057             * equations</a>", Survey Review, vol XXII no 176, 1975.
058             *
059             * @param location
060             *            the destination location
061             */
062            public static double getGeodesicFinalBearing(GeoLocation location, GeoLocation destination) {
063                    return vincentyFormula(location, destination, FINAL_BEARING);
064            }
065    
066            /**
067             * Calculate <a
068             * href="http://en.wikipedia.org/wiki/Great-circle_distance">geodesic
069             * distance</a> in Meters between this Object and a second Object passed to
070             * this method using <a
071             * href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a>
072             * inverse formula See T Vincenty, "<a
073             * href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse
074             * Solutions of Geodesics on the Ellipsoid with application of nested
075             * equations</a>", Survey Review, vol XXII no 176, 1975.
076             *
077             * @param location
078             *            the destination location
079             */
080            public static double getGeodesicDistance(GeoLocation location, GeoLocation destination) {
081                    return vincentyFormula(location, destination, DISTANCE);
082            }
083    
084            /**
085             * Calculate <a
086             * href="http://en.wikipedia.org/wiki/Great-circle_distance">geodesic
087             * distance</a> in Meters between this Object and a second Object passed to
088             * this method using <a
089             * href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a>
090             * inverse formula See T Vincenty, "<a
091             * href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse
092             * Solutions of Geodesics on the Ellipsoid with application of nested
093             * equations</a>", Survey Review, vol XXII no 176, 1975.
094             *
095             * @param location
096             *            the destination location
097             * @param formula
098             *            This formula calculates initial bearing ({@link #INITIAL_BEARING}),
099             *            final bearing ({@link #FINAL_BEARING}) and distance ({@link #DISTANCE}).
100             */
101            private static double vincentyFormula(GeoLocation location, GeoLocation destination, int formula) {
102                    double a = 6378137;
103                    double b = 6356752.3142;
104                    double f = 1 / 298.257223563; // WGS-84 ellipsiod
105                    double L = Math.toRadians(destination.getLongitude() - location.getLongitude());
106                    double U1 = Math
107                                    .atan((1 - f) * Math.tan(Math.toRadians(location.getLatitude())));
108                    double U2 = Math.atan((1 - f)
109                                    * Math.tan(Math.toRadians(destination.getLatitude())));
110                    double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
111                    double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
112    
113                    double lambda = L;
114                    double lambdaP = 2 * Math.PI;
115                    double iterLimit = 20;
116                    double sinLambda = 0;
117                    double cosLambda = 0;
118                    double sinSigma = 0;
119                    double cosSigma = 0;
120                    double sigma = 0;
121                    double sinAlpha = 0;
122                    double cosSqAlpha = 0;
123                    double cos2SigmaM = 0;
124                    double C;
125                    while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0) {
126                            sinLambda = Math.sin(lambda);
127                            cosLambda = Math.cos(lambda);
128                            sinSigma = Math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda)
129                                            + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda)
130                                            * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda));
131                            if (sinSigma == 0)
132                                    return 0; // co-incident points
133                            cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
134                            sigma = Math.atan2(sinSigma, cosSigma);
135                            sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
136                            cosSqAlpha = 1 - sinAlpha * sinAlpha;
137                            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
138                            if (Double.isNaN(cos2SigmaM))
139                                    cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6)
140                            C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
141                            lambdaP = lambda;
142                            lambda = L
143                                            + (1 - C)
144                                            * f
145                                            * sinAlpha
146                                            * (sigma + C
147                                                            * sinSigma
148                                                            * (cos2SigmaM + C * cosSigma
149                                                                            * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
150                    }
151                    if (iterLimit == 0)
152                            return Double.NaN; // formula failed to converge
153    
154                    double uSq = cosSqAlpha * (a * a - b * b) / (b * b);
155                    double A = 1 + uSq / 16384
156                                    * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
157                    double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
158                    double deltaSigma = B
159                                    * sinSigma
160                                    * (cos2SigmaM + B
161                                                    / 4
162                                                    * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B
163                                                                    / 6 * cos2SigmaM
164                                                                    * (-3 + 4 * sinSigma * sinSigma)
165                                                                    * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
166                    double distance = b * A * (sigma - deltaSigma);
167    
168                    // initial bearing
169                    double fwdAz = Math.toDegrees(Math.atan2(cosU2 * sinLambda, cosU1
170                                    * sinU2 - sinU1 * cosU2 * cosLambda));
171                    // final bearing
172                    double revAz = Math.toDegrees(Math.atan2(cosU1 * sinLambda, -sinU1
173                                    * cosU2 + cosU1 * sinU2 * cosLambda));
174                    if (formula == DISTANCE) {
175                            return distance;
176                    } else if (formula == INITIAL_BEARING) {
177                            return fwdAz;
178                    } else if (formula == FINAL_BEARING) {
179                            return revAz;
180                    } else { // should never happpen
181                            return Double.NaN;
182                    }
183            }
184    
185            /**
186             * Returns the <a href="http://en.wikipedia.org/wiki/Rhumb_line">rhumb line</a>
187             * bearing from the current location to the GeoLocation passed in.
188             *
189             * @param location
190             *            destination location
191             * @return the bearing in degrees
192             */
193            public static double getRhumbLineBearing(GeoLocation location, GeoLocation destination) {
194                    double dLon = Math.toRadians(destination.getLongitude() - location.getLongitude());
195                    double dPhi = Math.log(Math.tan(Math.toRadians(destination.getLatitude())
196                                    / 2 + Math.PI / 4)
197                                    / Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4));
198                    if (Math.abs(dLon) > Math.PI)
199                            dLon = dLon > 0 ? -(2 * Math.PI - dLon) : (2 * Math.PI + dLon);
200                    return Math.toDegrees(Math.atan2(dLon, dPhi));
201            }
202    
203            /**
204             * Returns the <a href="http://en.wikipedia.org/wiki/Rhumb_line">rhumb line</a>
205             * distance from the current location to the GeoLocation passed in.
206             * Ported from <a href="http://www.movable-type.co.uk/">Chris Veness'</a> Javascript Implementation
207             *
208             * @param location
209             *            the destination location
210             * @return the distance in Meters
211             */
212            public static double getRhumbLineDistance(GeoLocation location, GeoLocation destination) {
213                    double R = 6371; // earth's mean radius in km
214                    double dLat = Math.toRadians(destination.getLatitude() - location.getLatitude());
215                    double dLon = Math.toRadians(Math.abs(destination.getLongitude()
216                                    - location.getLongitude()));
217                    double dPhi = Math.log(Math.tan(Math.toRadians(destination.getLongitude())
218                                    / 2 + Math.PI / 4)
219                                    / Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4));
220                    double q = (Math.abs(dLat) > 1e-10) ? dLat / dPhi : Math.cos(Math
221                                    .toRadians(location.getLatitude()));
222                    // if dLon over 180° take shorter rhumb across 180° meridian:
223                    if (dLon > Math.PI)
224                            dLon = 2 * Math.PI - dLon;
225                    double d = Math.sqrt(dLat * dLat + q * q * dLon * dLon);
226                    return d * R;
227            }
228    
229    }