GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
area_ellipse.c
Go to the documentation of this file.
1 /*!
2  * \file lib/gis/area_ellipse.c
3  *
4  * \brief GIS Library - Ellipse area 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 
14 #include <math.h>
15 #include <grass/gis.h>
16 #include "pi.h"
17 
18 static struct state {
19  double E;
20  double M;
21 } state;
22 
23 static struct state *st = &state;
24 
25 /*
26  * a is semi-major axis, e2 is eccentricity squared, s is a scale factor
27  * code will fail if e2==0 (sphere)
28  */
29 
30 /*!
31  * \brief Begin area calculations for ellipsoid.
32  *
33  * Initializes raster area calculations for an ellipsoid, where <i>a</i>
34  * is the semi-major axis of the ellipse (in meters), <i>e2</i> is the
35  * ellipsoid eccentricity squared, and <i>s</i> is a scale factor to
36  * allow for calculations of part of the zone (<i>s</i>=1.0 is full
37  * zone, <i>s</i>=0.5 is half the zone, and <i>s</i>=360/ew_res is for a
38  * single grid cell).
39  *
40  * <b>Note:</b> <i>e2</i> must be positive. A negative value makes no
41  * sense, and zero implies a sphere.
42  *
43  * \param a semi-major axis
44  * \param e2 ellipsoid eccentricity
45  * \param s scale factor
46  */
47 void G_begin_zone_area_on_ellipsoid(double a, double e2, double s)
48 {
49  st->E = sqrt(e2);
50  st->M = s * a * a * M_PI * (1 - e2) / st->E;
51 }
52 
53 /*!
54  * \brief Calculate integral for area between two latitudes.
55  *
56  * This routine is part of the integral for the area between two
57  * latitudes.
58  *
59  * \param lat latitude
60  *
61  * \return cell area
62  */
63 double G_darea0_on_ellipsoid(double lat)
64 {
65  double x;
66 
67  x = st->E * sin(Radians(lat));
68 
69  return (st->M * (x / (1.0 - x * x) + 0.5 * log((1.0 + x) / (1.0 - x))));
70 }
71 
72 /*!
73  * \brief Calculates area between latitudes.
74  *
75  * This routine shows how to calculate area between two lats, but
76  * isn't efficient for row by row since G_darea0_on_ellipsoid()
77  * will be called twice for the same lat, once as a <i>south</i> then
78  * again as a <i>north</i>.
79  *
80  * Returns the area between latitudes <i>north</i> and <i>south</i>
81  * scaled by the factor <i>s</i> passed to
82  * G_begin_zone_area_on_ellipsoid().
83  *
84  * \param north north coordinate
85  * \param south south coordinate
86  *
87  * \return cell area
88  */
89 double G_area_for_zone_on_ellipsoid(double north, double south)
90 {
91  return (G_darea0_on_ellipsoid(north) - G_darea0_on_ellipsoid(south));
92 }
void G_begin_zone_area_on_ellipsoid(double a, double e2, double s)
Begin area calculations for ellipsoid.
Definition: area_ellipse.c:47
double G_area_for_zone_on_ellipsoid(double north, double south)
Calculates area between latitudes.
Definition: area_ellipse.c:89
#define M_PI
Definition: gis.h:134
#define x
double G_darea0_on_ellipsoid(double lat)
Calculate integral for area between two latitudes.
Definition: area_ellipse.c:63
struct state * st
Definition: parser.c:104
#define Radians(x)
Definition: pi.h:6
#define M(row, col)
Definition: georef.c:49
struct state state
Definition: parser.c:103