GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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
29static void adjust_lat(double *);
30static void adjust_lon(double *);
31
32static struct state {
33 double A, B;
34} state;
35
36static struct state *st = &state;
37
38int 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
78double 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
86static 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
94static 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
#define Degrees(x)
Definition pi.h:7
#define Radians(x)
Definition pi.h:6