1 /*
2 * geography.c --
3 *
4 * This file provides general purpose geography functions.
5 * See the user documentation for more information.
6 *
7 * Copyright (c) 2004 Gordon D. Carrie. All rights reserved.
8 *
9 * Licensed under the Open Software License version 2.1
10 *
11 * Please address questions and feedback to user0@tkgeomap.org
12 *
13 * @(#) $Id: geography.c,v 1.9 2007/06/26 21:49:30 tkgeomap Exp $
14 *
15 ********************************************
16 *
17 */
18
19 #include <limits.h>
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <string.h>
23 #include <assert.h>
24 #include "geography.h"
25
26 /*
27 * The following constants enable conversion between floating point degree
28 * and radian measurements and the Angle data type.
29 *
30 * Angles are stored as integer microdegrees. This helps reduce the effects
31 * of round off. A microdegree is about 0.11 meters on the spherical Earth,
32 * so points closer than this are assumed to be at the same location.
33 */
34
35 #define UNITPERDEG 1.0e+6
36 #define DEGPERUNIT 1.0e-6
37 #define RADPERUNIT 0.017453292519943294892e-6
38 #define UNITPERRAD 57.29577951308232088e+6
39 #define MAXDEG (INT_MAX * DEGPERUNIT)
40 #define MINDEG (INT_MIN * DEGPERUNIT)
41 #define MAXRAD (INT_MAX * RADPERUNIT)
42 #define MINRAD (INT_MIN * RADPERUNIT)
43
44 /*
45 * The following constants give Angle representations of some frequently used
46 * angles.
47 */
48
49 #define D360 360000000
50 #define D270 270000000
51 #define D180 180000000
52 #define D90 90000000
53
54 static double rearth = 6366707.01949;
55
56 /*
57 *----------------------------------------------------------------------
58 *
59 * REarth --
60 *
61 * This procedure retrieves the radius of the Earth.
62 *
63 * Results:
64 * See the user documentation.
65 *
66 * Side effects:
67 * None.
68 *
69 *----------------------------------------------------------------------
70 */
71
72 double
REarth(void)73 REarth(void)
74 {
75 return rearth;
76 }
77
78 /*
79 *----------------------------------------------------------------------
80 *
81 * SetREarth --
82 *
83 * This procedure sets the assumed radius of the Earth.
84 *
85 * Results:
86 * None.
87 *
88 * Side effects:
89 * A static variable is set.
90 *
91 *----------------------------------------------------------------------
92 */
93
94 void
SetREarth(r)95 SetREarth(r)
96 double r; /* New Earth radius */
97 {
98 rearth = r;
99 }
100
101 /*
102 *----------------------------------------------------------------------
103 *
104 * BadAngle --
105 *
106 * This procedure returns a bogus angle. It is used to indicate
107 * error conditions.
108 *
109 * Results:
110 * See the user documentation.
111 *
112 * Side effects:
113 * None.
114 *
115 *----------------------------------------------------------------------
116 */
117
118 Angle
BadAngle(void)119 BadAngle(void)
120 {
121 return -INT_MAX;
122 }
123
124 /*
125 *----------------------------------------------------------------------
126 *
127 * AngleIsOK --
128 *
129 * This procedure determines whether an angle value is valid.
130 *
131 * Results:
132 * See the user documentation.
133 *
134 * Side effects:
135 * None.
136 *
137 *----------------------------------------------------------------------
138 */
139
140 int
AngleIsOK(a)141 AngleIsOK(a)
142 Angle a; /* Angle to evaluate */
143 {
144 return a != BadAngle();
145 }
146
147 /*
148 *----------------------------------------------------------------------
149 *
150 * AngleIsBad --
151 *
152 * This procedure determines whether an angle value is invalid.
153 *
154 * Results:
155 * See the user documentation.
156 *
157 * Side effects:
158 * None.
159 *
160 *----------------------------------------------------------------------
161 */
162
163 int
AngleIsBad(a)164 AngleIsBad(a)
165 Angle a; /* Angle to evaluate */
166 {
167 return a == BadAngle();
168 }
169
170 /*
171 *----------------------------------------------------------------------
172 *
173 * AngleFmDeg --
174 *
175 * This procedure converts a degree value to an Angle value.
176 *
177 * Results:
178 * See the user documentation.
179 *
180 * Side effects:
181 * None.
182 *
183 *----------------------------------------------------------------------
184 */
185
186 Angle
AngleFmDeg(deg)187 AngleFmDeg(deg)
188 double deg; /* Degree value to convert */
189 {
190 return (deg > MAXDEG || deg < MINDEG)
191 ? BadAngle() : UNITPERDEG * deg + (deg > 0.0 ? 0.5 : -0.5);
192 }
193
194 /*
195 *----------------------------------------------------------------------
196 *
197 * AngleToDeg --
198 *
199 * This procedure converts an Angle value to degrees.
200 *
201 * Results:
202 * See the user documentation.
203 *
204 * Side effects:
205 * None.
206 *
207 *----------------------------------------------------------------------
208 */
209
210 double
AngleToDeg(a)211 AngleToDeg(a)
212 Angle a; /* Angle value to convert */
213 {
214 return a * DEGPERUNIT;
215 }
216
217 /*
218 *----------------------------------------------------------------------
219 *
220 * AngleFmRad --
221 *
222 * This procedure converts a radian value to an Angle value.
223 *
224 * Results:
225 * See the user documentation.
226 *
227 * Side effects:
228 * None.
229 *
230 *----------------------------------------------------------------------
231 */
232
233 Angle
AngleFmRad(rad)234 AngleFmRad(rad)
235 double rad; /* Radian value to convert */
236 {
237 return (rad > MAXRAD || rad < MINRAD)
238 ? BadAngle() : UNITPERRAD * rad + (rad > 0.0 ? 0.5 : -0.5);
239 }
240
241 /*
242 *----------------------------------------------------------------------
243 *
244 * AngleToRad --
245 *
246 * This procedure converts an Angle value to radians.
247 *
248 * Results:
249 * See the user documentation.
250 *
251 * Side effects:
252 * None.
253 *
254 *----------------------------------------------------------------------
255 */
256
257 double
AngleToRad(a)258 AngleToRad(a)
259 Angle a; /* Angle value to convert */
260 {
261 return a * RADPERUNIT;
262 }
263
264 /*
265 *----------------------------------------------------------------------
266 *
267 * ISin --
268 *
269 * This procedure computes the sine of an Angle value.
270 *
271 * Results:
272 * See the user documentation.
273 *
274 * Side effects:
275 * None.
276 *
277 *----------------------------------------------------------------------
278 */
279
280 double
ISin(a)281 ISin(a)
282 Angle a; /* Angle value */
283 {
284 return sin(a * RADPERUNIT);
285 }
286
287 /*
288 *----------------------------------------------------------------------
289 *
290 * ICos --
291 *
292 * This procedure computes the cosine of an Angle value.
293 *
294 * Results:
295 * See the user documentation.
296 *
297 * Side effects:
298 * None.
299 *
300 *----------------------------------------------------------------------
301 */
302
303 double
ICos(a)304 ICos(a)
305 Angle a; /* Angle value */
306 {
307 return cos(a * RADPERUNIT);
308 }
309
310 /*
311 *----------------------------------------------------------------------
312 *
313 * GeoPtFmDeg --
314 *
315 * This procedure sets a geopoint given latitude and longitude values
316 * in degrees.
317 *
318 * Results:
319 * See the user documentation.
320 *
321 * Side effects:
322 * None.
323 *
324 *----------------------------------------------------------------------
325 */
326
327 GeoPt
GeoPtFmDeg(dLat,dLon)328 GeoPtFmDeg(dLat, dLon)
329 double dLat; /* Latitude, degrees */
330 double dLon; /* Longitude, degrees */
331 {
332 GeoPt geoPt; /* Return value */
333 geoPt.lat = AngleFmDeg(dLat);
334 geoPt.lon = AngleFmDeg(dLon);
335 return (AngleIsBad(geoPt.lat) || AngleIsBad(geoPt.lon))
336 ? GeoPtNowhere() : geoPt;
337 }
338
339 /*
340 *----------------------------------------------------------------------
341 *
342 * GeoPtFmRad --
343 *
344 * This procedure sets a geopoint given latitude and longitude values
345 * in radians.
346 *
347 * Results:
348 * See the user documentation.
349 *
350 * Side effects:
351 * None.
352 *
353 *----------------------------------------------------------------------
354 */
355
356 GeoPt
GeoPtFmRad(dLat,dLon)357 GeoPtFmRad(dLat, dLon)
358 double dLat; /* Latitude, radians */
359 double dLon; /* Longitude, radians */
360 {
361 GeoPt geoPt;
362 geoPt.lat = AngleFmRad(dLat);
363 geoPt.lon = AngleFmRad(dLon);
364 return (AngleIsBad(geoPt.lat) || AngleIsBad(geoPt.lon))
365 ? GeoPtNowhere() : geoPt;
366 }
367
368 /*
369 *----------------------------------------------------------------------
370 *
371 * GeoPtGetDeg --
372 *
373 * This procedure retrieves latitude and longitude values from a
374 * geopoint in degrees.
375 *
376 * Results:
377 * See the user documentation.
378 *
379 * Side effects:
380 * None.
381 *
382 *----------------------------------------------------------------------
383 */
384
385 void
GeoPtGetDeg(geoPt,dLatPtr,dLonPtr)386 GeoPtGetDeg(geoPt, dLatPtr, dLonPtr)
387 GeoPt geoPt; /* Geographic point */
388 double *dLatPtr; /* Recipient of latitude */
389 double *dLonPtr; /* Recipient of longitude */
390 {
391 *dLatPtr = AngleToDeg(geoPt.lat);
392 *dLonPtr = AngleToDeg(geoPt.lon);
393 }
394
395 /*
396 *----------------------------------------------------------------------
397 *
398 * GeoPtGetRad --
399 *
400 * This procedure retrieves latitude and longitude values from a
401 * geopoint in degrees.
402 *
403 * Results:
404 * See the user documentation.
405 *
406 * Side effects:
407 * None.
408 *
409 *----------------------------------------------------------------------
410 */
411
412 void
GeoPtGetRad(geoPt,dLatPtr,dLonPtr)413 GeoPtGetRad(geoPt, dLatPtr, dLonPtr)
414 GeoPt geoPt; /* Geographic point */
415 double *dLatPtr; /* Recipient of latitude */
416 double *dLonPtr; /* Recipient of longitude */
417 {
418 *dLatPtr = AngleToRad(geoPt.lat);
419 *dLonPtr = AngleToRad(geoPt.lon);
420 }
421
422 /*
423 *----------------------------------------------------------------------
424 *
425 * GeoPtIsSomewhere --
426 *
427 * This procedure determines whether a geopoint is valid.
428 *
429 * Results:
430 * See the user documentation.
431 *
432 * Side effects:
433 * None.
434 *
435 *----------------------------------------------------------------------
436 */
437
438 int
GeoPtIsSomewhere(geoPt)439 GeoPtIsSomewhere(geoPt)
440 GeoPt geoPt; /* Geographic point to evaluate */
441 {
442 return AngleIsOK(geoPt.lat) && AngleIsOK(geoPt.lon);
443 }
444 /*
445 *----------------------------------------------------------------------
446 *
447 * GeoPtIsNowhere --
448 *
449 * This procedure determines whether a geopoint is invalid.
450 *
451 * Results:
452 * See the user documentation.
453 *
454 * Side effects:
455 * None.
456 *
457 *----------------------------------------------------------------------
458 */
459
460 int
GeoPtIsNowhere(geoPt)461 GeoPtIsNowhere(geoPt)
462 GeoPt geoPt; /* Geographic point to evaluate */
463 {
464 return AngleIsBad(geoPt.lat) || AngleIsBad(geoPt.lon);
465 }
466
467 /*
468 *----------------------------------------------------------------------
469 *
470 * GeoPtNowhere --
471 *
472 * This procedure returns an invalid geopoint.
473 *
474 * Results:
475 * See the user documentation.
476 *
477 * Side effects:
478 * None.
479 *
480 *----------------------------------------------------------------------
481 */
482
483 GeoPt
GeoPtNowhere(void)484 GeoPtNowhere(void)
485 {
486 GeoPt nowhere; /* Return value */
487 nowhere.lat = nowhere.lon = BadAngle();
488 return nowhere;
489 }
490 /*
491 *----------------------------------------------------------------------
492 *
493 * MapPtIsSomewhere --
494 *
495 * This procedure determines whether a map-point is valid.
496 *
497 * Results:
498 * See the user documentation.
499 *
500 * Side effects:
501 * None.
502 *
503 *----------------------------------------------------------------------
504 */
505
506
507 int
MapPtIsSomewhere(mapPt)508 MapPtIsSomewhere(mapPt)
509 MapPt mapPt; /* Map point to evaluate */
510 {
511 return mapPt.abs != FLT_MAX && mapPt.ord != FLT_MAX;
512 }
513
514 /*
515 *----------------------------------------------------------------------
516 *
517 * MapPtIsNowhere --
518 *
519 * This procedure determines whether a map-point is invalid.
520 *
521 * Results:
522 * See the user documentation.
523 *
524 * Side effects:
525 * None.
526 *
527 *----------------------------------------------------------------------
528 */
529
530 int
MapPtIsNowhere(mapPt)531 MapPtIsNowhere(mapPt)
532 MapPt mapPt; /* Map point to evaluate */
533 {
534 return mapPt.abs == FLT_MAX || mapPt.ord == FLT_MAX;
535 }
536
537 /*
538 *----------------------------------------------------------------------
539 *
540 * MapPtNowhere --
541 *
542 * This procedure returns an invalid map-point.
543 *
544 * Results:
545 * See the user documentation.
546 *
547 * Side effects:
548 * None.
549 *
550 *----------------------------------------------------------------------
551 */
552
553 MapPt
MapPtNowhere(void)554 MapPtNowhere(void)
555 {
556 MapPt nowhere = {FLT_MAX, FLT_MAX};
557 return nowhere;
558 }
559
560 /*
561 *----------------------------------------------------------------------
562 *
563 * LatLonToCart --
564 *
565 * This procedure convert a geopoint to 3D Geocentric Cartesian
566 * coordinates where Earth has unit radius
567 *
568 * Results:
569 * See the user documentation.
570 *
571 * Side effects:
572 * None.
573 *
574 *----------------------------------------------------------------------
575 */
576
577 CartPt
LatLonToCart(geoPt)578 LatLonToCart(geoPt)
579 GeoPt geoPt; /* Geographic point value to convert */
580 {
581 double lat, lon; /* Latitude, longitude in radians */
582 double coslat;
583 CartPt cpt; /* Return value */
584
585 GeoPtGetRad(geoPt, &lat, &lon);
586 coslat = cos(lat);
587 cpt.x = coslat * cos(lon);
588 cpt.y = coslat * sin(lon);
589 cpt.z = sin(lat);
590 return cpt;
591 }
592
593 /*
594 *----------------------------------------------------------------------
595 *
596 * ScaleMapPt --
597 *
598 * This procedure scales a map-point.
599 *
600 * Results:
601 * See the user documentation.
602 *
603 * Side effects:
604 * None.
605 *
606 *----------------------------------------------------------------------
607 */
608
609 MapPt
ScaleMapPt(mapPt,scale)610 ScaleMapPt(mapPt, scale)
611 MapPt mapPt; /* Map point to scale */
612 double scale; /* Amount to scale it by */
613 {
614 if (MapPtIsNowhere(mapPt) || scale <= 0.0) {
615 return MapPtNowhere();
616 }
617 mapPt.abs *= scale;
618 mapPt.ord *= scale;
619 return mapPt;
620 }
621
622 /*
623 *----------------------------------------------------------------------
624 *
625 * GeoDistance --
626 *
627 * This procedure computes the distance between two geopoints.
628 *
629 * Results:
630 * See the user documentation.
631 *
632 * Side effects:
633 * None.
634 *
635 *----------------------------------------------------------------------
636 */
637
638 Angle
GeoDistance(p1,p2)639 GeoDistance(p1, p2)
640 GeoPt p1; /* First point */
641 GeoPt p2; /* Second point */
642 {
643 double lat1, lon1, lat2, lon2;
644 double sin_dlon_2, sin_dlat_2;
645 double a;
646
647 /*
648 * Reference:
649 *
650 * R.W. Sinnott, "Virtues of the Haversine", Sky and Telescope,
651 * vol. 68, no. 2, 1984, p. 159
652 *
653 * cited in:
654 * http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
655 */
656
657 GeoPtGetRad(p1, &lat1, &lon1);
658 GeoPtGetRad(p2, &lat2, &lon2);
659 sin_dlon_2 = sin(0.5 * (lon2 - lon1));
660 sin_dlat_2 = sin(0.5 * (lat2 - lat1));
661 a = sqrt(sin_dlat_2 * sin_dlat_2
662 + cos(lat1) * cos(lat2) * sin_dlon_2 * sin_dlon_2);
663 return AngleFmRad(a > 1.0 ? M_PI : 2.0 * asin(a));
664 }
665 /*
666 *----------------------------------------------------------------------
667 *
668 * GeoQuickDistance --
669 *
670 * This procedure computes a measure of the distance between two geopoints.
671 *
672 * Results:
673 * See the user documentation.
674 *
675 * Side effects:
676 * None.
677 *
678 *----------------------------------------------------------------------
679 */
680
681 double
GeoQuickDistance(p1,p2)682 GeoQuickDistance(p1, p2)
683 GeoPt p1; /* First point */
684 GeoPt p2; /* Second point */
685 {
686 double lat1, lon1, lat2, lon2;
687 double sin_dlon_2, sin_dlat_2;
688
689 /*
690 * Reference:
691 *
692 * R.W. Sinnott, "Virtues of the Haversine", Sky and Telescope,
693 * vol. 68, no. 2, 1984, p. 159
694 *
695 * cited in:
696 * http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
697 */
698
699 GeoPtGetRad(p1, &lat1, &lon1);
700 GeoPtGetRad(p2, &lat2, &lon2);
701 sin_dlon_2 = sin(0.5 * (lon2 - lon1));
702 sin_dlat_2 = sin(0.5 * (lat2 - lat1));
703 return sin_dlat_2 * sin_dlat_2
704 + cos(lat1) * cos(lat2) * sin_dlon_2 * sin_dlon_2;
705 }
706
707 /*
708 *----------------------------------------------------------------------
709 *
710 * Azimuth --
711 *
712 * This procedure computes the azimuth (bearing) from one point to another.
713 *
714 * Results:
715 * See the user documentation.
716 *
717 * Side effects:
718 * None.
719 *
720 *----------------------------------------------------------------------
721 */
722
723 Angle
Azimuth(p1,p2)724 Azimuth(p1, p2)
725 GeoPt p1; /* First point */
726 GeoPt p2; /* Second point */
727 {
728 double lat1, lon1, lat2, lon2;
729 double cosDLon, sinDLon, sinDLat, sinMLat;
730
731 GeoPtGetRad(p1, &lat1, &lon1);
732 GeoPtGetRad(p2, &lat2, &lon2);
733 cosDLon = cos(lon2 - lon1);
734 sinDLon = sin(lon2 - lon1);
735 sinDLat = sin(lat1 - lat2);
736 sinMLat = sin(lat2 + lat1);
737 return AngleFmRad(atan2(cos(lat2) * sinDLon,
738 0.5 * (sinMLat - sinDLat - (sinMLat + sinDLat) * cosDLon)));
739 }
740
741 /*
742 *----------------------------------------------------------------------
743 *
744 * GeoStep --
745 *
746 * This procedure finds the point a given distance and direction from
747 * a another point.
748 *
749 * Results:
750 * See the user documentation.
751 *
752 * Side effects:
753 * None.
754 *
755 *----------------------------------------------------------------------
756 */
757
758 GeoPt
GeoStep(geoPt,iDir,iDist)759 GeoStep(geoPt, iDir, iDist)
760 GeoPt geoPt; /* Starting point */
761 Angle iDir; /* Which way to step */
762 Angle iDist; /* How far to step */
763 {
764 double dist, dir;
765 double lat, lon;
766 double cos_dist, sin_dist, cos_dir, sin_dir;
767 double cos_lat, cos_lon, sin_lon, sin_lat;
768 double x, y, z, h_1, h_2, dh;
769
770 dist = AngleToRad(iDist);
771 cos_dist = cos(dist);
772 sin_dist = sin(dist);
773
774 dir = AngleToRad(iDir);
775 cos_dir = cos(dir);
776 sin_dir = sin(dir);
777
778 GeoPtGetRad(geoPt, &lat, &lon);
779 cos_lat = cos(lat);
780 cos_lon = cos(lon);
781 sin_lon = sin(lon);
782 sin_lat = sin(lat);
783 x = cos_dist * cos_lon * cos_lat - sin_dir * sin_dist * sin_lon
784 - cos_lon * cos_dir * sin_dist * sin_lat;
785 y = sin_dir * cos_lon * sin_dist + cos_dist * cos_lat * sin_lon
786 - cos_dir * sin_dist * sin_lon * sin_lat;
787 z = cos_lat * cos_dir * sin_dist + cos_dist * sin_lat;
788
789 /*
790 * Compute x^2 + y^2
791 */
792
793 h_1 = cos_dist * cos_lat - cos_dir * sin_dist * sin_lat;
794 h_2 = (sin_dir) * (sin_dist);
795 dh = h_1 * h_1 + h_2 * h_2;
796
797 lat = (dh == 0.0) ? (z > 0 ? M_PI_2 : -M_PI_2) : atan(z / sqrt(dh));
798 lon = atan2(y, x);
799
800 return GeoPtFmRad(lat, lon);
801 }
802
803 /*
804 *----------------------------------------------------------------------
805 *
806 * GCircleX --
807 *
808 * This procedure computes the intersection of two great circles.
809 *
810 * Results:
811 * See the user documentation.
812 *
813 * Side effects:
814 * None.
815 *
816 *----------------------------------------------------------------------
817 */
818
819 GeoPt
GCircleX(ln1pt1,ln1pt2,ln2pt1,ln2pt2)820 GCircleX(ln1pt1, ln1pt2, ln2pt1, ln2pt2)
821 GeoPt ln1pt1; /* First point on first great circle */
822 GeoPt ln1pt2; /* Second point on first great circle */
823 GeoPt ln2pt1; /* First point on second great circle */
824 GeoPt ln2pt2; /* Second point on second great circle */
825 {
826 CartPt ln1pt1c, ln1pt2c, ln2pt1c, ln2pt2c;
827 /* Points on the great circles in
828 * Cartesian Geocentric coordinates */
829 CartPt p1; /* Normal to plane containing ln1pt1 ln1pt2 */
830 CartPt p2; /* Normal to plane containing ln2pt2 ln2pt2 */
831 CartPt c; /* Vector perpendicular to p1 and p2 */
832 double cs; /* Normalizing factor */
833 CartPt m; /* Mean of the four input points */
834 double dp, dm; /* Distance from mean to +-(cx, cy, cz) */
835 double lat, lon; /* Lat and lon of return value */
836
837 /*
838 * Points in 3D Geocentric Cartesian Coordinates
839 */
840
841 ln1pt1c = LatLonToCart(ln1pt1);
842 ln1pt2c = LatLonToCart(ln1pt2);
843 ln2pt1c = LatLonToCart(ln2pt1);
844 ln2pt2c = LatLonToCart(ln2pt2);
845
846 /*
847 * Normal to plane containing ln1pt1 and ln1pt2 (given by cross product)
848 */
849
850 p1.x = ln1pt1c.y * ln1pt2c.z - ln1pt1c.z * ln1pt2c.y;
851 p1.y = ln1pt1c.z * ln1pt2c.x - ln1pt1c.x * ln1pt2c.z;
852 p1.z = ln1pt1c.x * ln1pt2c.y - ln1pt1c.y * ln1pt2c.x;
853
854 /*
855 * Normal to plane containing ln2pt1 and ln2pt2 (given by cross product)
856 */
857
858 p2.x = ln2pt1c.y * ln2pt2c.z - ln2pt1c.z * ln2pt2c.y;
859 p2.y = ln2pt1c.z * ln2pt2c.x - ln2pt1c.x * ln2pt2c.z;
860 p2.z = ln2pt1c.x * ln2pt2c.y - ln2pt1c.y * ln2pt2c.x;
861
862 /*
863 * Perpendicular to p1 and p2 (another cross product)
864 */
865
866 c.x = p1.y * p2.z - p1.z * p2.y;
867 c.y = p1.z * p2.x - p1.x * p2.z;
868 c.z = p1.x * p2.y - p1.y * p2.x;
869
870 /*
871 * If all pts on one line, undefined intersection
872 */
873
874 if (c.x == 0.0 && c.y == 0.0 && c.z == 0.0) {
875 return GeoPtNowhere();
876 }
877
878 /*
879 * Normalize
880 */
881
882 cs = 1.0 / sqrt(c.x * c.x + c.y * c.y + c.z * c.z);
883 c.x *= cs;
884 c.y *= cs;
885 c.z *= cs;
886
887 /*
888 * Intersection is at a pole
889 */
890
891 if (c.x == 0.0 && c.y == 0.0) {
892 return GeoPtFmDeg((c.z > 0) ? 90.0 : -90.0, 0.0);
893 }
894
895 /*
896 * Find which of the two intersections of the two great circles is
897 * closer to the mean of the four input points
898 */
899
900 m.x = 0.25 * (ln1pt1c.x + ln1pt2c.x + ln2pt2c.x + ln2pt2c.x);
901 m.y = 0.25 * (ln1pt1c.y + ln1pt2c.y + ln2pt2c.y + ln2pt2c.y);
902 m.z = 0.25 * (ln1pt1c.z + ln1pt2c.z + ln2pt2c.z + ln2pt2c.z);
903 dp = (m.x-c.x) * (m.x-c.x) + (m.y-c.y) * (m.y-c.y) + (m.z-c.z) * (m.z-c.z);
904 dm = (m.x+c.x) * (m.x+c.x) + (m.y+c.y) * (m.y+c.y) + (m.z+c.z) * (m.z+c.z);
905 if (dm < dp) {
906 /*
907 * Mean point is closer to (-c.x, -c.y, -c.z)
908 */
909
910 c.x = -c.x;
911 c.y = -c.y;
912 c.z = -c.z;
913 }
914
915 lat = atan(c.z / sqrt(c.x * c.x + c.y * c.y));
916 lon = atan2(c.y, c.x);
917 return GeoPtFmRad(lat, lon);
918 }
919
920 /*
921 *----------------------------------------------------------------------
922 *
923 * DomainLat --
924 *
925 * This procedure puts an angle into the range -90.0 <= angle <= 90.0
926 * without changing its sine.
927 *
928 * Results:
929 * See the user documentation.
930 *
931 * Side effects:
932 * None.
933 *
934 *----------------------------------------------------------------------
935 */
936
937 Angle
DomainLat(lat)938 DomainLat(lat)
939 Angle lat; /* Latitude value */
940 {
941 if (lat > D360) {
942 lat = lat - D360 * (lat / D360);
943 } else if (lat < 0) {
944 lat = lat + D360 * (-lat / D360 + 1);
945 }
946 if (lat > D90 && lat < D270) {
947 lat = D180 - lat;
948 } else if (lat >= D270) {
949 lat = lat - D360;
950 }
951 return lat;
952 }
953
954 /*
955 *----------------------------------------------------------------------
956 *
957 * DomainLon --
958 *
959 * This procedure puts an angle within 180 degrees of a reference angle.
960 *
961 * Results:
962 * See the user documentation.
963 *
964 * Side effects:
965 * None.
966 *
967 *----------------------------------------------------------------------
968 */
969
970 Angle
DomainLon(lon,refLon)971 DomainLon(lon, refLon)
972 Angle lon; /* Longitude value */
973 Angle refLon; /* Reference longitude */
974 {
975 if (lon == refLon) {
976 return lon;
977 }
978 if (lon > refLon + D360) {
979 lon -= D360 * ((lon - refLon) / D360);
980 } else if (lon < refLon - D360) {
981 lon += D360 * ((refLon - lon) / D360);
982 }
983 if (lon > refLon + D180) {
984 return lon - D360;
985 } else if (lon < refLon - D180) {
986 return lon + D360;
987 } else {
988 return lon;
989 }
990 }
991
992 /*
993 *----------------------------------------------------------------------
994 *
995 * GwchLon --
996 *
997 * This procedure is equivalent to DomainLon(lon, 0.0)
998 *
999 * Results:
1000 * See the user documentation.
1001 *
1002 * Side effects:
1003 * None.
1004 *
1005 *----------------------------------------------------------------------
1006 */
1007
1008 Angle
GwchLon(lon)1009 GwchLon(lon)
1010 Angle lon; /* Longitude value */
1011 {
1012 if (lon == 0) {
1013 return lon;
1014 }
1015 if (lon > D360) {
1016 lon -= D360 * ((lon ) / D360);
1017 } else if (lon < -D360) {
1018 lon += D360 * ((-lon) / D360);
1019 }
1020 if (lon > D180) {
1021 return lon - D360;
1022 } else if (lon < -D180) {
1023 return lon + D360;
1024 } else {
1025 return lon;
1026 }
1027 }
1028
1029 /*
1030 *----------------------------------------------------------------------
1031 *
1032 * DomainLonPt --
1033 *
1034 * This procedure returns a geographic point whose longitude is within
1035 * 180 degrees of a reference longitude.
1036 *
1037 * Results:
1038 * See the user documentation.
1039 *
1040 * Side effects:
1041 * None.
1042 *
1043 *----------------------------------------------------------------------
1044 */
1045
1046 GeoPt
DomainLonPt(geoPt,refLon)1047 DomainLonPt(geoPt, refLon)
1048 GeoPt geoPt; /* Geographic point */
1049 Angle refLon; /* Reference longitude */
1050 {
1051 geoPt.lat = DomainLat(geoPt.lat);
1052 geoPt.lon = DomainLon(geoPt.lon, refLon);
1053 return geoPt;
1054 }
1055
1056 /*
1057 *----------------------------------------------------------------------
1058 *
1059 * GwchLonPt --
1060 *
1061 * This procedure is equivalent to DomainLonPt(geoPt, 0.0).
1062 *
1063 * Results:
1064 * See the user documentation.
1065 *
1066 * Side effects:
1067 * None.
1068 *
1069 *----------------------------------------------------------------------
1070 */
1071
1072 GeoPt
GwchLonPt(geoPt)1073 GwchLonPt(geoPt)
1074 GeoPt geoPt; /* Geographic point */
1075 {
1076 return DomainLonPt(geoPt, 0);
1077 }
1078
1079 /*
1080 *----------------------------------------------------------------------
1081 *
1082 * LonCmp --
1083 *
1084 * This procedure compares two longitudes.
1085 *
1086 * Results:
1087 * See the user documentation.
1088 *
1089 * Side effects:
1090 * None.
1091 *
1092 *----------------------------------------------------------------------
1093 */
1094
1095 enum LonSgn
LonCmp(lon0,lon1)1096 LonCmp(lon0, lon1)
1097 Angle lon0; /* First longitude */
1098 Angle lon1; /* Second longitude */
1099 {
1100 lon0 = DomainLon(lon0, lon1);
1101 return (lon0 < lon1) ? West : (lon0 > lon1) ? East : PrMd;
1102 }
1103
1104 /*
1105 *----------------------------------------------------------------------
1106 *
1107 * LonBtwn --
1108 *
1109 * This procedure determines whether a meridian lies between two others.
1110 *
1111 * Results:
1112 * See the user documentation.
1113 *
1114 * Side effects:
1115 * None.
1116 *
1117 *----------------------------------------------------------------------
1118 */
1119
1120 int
LonBtwn(lon,lon0,lon1)1121 LonBtwn(lon, lon0, lon1)
1122 Angle lon; /* Longitude to evaluate */
1123 Angle lon0; /* Longitude at one end of interval */
1124 Angle lon1; /* Longitude at other end of interval */
1125 {
1126 lon0 = DomainLon(lon0, lon);
1127 lon1 = DomainLon(lon1, lon);
1128 return ((lon0 > lon1) ? lon0 - lon1 : lon1 - lon0) < D180
1129 && ((lon0 < lon && lon < lon1) || (lon1 < lon && lon < lon0));
1130 }
1131
1132 /*
1133 *----------------------------------------------------------------------
1134 *
1135 * LonBtwn1 --
1136 *
1137 * This procedure determines whether a meridian lies between or coincides
1138 * with two others,
1139 *
1140 * Results:
1141 * See the user documentation.
1142 *
1143 * Side effects:
1144 * None.
1145 *
1146 *----------------------------------------------------------------------
1147 */
1148
1149 int
LonBtwn1(lon,lon0,lon1)1150 LonBtwn1(lon, lon0, lon1)
1151 Angle lon; /* Longitude to evaluate */
1152 Angle lon0; /* Longitude at one end of interval */
1153 Angle lon1; /* Longitude at other end of interval */
1154 {
1155 lon0 = DomainLon(lon0, lon);
1156 lon1 = DomainLon(lon1, lon);
1157 return ((lon0 > lon1) ? lon0 - lon1 : lon1 - lon0) < D180
1158 && ((lon0 < lon && lon <= lon1) || (lon1 < lon && lon <= lon0));
1159 }
1160
1161 /*
1162 *----------------------------------------------------------------------
1163 *
1164 * LatCmp --
1165 *
1166 * This procedure compares two latitudes.
1167 *
1168 * Results:
1169 * See the user documentation.
1170 *
1171 * Side effects:
1172 * None.
1173 *
1174 *----------------------------------------------------------------------
1175 */
1176
1177 enum LatSgn
LatCmp(lat0,lat1)1178 LatCmp(lat0, lat1)
1179 Angle lat0; /* First latitude */
1180 Angle lat1; /* Second latitude */
1181 {
1182 return (lat0 < lat1) ? South : (lat0 > lat1) ? North : Eq;
1183 }
1184
1185 /*
1186 *----------------------------------------------------------------------
1187 *
1188 * AngleCmp --
1189 *
1190 * This procedure compares to angles.
1191 *
1192 * Results:
1193 * See the user documentation.
1194 *
1195 * Side effects:
1196 * None.
1197 *
1198 *----------------------------------------------------------------------
1199 */
1200
1201 int
AngleCmp(a0,a1)1202 AngleCmp(a0, a1)
1203 Angle a0; /* First angle to compare */
1204 Angle a1; /* Second angle to compare */
1205 {
1206 return (a0 < a1) ? -1 : (a0 > a1) ? 1 : 0;
1207 }
1208