00001
00017 #include <math.h>
00018 #include <grass/gis.h>
00019 #include "pi.h"
00020
00021
00022 #define TWOPI M_PI + M_PI
00023
00024 static double QA, QB, QC;
00025 static double QbarA, QbarB, QbarC, QbarD;
00026
00027 static double AE;
00029 static double Qp;
00031 static double E;
00034 static double Q(double x)
00035 {
00036 double sinx, sinx2;
00037
00038 sinx = sin(x);
00039 sinx2 = sinx * sinx;
00040
00041 return sinx * (1 + sinx2 * (QA + sinx2 * (QB + sinx2 * QC)));
00042 }
00043
00044 static double Qbar(double x)
00045 {
00046 double cosx, cosx2;
00047
00048 cosx = cos(x);
00049 cosx2 = cosx * cosx;
00050
00051 return cosx * (QbarA + cosx2 * (QbarB + cosx2 * (QbarC + cosx2 * QbarD)));
00052 }
00053
00054
00067 int G_begin_ellipsoid_polygon_area(double a, double e2)
00068 {
00069 double e4, e6;
00070
00071 e4 = e2 * e2;
00072 e6 = e4 * e2;
00073
00074 AE = a * a * (1 - e2);
00075
00076 QA = (2.0 / 3.0) * e2;
00077 QB = (3.0 / 5.0) * e4;
00078 QC = (4.0 / 7.0) * e6;
00079
00080 QbarA = -1.0 - (2.0 / 3.0) * e2 - (3.0 / 5.0) * e4 - (4.0 / 7.0) * e6;
00081 QbarB = (2.0 / 9.0) * e2 + (2.0 / 5.0) * e4 + (4.0 / 7.0) * e6;
00082 QbarC = -(3.0 / 25.0) * e4 - (12.0 / 35.0) * e6;
00083 QbarD = (4.0 / 49.0) * e6;
00084
00085 Qp = Q(M_PI_2);
00086 E = 4 * M_PI * Qp * AE;
00087 if (E < 0.0)
00088 E = -E;
00089
00090 return 0;
00091 }
00092
00093
00110 double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
00111 {
00112 double x1, y1, x2, y2, dx, dy;
00113 double Qbar1, Qbar2;
00114 double area;
00115
00116 x2 = Radians(lon[n - 1]);
00117 y2 = Radians(lat[n - 1]);
00118 Qbar2 = Qbar(y2);
00119
00120 area = 0.0;
00121
00122 while (--n >= 0) {
00123 x1 = x2;
00124 y1 = y2;
00125 Qbar1 = Qbar2;
00126
00127 x2 = Radians(*lon++);
00128 y2 = Radians(*lat++);
00129 Qbar2 = Qbar(y2);
00130
00131 if (x1 > x2)
00132 while (x1 - x2 > M_PI)
00133 x2 += TWOPI;
00134 else if (x2 > x1)
00135 while (x2 - x1 > M_PI)
00136 x1 += TWOPI;
00137
00138 dx = x2 - x1;
00139 area += dx * (Qp - Q(y2));
00140
00141 if ((dy = y2 - y1) != 0.0)
00142 area += dx * Q(y2) - (dx / dy) * (Qbar2 - Qbar1);
00143 }
00144 if ((area *= AE) < 0.0)
00145 area = -area;
00146
00147
00148
00149
00150
00151
00152 if (area > E)
00153 area = E;
00154 if (area > E / 2)
00155 area = E - area;
00156
00157 return area;
00158 }