geodesic.c 2.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102
  1. #include <math.h>
  2. #include <grass/gis.h>
  3. #include "pi.h"
  4. /*
  5. * This code is preliminary. I don't know if it is even
  6. * correct.
  7. */
  8. /*
  9. * From "Map Projections" by Peter Richardus and Ron K. Alder, 1972
  10. * (526.8 R39m in Map & Geography Library)
  11. * page 19, formula 2.17
  12. *
  13. * Formula is the equation of a geodesic from (lat1,lon1) to (lat2,lon2)
  14. * Input is lon, output is lat (all in degrees)
  15. *
  16. * Note formula only works if 0 < abs(lon2-lon1) < 180
  17. * If lon1 == lon2 then geodesic is the merdian lon1
  18. * (and the formula will fail)
  19. * if lon2-lon1=180 then the geodesic is either meridian lon1 or lon2
  20. */
  21. /* TODO:
  22. *
  23. * integrate code from raster/r.surf.idw/ll.c
  24. */
  25. static int adjust_lat(double *);
  26. static int adjust_lon(double *);
  27. static struct state {
  28. double A, B;
  29. } state;
  30. static struct state *st = &state;
  31. int G_begin_geodesic_equation(double lon1, double lat1, double lon2,
  32. double lat2)
  33. {
  34. double sin21, tan1, tan2;
  35. adjust_lon(&lon1);
  36. adjust_lon(&lon2);
  37. adjust_lat(&lat1);
  38. adjust_lat(&lat2);
  39. if (lon1 > lon2) {
  40. double temp;
  41. temp = lon1; lon1 = lon2; lon2 = temp;
  42. temp = lat1; lat1 = lat2; lat2 = temp;
  43. }
  44. if (lon1 == lon2) {
  45. st->A = st->B = 0.0;
  46. return 0;
  47. }
  48. lon1 = Radians(lon1);
  49. lon2 = Radians(lon2);
  50. lat1 = Radians(lat1);
  51. lat2 = Radians(lat2);
  52. sin21 = sin(lon2 - lon1);
  53. tan1 = tan(lat1);
  54. tan2 = tan(lat2);
  55. st->A = (tan2 * cos(lon1) - tan1 * cos(lon2)) / sin21;
  56. st->B = (tan2 * sin(lon1) - tan1 * sin(lon2)) / sin21;
  57. return 1;
  58. }
  59. /* only works if lon1 < lon < lon2 */
  60. double G_geodesic_lat_from_lon(double lon)
  61. {
  62. adjust_lon(&lon);
  63. lon = Radians(lon);
  64. return Degrees(atan(st->A * sin(lon) - st->B * cos(lon)));
  65. }
  66. static int adjust_lon(double *lon)
  67. {
  68. while (*lon > 180.0)
  69. *lon -= 360.0;
  70. while (*lon < -180.0)
  71. *lon += 360.0;
  72. return 0;
  73. }
  74. static int adjust_lat(double *lat)
  75. {
  76. if (*lat > 90.0)
  77. *lat = 90.0;
  78. if (*lat < -90.0)
  79. *lat = -90.0;
  80. return 0;
  81. }