2 * This file is part of mapper
4 * Copyright (C) 2007 Kaj-Michael Lang
6 * Original Algorithms from misc web sites
8 * Default map data provided by http://www.openstreetmap.org/
10 * This program is free software; you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation; either version 2 of the License, or
13 * (at your option) any later version.
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * Different Lat/Lon and related functions
27 * (conversion, distances, etc)
33 #include <glib/gi18n.h>
38 /* Int ranges for integerized lat/lon */
39 #define LATLON_MAX 2147483646
40 #define LATLON_MIN -2147483646
42 #define EARTH_RADIUS (3440.06479f)
47 DEG_FORMAT_TEXT[DDPDDDDD] = "-dd.ddddd°";
48 DEG_FORMAT_TEXT[DD_MMPMMM] = "-dd°mm.mmm'";
49 DEG_FORMAT_TEXT[DD_MM_SSPS] = "-dd°mm'ss.s\"";
50 DEG_FORMAT_TEXT[DDPDDDDD_NSEW] = "dd.ddddd° S";
51 DEG_FORMAT_TEXT[DD_MMPMMM_NSEW] = "dd°mm.mmm' S";
52 DEG_FORMAT_TEXT[DD_MM_SSPS_NSEW] = "dd°mm'ss.s\" S";
54 UNITS_TEXT[UNITS_KM] = _("km");
55 UNITS_TEXT[UNITS_MI] = _("mi.");
56 UNITS_TEXT[UNITS_NM] = _("n.m.");
60 magnitude(gdouble x1, gdouble y1, gdouble x2, gdouble y2)
66 return sqrt((x*x)+(y*y));
70 distance_point_to_line(gdouble x, gdouble y, gdouble x1, gdouble y1, gdouble x2, gdouble y2, gdouble *d)
75 lm=magnitude(x1,y1,x2,y2);
79 tmp=((x-x1)*(x2-x1))+((y-y1)*(y2-y1));
88 *d=magnitude(x,y, ix, iy);
94 * QuadTile calculation functions, from:
95 * http://trac.openstreetmap.org/browser/sites/rails_port/lib/quad_tile/quad_tile.h
98 xy2tile(guint x, guint y)
103 for (i=15;i>=0;i--) {
104 tile=(tile << 1) | ((x >> i) & 1);
105 tile=(tile << 1) | ((y >> i) & 1);
113 return round((lon + 180.0) * 65535.0 / 360.0);
119 return round((lat + 90.0) * 65535.0 / 180.0);
123 * Convert latitude to integerized+mercator projected value
126 lat2mp_int(gdouble lat)
128 return lat > 85.051128779 ? LATLON_MAX : lat < -85.051128779 ? LATLON_MIN :
129 lrint(log(tan(M_PI_4l+lat*M_PIl/360))/M_PIl*LATLON_MAX);
133 * Convert longitude to integerized+mercator projected value
136 lon2mp_int(gdouble lon)
138 return lrint(lon/180*LATLON_MAX);
144 return (gdouble)lon/LATLON_MAX*180;
154 * Quick distance calculation, (does not handle earth curvature, so use for quick non exact needs only)
157 calculate_idistance(gint lat1, gint lon1, gint lat2, gint lon2)
159 return lrint(sqrt((double)((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2))));
163 calculate_ddistance(gint lat1, gint lon1, gint lat2, gint lon2)
165 return sqrt((double)((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2)));
169 * Quick distance for comparing, skips the final square root
172 calculate_idistance_cmp(gint lat1, gint lon1, gint lat2, gint lon2)
174 return ((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2));
178 * Calculate the distance between two lat/lon pairs.
179 * The distance is returned in kilometers.
183 calculate_distance(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
185 gdouble dlat, dlon, slat, slon, a;
187 /* Convert to radians. */
188 lat1*=(M_PIl / 180.f);
189 lon1*=(M_PIl / 180.f);
190 lat2*=(M_PIl / 180.f);
191 lon2*=(M_PIl / 180.f);
196 slat=sin(dlat / 2.f);
197 slon=sin(dlon / 2.f);
199 a=(slat * slat) + (cos(lat1) * cos(lat2) * slon * slon);
201 return ((2.f * atan2(sqrt(a), sqrt(1.f - a))) * EARTH_RADIUS);
205 * Calculate course from lat1,lon1 to lat2,lon2
209 calculate_course_rad(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
213 lat1*=(M_PIl / 180.f);
214 lon1*=(M_PIl / 180.f);
215 lat2*=(M_PIl / 180.f);
216 lon2*=(M_PIl / 180.f);
219 y=sin(dlon)*cos(lat2);
220 x=cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dlon);
227 calculate_course(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
231 c=calculate_course_rad(lat1, lon1, lat2, lon2);
239 deg_format(DegFormat degformat, gdouble coor, gchar * scoor, gchar neg_char, gchar pos_char)
242 gdouble acoor=fabs(coor);
246 g_sprintf(scoor, "%.5f°", coor);
249 g_sprintf(scoor, "%.5f° %c", acoor,
250 coor < 0.f ? neg_char : pos_char);
254 g_sprintf(scoor, "%d°%06.3f'",
255 (int)coor, (acoor - (int)acoor) * 60.0);
258 g_sprintf(scoor, "%d°%06.3f' %c",
259 (int)acoor, (acoor - (int)acoor) * 60.0,
260 coor < 0.f ? neg_char : pos_char);
264 min = (acoor - (int)acoor) * 60.0;
265 g_sprintf(scoor, "%d°%02d'%04.1f\"", (int)coor, (int)min,
266 ((min - (int)min) * 60.0));
268 case DD_MM_SSPS_NSEW:
269 min = (acoor - (int)acoor) * 60.0;
270 g_sprintf(scoor, "%d°%02d'%04.1f\" %c", (int)acoor, (int)min,
271 ((min - (int)min) * 60.0),
272 coor < 0.f ? neg_char : pos_char);