GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-a16a02f7eb
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 GIS Development Team
25  *
26  * \date 1999-2014
27  */
28 
29 #include <math.h>
30 #include <grass/gis.h>
31 #include "pi.h"
32 
33 static void adjust_lat(double *);
34 
35 #if 0
36 static void adjust_lon(double *);
37 #endif /* unused */
38 
39 static struct state {
40  double TAN_A, TAN1, TAN2, L;
41  int parallel;
42 } state;
43 
44 static 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  */
57 
58 int G_begin_rhumbline_equation(double lon1, double lat1, double lon2,
59  double lat2)
60 {
61  adjust_lat(&lat1);
62  adjust_lat(&lat2);
63 
64  if (lon1 == lon2) {
65  st->parallel = 1; /* a lie */
66  st->L = lat1;
67  return 0;
68  }
69  if (lat1 == lat2) {
70  st->parallel = 1;
71  st->L = lat1;
72  return 1;
73  }
74  st->parallel = 0;
75  lon1 = Radians(lon1);
76  lon2 = Radians(lon2);
77  lat1 = Radians(lat1);
78  lat2 = Radians(lat2);
79 
80  st->TAN1 = tan(M_PI_4 + lat1 / 2.0);
81  st->TAN2 = tan(M_PI_4 + lat2 / 2.0);
82  st->TAN_A = (lon2 - lon1) / (log(st->TAN2) - log(st->TAN1));
83  st->L = lon1;
84 
85  return 1;
86 }
87 
88 /**
89  * \brief Calculates rhumbline latitude.
90  *
91  * <b>Note:</b> Function only works if lon1 < lon < lon2.
92  *
93  * \param[in] lon longitude
94  * \return double latitude in degrees
95  */
96 
97 double G_rhumbline_lat_from_lon(double lon)
98 {
99  if (st->parallel)
100  return st->L;
101 
102  lon = Radians(lon);
103 
104  return Degrees(2 * atan(exp((lon - st->L) / st->TAN_A) * st->TAN1) -
105  M_PI_2);
106 }
107 
108 #if 0
109 static void adjust_lon(double *lon)
110 {
111  while (*lon > 180.0)
112  *lon -= 360.0;
113  while (*lon < -180.0)
114  *lon += 360.0;
115 }
116 #endif /* unused */
117 
118 static void adjust_lat(double *lat)
119 {
120  if (*lat > 90.0)
121  *lat = 90.0;
122  if (*lat < -90.0)
123  *lat = -90.0;
124 }
#define M_PI_2
Definition: gis.h:161
#define M_PI_4
Definition: gis.h:164
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:58
double G_rhumbline_lat_from_lon(double lon)
Calculates rhumbline latitude.
Definition: rhumbline.c:97