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