00001 #include <math.h>
00002 #include <grass/gis.h>
00003 #include "pi.h"
00004
00005 static double QA, QB, QC;
00006 static double QbarA, QbarB, QbarC, QbarD;
00007 static double AE;
00008 static double Qp;
00009 static double E;
00010 static double TwoPI;
00011
00012 static double Q(double x)
00013 {
00014 double sinx, sinx2;
00015
00016 sinx = sin(x);
00017 sinx2 = sinx * sinx;
00018
00019 return sinx * (1 + sinx2 * (QA + sinx2 * (QB + sinx2 * QC)));
00020 }
00021
00022 static double Qbar(double x)
00023 {
00024 double cosx, cosx2;
00025
00026 cosx = cos(x);
00027 cosx2 = cosx * cosx;
00028
00029 return cosx * (QbarA + cosx2 * (QbarB + cosx2 * (QbarC + cosx2 * QbarD)));
00030 }
00031
00032
00045 int G_begin_ellipsoid_polygon_area (double a,double e2)
00046 {
00047 double e4, e6;
00048
00049 TwoPI = PI+PI;
00050
00051 e4 = e2 * e2;
00052 e6 = e4 * e2;
00053
00054 AE = a * a * (1 - e2);
00055
00056 QA = (2.0/3.0)*e2;
00057 QB = (3.0/5.0)*e4;
00058 QC = (4.0/7.0)*e6;
00059
00060 QbarA = -1.0 - (2.0/3.0)*e2 - (3.0/5.0)*e4 - (4.0/7.0)*e6;
00061 QbarB = (2.0/9.0)*e2 + (2.0/5.0)*e4 + (4.0/7.0)*e6;
00062 QbarC = - (3.0/25.0)*e4 - (12.0/35.0)*e6;
00063 QbarD = (4.0/49.0)*e6;
00064
00065 Qp = Q(PI/2);
00066 E = 4 * PI * Qp * AE;
00067 if (E < 0.0) E = -E;
00068
00069 return 0;
00070 }
00071
00072
00088 double G_ellipsoid_polygon_area (double *lon,double *lat,int n)
00089 {
00090 double x1,y1,x2,y2,dx,dy;
00091 double Qbar1, Qbar2;
00092 double area;
00093
00094 x2 = Radians (lon[n-1]);
00095 y2 = Radians (lat[n-1]);
00096 Qbar2 = Qbar(y2);
00097
00098 area = 0.0;
00099
00100 while (--n >= 0)
00101 {
00102 x1 = x2;
00103 y1 = y2;
00104 Qbar1 = Qbar2;
00105
00106 x2 = Radians (*lon++);
00107 y2 = Radians (*lat++);
00108 Qbar2 = Qbar(y2);
00109
00110 if (x1 > x2)
00111 while (x1 - x2 > PI)
00112 x2 += TwoPI;
00113 else if (x2 > x1)
00114 while (x2 - x1 > PI)
00115 x1 += TwoPI;
00116
00117 dx = x2 - x1;
00118 area += dx * (Qp - Q(y2));
00119
00120 if ((dy = y2 - y1) != 0.0)
00121 area += dx * Q(y2) - (dx/dy)*(Qbar2-Qbar1);
00122 }
00123 if((area *= AE) < 0.0)
00124 area = -area;
00125
00126
00127
00128
00129
00130
00131 if (area > E) area = E;
00132 if (area > E/2) area = E - area;
00133
00134 return area;
00135 }