area_poly1.c

Go to the documentation of this file.
00001 #include <math.h>
00002 #include "gis.h"
00003 #include "pi.h"
00004 
00005 static double QA, QB, QC;
00006 static double QbarA, QbarB, QbarC, QbarD;
00007 static double AE;  /* a^2(1-e^2) */
00008 static double Qp;  /* Q at the north pole */
00009 static double E;   /* area of the earth */
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     /* kludge - if polygon circles the south pole the area will be
00127      * computed as if it cirlced the north pole. The correction is
00128      * the difference between total surface area of the earth and
00129      * the "north pole" area.
00130      */
00131     if (area > E) area = E;
00132     if (area > E/2) area = E - area;
00133 
00134     return area;
00135 }

Generated on Mon Jan 1 19:49:25 2007 for GRASS by  doxygen 1.5.1