area_sphere.c 1.7 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374
  1. /*!
  2. * \file lib/gis/area_sphere.c
  3. *
  4. * \brief GIS Library - Sphereical area calculation routines.
  5. *
  6. * (C) 2001-2009 by the GRASS Development Team
  7. *
  8. * This program is free software under the GNU General Public License
  9. * (>=v2). Read the file COPYING that comes with GRASS for details.
  10. *
  11. * \author Original author CERL
  12. */
  13. #include <math.h>
  14. #include <grass/gis.h>
  15. #include "pi.h"
  16. static struct state {
  17. double M;
  18. } state;
  19. static struct state *st = &state;
  20. /*!
  21. * \brief Initialize calculations for sphere.
  22. *
  23. * Initializes raster area calculations for a sphere.
  24. * The radius of the sphere is <i>r</i> and <i>s</i> is a scale factor to
  25. * allow for calculations of a part of the zone (see
  26. * G_begin_zone_area_on_ellipsoid()).
  27. *
  28. * \param r radius of sphere
  29. * \param s scale factor
  30. */
  31. void G_begin_zone_area_on_sphere(double r, double s)
  32. {
  33. st->M = s * 2.0 * r * r * M_PI;
  34. }
  35. /*!
  36. * \brief Calculates integral for area between two latitudes.
  37. *
  38. * \param lat latitude
  39. *
  40. * \return area value
  41. */
  42. double G_darea0_on_sphere(double lat)
  43. {
  44. return (st->M * sin(Radians(lat)));
  45. }
  46. /*!
  47. * \brief Calculates area between latitudes.
  48. *
  49. * This routine shows how to calculate area between two lats, but
  50. * isn't efficient for row by row since G_darea0_on_sphere() will
  51. * be called twice for the same lat, once as a <i>south</i> then
  52. * again as a <i>north</i>.
  53. *
  54. * Returns the area between latitudes <i>north</i> and <i>south</i>
  55. * scaled by the factor <i>s</i> passed to
  56. * G_begin_zone_area_on_sphere().
  57. *
  58. * \param north
  59. * \param[in] south
  60. * \return double
  61. */
  62. double G_area_for_zone_on_sphere(double north, double south)
  63. {
  64. return (G_darea0_on_sphere(north) - G_darea0_on_sphere(south));
  65. }