locator.c
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 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 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 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 appropriately. 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 */ 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 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 */ 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 */ 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 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 594 double azimuth_long_path (double azimuth) { 595 return 360.0 - azimuth; 596 }