Application: gvSIG desktop

package com.iver.cit.gvsig.fmap.tools.geo;

/**

 * GeoC is a convenience class used for a couple of things. It can

 * be used to specifiy constants, and can be used for calculate the geographical area.

 *

 * @author Vicente Caballero Navarro

 */

public class Geo {

    public static double HalfPi = 1.5707963267948966192313;

    public static double Degree = 57.295779513082320876798; /* degrees per radian */

    public static double SqMi = 273218.4; /* Square mi per spherical degree. */

    public static double SqKm = 707632.4; /* Square km per spherical degree. */

    public static double SqM = 707632400000.0; /* Square M per spherical degree. */

    public static double getDecimalDegrees(double m) {

            ///(m*180)/ (6378137.0 * Math.PI)

            return (m*8.983152841195214E-6);

    }

    /*  Haversine function: hav(x)= (1-cos(x))/2.  */

    private static double hav(double X){

        return (1.0 - Math.cos(X)) / 2.0;

    }

    /*  Returns the area of a spherical polygon in spherical degrees,

    given the latitudes and longitudes in Lat and Lon, respectively.

    The N data points have indices which range from 0 to N-1. */

    public static double sphericalPolyArea(double[] lat, double[] lon, int n) {

        int j;

        int k;

        double lam1;

        double lam2 = 0;

        double beta1;

        double beta2 = 0;

        double cosB1;

        double cosB2 = 0;

        double havA;

        double t;

        double a;

        double b;

        double c;

        double s;

        double sum;

        double excess;

        sum = 0;

        for (j = 0; j <= n; j++) {

            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) % (n + 1);

                lam1 = lam2;

                beta1 = beta2;

                lam2 = lon[k];

                beta2 = lat[k];

                cosB1 = cosB2;

                cosB2 = Math.cos(beta2);

            }

            if (lam1 != lam2) {

                havA = hav(beta2 - beta1) + (cosB1 * cosB2 * hav(lam2 - lam1));

                a = 2 * Math.asin(Math.sqrt(havA));

                b = HalfPi - beta2;

                c = HalfPi - beta1;

                s = 0.5 * (a + b + c);

                t = Math.tan(s / 2) * Math.tan((s - a) / 2) * Math.tan((s - b) / 2) * Math.tan((s -

                        c) / 2);

                excess = Math.abs(4 * Math.atan(Math.sqrt(Math.abs(t)))) * Degree;

                if (lam2 < lam1) {

                    excess = -excess;

                }

                sum = sum + excess;

            }

        }

        return Math.abs(sum);

    } /*        SphericalPolyArea. */

}