daoas/fulddaekning/GeoPoint.java

package dk.daoas.fulddaekning;

public class GeoPoint {

        public double latitude;
        public double longitude;

        public GeoPoint() { // Default
        }

        public GeoPoint(double lat, double lng) { // Default
                latitude = lat;
                longitude = lng;
        }

        // Denne er alt for upræcis
        @Deprecated
        public static double beregnAfstand_old(GeoPoint point1, GeoPoint point2) {
                // (62.8*sqrt(3.1*(Power(a.Latitude-x.Latitude,2)+Power(a.Longitude-x.Longitude,2))))
                // as Afstand,

                double pwrLat = Math.pow(point1.latitude - point2.latitude, 2);
                double pwrLng = Math.pow(point1.longitude - point2.longitude, 2);

                return 62.8 * Math.sqrt(3.1 * (pwrLat + pwrLng));
        }
        
        
        // Latitude (horizonal), longitude(vertical) so
        // 1 degree latitude is ~ 111320 meters, since the distance between the
        // horizonal lines is always the same
        // 1 degree longitude is ~111320 meters at equator but gets shorter as we
        // get closer to the poles.
        // so 1 degree longitude is 64.5 km at denmarks southern point
        // (gedser=54.55,11.95)
        // and is 59.4km at northern point (skagen = 57.75,10.65)

        public static double kmToLatitude(double km) {
                return km / 111.320;
        }

        public static double kmToLongitude(double km) {// denne er kun ca
                return km / 62.0;
        }
        
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

