Formula for Calculating the area of a spherical polygon under latitude and longitude coordinates

// Calculate Area
Public static double calcarea (arraylist pointx, arraylist pointy, string mapunits)
{

Int COUNT = pointx. count;
If (count> 2)
{
Double mtotalarea = 0;

If (mapunits = "Degrees") // spherical polygon under the coordinate of latitude and longitude
{
Double lowx = 0.0;
Double Lowy = 0.0;
Double middlex = 0.0;
Double middley = 0.0;
Double highx = 0.0;
Double highy = 0.0;

Double AM = 0.0;
Double Bm = 0.0;
Double CM = 0.0;

Double Al = 0.0;
Double BL = 0.0;
Double Cl = 0.0;

Double Ah = 0.0;
Double bH = 0.0;
Double CH = 0.0;

Double coefficientl = 0.0;
Double coefficienth = 0.0;

Double altangent = 0.0;
Double bltangent = 0.0;
Double cltangent = 0.0;

Double ahtangent = 0.0;
Double bhtangent = 0.0;
Double chtangent = 0.0;

Double anormalline = 0.0;
Double bnormalline= 0.0;
Double cnormalline= 0.0;

Double orientationvalue = 0.0;

Double anglecos = 0.0;

Double sum1 = 0.0;
Double sum2 = 0.0;
Double count2 = 0;
Double count1 = 0;

Double sum = 0.0;
Double Radius = 6378000;

For (INT I = 0; I <count; I ++)
{
If (I = 0)
{
Lowx = (double) pointx [count-1] * Math. PI/180;
Lowy = (double) pointy [count-1] * Math. PI/180;
Middlex = (double) pointx [0] * Math. PI/180;
Middley = (double) pointy [0] * Math. PI/180;
Highx = (double) pointx [1] * Math. PI/180;
Highy = (double) pointy [1] * Math. PI/180;
}
Else if (I = count-1)
{
Lowx = (double) pointx [count-2] * Math. PI/180;
Lowy = (double) pointy [count-2] * Math. PI/180;
Middlex = (double) pointx [count-1] * Math. PI/180;
Middley = (double) pointy [count-1] * Math. PI/180;
Highx = (double) pointx [0] * Math. PI/180;
Highy = (double) pointy [0] * Math. PI/180;
}
Else
{
Lowx = (double) pointx [I-1] * Math. PI/180;
Lowy = (double) pointy [I-1] * Math. PI/180;
Middlex = (double) pointx [I] * Math. PI/180;
Middley = (double) pointy [I] * Math. PI/180;
Highx = (double) pointx [I + 1] * Math. PI/180;
Highy = (double) pointy [I + 1] * Math. PI/180;
}

AM = math. Cos (middley) * Math. Cos (middlex );
Bm = math. Cos (middley) * Math. Sin (middlex );
Cm = math. Sin (middley );
Al = math. Cos (Lowy) * Math. Cos (lowx );
BL = math. Cos (Lowy) * Math. Sin (lowx );
CL = math. Sin (Lowy );
Ah = math. Cos (highy) * Math. Cos (highx );
BH = math. Cos (highy) * Math. Sin (highx );
Ch = math. Sin (highy );

Coefficientl = (AM * Am + BM * BM + cm * cm)/(AM * Al + BM * BL + cm * Cl );
Coefficienth = (AM * Am + BM * BM + cm * cm)/(AM * Ah + BM * BH + cm * Ch );

Altangent = coefficientl * al-am;
Bltangent = coefficientl * bl-bm;
Cltangent = coefficientl * Cl-cm;
Ahtangent = coefficienth * Ah-am;
Bhtangent = coefficienth * bh-bm;
Chtangent = coefficienth * Ch-cm;

Anglecos = (ahtangent * altangent + bhtangent * bltangent + chtangent * cltangent)/(math. SQRT (ahtangent * ahtangent + bhtangent * bhtangent + chtangent * chtangent) * Math. SQRT (altangent * altangent + bltangent * bltangent + cltangent * cltangent ));

Anglecos = math. ACOs (anglecos );

Anormalline = bhtangent * cltangent-chtangent * bltangent;
Bnormalline = 0-(ahtangent * cltangent-chtangent * altangent );
Cnormalline = ahtangent * bltangent-bhtangent * altangent;

If (AM! = 0)
Orientationvalue = anormalline/am;
Else if (BM! = 0)
Orientationvalue = bnormalline/BM;
Else
Orientationvalue = cnormalline/cm;

If (orientationvalue> 0)
{
Sum1 + = anglecos;
Count1 ++;

}
Else
{
Sum2 + = anglecos;
Count2 ++;
// Sum + = 2 * Math. Pi-anglecos;
}

}

If (sum1> sum2)
{
Sum = sum1 + (2 * Math. Pi * count2-sum2 );
}
Else
{
Sum = (2 * Math. Pi * count1-sum1) + sum2;
}

// Square meter
Mtotalarea = (Sum-(count-2) * Math. Pi) * radius;
}
Else
{// Plane polygon under non-latitude and longitude coordinates

Int I, J;
// Double J;
Double p1x, p1y;
Double P2x, p2y;
For (I = count-1, j = 0; j <count; I = J, J ++)
{

P1x = (double) pointx [I];
P1y = (double) pointy [I];

P2x = (double) pointx [J];
P2y = (double) pointy [J];

Mtotalarea + = p1x * p2y-P2x * p1y;
}
Mtotalarea/= 2.0;
}
Return mtotalarea;
}
Return 0;
}