GRASS GIS 7 Programmer's Manual  7.0.2(2015)-r00000
geodist.c
Go to the documentation of this file.
1 
23 #include <math.h>
24 #include <grass/gis.h>
25 #include "pi.h"
26 
27 static struct state {
28  double boa;
29  double f;
30  double ff64;
31  double al;
32  double t1, t2, t3, t4, t1r, t2r;
33 } state;
34 
35 static struct state *st = &state;
36 
50 void G_begin_geodesic_distance(double a, double e2)
51 {
52  st->al = a;
53  st->boa = sqrt(1 - e2);
54  st->f = 1 - st->boa;
55  st->ff64 = st->f * st->f / 64;
56 }
57 
70 {
71  st->t1r = atan(st->boa * tan(Radians(lat1)));
72 }
73 
84 {
85  double stm, ctm, sdtm, cdtm;
86  double tm, dtm;
87 
88  st->t2r = atan(st->boa * tan(Radians(lat2)));
89 
90  tm = (st->t1r + st->t2r) / 2;
91  dtm = (st->t2r - st->t1r) / 2;
92 
93  stm = sin(tm);
94  ctm = cos(tm);
95  sdtm = sin(dtm);
96  cdtm = cos(dtm);
97 
98  st->t1 = stm * cdtm;
99  st->t1 = st->t1 * st->t1 * 2;
100 
101  st->t2 = sdtm * ctm;
102  st->t2 = st->t2 * st->t2 * 2;
103 
104  st->t3 = sdtm * sdtm;
105  st->t4 = cdtm * cdtm - stm * stm;
106 }
107 
121 double G_geodesic_distance_lon_to_lon(double lon1, double lon2)
122 {
123  double a, cd, d, e, /*dl, */
124  q, sd, sdlmr, t, u, v, x, y;
125 
126 
127  sdlmr = sin(Radians(lon2 - lon1) / 2);
128 
129  /* special case - shapiro */
130  if (sdlmr == 0.0 && st->t1r == st->t2r)
131  return 0.0;
132 
133  q = st->t3 + sdlmr * sdlmr * st->t4;
134 
135  /* special case - shapiro */
136  if (q == 1.0)
137  return M_PI * st->al;
138 
139  /* Mod: shapiro
140  * cd=1-2q is ill-conditioned if q is small O(10**-23)
141  * (for high lats? with lon1-lon2 < .25 degrees?)
142  * the computation of cd = 1-2*q will give cd==1.0.
143  * However, note that t=dl/sd is dl/sin(dl) which approaches 1 as dl->0.
144  * So the first step is to compute a good value for sd without using sin()
145  * and then check cd && q to see if we got cd==1.0 when we shouldn't.
146  * Note that dl isn't used except to get t,
147  * but both cd and sd are used later
148  */
149 
150  /* original code
151  cd=1-2*q;
152  dl=acos(cd);
153  sd=sin(dl);
154  t=dl/sd;
155  */
156 
157  cd = 1 - 2 * q; /* ill-conditioned subtraction for small q */
158  /* mod starts here */
159  sd = 2 * sqrt(q - q * q); /* sd^2 = 1 - cd^2 */
160  if (q != 0.0 && cd == 1.0) /* test for small q */
161  t = 1.0;
162  else if (sd == 0.0)
163  t = 1.0;
164  else
165  t = acos(cd) / sd; /* don't know how to fix acos(1-2*q) yet */
166  /* mod ends here */
167 
168  u = st->t1 / (1 - q);
169  v = st->t2 / q;
170  d = 4 * t * t;
171  x = u + v;
172  e = -2 * cd;
173  y = u - v;
174  a = -d * e;
175 
176  return st->al * sd * (t
177  - st->f / 4 * (t * x - y)
178  + st->ff64 * (x * (a + (t - (a + e) / 2) * x)
179  + y * (-2 * d + e * y) + d * x * y));
180 }
181 
196 double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
197 {
200  return G_geodesic_distance_lon_to_lon(lon1, lon2);
201 }
double G_geodesic_distance_lon_to_lon(double lon1, double lon2)
Calculates geodesic distance.
Definition: geodist.c:121
void G_begin_geodesic_distance(double a, double e2)
Begin geodesic distance.
Definition: geodist.c:50
struct state * st
Definition: parser.c:101
void G_set_geodesic_distance_lat2(double lat2)
Sets geodesic distance lat2.
Definition: geodist.c:83
double t
#define Radians(x)
Definition: pi.h:6
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition: geodist.c:196
void G_set_geodesic_distance_lat1(double lat1)
Sets geodesic distance lat1.
Definition: geodist.c:69
struct state state
Definition: parser.c:100