GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
geodesic.c
Go to the documentation of this file.
1 #include <math.h>
2 #include <grass/gis.h>
3 #include "pi.h"
4 
5 
6 /*
7  * This code is preliminary. I don't know if it is even
8  * correct.
9  */
10 
11 /*
12  * From "Map Projections" by Peter Richardus and Ron K. Alder, 1972
13  * (526.8 R39m in Map & Geography Library)
14  * page 19, formula 2.17
15  *
16  * Formula is the equation of a geodesic from (lat1,lon1) to (lat2,lon2)
17  * Input is lon, output is lat (all in degrees)
18  *
19  * Note formula only works if 0 < abs(lon2-lon1) < 180
20  * If lon1 == lon2 then geodesic is the merdian lon1
21  * (and the formula will fail)
22  * if lon2-lon1=180 then the geodesic is either meridian lon1 or lon2
23  */
24 
25 /* TODO:
26  *
27  * integrate code from raster/r.surf.idw/ll.c
28  */
29 
30 
31 static void adjust_lat(double *);
32 static void adjust_lon(double *);
33 
34 static struct state {
35  double A, B;
36 } state;
37 
38 static struct state *st = &state;
39 
40 int G_begin_geodesic_equation(double lon1, double lat1, double lon2,
41  double lat2)
42 {
43  double sin21, tan1, tan2;
44 
45  adjust_lon(&lon1);
46  adjust_lon(&lon2);
47  adjust_lat(&lat1);
48  adjust_lat(&lat2);
49  if (lon1 > lon2) {
50  double temp;
51  temp = lon1; lon1 = lon2; lon2 = temp;
52  temp = lat1; lat1 = lat2; lat2 = temp;
53  }
54  if (lon1 == lon2) {
55  st->A = st->B = 0.0;
56  return 0;
57  }
58  lon1 = Radians(lon1);
59  lon2 = Radians(lon2);
60  lat1 = Radians(lat1);
61  lat2 = Radians(lat2);
62 
63  sin21 = sin(lon2 - lon1);
64  tan1 = tan(lat1);
65  tan2 = tan(lat2);
66 
67  st->A = (tan2 * cos(lon1) - tan1 * cos(lon2)) / sin21;
68  st->B = (tan2 * sin(lon1) - tan1 * sin(lon2)) / sin21;
69 
70  return 1;
71 }
72 
73 /* only works if lon1 < lon < lon2 */
74 
75 double G_geodesic_lat_from_lon(double lon)
76 {
77  adjust_lon(&lon);
78  lon = Radians(lon);
79 
80  return Degrees(atan(st->A * sin(lon) - st->B * cos(lon)));
81 }
82 
83 static void adjust_lon(double *lon)
84 {
85  while (*lon > 180.0)
86  *lon -= 360.0;
87  while (*lon < -180.0)
88  *lon += 360.0;
89 }
90 
91 static void adjust_lat(double *lat)
92 {
93  if (*lat > 90.0)
94  *lat = 90.0;
95  if (*lat < -90.0)
96  *lat = -90.0;
97 }
int G_begin_geodesic_equation(double lon1, double lat1, double lon2, double lat2)
Definition: geodesic.c:40
struct state * st
Definition: parser.c:104
#define Radians(x)
Definition: pi.h:6
#define Degrees(x)
Definition: pi.h:7
struct state state
Definition: parser.c:103
double G_geodesic_lat_from_lon(double lon)
Definition: geodesic.c:75