        /* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
        /* :: This function converts decimal degrees to radians : */
        /* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
        private static double deg2rad(double deg) {
                return (deg * Math.PI / 180.0);
        }

        /* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
        /* :: This function converts radians to decimal degrees : */
        /* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
        private static double rad2deg(double rad) {
                return (rad * 180 / Math.PI);
        }

        // http://www.geodatasource.com/developers/java
        private static double distanceHaversine(double lat1, double lon1,
                        double lat2, double lon2) {
                double theta = lon1 - lon2;
                double dist = Math.sin(deg2rad(lat1)) * Math.sin(deg2rad(lat2))
                                + Math.cos(deg2rad(lat1)) * Math.cos(deg2rad(lat2))
                                * Math.cos(deg2rad(theta));
                dist = Math.acos(dist);
                dist = rad2deg(dist);
                dist = dist * 60 * 1.1515;

                // indtil nu er dist i miles - så vi omregner lige til km
                dist = dist * 1.609344;
                return (dist);
        }

        public static double beregnAfstand(GeoPoint p1, GeoPoint p2) {
                return distanceHaversine(p1.latitude, p1.longitude, p2.latitude,
                                p2.longitude);
        }
        
        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

        // denne er nok den mest præcise - men er også den langsomste
        public static float beregnAfstand_google(GeoPoint p1, GeoPoint p2) {
                float[] result = new float[1];

                computeDistanceAndBearing(p1.latitude, p1.longitude, p2.latitude,
                                p2.longitude, result);

                return (result[0] / 1000.0f);
        }


        // Kopieret fra android.location.Location
        private static void computeDistanceAndBearing(double lat1, double lon1,
                        double lat2, double lon2, float[] results) {
                // Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
                // using the "Inverse Formula" (section 4)

                int MAXITERS = 20;
                // Convert lat/long to radians
                lat1 *= Math.PI / 180.0;
                lat2 *= Math.PI / 180.0;
                lon1 *= Math.PI / 180.0;
                lon2 *= Math.PI / 180.0;

                double a = 6378137.0; // WGS84 major axis
                double b = 6356752.3142; // WGS84 semi-major axis
                double f = (a - b) / a;
                double aSqMinusBSqOverBSq = (a * a - b * b) / (b * b);

                double L = lon2 - lon1;
                double A = 0.0;
                double U1 = Math.atan((1.0 - f) * Math.tan(lat1));
                double U2 = Math.atan((1.0 - f) * Math.tan(lat2));

                double cosU1 = Math.cos(U1);
                double cosU2 = Math.cos(U2);
                double sinU1 = Math.sin(U1);
                double sinU2 = Math.sin(U2);
                double cosU1cosU2 = cosU1 * cosU2;
                double sinU1sinU2 = sinU1 * sinU2;

                double sigma = 0.0;
                double deltaSigma = 0.0;
                double cosSqAlpha = 0.0;
                double cos2SM = 0.0;
                double cosSigma = 0.0;
                double sinSigma = 0.0;
                double cosLambda = 0.0;
                double sinLambda = 0.0;

                double lambda = L; // initial guess
                for (int iter = 0; iter < MAXITERS; iter++) {
                        double lambdaOrig = lambda;
                        cosLambda = Math.cos(lambda);
                        sinLambda = Math.sin(lambda);
                        double t1 = cosU2 * sinLambda;
                        double t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
                        double sinSqSigma = t1 * t1 + t2 * t2; // (14)
                        sinSigma = Math.sqrt(sinSqSigma);
                        cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda; // (15)
                        sigma = Math.atan2(sinSigma, cosSigma); // (16)
                        double sinAlpha = (sinSigma == 0) ? 0.0 : cosU1cosU2 * sinLambda
                                        / sinSigma; // (17)
                        cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
                        cos2SM = (cosSqAlpha == 0) ? 0.0 : cosSigma - 2.0 * sinU1sinU2
                                        / cosSqAlpha; // (18)

                        double uSquared = cosSqAlpha * aSqMinusBSqOverBSq; // defn
                        A = 1
                                        + (uSquared / 16384.0)
                                        * // (3)
                                        (4096.0 + uSquared
                                                        * (-768 + uSquared * (320.0 - 175.0 * uSquared)));
                        double B = (uSquared / 1024.0) * // (4)
                                        (256.0 + uSquared
                                                        * (-128.0 + uSquared * (74.0 - 47.0 * uSquared)));
                        double C = (f / 16.0) * cosSqAlpha
                                        * (4.0 + f * (4.0 - 3.0 * cosSqAlpha)); // (10)
                        double cos2SMSq = cos2SM * cos2SM;
                        deltaSigma = B
                                        * sinSigma
                                        * // (6)
                                        (cos2SM + (B / 4.0)
                                                        * (cosSigma * (-1.0 + 2.0 * cos2SMSq) - (B / 6.0)
                                                                        * cos2SM
                                                                        * (-3.0 + 4.0 * sinSigma * sinSigma)
                                                                        * (-3.0 + 4.0 * cos2SMSq)));

                        lambda = L
                                        + (1.0 - C)
                                        * f
                                        * sinAlpha
                                        * (sigma + C
                                                        * sinSigma
                                                        * (cos2SM + C * cosSigma
                                                                        * (-1.0 + 2.0 * cos2SM * cos2SM))); // (11)

                        double delta = (lambda - lambdaOrig) / lambda;
                        if (Math.abs(delta) < 1.0e-12) {
                                break;
                        }
                }

                float distance = (float) (b * A * (sigma - deltaSigma));
                results[0] = distance;
                if (results.length > 1) {
                        float initialBearing = (float) Math.atan2(cosU2 * sinLambda, cosU1
                                        * sinU2 - sinU1 * cosU2 * cosLambda);
                        initialBearing *= 180.0 / Math.PI;
                        results[1] = initialBearing;
                        if (results.length > 2) {
                                float finalBearing = (float) Math.atan2(cosU1 * sinLambda,
                                                -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda);
                                finalBearing *= 180.0 / Math.PI;
                                results[2] = finalBearing;
                        }
                }
        }

        /**
         * Computes the approximate distance in meters between two locations, and
         * optionally the initial and final bearings of the shortest path between
         * them. Distance and bearing are defined using the WGS84 ellipsoid.
         *
         * <p>
         * The computed distance is stored in results[0]. If results has length 2 or
         * greater, the initial bearing is stored in results[1]. If results has
         * length 3 or greater, the final bearing is stored in results[2].
         *
         * @param startLatitude
         *            the starting latitude
         * @param startLongitude
         *            the starting longitude
         * @param endLatitude
         *            the ending latitude
         * @param endLongitude
         *            the ending longitude
         * @param results
         *            an array of floats to hold the results
         *
         * @throws IllegalArgumentException
         *             if results is null or has length < 1
         */
        public static void distanceBetween(double startLatitude,
                        double startLongitude, double endLatitude, double endLongitude,
                        float[] results) {
                if (results == null || results.length < 1) {
                        throw new IllegalArgumentException(
                                        "results is null or has length < 1");
                }
                computeDistanceAndBearing(startLatitude, startLongitude, endLatitude,
                                endLongitude, results);
        }

}
1	torben	2211	package dk.daoas.fulddaekning;
2
3			public class GeoPoint {
4	torben	2580
5	torben	2211	public double latitude;
6			public double longitude;
7	torben	2580
8			public GeoPoint() { // Default
9			}
10
11			public GeoPoint(double lat, double lng) { // Default
12			latitude = lat;
13			longitude = lng;
14			}
15
16			// Denne er alt for upræcis
17			@Deprecated
18			public static double beregnAfstand_old(GeoPoint point1, GeoPoint point2) {
19			// (62.8sqrt(3.1(Power(a.Latitude-x.Latitude,2)+Power(a.Longitude-x.Longitude,2))))
20			// as Afstand,
21
22			double pwrLat = Math.pow(point1.latitude - point2.latitude, 2);
23			double pwrLng = Math.pow(point1.longitude - point2.longitude, 2);
24
25			return 62.8 * Math.sqrt(3.1 * (pwrLat + pwrLng));
26			}
27	torben	2211
28	torben	2580
29			// Latitude (horizonal), longitude(vertical) so
30			// 1 degree latitude is ~ 111320 meters, since the distance between the
31			// horizonal lines is always the same
32			// 1 degree longitude is ~111320 meters at equator but gets shorter as we
33			// get closer to the poles.
34			// so 1 degree longitude is 64.5 km at denmarks southern point
35			// (gedser=54.55,11.95)
36			// and is 59.4km at northern point (skagen = 57.75,10.65)
37
38			public static double kmToLatitude(double km) {
39			return km / 111.320;
40	torben	2211	}
41	torben	2580
42			public static double kmToLongitude(double km) {// denne er kun ca
43			return km / 62.0;
44			}
45	torben	2211
46	torben	2580	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
47
48			/* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
49			/* :: This function converts decimal degrees to radians : */
50			/* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
51			private static double deg2rad(double deg) {
52			return (deg * Math.PI / 180.0);
53	torben	2211	}
54	torben	2578
55	torben	2580	/* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
56			/* :: This function converts radians to decimal degrees : */
57			/* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
58			private static double rad2deg(double rad) {
59			return (rad * 180 / Math.PI);
60			}
61	torben	2577
62	torben	2580	// http://www.geodatasource.com/developers/java
63			private static double distanceHaversine(double lat1, double lon1,
64			double lat2, double lon2) {
65			double theta = lon1 - lon2;
66			double dist = Math.sin(deg2rad(lat1)) * Math.sin(deg2rad(lat2))
67			+ Math.cos(deg2rad(lat1)) * Math.cos(deg2rad(lat2))
68			* Math.cos(deg2rad(theta));
69			dist = Math.acos(dist);
70			dist = rad2deg(dist);
71			dist = dist * 60 * 1.1515;
72	torben	2579
73	torben	2580	// indtil nu er dist i miles - så vi omregner lige til km
74			dist = dist * 1.609344;
75			return (dist);
76			}
77	torben	2579
78			public static double beregnAfstand(GeoPoint p1, GeoPoint p2) {
79	torben	2580	return distanceHaversine(p1.latitude, p1.longitude, p2.latitude,
80			p2.longitude);
81	torben	2579	}
82	torben	2580
83			////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
84	torben	2579
85			// denne er nok den mest præcise - men er også den langsomste
86			public static float beregnAfstand_google(GeoPoint p1, GeoPoint p2) {
87	torben	2577	float[] result = new float[1];
88
89	torben	2580	computeDistanceAndBearing(p1.latitude, p1.longitude, p2.latitude,
90			p2.longitude, result);
91
92	torben	2578	return (result[0] / 1000.0f);
93	torben	2577	}
94
95
96	torben	2580	// Kopieret fra android.location.Location
97			private static void computeDistanceAndBearing(double lat1, double lon1,
98			double lat2, double lon2, float[] results) {
99			// Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
100			// using the "Inverse Formula" (section 4)
101	torben	2577
102	torben	2580	int MAXITERS = 20;
103			// Convert lat/long to radians
104			lat1 *= Math.PI / 180.0;
105			lat2 *= Math.PI / 180.0;
106			lon1 *= Math.PI / 180.0;
107			lon2 *= Math.PI / 180.0;
108	torben	2577
109	torben	2580	double a = 6378137.0; // WGS84 major axis
110			double b = 6356752.3142; // WGS84 semi-major axis
111			double f = (a - b) / a;
112			double aSqMinusBSqOverBSq = (a * a - b * b) / (b * b);
113	torben	2577
114	torben	2580	double L = lon2 - lon1;
115			double A = 0.0;
116			double U1 = Math.atan((1.0 - f) * Math.tan(lat1));
117			double U2 = Math.atan((1.0 - f) * Math.tan(lat2));
118	torben	2577
119	torben	2580	double cosU1 = Math.cos(U1);
120			double cosU2 = Math.cos(U2);
121			double sinU1 = Math.sin(U1);
122			double sinU2 = Math.sin(U2);
123			double cosU1cosU2 = cosU1 * cosU2;
124			double sinU1sinU2 = sinU1 * sinU2;
125	torben	2577
126	torben	2580	double sigma = 0.0;
127			double deltaSigma = 0.0;
128			double cosSqAlpha = 0.0;
129			double cos2SM = 0.0;
130			double cosSigma = 0.0;
131			double sinSigma = 0.0;
132			double cosLambda = 0.0;
133			double sinLambda = 0.0;
134	torben	2577
135	torben	2580	double lambda = L; // initial guess
136			for (int iter = 0; iter < MAXITERS; iter++) {
137			double lambdaOrig = lambda;
138			cosLambda = Math.cos(lambda);
139			sinLambda = Math.sin(lambda);
140			double t1 = cosU2 * sinLambda;
141			double t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
142			double sinSqSigma = t1 * t1 + t2 * t2; // (14)
143			sinSigma = Math.sqrt(sinSqSigma);
144			cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda; // (15)
145			sigma = Math.atan2(sinSigma, cosSigma); // (16)
146			double sinAlpha = (sinSigma == 0) ? 0.0 : cosU1cosU2 * sinLambda
147			/ sinSigma; // (17)
148			cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
149			cos2SM = (cosSqAlpha == 0) ? 0.0 : cosSigma - 2.0 * sinU1sinU2
150			/ cosSqAlpha; // (18)
151	torben	2577
152	torben	2580	double uSquared = cosSqAlpha * aSqMinusBSqOverBSq; // defn
153			A = 1
154			+ (uSquared / 16384.0)
155			* // (3)
156			(4096.0 + uSquared
157			* (-768 + uSquared * (320.0 - 175.0 * uSquared)));
158			double B = (uSquared / 1024.0) * // (4)
159			(256.0 + uSquared
160			* (-128.0 + uSquared * (74.0 - 47.0 * uSquared)));
161			double C = (f / 16.0) * cosSqAlpha
162			* (4.0 + f * (4.0 - 3.0 * cosSqAlpha)); // (10)
163			double cos2SMSq = cos2SM * cos2SM;
164			deltaSigma = B
165			* sinSigma
166			* // (6)
167			(cos2SM + (B / 4.0)
168			* (cosSigma * (-1.0 + 2.0 * cos2SMSq) - (B / 6.0)
169			* cos2SM
170			* (-3.0 + 4.0 * sinSigma * sinSigma)
171			* (-3.0 + 4.0 * cos2SMSq)));
172	torben	2577
173	torben	2580	lambda = L
174			+ (1.0 - C)
175			* f
176			* sinAlpha
177			* (sigma + C
178			* sinSigma
179			* (cos2SM + C * cosSigma
180			* (-1.0 + 2.0 * cos2SM * cos2SM))); // (11)
181	torben	2577
182	torben	2580	double delta = (lambda - lambdaOrig) / lambda;
183			if (Math.abs(delta) < 1.0e-12) {
184			break;
185			}
186			}
187	torben	2577
188	torben	2580	float distance = (float) (b * A * (sigma - deltaSigma));
189			results[0] = distance;
190			if (results.length > 1) {
191			float initialBearing = (float) Math.atan2(cosU2 * sinLambda, cosU1
192			* sinU2 - sinU1 * cosU2 * cosLambda);
193			initialBearing *= 180.0 / Math.PI;
194			results[1] = initialBearing;
195			if (results.length > 2) {
196			float finalBearing = (float) Math.atan2(cosU1 * sinLambda,
197			-sinU1 * cosU2 + cosU1 * sinU2 * cosLambda);
198			finalBearing *= 180.0 / Math.PI;
199			results[2] = finalBearing;
200			}
201			}
202			}
203	torben	2577
204	torben	2580	/**
205			* Computes the approximate distance in meters between two locations, and
206			* optionally the initial and final bearings of the shortest path between
207			* them. Distance and bearing are defined using the WGS84 ellipsoid.
208			*
209			* <p>
210			* The computed distance is stored in results[0]. If results has length 2 or
211			* greater, the initial bearing is stored in results[1]. If results has
212			* length 3 or greater, the final bearing is stored in results[2].
213			*
214			* @param startLatitude
215			* the starting latitude
216			* @param startLongitude
217			* the starting longitude
218			* @param endLatitude
219			* the ending latitude
220			* @param endLongitude
221			* the ending longitude
222			* @param results
223			* an array of floats to hold the results
224			*
225			* @throws IllegalArgumentException
226			* if results is null or has length < 1
227			*/
228			public static void distanceBetween(double startLatitude,
229			double startLongitude, double endLatitude, double endLongitude,
230			float[] results) {
231			if (results == null \|\| results.length < 1) {
232			throw new IllegalArgumentException(
233			"results is null or has length < 1");
234			}
235			computeDistanceAndBearing(startLatitude, startLongitude, endLatitude,
236			endLongitude, results);
237			}
238	torben	2577
239	torben	2211	}