15 #include <grass/gis.h>
18 #define TWOPI M_PI + M_PI
22 double QbarA, QbarB, QbarC, QbarD;
30 static double Q(
double x)
37 return sinx * (1 + sinx2 * (
st->QA + sinx2 * (
st->QB + sinx2 *
st->QC)));
40 static double Qbar(
double x)
47 return cosx * (
st->QbarA + cosx2 * (
st->QbarB + cosx2 * (
st->QbarC + cosx2 *
st->QbarD)));
68 st->AE = a * a * (1 - e2);
70 st->QA = (2.0 / 3.0) * e2;
71 st->QB = (3.0 / 5.0) * e4;
72 st->QC = (4.0 / 7.0) * e6;
74 st->QbarA = -1.0 - (2.0 / 3.0) * e2 - (3.0 / 5.0) * e4 - (4.0 / 7.0) * e6;
75 st->QbarB = (2.0 / 9.0) * e2 + (2.0 / 5.0) * e4 + (4.0 / 7.0) * e6;
76 st->QbarC = -(3.0 / 25.0) * e4 - (12.0 / 35.0) * e6;
77 st->QbarD = (4.0 / 49.0) * e6;
80 st->E = 4 * M_PI *
st->Qp *
st->AE;
131 double x1, y1, x2, y2, dx, dy;
134 double thresh = 1e-6;
152 while (x1 - x2 > M_PI)
155 while (x2 - x1 > M_PI)
161 if (fabs(dy) > thresh) {
163 area += dx * (
st->Qp - (Qbar2 - Qbar1) / dy);
176 area += dx * (
st->Qp - Q((y1 + y2) / 2));
179 if ((area *=
st->AE) < 0.0)
189 if (area >
st->E / 2)