1
/* -*- Mode: C; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
5
* \brief locator and bearing conversion interface
6
* \author Stephane Fillod and the Hamlib Group
9
* Hamlib Interface - locator, bearing, and conversion calls
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
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
23
* Code to calculate distance and azimuth between two Maidenhead locators,
24
* taken from wwl, by IK0ZSN Mirko Caserta.
26
* New bearing code added by N0NB was found at:
27
* http://williams.best.vwh.net/avform.htm#Crs
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.
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.
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.
47
/*! \page hamlib Hamlib general purpose API
49
* Here are grouped some often used functions, like locator conversion
54
#include "build-config.h"
69
#define RADIAN (180.0 / M_PI)
71
/* arc length for 1 degree, 60 Nautical Miles */
72
#define ARC_IN_KM 111.2
74
/* The following is contributed by Dave Hines M1CXW
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
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.
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 };
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 };
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
108
const static int loc_char_range[] = { 18, 10, 24, 10, 24, 10 };
110
#define MAX_LOCATOR_PAIRS 6
111
#define MIN_LOCATOR_PAIRS 1
115
#endif /* !DOC_HIDDEN */
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
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.
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
133
* \return The angle in decimal degrees.
138
double dms2dec (int degrees, int minutes, double seconds, int sw) {
142
degrees = abs(degrees);
144
minutes = abs(minutes);
146
seconds = fabs(seconds);
148
st = (double)degrees + (double)minutes / 60. + seconds / 3600.;
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
162
* Convert a degrees, decimal minutes notation common on
163
* many GPS units to its decimal degrees value.
165
* \a degrees > 360, \a minutes > 60.0 are allowed, but
166
* resulting angle won't be normalized.
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
173
* \return The angle in decimal degrees.
178
double dmmm2dec (int degrees, double minutes, int sw) {
182
degrees = abs(degrees);
184
minutes = fabs(minutes);
186
st = (double)degrees + minutes / 60.;
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
202
* Convert decimal degrees angle into its degree/minute/second
205
* When \a dec < -180 or \a dec > 180, the angle will be normalized
206
* within these limits and the sign set appropriately.
208
* Upon return dec2dms guarantees 0 >= \a degrees <= 180,
209
* 0 >= \a minutes < 60, and 0.0 >= \a seconds < 60.0.
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
216
* \retval -RIG_EINVAL if any of the pointers are NULL.
217
* \retval RIG_OK if conversion went OK.
222
int dec2dms (double dec, int *degrees, int *minutes, double *seconds, int *sw) {
226
/* bail if NULL pointers passed */
227
if (!degrees || !minutes || !seconds || !sw)
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.
238
st = fmod(dec + 180, 360) - 180;
240
st = fmod(dec - 180, 360) + 180;
242
/* if after all of that st is negative, we want deg
243
* to be negative as well except for 180 which we want
246
if (st < 0.0 && st != -180)
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.
259
deg = (int)floor(st);
260
st = 60. * (st - (double)deg);
261
min = (int)floor(st);
262
st = 60. * (st - (double)min);
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
278
* Convert a decimal angle into its degree, decimal minute
279
* notation common on many GPS units.
281
* When passed a value < -180 or > 180, the value will be normalized
282
* within these limits and the sign set apropriately.
284
* Upon return dec2dmmm guarantees 0 >= \a degrees <= 180,
285
* 0.0 >= \a minutes < 60.0.
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
292
* \retval -RIG_EINVAL if any of the pointers are NULL.
293
* \retval RIG_OK if conversion went OK.
297
int dec2dmmm (double dec, int *degrees, double *minutes, int *sw) {
301
/* bail if NULL pointers passed */
302
if (!degrees || !minutes || !sw)
305
r = dec2dms(dec, degrees, &min, &sec, sw);
309
*minutes = (double)min + sec / 60;
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
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
325
* Decimal long/lat is computed to center of grid square, i.e. given
326
* EM19 will return coordinates equivalent to the southwest corner
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.
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
337
* \sa longlat2locator()
342
int locator2longlat (double *longitude, double *latitude, const char *locator) {
343
int x_or_y, paircount;
346
double xy[2], ordinate;
348
/* bail if NULL pointers passed */
349
if (!longitude || !latitude)
352
paircount = strlen(locator) / 2;
354
/* verify paircount is within limits */
355
if (paircount > MAX_LOCATOR_PAIRS)
356
paircount = MAX_LOCATOR_PAIRS;
357
else if (paircount < MIN_LOCATOR_PAIRS)
360
/* For x(=longitude) and y(=latitude) */
361
for (x_or_y = 0; x_or_y < 2; ++x_or_y) {
365
for (pair = 0; pair < paircount; ++pair) {
366
locvalue = locator[pair*2 + x_or_y];
368
/* Value of digit or letter */
369
locvalue -= (loc_char_range[pair] == 10) ? '0' :
370
(isupper(locvalue)) ? 'A' : 'a';
372
/* Check range for non-letter/digit or out of range */
373
if ((locvalue < 0) || (locvalue >= loc_char_range[pair]))
376
divisions *= loc_char_range[pair];
377
ordinate += locvalue * 180.0 / divisions;
379
/* Center ordinate in the Maidenhead "square" or "subsquare" */
380
ordinate += 90.0 / divisions;
382
xy[x_or_y] = ordinate;
385
*longitude = xy[0] * 2.0;
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
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'.
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.
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
411
* \sa locator2longlat()
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;
422
if (pair_count < MIN_LOCATOR_PAIRS || pair_count > MAX_LOCATOR_PAIRS)
425
for (x_or_y = 0; x_or_y < 2; ++x_or_y) {
426
ordinate = (x_or_y == 0) ? longitude / 2.0 : latitude;
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;
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;
441
locator[pair_count * 2] = '\0';
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
457
* Calculate the QRB between \a lon1, \a lat1 and \a lon2, \a lat2.
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).
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.
468
* \return The distance in kilometers and azimuth in decimal degrees
469
* for the short path are stored in \a distance and \a azimuth.
471
* \sa distance_long_path(), azimuth_long_path()
473
int qrb (double lon1, double lat1, double lon2, double lat2,
474
double *distance, double *azimuth)
477
double delta_long, tmp, arc, az;
479
/* bail if NULL pointers passed */
480
if (!distance || !azimuth)
483
if ((lat1 > 90.0 || lat1 < -90.0) || (lat2 > 90.0 || lat2 < -90.0))
486
if ((lon1 > 180.0 || lon1 < -180.0) || (lon2 > 180.0 || lon2 < -180.0))
489
/* Prevent ACOS() Domain Error */
492
else if (lat1 == -90.0)
493
lat1 = -89.999999999;
497
else if (lat2 == -90.0)
498
lat2 = -89.999999999;
500
/* Convert variables to Radians */
506
delta_long = lon2 - lon1;
508
tmp = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(delta_long);
510
if (tmp > .999999999999999) {
511
/* Station points coincide, use an Omni! */
517
if (tmp < -.999999) {
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...
524
*distance = 180.0 * ARC_IN_KM;
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
538
*distance = ARC_IN_KM * RADIAN * arc;
540
/* This formula seems to work with very small distances
542
* I found it on the Web at:
543
* http://williams.best.vwh.net/avform.htm#Crs
545
* Strangely, all the computed values were negative thus the
546
* sign reversal below.
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);
567
* \brief Calculate the long path distance between two points.
568
* \param distance The shortpath distance
570
* Calculate the long path (respective of the short path)
571
* of a given distance.
573
* \return the distance in kilometers for the opposite path.
578
double distance_long_path (double distance) {
579
return (ARC_IN_KM * 360.0) - distance;
583
* \brief Calculate the long path bearing between two points.
584
* \param azimuth The shortpath bearing
586
* Calculate the long path (respective of the short path)
587
* of a given bearing.
589
* \return the azimuth in decimal degrees for the opposite path.
594
double azimuth_long_path (double azimuth) {
595
return 360.0 - azimuth;