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