tim@185: /**
tim@185:  *
tim@185:  */
tim@185: package de.intevation.gnv.utils;
tim@185: 
tim@185: import com.vividsolutions.jts.geom.Point;
tim@185: 
tim@185: /**
tim@185:  * @author Tim Englich <tim.englich@intevation.de>
tim@185:  *
tim@185:  */
tim@185: public class DistanceCalculator {
tim@185: 
tim@185:     private final static double flattening = 1.0 / 298.257233563;
tim@185:     
tim@185:     private final static double earthRadius = 6378137.0 / 1000.0 ;
tim@185:     
tim@185:     /**
tim@185:      * Constructor
tim@185:      */
tim@185:     public DistanceCalculator() {
tim@185:     }
tim@185:     
ingo@296:     public static double calculateDistance(Point p1, Point p2){
tim@185:         double resultValue = 0.0;
tim@185:         
tim@185:         double b1 = p1.getY();
tim@185:         double b2 = p2.getY();
tim@185:         
tim@185:         double l1 = p1.getX();
tim@185:         double l2 = p2.getX();
tim@185:         
tim@185:         
tim@185:         double F = (b1 + b2) / 2.0;
tim@185:         double G = (b1 - b2) / 2.0;
tim@185:         double l = (l1 - l2) / 2.0;
tim@185:         
tim@185:         F = (Math.PI / 180.0) * F;
tim@185:         G = (Math.PI / 180.0) * G;
tim@185:         l = (Math.PI / 180.0) * l;
tim@185:         
tim@185:         double S = ((Math.sin(G) * Math.sin(G)) * ((Math.cos(l) * Math.cos(l))))+
tim@185:                    ((Math.cos(F) * Math.cos(F)) * ((Math.sin(l) * Math.sin(l))));
tim@185:         
tim@185:         double C = ((Math.cos(G) * Math.cos(G)) * ((Math.cos(l) * Math.cos(l))))+
tim@185:                    ((Math.sin(F) * Math.sin(F)) * ((Math.sin(l) * Math.sin(l))));
tim@185:         
tim@185:         double w = Math.atan(Math.sqrt((S/C)));
tim@185:         
tim@185:         double D = 2.0 * w * earthRadius;
tim@185:         
tim@185:         double R = Math.sqrt((S*C)) / w;
tim@185:         
tim@185:         double H1 = (3.0 * R - 1.0 ) / (2.0 * C);
tim@185:         double H2 = (3.0 * R + 1.0 ) / (2.0 * S);
tim@185:         
tim@185:         resultValue = D * (1 + (flattening * H1 * (Math.sin(F) * Math.sin(F)) * 
tim@185:                                                   (Math.cos(G) * Math.cos(G))) - 
tim@185:                            (flattening * H2 * (Math.cos(F) * Math.cos(F)) * 
tim@185:                                               (Math.sin(G) * Math.sin(G))));
tim@185:         
tim@185:         return resultValue;
tim@185:     }
tim@185: 
tim@185: }