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	package dk.daoas.fulddaekning;
2
3	public class GeoPoint {
4
5	public double latitude;
6	public double longitude;
7
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
28
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	}
41
42	public static double kmToLongitude(double km) {// denne er kun ca
43	return km / 62.0;
44	}
45
46	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
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	}
54
55	/* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
56	/* :: This function converts radians to decimal degrees : */
57	/* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
58	private static double rad2deg(double rad) {
59	return (rad * 180 / Math.PI);
60	}
61
62	// 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
73	// indtil nu er dist i miles - så vi omregner lige til km
74	dist = dist * 1.609344;
75	return (dist);
76	}
77
78	public static double beregnAfstand(GeoPoint p1, GeoPoint p2) {
79	return distanceHaversine(p1.latitude, p1.longitude, p2.latitude,
80	p2.longitude);
81	}
82
83	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
84
85	// denne er nok den mest præcise - men er også den langsomste
86	public static float beregnAfstand_google(GeoPoint p1, GeoPoint p2) {
87	float[] result = new float[1];
88
89	computeDistanceAndBearing(p1.latitude, p1.longitude, p2.latitude,
90	p2.longitude, result);
91
92	return (result[0] / 1000.0f);
93	}
94
95
96	// 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
102	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
109	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
114	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
119	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
126	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
135	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
152	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
173	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
182	double delta = (lambda - lambdaOrig) / lambda;
183	if (Math.abs(delta) < 1.0e-12) {
184	break;
185	}
186	}
187
188	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
204	/**
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
239	}