GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
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
31static void adjust_lat(double *);
32static void adjust_lon(double *);
33
34static struct state {
35 double A, B;
36} state;
37
38static struct state *st = &state;
39
40int 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
75double 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
83static 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
91static 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
double G_geodesic_lat_from_lon(double lon)
Definition: geodesic.c:75
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