-
Notifications
You must be signed in to change notification settings - Fork 0
/
geom2.c
326 lines (300 loc) · 11.7 KB
/
geom2.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
// src/aksl/geom2.c 2017-10-25 Alan U. Kennington.
/*-----------------------------------------------------------------------------
Copyright (C) 1989-2018, Alan U. Kennington.
You may distribute this software under the terms of Alan U. Kennington's
modified Artistic Licence, as specified in the accompanying LICENCE file.
-----------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------
Functions in this file:
dms2deg
E_dist2
// E_dist
E_dist2_lin
E_dist_lin
// translate
int_rect::
print
real_rect::
positivize
view
rect_int16::
------------------------------------------------------------------------------*/
#include "aksl/geom2.h"
/*------------------------------------------------------------------------------
dms2deg() converts an integer in the form d*10000 + m*100 + s to degrees.
------------------------------------------------------------------------------*/
//----------------------//
// dms2deg //
//----------------------//
double dms2deg(long i) {
double r = (i%100);
i /= 100;
r /= 60;
r += (i%100);
i /= 100;
r /= 60;
return r + i;
} // End of function dms2deg.
/*------------------------------------------------------------------------------
Function E_dist2() calculates distance between two points on the earth,
given the coordinates of the two points, and the radius of the earth
(a "global" constant). Distances in metres, angles in radians!
The formula is (probably):
d = radius_E * sqrt(2*(1 - cos(theta1 - theta0)
+ (1 - cos(phi1 - phi0))*cos(theta0)*cos(theta1)));
------------------------------------------------------------------------------*/
//----------------------//
// E_dist2 //
//----------------------//
double E_dist2(double phi0, double theta0, double phi1, double theta1) {
register double t0 = (1 - cos(phi1 - phi0))*cos(theta1)*cos(theta0);
register double t1 = 1 - cos(theta1 - theta0);
t0 += t1;
t0 += t0;
return R_earth_mean_2 * t0;
} // End of function E_dist2.
/*------------------------------------------------------------------------------
Function E_dist() calculates distance between two points on the earth,
given the coordinates of the two points, and the radius of the earth
(a "global" constant). Distances in metres, angles in radians!
The formula is (probably):
d = radius_E * sqrt(2*(1 - cos(theta1 - theta0)
+ (1 - cos(phi1 - phi0))*cos(theta0)*cos(theta1)));
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
This is commented out because there is an inline function E_dist() which simply
invokes E_dist2, and which therefore avoids duplication of code.
------------------------------------------------------------------------------*/
//----------------------//
// E_dist //
//----------------------//
/*------------------------------------------------------------------------------
double E_dist(double phi0, double theta0, double phi1, double theta1) {
register double t0 = (1 - cos(phi1 - phi0))*cos(theta1)*cos(theta0);
register double t1 = 1 - cos(theta1 - theta0);
t0 += t1;
t0 += t0;
return R_earth_mean * sqrt(t0);
} // End of function E_dist.
------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------
Function E_dist2_lin() calculates the approximate distance between two points
on the earth, given the coordinates of the two points, and the radius of the
earth (a "global" constant). Distances in metres, angles in radians!
The formula for the linearised distance is (probably):
d = radius_E *
sqrt((cos((theta0 + theta1)/2)*(phi1 - phi0))^2 + (theta1 - theta0)^2);
------------------------------------------------------------------------------*/
//----------------------//
// E_dist2_lin //
//----------------------//
double E_dist2_lin(double phi0, double theta0, double phi1, double theta1) {
register double t0 = (phi1 - phi0)*cos((theta0 + theta1)/2);
register double t1 = theta1 - theta0;
return R_earth_mean_2 * (t0*t0 + t1*t1);
} // End of function E_dist2_lin.
//----------------------//
// E_dist_lin //
//----------------------//
double E_dist_lin(double phi0, double theta0, double phi1, double theta1) {
register double t0 = (phi1 - phi0)*cos((theta0 + theta1)/2);
register double t1 = theta1 - theta0;
return R_earth_mean * hypot(t0, t1);
} // End of function E_dist_lin.
/*------------------------------------------------------------------------------
Function translate() calculates the new coordinates of a point, given that
it has moved a given distance in a given direction on the earth.
Distances in metres, angles in radians! The bearing alpha is measured with east
as 0 radians, and north as PI/2 radians.
------------------------------------------------------------------------------*/
/*
//----------------------//
// translate //
//----------------------//
void translate(double phi0, double theta0, double delta, double alpha,
double& phi1, double& theta1) {
} // End of function translate.
*/
//----------------------//
// int_rect::print //
//----------------------//
void int_rect::print(ostream& os) {
os << "(x, y, w, h) = (" << x << ", " << y
<< ", " << w << ", " << h << ")\n";
} // End of function int_rect::print.
//--------------------------//
// real_rect::positivize //
//--------------------------//
void real_rect::positivize() {
if (w < 0) {
x += w;
w = -w;
}
if (h < 0) {
y += h;
h = -h;
}
} // End of function real_rect::positivize.
/*------------------------------------------------------------------------------
This function returns a "position2" structure if the given point lies within the
real rectangle. If the rectangle has zero width or height and the point lies in
the zero-area region, then it is still considered to be in the region.
If w or h is negative, then the region is considered to be empty.
If the rectangle parameter is non-null, the returned point is scaled to lie in
that rectangle. Otherwise it lies in [0,1] x [0,1].
------------------------------------------------------------------------------*/
//----------------------//
// real_rect::view //
//----------------------//
position2* real_rect::view(position2* pp, real_rect* rr) {
if (!pp)
return 0;
if (pp->x < x || pp->x > x + w || pp->y < y || pp->y > y + h)
return 0;
double x_view = pp->x - x;
double y_view = pp->y - y;
if (w == 0)
x_view = 0;
else
x_view /= w;
if (h == 0)
y_view = 0;
else
y_view /= h;
if (rr) {
x_view *= rr->w;
x_view += rr->x;
y_view *= rr->h;
y_view += rr->y;
}
position2* pp2 = new position2(x_view, y_view);
return pp2;
} // End of function real_rect::view.
/*------------------------------------------------------------------------------
This returns the coords of the intersections of a given line with
the rectangle. This is intended for use in drawing the line segment which is
created by the intersection of the line with the rectangle.
If there is no intersection, then return a negative number.
If there is an intersection, then 0 is returned, and the intersection
points are returned in pt1 and pt2.
This function tries to avoid many kinds of integer overflow conditions.
The line pt1 to pt2 always goes from left to right, except that a vertical line
goes from bottom to top.
------------------------------------------------------------------------------*/
//----------------------//
// rect_int16::cuts //
//----------------------//
int rect_int16::cuts(const pt_int16& pt0, int16 kx, int16 ky,
pt_int16& pt1, pt_int16& pt2) const {
// The case of a horizontal line.
if (kx == 0) {
if (pt0.y < ymin)
return -1;
if (pt0.y > ymax)
return -2;
pt1.x = xmin; pt1.y = pt0.y;
pt2.x = xmax; pt2.y = pt0.y;
return 0;
}
// The case of a vertical line.
if (ky == 0) {
if (pt0.x < xmin)
return -3;
if (pt0.x > xmax)
return -4;
pt1.x = pt0.x; pt1.y = ymin;
pt2.x = pt0.x; pt2.y = ymax;
return 0;
}
// Normalize the slope vector.
if (kx < 0) {
kx = -kx;
ky = -ky;
}
// The case of positive slope.
if (ky > 0) {
// Find intersection with bottom edge.
// (Note excessively zealous casting.)
int32 x1 = (int32)pt0.x
+ (((int32)ymin - pt0.y) * (int32)ky)/(int32)kx;
if (x1 > (int32)xmax)
return -5;
// Find intersection with top edge.
int32 x2 = (int32)pt0.x
+ (((int32)ymax - pt0.y) * (int32)ky)/(int32)kx;
if (x2 < (int32)xmin)
return -6;
// At this point, the line definitely intersects the rectangle.
if (x1 >= (int32)xmin) {
// Intersects the bottom edge.
if (x2 <= (int32)xmax) {
// Intersects the top edge too.
pt1.x = (int16)x1; pt1.y = ymin;
pt2.x = (int16)x2; pt2.y = ymax;
return 0;
}
// Intersects the right edge.
int32 y2 = (int32)pt0.y
+ (((int32)xmax - pt0.x) * (int32)kx)/(int32)ky;
pt1.x = (int16)x1; pt1.y = ymin;
pt2.x = xmax; pt2.y = (int16)y2;
return 0;
}
// Intersects the left edge. (x1 < (int32)xmin.)
int32 y1 = (int32)pt0.y
+ (((int32)xmin - pt0.x) * (int32)kx)/(int32)ky;
if (x2 <= (int32)xmax) {
// Intersects the top edge.
pt1.x = xmin; pt1.y = (int16)y1;
pt2.x = (int16)x2; pt2.y = ymax;
return 0;
}
// Intersects the right edge.
int32 y2 = (int32)pt0.y
+ (((int32)xmax - pt0.x) * (int32)kx)/(int32)ky;
pt1.x = xmin; pt1.y = (int16)y1;
pt2.x = xmax; pt2.y = (int16)y2;
return 0;
}
// The case of negative slope.
// Find intersection with bottom edge.
int32 x1 = (int32)pt0.x
+ (((int32)ymin - pt0.y) * (int32)ky)/(int32)kx;
if (x1 < (int32)xmin)
return -7;
// Find intersection with top edge.
int32 x2 = (int32)pt0.x
+ (((int32)ymax - pt0.y) * (int32)ky)/(int32)kx;
if (x2 > (int32)xmax)
return -8;
// At this point, the line definitely intersects the rectangle.
if (x1 <= (int32)xmax) {
// Intersects the bottom edge.
if (x2 >= (int32)xmin) {
// Intersects the top edge too.
pt1.x = (int16)x2; pt1.y = ymax;
pt2.x = (int16)x1; pt2.y = ymin;
return 0;
}
// Intersects the left edge.
int32 y2 = (int32)pt0.y
+ (((int32)xmin - pt0.x) * (int32)kx)/(int32)ky;
pt1.x = xmin; pt1.y = (int16)y2;
pt2.x = (int16)x1; pt2.y = ymin;
return 0;
}
// Intersects the right edge. (x1 > (int32)xmax.)
int32 y1 = (int32)pt0.y
+ (((int32)xmax - pt0.x) * (int32)kx)/(int32)ky;
if (x2 >= (int32)xmin) {
// Intersects the top edge.
pt1.x = (int16)x2; pt1.y = ymax;
pt2.x = xmax; pt2.y = (int16)y1;
return 0;
}
// Intersects the left edge.
int32 y2 = (int32)pt0.y
+ (((int32)xmin - pt0.x) * (int32)kx)/(int32)ky;
pt1.x = xmin; pt1.y = (int16)y2;
pt2.x = xmax; pt2.y = (int16)y1;
return 0;
} // End of function rect_int16::cuts.