/*
- * Different Lat/Lon related function (conversion, distances, etc)
+ * This file is part of mapper
+ *
+ * Copyright (C) 2007 Kaj-Michael Lang
+ *
+ * Original Algorithms from misc web sites
+ *
+ * Default map data provided by http://www.openstreetmap.org/
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ */
+
+/*
+ * Different Lat/Lon and related functions
+ * (conversion, distances, etc)
*/
-#define _GNU_SOURCE
+#include "config.h"
#include <math.h>
#include <glib.h>
+#include <glib/gi18n.h>
+#include <glib/gstdio.h>
+#include "latlon.h"
#include "osm.h"
-/* Int ranges for integerized lat/lon */
-#define LATLON_MAX 2147483646
-#define LATLON_MIN -2147483646
+void
+latlon_init(void)
+{
+DEG_FORMAT_TEXT[DDPDDDDD] = "-dd.ddddd°";
+DEG_FORMAT_TEXT[DD_MMPMMM] = "-dd°mm.mmm'";
+DEG_FORMAT_TEXT[DD_MM_SSPS] = "-dd°mm'ss.s\"";
+DEG_FORMAT_TEXT[DDPDDDDD_NSEW] = "dd.ddddd° S";
+DEG_FORMAT_TEXT[DD_MMPMMM_NSEW] = "dd°mm.mmm' S";
+DEG_FORMAT_TEXT[DD_MM_SSPS_NSEW] = "dd°mm'ss.s\" S";
+
+UNITS_TEXT[UNITS_KM] = _("km");
+UNITS_TEXT[UNITS_MI] = _("mi.");
+UNITS_TEXT[UNITS_NM] = _("n.m.");
+}
+
+gfloat
+angle_diff(gfloat a, gfloat b)
+{
+gfloat tmax, tmin, t;
-#define EARTH_RADIUS (3440.06479f)
+tmax=MAX(a,b);
+tmin=MIN(a,b);
-/* Convert latitude to integerized+mercator projected value */
+t=tmax-tmin;
+if (t>180) t-=360;
+
+return t;
+}
+
+static inline gdouble
+magnitude(gdouble x1, gdouble y1, gdouble x2, gdouble y2)
+{
+gdouble x,y;
+x=x2-x1;
+y=y2-y1;
+
+return sqrt((x*x)+(y*y));
+}
+
+gboolean
+distance_point_to_line(gdouble x, gdouble y, gdouble x1, gdouble y1, gdouble x2, gdouble y2, gdouble *d)
+{
+gdouble lm,u,tmp;
+gdouble ix,iy;
+
+lm=magnitude(x1,y1,x2,y2);
+if (lm==0.0f)
+ return FALSE;
+
+tmp=((x-x1)*(x2-x1))+((y-y1)*(y2-y1));
+u=tmp/(lm*lm);
+
+if (u<0.0f || u>1.0f)
+ return FALSE;
+
+ix=x1+u*(x2-x1);
+iy=y1+u*(y2-y1);
+
+*d=magnitude(x,y, ix, iy);
+
+return TRUE;
+}
+
+/*
+ * QuadTile calculation functions, from:
+ * http://trac.openstreetmap.org/browser/sites/rails_port/lib/quad_tile/quad_tile.h
+ */
+guint32
+xy2tile(guint x, guint y)
+{
+guint32 tile=0;
+gint i;
+
+for (i=15;i>=0;i--) {
+ tile=(tile << 1) | ((x >> i) & 1);
+ tile=(tile << 1) | ((y >> i) & 1);
+}
+return tile;
+}
+
+guint32
+lon2x(gdouble lon)
+{
+return round((lon + 180.0) * 65535.0 / 360.0);
+}
+
+guint32
+lat2y(gdouble lat)
+{
+return round((lat + 90.0) * 65535.0 / 180.0);
+}
+
+/*
+ * Convert latitude to integerized+mercator projected value
+ */
gint32
lat2mp_int(gdouble lat)
{
lrint(log(tan(M_PI_4l+lat*M_PIl/360))/M_PIl*LATLON_MAX);
}
-/* Convert longitude to integerized+mercator projected value */
+/*
+ * Convert longitude to integerized+mercator projected value
+ */
gint32
lon2mp_int(gdouble lon)
{
}
gdouble
-mp_int2lon(gint lon)
+mp_int2lon(gint32 lon)
{
return (gdouble)lon/LATLON_MAX*180;
}
gdouble
-mp_int2lat(gint lat)
+mp_int2lat(gint32 lat)
{
return 0;
}
+/**
+ * Quick distance calculation, (does not handle earth curvature, so use for quick non exact needs only)
+ */
gulong
calculate_idistance(gint lat1, gint lon1, gint lat2, gint lon2)
{
return lrint(sqrt((double)((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2))));
}
+gdouble
+calculate_ddistance(gint lat1, gint lon1, gint lat2, gint lon2)
+{
+return sqrt((double)((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2)));
+}
+
/**
* Quick distance for comparing, skips the final square root
*/
gdouble dlat, dlon, slat, slon, a;
/* Convert to radians. */
-lat1*=(M_PI_4l / 180.f);
-lon1*=(M_PI_4l / 180.f);
-lat2*=(M_PI_4l / 180.f);
-lon2*=(M_PI_4l / 180.f);
+lat1*=(M_PIl / 180.f);
+lon1*=(M_PIl / 180.f);
+lat2*=(M_PIl / 180.f);
+lon2*=(M_PIl / 180.f);
dlat=lat2 - lat1;
dlon=lon2 - lon1;
return ((2.f * atan2(sqrt(a), sqrt(1.f - a))) * EARTH_RADIUS);
}
+/**
+ * Calculate course from lat1,lon1 to lat2,lon2
+ *
+ */
+gdouble
+calculate_course_rad(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
+{
+gdouble c,dlon,x,y;
+
+lat1*=(M_PIl / 180.f);
+lon1*=(M_PIl / 180.f);
+lat2*=(M_PIl / 180.f);
+lon2*=(M_PIl / 180.f);
+
+dlon=lon2-lon1;
+y=sin(dlon)*cos(lat2);
+x=cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dlon);
+c=atan2(y,x);
+
+return c;
+}
+
+gdouble
+calculate_course(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
+{
+gdouble c;
+
+c=calculate_course_rad(lat1, lon1, lat2, lon2);
+
+c*=180.f/M_PIl;
+c=c+360 % 360;
+return c;
+}
+
+void
+deg_format(DegFormat degformat, gdouble coor, gchar * scoor, gchar neg_char, gchar pos_char)
+{
+gdouble min;
+gdouble acoor=fabs(coor);
+
+switch (degformat) {
+ case DDPDDDDD:
+ g_sprintf(scoor, "%.5f°", coor);
+ break;
+ case DDPDDDDD_NSEW:
+ g_sprintf(scoor, "%.5f° %c", acoor, coor < 0.f ? neg_char : pos_char);
+ break;
+ case DD_MMPMMM:
+ g_sprintf(scoor, "%d°%06.3f'", (int)coor, (acoor - (int)acoor) * 60.0);
+ break;
+ case DD_MMPMMM_NSEW:
+ g_sprintf(scoor, "%d°%06.3f' %c", (int)acoor, (acoor - (int)acoor) * 60.0, coor < 0.f ? neg_char : pos_char);
+ break;
+ case DD_MM_SSPS:
+ min = (acoor - (int)acoor) * 60.0;
+ g_sprintf(scoor, "%d°%02d'%04.1f\"", (int)coor, (int)min, ((min - (int)min) * 60.0));
+ break;
+ case DD_MM_SSPS_NSEW:
+ min = (acoor - (int)acoor) * 60.0;
+ g_sprintf(scoor, "%d°%02d'%04.1f\" %c", (int)acoor, (int)min, ((min - (int)min) * 60.0), coor < 0.f ? neg_char : pos_char);
+ break;
+ default:
+ g_return_if_reached();
+ break;
+}
+}
+
projects / mapper / blobdiff