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