Loading...
Searching...
No Matches
rhumbline.c
Go to the documentation of this file.
/*!
* \file lib/gis/rhumbline.c
*
* \brief GIS Library - Rhumbline calculation routines.
*
* From "Map Projections" by Peter Richardus and Ron K. Alder, 1972<br>
* (526.8 R39m in Map & Geography Library)<br>
* Page 20,21, formulas 2.21, 2.22
*
* Formula is the equation of a rhumbline from (lat1,lon1) to
* (lat2,lon2). Input is lon, output is lat (all in degrees).
*
* <b>Note:</b> Formula only works if 0 < abs(lon2-lon1) < 180.
* If lon1 == lon2 then rhumbline is the merdian lon1 (and the formula
* will fail).
* <br>
* <b>WARNING:</b> This code is preliminary. It may not even be correct.
*
* (C) 2001-2014 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 GRASS Development Team
*
* \date 1999-2014
*/
#include <math.h>
#include <grass/gis.h>
#include "pi.h"
static void adjust_lat(double *);
#if 0
static void adjust_lon(double *);
#endif /* unused */
static struct state {
double TAN_A, TAN1, TAN2, L;
int parallel;
} state;
static struct state *st = &state;
/**
* \brief Start rhumbline calculations.
*
* <b>Note:</b> This function must be called before other rhumbline
* functions to initialize parameters.
*
* \param[in] lon1,lat1 longitude, latitude of first point
* \param[in] lon2,lat2 longitude, latitude of second point
* \return 1 on success
* \return 0 on error
*/
int G_begin_rhumbline_equation(double lon1, double lat1, double lon2,
double lat2)
{
adjust_lat(&lat1);
adjust_lat(&lat2);
if (lon1 == lon2) {
st->parallel = 1; /* a lie */
st->L = lat1;
return 0;
}
if (lat1 == lat2) {
st->parallel = 1;
st->L = lat1;
return 1;
}
st->parallel = 0;
lon1 = Radians(lon1);
lon2 = Radians(lon2);
lat1 = Radians(lat1);
lat2 = Radians(lat2);
st->TAN1 = tan(M_PI_4 + lat1 / 2.0);
st->TAN2 = tan(M_PI_4 + lat2 / 2.0);
st->TAN_A = (lon2 - lon1) / (log(st->TAN2) - log(st->TAN1));
st->L = lon1;
return 1;
}
/**
* \brief Calculates rhumbline latitude.
*
* <b>Note:</b> Function only works if lon1 < lon < lon2.
*
* \param[in] lon longitude
* \return double latitude in degrees
*/
double G_rhumbline_lat_from_lon(double lon)
{
if (st->parallel)
return st->L;
lon = Radians(lon);
return Degrees(2 * atan(exp((lon - st->L) / st->TAN_A) * st->TAN1) -
M_PI_2);
}
#if 0
static void adjust_lon(double *lon)
{
while (*lon > 180.0)
*lon -= 360.0;
while (*lon < -180.0)
*lon += 360.0;
}
#endif /* unused */
static void adjust_lat(double *lat)
{
if (*lat > 90.0)
*lat = 90.0;
if (*lat < -90.0)
*lat = -90.0;
}
pi.h
Degrees
#define Degrees(x)
Definition pi.h:7
Radians
#define Radians(x)
Definition pi.h:6
G_begin_rhumbline_equation
int G_begin_rhumbline_equation(double lon1, double lat1, double lon2, double lat2)
Start rhumbline calculations.
Definition rhumbline.c:57
G_rhumbline_lat_from_lon
double G_rhumbline_lat_from_lon(double lon)
Calculates rhumbline latitude.
Definition rhumbline.c:95