Given a spherical polygon described by its vertices,find its area.

Code:

/// <summary>
/// Haversine function : hav(x) = (1-cos(x))/2
/// </summary>
/// <param name="x"></param>
/// <returns>Returns the value of Haversine function</returns>
public double Haversine( double x )
{
return ( 1.0 - Math.Cos( x ) ) / 2.0;
}
/// <summary>
/// Compute the Area of a Spherical Polygon
/// </summary>
/// <param name="lat">the latitudes of all vertices(in radian)</param>
/// <param name="lon">the longitudes of all vertices(in radian)</param>
/// <param name="r">spherical radius</param>
/// <returns>Returns the area of a spherical polygon</returns>
public double SphericalPolygonArea( double[ ] lat , double[ ] lon , double r )
{
double lam1 = 0, lam2 = 0, beta1 =0, beta2 = 0, cosB1 =0, cosB2 = 0;
double hav = 0;
double sum = 0;
for( int j = 0 ; j < lat.Length ; j++ )
{
int k = j + 1;
if( j == 0 )
{
lam1 = lon[j];
beta1 = lat[j];
lam2 = lon[j + 1];
beta2 = lat[j + 1];
cosB1 = Math.Cos( beta1 );
cosB2 = Math.Cos( beta2 );
}
else
{
k = ( j + 1 ) % lat.Length;
lam1 = lam2;
beta1 = beta2;
lam2 = lon[k];
beta2 = lat[k];
cosB1 = cosB2;
cosB2 = Math.Cos( beta2 );
}
if( lam1 != lam2 )
{
hav = Haversine( beta2 - beta1 ) +
cosB1 * cosB2 * Haversine( lam2 - lam1 );
double a = 2 * Math.Asin( Math.Sqrt( hav ) );
double b = Math.PI / 2 - beta2;
double c = Math.PI / 2 - beta1;
double s = 0.5 * ( a + b + c );
double t = Math.Tan( s / 2 ) * Math.Tan( ( s - a ) / 2 ) *
Math.Tan( ( s - b ) / 2 ) * Math.Tan( ( s - c ) / 2 );
double excess = Math.Abs( 4 * Math.Atan( Math.Sqrt(
Math.Abs( t ) ) ) );
if( lam2 < lam1 )
{
excess = -excess;
}
sum += excess;
}
}
return Math.Abs( sum ) * r * r;
}

The lat and lon parameters is a one-to-one relation .

For example

(lat[0],lon[0]) is a coordinate pair,and so on.

This theory for this algorithm is as follows:

http://mathworld.wolfram.com/SphericalPolygon.html