GRASS 8 Programmer's Manual 8.6.0dev(2026)-f6f2c534ea
Loading...
Searching...
No Matches
rhumbline.c
Go to the documentation of this file.
1/*!
2 * \file lib/gis/rhumbline.c
3 *
4 * \brief GIS Library - Rhumbline calculation routines.
5 *
6 * From "Map Projections" by Peter Richardus and Ron K. Alder, 1972<br>
7 * (526.8 R39m in Map & Geography Library)<br>
8 * Page 20,21, formulas 2.21, 2.22
9 *
10 * Formula is the equation of a rhumbline from (lat1,lon1) to
11 * (lat2,lon2). Input is lon, output is lat (all in degrees).
12 *
13 * <b>Note:</b> Formula only works if 0 < abs(lon2-lon1) < 180.
14 * If lon1 == lon2 then rhumbline is the merdian lon1 (and the formula
15 * will fail).
16 * <br>
17 * <b>WARNING:</b> This code is preliminary. It may not even be correct.
18 *
19 * (C) 2001-2014 by the GRASS Development Team
20 *
21 * This program is free software under the GNU General Public License
22 * (>=v2). Read the file COPYING that comes with GRASS for details.
23 *
24 * \author GRASS Development Team
25 *
26 * \date 1999-2014
27 */
28
29#include <math.h>
30#include <grass/gis.h>
31#include "pi.h"
32
33static void adjust_lat(double *);
34
35#if 0
36static void adjust_lon(double *);
37#endif /* unused */
38
39static struct state {
40 double TAN_A, TAN1, TAN2, L;
41 int parallel;
42} state;
43
44static struct state *st = &state;
45
46/**
47 * \brief Start rhumbline calculations.
48 *
49 * <b>Note:</b> This function must be called before other rhumbline
50 * functions to initialize parameters.
51 *
52 * \param[in] lon1,lat1 longitude, latitude of first point
53 * \param[in] lon2,lat2 longitude, latitude of second point
54 * \return 1 on success
55 * \return 0 on error
56 */
57int G_begin_rhumbline_equation(double lon1, double lat1, double lon2,
58 double lat2)
59{
60 adjust_lat(&lat1);
61 adjust_lat(&lat2);
62
63 if (lon1 == lon2) {
64 st->parallel = 1; /* a lie */
65 st->L = lat1;
66 return 0;
67 }
68 if (lat1 == lat2) {
69 st->parallel = 1;
70 st->L = lat1;
71 return 1;
72 }
73 st->parallel = 0;
74 lon1 = Radians(lon1);
75 lon2 = Radians(lon2);
76 lat1 = Radians(lat1);
77 lat2 = Radians(lat2);
78
79 st->TAN1 = tan(M_PI_4 + lat1 / 2.0);
80 st->TAN2 = tan(M_PI_4 + lat2 / 2.0);
81 st->TAN_A = (lon2 - lon1) / (log(st->TAN2) - log(st->TAN1));
82 st->L = lon1;
83
84 return 1;
85}
86
87/**
88 * \brief Calculates rhumbline latitude.
89 *
90 * <b>Note:</b> Function only works if lon1 < lon < lon2.
91 *
92 * \param[in] lon longitude
93 * \return double latitude in degrees
94 */
96{
97 if (st->parallel)
98 return st->L;
99
100 lon = Radians(lon);
101
102 return Degrees(2 * atan(exp((lon - st->L) / st->TAN_A) * st->TAN1) -
103 M_PI_2);
104}
105
106#if 0
107static void adjust_lon(double *lon)
108{
109 while (*lon > 180.0)
110 *lon -= 360.0;
111 while (*lon < -180.0)
112 *lon += 360.0;
113}
114#endif /* unused */
115
116static void adjust_lat(double *lat)
117{
118 if (*lat > 90.0)
119 *lat = 90.0;
120 if (*lat < -90.0)
121 *lat = -90.0;
122}
#define M_PI_2
Definition gis.h:160
#define M_PI_4
Definition gis.h:168
struct state state
Definition parser.c:103
struct state * st
Definition parser.c:104
#define Degrees(x)
Definition pi.h:7
#define Radians(x)
Definition pi.h:6
int G_begin_rhumbline_equation(double lon1, double lat1, double lon2, double lat2)
Start rhumbline calculations.
Definition rhumbline.c:57
double G_rhumbline_lat_from_lon(double lon)
Calculates rhumbline latitude.
Definition rhumbline.c:95