ingo@1115: /*
ingo@1115:  * Copyright (c) 2010 by Intevation GmbH
ingo@1115:  *
ingo@1115:  * This program is free software under the LGPL (>=v2.1)
ingo@1115:  * Read the file LGPL.txt coming with the software for details
ingo@1115:  * or visit http://www.gnu.org/licenses/ if it does not exist.
ingo@1115:  */
ingo@1115: 
tim@185: package de.intevation.gnv.utils;
tim@185: 
sascha@779: import com.vividsolutions.jts.geom.Coordinate;
tim@185: import com.vividsolutions.jts.geom.Point;
sascha@362: 
sascha@362: import java.util.List;
tim@185: 
tim@185: /**
ingo@806:  * A helper class to calculate distances between points and coordinates.
sascha@807:  *
sascha@780:  * @author <a href="mailto:tim.englich@intevation.de">Tim Englich</a>
tim@185:  *
tim@185:  */
tim@185: public class DistanceCalculator {
tim@185: 
tim@185:     private final static double flattening = 1.0 / 298.257233563;
sascha@778: 
tim@185:     private final static double earthRadius = 6378137.0 / 1000.0 ;
sascha@778: 
tim@185:     /**
tim@185:      * Constructor
tim@185:      */
tim@185:     public DistanceCalculator() {
tim@185:     }
sascha@778: 
ingo@806: 
ingo@806:     /**
ingo@806:      * Calculates the distance between two points.
ingo@806:      *
ingo@806:      * @param p1 First point.
ingo@806:      * @param p2 Second point.
ingo@806:      * @return the distance.
ingo@806:      */
ingo@296:     public static double calculateDistance(Point p1, Point p2){
sascha@362:         return calculateDistance(p1.getCoordinate(), p2.getCoordinate());
sascha@362:     }
sascha@362: 
ingo@806: 
ingo@806:     /**
ingo@806:      * Calculates the distance between two coordinates.
ingo@806:      *
ingo@806:      * @param p1 First coordinate.
ingo@806:      * @param p2 Second coordinate.
ingo@806:      * @return the distance.
ingo@806:      */
sascha@362:     public static double calculateDistance(Coordinate p1, Coordinate p2){
tim@185:         double resultValue = 0.0;
sascha@778: 
sascha@362:         double b1 = p1.y;
sascha@362:         double b2 = p2.y;
sascha@778: 
sascha@362:         double l1 = p1.x;
sascha@362:         double l2 = p2.x;
sascha@778: 
sascha@778: 
tim@185:         double F = (b1 + b2) / 2.0;
tim@185:         double G = (b1 - b2) / 2.0;
tim@185:         double l = (l1 - l2) / 2.0;
sascha@778: 
tim@185:         F = (Math.PI / 180.0) * F;
tim@185:         G = (Math.PI / 180.0) * G;
tim@185:         l = (Math.PI / 180.0) * l;
sascha@778: 
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))));
sascha@778: 
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))));
sascha@778: 
tim@185:         double w = Math.atan(Math.sqrt((S/C)));
sascha@778: 
tim@185:         double D = 2.0 * w * earthRadius;
sascha@778: 
tim@185:         double R = Math.sqrt((S*C)) / w;
sascha@778: 
tim@185:         double H1 = (3.0 * R - 1.0 ) / (2.0 * C);
tim@185:         double H2 = (3.0 * R + 1.0 ) / (2.0 * S);
sascha@778: 
sascha@778:         resultValue = D * (1 + (flattening * H1 * (Math.sin(F) * Math.sin(F)) *
sascha@778:                                                   (Math.cos(G) * Math.cos(G))) -
sascha@778:                            (flattening * H2 * (Math.cos(F) * Math.cos(F)) *
tim@185:                                               (Math.sin(G) * Math.sin(G))));
sascha@778: 
tim@185:         return resultValue;
tim@185:     }
tim@185: 
ingo@806:     /**
ingo@806:      * Calculates the length of a path specified by coordinates in <i>path</i>.
ingo@806:      *
ingo@806:      * @param path A list of coordinates.
ingo@806:      * @return the length of the given path.
ingo@806:      */
sascha@362:     public static final double calculateDistance(List<Coordinate> path) {
sascha@362:         int N = path.size();
sascha@362:         if (N < 2) {
sascha@362:             return 0d;
sascha@362:         }
sascha@362:         double sum = 0d;
sascha@362:         for (int i = 1; i < N; ++i) {
sascha@362:             sum += calculateDistance(path.get(i-1), path.get(i));
sascha@362:         }
sascha@362:         return sum;
sascha@362:     }
sascha@362: 
tim@185: }
sascha@836: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :