Loading...
Searching...
No Matches
geodist.c
Go to the documentation of this file.
/*!
* \file lib/gis/geodist.c
*
* \brief GIS Library - Geodesic distance routines.
*
* Distance from point to point along a geodesic code from Paul
* D. Thomas, 1970 "Spheroidal Geodesics, Reference Systems, and Local
* Geometry" U.S. Naval Oceanographic Office, p. 162 Engineering
* Library 526.3 T36s
* http://stinet.dtic.mil/oai/oai?&verb=getRecord&metadataPrefix=html&identifier=AD0703541
*
* <b>WARNING:</b> this code is preliminary and may be changed,
* including calling sequences to any of the functions defined here.
*
* (C) 2001-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public License
* (>=v2). Read the file COPYING that comes with GRASS for details.
*
* \author Original author CERL
*/
#include <math.h>
#include <grass/gis.h>
#include "pi.h"
static struct state {
double boa;
double f;
double ff64;
double al;
double t1, t2, t3, t4, t1r, t2r;
} state;
static struct state *st = &state;
/*!
* \brief Begin geodesic distance.
*
* Initializes the distance calculations for the ellipsoid with
* semi-major axis <i>a</i> (in meters) and ellipsoid eccentricity squared
* <i>e2</i>. It is used only for the latitude-longitude projection.
*
* <b>Note:</b> Must be called once to establish the ellipsoid.
*
* \param a semi-major axis in meters
* \param e2 ellipsoid eccentricity
*/
void G_begin_geodesic_distance(double a, double e2)
{
st->al = a;
st->boa = sqrt(1 - e2);
st->f = 1 - st->boa;
st->ff64 = st->f * st->f / 64;
}
/*!
* \brief Sets geodesic distance lat1.
*
* Set the first latitude.
*
* <b>Note:</b> Must be called first.
*
* \param[in] lat1 first latitude
* \return
*/
void G_set_geodesic_distance_lat1(double lat1)
{
st->t1r = atan(st->boa * tan(Radians(lat1)));
}
/*!
* \brief Sets geodesic distance lat2.
*
* Set the second latitude.
*
* <b>Note:</b> Must be called second.
*
* \param[in] lat2 second latitude
*/
void G_set_geodesic_distance_lat2(double lat2)
{
double stm, ctm, sdtm, cdtm;
double tm, dtm;
st->t2r = atan(st->boa * tan(Radians(lat2)));
tm = (st->t1r + st->t2r) / 2;
dtm = (st->t2r - st->t1r) / 2;
stm = sin(tm);
ctm = cos(tm);
sdtm = sin(dtm);
cdtm = cos(dtm);
st->t1 = stm * cdtm;
st->t1 = st->t1 * st->t1 * 2;
st->t2 = sdtm * ctm;
st->t2 = st->t2 * st->t2 * 2;
st->t3 = sdtm * sdtm;
st->t4 = cdtm * cdtm - stm * stm;
}
/*!
* \brief Calculates geodesic distance.
*
* Calculates the geodesic distance from <i>lon1,lat1</i> to
* <i>lon2,lat2</i> in meters where <i>lat1</i> was the latitude
* passed to G_set_geodesic_distance_latl() and <i>lat2</i> was the
* latitude passed to G_set_geodesic_distance_lat2().
*
* \param[in] lon1 first longitude
* \param[in] lon2 second longitude
*
* \return double distance in meters
*/
double G_geodesic_distance_lon_to_lon(double lon1, double lon2)
{
double a, cd, d, e, /*dl, */
q, sd, sdlmr, t, u, v, x, y;
sdlmr = sin(Radians(lon2 - lon1) / 2);
/* special case - shapiro */
if (sdlmr == 0.0 && st->t1r == st->t2r)
return 0.0;
q = st->t3 + sdlmr * sdlmr * st->t4;
/* special case - shapiro */
if (q == 1.0)
return M_PI * st->al;
/* Mod: shapiro
* cd=1-2q is ill-conditioned if q is small O(10**-23)
* (for high lats? with lon1-lon2 < .25 degrees?)
* the computation of cd = 1-2*q will give cd==1.0.
* However, note that t=dl/sd is dl/sin(dl) which approaches 1 as dl->0.
* So the first step is to compute a good value for sd without using sin()
* and then check cd && q to see if we got cd==1.0 when we shouldn't.
* Note that dl isn't used except to get t,
* but both cd and sd are used later
*/
/* original code
cd=1-2*q;
dl=acos(cd);
sd=sin(dl);
t=dl/sd;
*/
cd = 1 - 2 * q; /* ill-conditioned subtraction for small q */
/* mod starts here */
sd = 2 * sqrt(q - q * q); /* sd^2 = 1 - cd^2 */
if (q != 0.0 && cd == 1.0) /* test for small q */
t = 1.0;
else if (sd == 0.0)
t = 1.0;
else
t = acos(cd) / sd; /* don't know how to fix acos(1-2*q) yet */
/* mod ends here */
u = st->t1 / (1 - q);
v = st->t2 / q;
d = 4 * t * t;
x = u + v;
e = -2 * cd;
y = u - v;
a = -d * e;
return st->al * sd *
(t - st->f / 4 * (t * x - y) +
st->ff64 * (x * (a + (t - (a + e) / 2) * x) + y * (-2 * d + e * y) +
d * x * y));
}
/*!
* \brief Calculates geodesic distance.
*
* Calculates the geodesic distance from <i>lon1,lat1</i> to
* <i>lon2,lat2</i> in meters.
*
* <b>Note:</b> The calculation of the geodesic distance is fairly
* costly.
*
* \param[in] lon1,lat1 longitude,latitude of first point
* \param[in] lon2,lat2 longitude,latitude of second point
*
* \return distance in meters
*/
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
{
G_set_geodesic_distance_lat1(lat1);
G_set_geodesic_distance_lat2(lat2);
return G_geodesic_distance_lon_to_lon(lon1, lon2);
}
t
double t
Definition driver/set_window.c:5
G_set_geodesic_distance_lat2
void G_set_geodesic_distance_lat2(double lat2)
Sets geodesic distance lat2.
Definition geodist.c:81
G_set_geodesic_distance_lat1
void G_set_geodesic_distance_lat1(double lat1)
Sets geodesic distance lat1.
Definition geodist.c:67
G_geodesic_distance
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition geodist.c:193
G_geodesic_distance_lon_to_lon
double G_geodesic_distance_lon_to_lon(double lon1, double lon2)
Calculates geodesic distance.
Definition geodist.c:119
G_begin_geodesic_distance
void G_begin_geodesic_distance(double a, double e2)
Begin geodesic distance.
Definition geodist.c:49
pi.h
Radians
#define Radians(x)
Definition pi.h:6
x
#define x