1 /**
2 * \addtogroup utilities
3 * @{
4 */
5
6 /**
7 * \file src/locator.cxx
8 * \brief locator and bearing conversion interface
9 * \author Stephane Fillod and the Hamlib Group
10 * \date 2000-2006
11 *
12 * Hamlib Interface - locator, bearing, and conversion calls
13 */
14
15 /*
16 * Hamlib Interface - locator and bearing conversion calls
17 * Copyright (c) 2001-2006 by Stephane Fillod
18 * Copyright (c) 2003 by Nate Bargmann
19 * Copyright (c) 2003 by Dave Hines
20 *
21 * $Id$
22 *
23 * Code to determine bearing and range was taken from the Great Circle,
24 * by S. R. Sampson, N5OWK.
25 * Ref: "Air Navigation", Air Force Manual 51-40, 1 February 1987
26 * Ref: "ARRL Satellite Experimenters Handbook", August 1990
27 *
28 * Code to calculate distance and azimuth between two Maidenhead locators,
29 * taken from wwl, by IK0ZSN Mirko Caserta.
30 *
31 * New bearing code added by N0NB was found at:
32 * http://williams.best.vwh.net/avform.htm#Crs
33 *
34 *
35 *
36 * This library is free software; you can redistribute it and/or modify
37 * it under the terms of the GNU Library General Public License as
38 * published by the Free Software Foundation; either version 2 of
39 * the License, or (at your option) any later version.
40 *
41 * This program is distributed in the hope that it will be useful,
42 * but WITHOUT ANY WARRANTY; without even the implied warranty of
43 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
44 * GNU Library General Public License for more details.
45 *
46 * You should have received a copy of the GNU Library General Public
47 * License along with this library; if not, write to the
48 *
49 * Free Software Foundation, Inc.
50 * 51 Franklin Street, Fifth Floor
51 * Boston, MA 02110-1301 USA.
52 *
53 */
54
55 /*! \page hamlib Hamlib general purpose API
56 *
57 * Here are grouped some often used functions, like locator conversion
58 * routines.
59 */
60
61 #ifdef HAVE_CONFIG_H
62 #include "config.h"
63 #endif
64
65 #include <stdio.h>
66 #include <stdlib.h>
67 #include <string.h>
68 #include <ctype.h>
69 #include <math.h>
70
71 #include "configuration.h"
72 #include "locator.h"
73
74 namespace QRB {
75
76 #define RADIAN (180.0 / M_PI)
77
78 /* arc length for 1 degree, 60 Nautical Miles
79 * 109.728 Kilometers
80 * 68.182 statute miles
81 * arc length for 1 radian, 6286.951 Km
82 * 3437.746 Nm
83 * 3906.41 Sm
84 */
85 #define ARC_IN_KM 6372.5639
86 #define ARC_IN_NM 3484.5603
87 #define ARC_IN_SM 3959.7276
88
89 /* The following is contributed by Dave Hines M1CXW
90 *
91 * begin dph
92 */
93 /*
94 * These are the constants used when converting between Maidenhead grid
95 * locators and longitude/latitude values. MAX_LOCATOR_PAIRS is the maximum
96 * number of locator character pairs to convert. This number MUST NOT exceed
97 * the number of pairs of values in loc_char_range[].
98 * Setting MAX_LOCATOR_PAIRS to 3 will convert the currently defined 6
99 * character locators. A value of 4 will convert the extended 8 character
100 * locators described in section 3L of "The IARU region 1 VHF managers
101 * handbook". Values of 5 and 6 will extent the format even more, to the
102 * longest definition I have seen for locators, see
103 * http://www.btinternet.com/~g8yoa/geog/non-ra.html
104 * Beware that there seems to be no universally accepted standard for 10 & 12
105 * character locators.
106 *
107 * The ranges of characters which will be accepted by locator2longlat, and
108 * generated by longlat2locator, are specified by the loc_char_range[] array.
109 * This array may be changed without requiring any other code changes.
110 *
111 * For the fifth pair to range from aa to xx use:
112 * const static int loc_char_range[] = { 18, 10, 24, 10, 24, 10 };
113 *
114 * For the fifth pair to range from aa to yy use:
115 * const static int loc_char_range[] = { 18, 10, 24, 10, 25, 10 };
116 *
117 * MAX_LOCATOR_PAIRS now sets the limit locator2longlat() will convert and
118 * sets the maximum length longlat2locator() will generate. Each function
119 * properly handles any value from 1 to 6 so MAX_LOCATOR_PAIRS should be
120 * left at 6. MIN_LOCATOR_PAIRS sets a floor on the shortest locator that
121 * should be handled. -N0NB
122 */
123 const static int loc_char_range[] = { 18, 10, 24, 10, 24, 10 };
124
125 #define MAX_LOCATOR_PAIRS 6
126 #define MIN_LOCATOR_PAIRS 1
127
128 /**
129 * \brief Convert DMS to decimal degrees
130 * \param degrees Degrees, whole degrees
131 * \param minutes Minutes, whole minutes
132 * \param seconds Seconds, decimal seconds
133 * \param sw South or West
134 *
135 * Convert degree/minute/second angle to decimal degrees angle.
136 * \a degrees >360, \a minutes > 60, and \a seconds > 60.0 are allowed,
137 * but resulting angle won't be normalized.
138 *
139 * When the variable sw is passed a value of 1, the returned decimal
140 * degrees value will be negative (south or west). When passed a
141 * value of 0 the returned decimal degrees value will be positive
142 * (north or east).
143 *
144 * \return The angle in decimal degrees.
145 *
146 * \sa dec2dms()
147 */
148
dms2dec(int degrees,int minutes,double seconds,int sw)149 double dms2dec(int degrees, int minutes, double seconds, int sw) {
150 double st;
151
152 if (degrees < 0)
153 degrees = abs(degrees);
154 if (minutes < 0)
155 minutes = abs(minutes);
156 if (seconds < 0)
157 seconds = fabs(seconds);
158
159 st = (double)degrees + (double)minutes / 60. + seconds / 3600.;
160
161 if (sw == 1)
162 return -st;
163 else
164 return st;
165 }
166
167 /**
168 * \brief Convert D M.MMM notation to decimal degrees
169 * \param degrees Degrees, whole degrees
170 * \param minutes Minutes, decimal minutes
171 * \param sw South or West
172 *
173 * Convert a degrees, decimal minutes notation common on
174 * many GPS units to its decimal degrees value.
175 *
176 * \a degrees > 360, \a minutes > 60.0 are allowed, but
177 * resulting angle won't be normalized.
178 *
179 * When the variable sw is passed a value of 1, the returned decimal
180 * degrees value will be negative (south or west). When passed a
181 * value of 0 the returned decimal degrees value will be positive
182 * (north or east).
183 *
184 * \return The angle in decimal degrees.
185 *
186 * \sa dec2dmmm()
187 */
188
dmmm2dec(int degrees,double minutes,int sw)189 double dmmm2dec(int degrees, double minutes, int sw) {
190 double st;
191
192 if (degrees < 0)
193 degrees = abs(degrees);
194 if (minutes < 0)
195 minutes = fabs(minutes);
196
197 st = (double)degrees + minutes / 60.;
198
199 if (sw == 1)
200 return -st;
201 else
202 return st;
203 }
204
205 /**
206 * \brief Convert decimal degrees angle into DMS notation
207 * \param dec Decimal degrees
208 * \param degrees Pointer for the calculated whole Degrees
209 * \param minutes Pointer for the calculated whole Minutes
210 * \param seconds Pointer for the calculated decimal Seconds
211 * \param sw Pointer for the calculated SW flag
212 *
213 * Convert decimal degrees angle into its degree/minute/second
214 * notation.
215 *
216 * When \a dec < -180 or \a dec > 180, the angle will be normalized
217 * within these limits and the sign set appropriately.
218 *
219 * Upon return dec2dms guarantees 0 >= \a degrees <= 180,
220 * 0 >= \a minutes < 60, and 0.0 >= \a seconds < 60.0.
221 *
222 * When \a dec is < 0.0 \a sw will be set to 1. When \a dec is
223 * >= 0.0 \a sw will be set to 0. This flag allows the application
224 * to determine whether the DMS angle should be treated as negative
225 * (south or west).
226 *
227 * \retval -QRB_EINVAL if any of the pointers are NULL.
228 * \retval QRB_OK if conversion went OK.
229 *
230 * \sa dms2dec()
231 */
232
dec2dms(double dec,int * degrees,int * minutes,double * seconds,int * sw)233 int dec2dms(double dec, int *degrees, int *minutes, double *seconds, int *sw) {
234 int deg, min;
235 double st;
236
237 /* bail if NULL pointers passed */
238 if (!degrees || !minutes || !seconds || !sw)
239 return -QRB_EINVAL;
240
241 /* reverse the sign if dec has a magnitude greater
242 * than 180 and factor out multiples of 360.
243 * e.g. when passed 270 st will be set to -90
244 * and when passed -270 st will be set to 90. If
245 * passed 361 st will be set to 1, etc. If passed
246 * a value > -180 || < 180, value will be unchanged.
247 */
248 if (dec >= 0.0)
249 st = fmod(dec + 180, 360) - 180;
250 else
251 st = fmod(dec - 180, 360) + 180;
252
253 /* if after all of that st is negative, we want deg
254 * to be negative as well except for 180 which we want
255 * to be positive.
256 */
257 if (st < 0.0 && st != -180)
258 *sw = 1;
259 else
260 *sw = 0;
261
262 /* work on st as a positive value to remove a
263 * bug introduced by the effect of floor() when
264 * passed a negative value. e.g. when passed
265 * -96.8333 floor() returns -95! Also avoids
266 * a rounding error introduced on negative values.
267 */
268 st = fabs(st);
269
270 deg = (int)floor(st);
271 st = 60. * (st - (double)deg);
272 min = (int)floor(st);
273 st = 60. * (st - (double)min);
274
275 *degrees = deg;
276 *minutes = min;
277 *seconds = st;
278
279 return QRB_OK;
280 }
281
282 /**
283 * \brief Convert a decimal angle into D M.MMM notation
284 * \param dec Decimal degrees
285 * \param degrees Pointer for the calculated whole Degrees
286 * \param minutes Pointer for the calculated decimal Minutes
287 * \param sw Pointer for the calculated SW flag
288 *
289 * Convert a decimal angle into its degree, decimal minute
290 * notation common on many GPS units.
291 *
292 * When passed a value < -180 or > 180, the value will be normalized
293 * within these limits and the sign set apropriately.
294 *
295 * Upon return dec2dmmm guarantees 0 >= \a degrees <= 180,
296 * 0.0 >= \a minutes < 60.0.
297 *
298 * When \a dec is < 0.0 \a sw will be set to 1. When \a dec is
299 * >= 0.0 \a sw will be set to 0. This flag allows the application
300 * to determine whether the D M.MMM angle should be treated as negative
301 * (south or west).
302 *
303 * \retval -QRB_EINVAL if any of the pointers are NULL.
304 * \retval QRB_OK if conversion went OK.
305 *
306 * \sa dmmm2dec()
307 */
308
dec2dmmm(double dec,int * degrees,double * minutes,int * sw)309 int dec2dmmm(double dec, int *degrees, double *minutes, int *sw) {
310 int r, min;
311 double sec;
312
313 /* bail if NULL pointers passed */
314 if (!degrees || !minutes || !sw)
315 return -QRB_EINVAL;
316
317 r = dec2dms(dec, degrees, &min, &sec, sw);
318 if (r != QRB_OK)
319 return r;
320
321 *minutes = (double)min + sec / 60;
322
323 return QRB_OK;
324 }
325
326 /**
327 * \brief Convert Maidenhead grid locator to Longitude/Latitude
328 * \param longitude Pointer for the calculated Longitude
329 * \param latitude Pointer for the calculated Latitude
330 * \param locator The Maidenhead grid locator--2 through 12 char + nul string
331 *
332 * Convert Maidenhead grid locator to Longitude/Latitude (decimal degrees).
333 * The locator should be in 2 through 12 chars long format.
334 * \a locator2longlat is case insensitive, however it checks for
335 * locator validity.
336 *
337 * Decimal long/lat is computed to center of grid square, i.e. given
338 * EM19 will return coordinates equivalent to the southwest corner
339 * of EM19mm.
340 *
341 * \retval -QRB_EINVAL if locator exceeds RR99xx99xx99 or exceeds length
342 * limit--currently 1 to 6 lon/lat pairs.
343 * \retval QRB_OK if conversion went OK.
344 *
345 * \bug The fifth pair ranges from aa to xx, there is another convention
346 * that ranges from aa to yy. At some point both conventions should be
347 * supported.
348 *
349 * \sa longlat2locator()
350 */
351
352 /* begin dph */
353
locator2longlat(double * longitude,double * latitude,const char * locator)354 int locator2longlat(double *longitude, double *latitude, const char *locator) {
355 int x_or_y, paircount;
356 int locvalue, pair;
357 int divisions;
358 double xy[2], ordinate;
359
360 /* bail if NULL pointers passed */
361 if (!longitude || !latitude)
362 return -QRB_EINVAL;
363
364 paircount = strlen(locator) / 2;
365
366 /* verify paircount is within limits */
367 if (paircount > MAX_LOCATOR_PAIRS)
368 paircount = MAX_LOCATOR_PAIRS;
369 else if (paircount < MIN_LOCATOR_PAIRS)
370 return -QRB_EINVAL;
371
372 /* For x(=longitude) and y(=latitude) */
373 for (x_or_y = 0; x_or_y < 2; ++x_or_y) {
374 ordinate = -90.0;
375 divisions = 1;
376
377 for (pair = 0; pair < paircount; ++pair) {
378 locvalue = locator[pair*2 + x_or_y];
379
380 /* Value of digit or letter */
381 locvalue -= (loc_char_range[pair] == 10) ? '0' :
382 (isupper(locvalue)) ? 'A' : 'a';
383
384 /* Check range for non-letter/digit or out of range */
385 if ((locvalue < 0) || (locvalue >= loc_char_range[pair]))
386 return -QRB_EINVAL;
387
388 divisions *= loc_char_range[pair];
389 ordinate += locvalue * 180.0 / divisions;
390 }
391 /* Center ordinate in the Maidenhead "square" or "subsquare" */
392 ordinate += 90.0 / divisions;
393
394 xy[x_or_y] = ordinate;
395 }
396
397 *longitude = xy[0] * 2.0;
398 *latitude = xy[1];
399
400 return QRB_OK;
401 }
402 /* end dph */
403
404 /**
405 * \brief Convert longitude/latitude to Maidenhead grid locator
406 * \param longitude Longitude, decimal degrees
407 * \param latitude Latitude, decimal degrees
408 * \param locator Pointer for the Maidenhead Locator
409 * \param pair_count Precision expressed as lon/lat pairs in the locator
410 *
411 * Convert longitude/latitude (decimal degrees) to Maidenhead grid locator.
412 * \a locator must point to an array at least \a pair_count * 2 char + '\\0'.
413 *
414 * \retval -QRB_EINVAL if \a locator is NULL or \a pair_count exceeds
415 * length limit. Currently 1 to 6 lon/lat pairs.
416 * \retval QRB_OK if conversion went OK.
417 *
418 * \bug \a locator is not tested for overflow.
419 * \bug The fifth pair ranges from aa to yy, there is another convention
420 * that ranges from aa to xx. At some point both conventions should be
421 * supported.
422 *
423 * \sa locator2longlat()
424 */
425
426 /* begin dph */
427
longlat2locator(double longitude,double latitude,char * locator,int pair_count)428 int longlat2locator(double longitude, double latitude, char *locator, int pair_count) {
429 int x_or_y, pair, locvalue, divisions;
430 double square_size, ordinate;
431
432 if (!locator)
433 return -QRB_EINVAL;
434
435 if (pair_count < MIN_LOCATOR_PAIRS || pair_count > MAX_LOCATOR_PAIRS)
436 return -QRB_EINVAL;
437
438 for (x_or_y = 0; x_or_y < 2; ++x_or_y) {
439 ordinate = (x_or_y == 0) ? longitude / 2.0 : latitude;
440 divisions = 1;
441
442 /* The 1e-6 here guards against floating point rounding errors */
443 ordinate = fmod(ordinate + 270.000001, 180.0);
444 for (pair = 0; pair < pair_count; ++pair) {
445 divisions *= loc_char_range[pair];
446 square_size = 180.0 / divisions;
447
448 locvalue = (int) (ordinate / square_size);
449 ordinate -= square_size * locvalue;
450 locvalue += (loc_char_range[pair] == 10) ? '0':'A';
451 locator[pair * 2 + x_or_y] = locvalue;
452 }
453 }
454 locator[pair_count * 2] = '\0';
455
456 return QRB_OK;
457 }
458
459 /* end dph */
460
461 /**
462 * \brief Calculate the distance and bearing between two points.
463 * \param lon1 The local Longitude, decimal degrees
464 * \param lat1 The local Latitude, decimal degrees
465 * \param lon2 The remote Longitude, decimal degrees
466 * \param lat2 The remote Latitude, decimal degrees
467 * \param distance Pointer for the distance, km
468 * \param azimuth Pointer for the bearing, decimal degrees
469 *
470 * Calculate the QRB between \a lon1, \a lat1 and \a lon2, \a lat2.
471 *
472 * This version will calculate the QRB to a precision sufficient
473 * for 12 character locators. Antipodal points, which are easily
474 * calculated, are considered equidistant and the bearing is
475 * simply resolved to be true north (0.0�).
476 *
477 * \retval -QRB_EINVAL if NULL pointer passed or lat and lon values
478 * exceed -90 to 90 or -180 to 180.
479 * \retval QRB_OK if calculations are successful.
480 *
481 * \return The distance in kilometers and azimuth in decimal degrees
482 * for the short path are stored in \a distance and \a azimuth.
483 *
484 * \sa distance_long_path(), azimuth_long_path()
485 */
486
qrb(double lon1,double lat1,double lon2,double lat2,double * distance,double * azimuth)487 int qrb(double lon1, double lat1, double lon2, double lat2, double *distance, double *azimuth) {
488 double delta_long, tmp, arc, az;
489
490 /* bail if NULL pointers passed */
491 if (!distance || !azimuth)
492 return -QRB_EINVAL;
493
494 if ((lat1 > 90.0 || lat1 < -90.0) || (lat2 > 90.0 || lat2 < -90.0))
495 return -QRB_EINVAL;
496
497 if ((lon1 > 180.0 || lon1 < -180.0) || (lon2 > 180.0 || lon2 < -180.0))
498 return -QRB_EINVAL;
499
500 /* Prevent ACOS() Domain Error */
501 if (lat1 == 90.0)
502 lat1 = 89.999999999;
503 else if (lat1 == -90.0)
504 lat1 = -89.999999999;
505
506 if (lat2 == 90.0)
507 lat2 = 89.999999999;
508 else if (lat2 == -90.0)
509 lat2 = -89.999999999;
510
511 /* Convert variables to Radians */
512 lat1 /= RADIAN;
513 lon1 /= RADIAN;
514 lat2 /= RADIAN;
515 lon2 /= RADIAN;
516
517 delta_long = lon2 - lon1;
518
519 tmp = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(delta_long);
520
521 if (tmp > .999999999999999) {
522 /* Station points coincide, use an Omni! */
523 *distance = 0.0;
524 *azimuth = 0.0;
525 return QRB_OK;
526 }
527
528 if (tmp < -.999999) {
529 /*
530 * points are antipodal, it's straight down.
531 * Station is equal distance in all Azimuths.
532 * So take 180 Degrees of arc times 60 nm,
533 * and you get 10800 nm, or whatever units...
534 */
535 *distance = M_PI * ARC_IN_KM;
536 *azimuth = 0.0;
537 return QRB_OK;
538 }
539
540 arc = acos(tmp);
541
542 /*
543 * One degree of arc is 60 Nautical miles
544 * at the surface of the earth, 111.2 km, or 69.1 sm
545 * This method is easier than the one in the handbook
546 */
547
548
549 *distance = arc * ARC_IN_KM;
550
551 /* Short Path */
552 /* Change to azimuth computation by Dave Freese, W1HKJ */
553
554 az = RADIAN * atan2(sin(lon2 - lon1) * cos(lat2),
555 (cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon2 - lon1)));
556
557 az = fmod(360.0 + az, 360.0);
558 if (az < 0.0)
559 az += 360.0;
560 else if (az >= 360.0)
561 az -= 360.0;
562
563 *azimuth = floor(az + 0.5);
564 return QRB_OK;
565 }
566
567 /**
568 * \brief Calculate the long path distance between two points.
569 * \param distance The shortpath distance
570 *
571 * Calculate the long path (respective of the short path)
572 * of a given distance.
573 *
574 * \return the distance in kilometers for the opposite path.
575 *
576 * \sa qrb()
577 */
578
distance_long_path(double distance)579 double distance_long_path(double distance) {
580 return (ARC_IN_KM * 2.0 * M_PI) - distance;
581 }
582
583 /**
584 * \brief Calculate the long path bearing between two points.
585 * \param azimuth The shortpath bearing
586 *
587 * Calculate the long path (respective of the short path)
588 * of a given bearing.
589 *
590 * \return the azimuth in decimal degrees for the opposite path.
591 *
592 * \sa qrb()
593 */
594
azimuth_long_path(double azimuth)595 double azimuth_long_path(double azimuth) {
596 return azimuth + (azimuth <= 180.0 ? 180.0 : -180.0);
597 }
598
599 } // namespace QRB
600
601 /*! @} */
602