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