/*
 * $Header: /cvs/CVS_LEBA/external/openmap/openmap/src/openmap/com/bbn/openmap/geo/RhumbCalculator.java,v 1.1 2004/04/18 14:38:00 pitek Exp $
 *
 * Copyright 2004 OBR Centrum Techniki Morskiej, All rights reserved.
 *
 */

package com.bbn.openmap.geo;
import java.lang.Math;
import com.bbn.openmap.LatLonPoint;

/**
 *
 * @version $Header: /cvs/CVS_LEBA/external/openmap/openmap/src/openmap/com/bbn/openmap/geo/RhumbCalculator.java,v 1.1 2004/04/18 14:38:00 pitek Exp $
 * @author pawkub
 *
 * Klasa zawiera metody pozwalające przeprowadzać obliczenia związane z
 * tzw. Rhumblinem. Klasa zawiera kilka ogólnych wzorów i poza OpenMapą nie jest 
 * z niczym związana.
 * Można ją wrzucić do OpenMapy i użyć zawartych w niej wzorów tam gdzie to będzie potrzebne
 * (a jest takich miejsc kilka).
 * @todo Tą klasę należy potraktować jako narzędzie i wrzucić do OpenMapy.
 */
public class RhumbCalculator {
    
    /** Konstruktor jest prywatny, bo (na razie) klasa zawiera tylko kilka metod statycznych. */
    private RhumbCalculator() {
    }

    /* Metoda pozwala pbliczyć położenie punktu na Rhumb Line, znając położenie punktu odniesienia,
     * azymut i odległość.
     * @param point Punkt odniesienia
     * @param azimuth azymut w stopniach
     * @oparam dist odległość w metrach
     * @return punkt (obiekt klasy LatLonPoint)
     */
    public static LatLonPoint calculatePointOnRhumbLine(LatLonPoint point, float azimuth, float dist) {
        double lat1 = (double)com.bbn.openmap.geo.Geo.radians(point.getLatitude());
        double lon1 = (double)com.bbn.openmap.geo.Geo.radians(point.getLongitude());
        double d = (double)dist / 1855.3 * Math.PI / 10800.0;
        double lat = 0.0;
        double lon = 0.0;
        lat = lat1 + d * Math.cos(azimuth);
        double dphi= Math.log((1+Math.sin(lat))/Math.cos(lat)) - Math.log((1+Math.sin(lat1))/Math.cos(lat1));        
        double dlon = 0.0;
        if (Math.abs(Math.cos(azimuth)) > Math.sqrt(0.00000000000001)) {
            dlon=dphi*Math.tan(azimuth);
        } else { // along parallel
            dlon=Math.sin(azimuth)*d/Math.cos(lat1);
        }
        lon = mod(lon1-dlon+Math.PI,2*Math.PI)-Math.PI;
        //System.out.println("calculatePointOnRhumbLine:  lat1 = "+lat1+"+  lon1 = "+lon1 + "    lat = "+lat+"+  lon = "+lon);
        return new LatLonPoint((float)lat, (float)lon, true);
    }
    
    /** Metoda pozwala znormować wartość do podanego zakresu
     * @param y normaowana wartość
     * @param x zakres
     * @return normowana wartość
     */
    private static double mod(double y, double x) {
        double ret;
        if (y>=0) {
            ret = y - x * (int)(y/x);
        } else {
            ret = y + x * ((int)(-y/x)+1);
        }
        return ret;
    }
    
    /** Metoda pozwala obliczyć odległość pomiędzy punktami liczoną po rhumbline.
     * @param p1 punkt geograficzny
     * @param p2 punkt geograficzny
     * @return odległość
     */
    public static float getDistanceBetweenPoints(LatLonPoint p1, LatLonPoint p2) {
        //float distWrong = (float) (com.bbn.openmap.geo.Geo.distanceKM(p1.getLatitude(), p1.getLongitude()
        //, p2.getLatitude(), p2.getLongitude()) * 1000.0);        
        double lat1 = (double)com.bbn.openmap.geo.Geo.radians(p1.getLatitude());
        double lon1 = (double)com.bbn.openmap.geo.Geo.radians(p1.getLongitude());
        double lat2 = (double)com.bbn.openmap.geo.Geo.radians(p2.getLatitude());
        double lon2 = (double)com.bbn.openmap.geo.Geo.radians(p2.getLongitude());        
        double d = 0.0;
        double tc = 0.0;        
        double dlon_W = mod(lon2-lon1,2*Math.PI);
        double dlon_E = mod(lon1-lon2, 2*Math.PI);
        double dphi = Math.log((1+Math.sin(lat2))/Math.cos(lat2)) -Math.log ((1+Math.sin(lat1))/Math.cos(lat1));
        tc = getAzimuthBetweenPoints(p1, p2);
        if (Math.abs(lat1-lat2) < Math.sqrt(0.00000000000001)) {
            d = Math.min(dlon_W,dlon_E) * Math.cos(lat1); // distance along parallel
        } else {
            d = Math.abs((lat2-lat1)/Math.cos(tc));            
        }
        float dist = (float)(d * 10800.0 / Math.PI * 1855.3);        
        //System.out.println("DWrong = " + distWrong + "   DOK = " + dist);
        return dist;
    }
 
    /** Metoda pozwala obliczyć azymut pomiędzy punktami (namiar od pierwszego do drugiego)
     * @param p1 pierwszy punkt geograficzny
     * @param p2 drugi punkt geograficzny
     * @return namiar w radianach
     */    
    public static float getAzimuthBetweenPoints(LatLonPoint p1, LatLonPoint p2) {
        double lat1 = (double)com.bbn.openmap.geo.Geo.radians(p1.getLatitude());
        double lon1 = (double)com.bbn.openmap.geo.Geo.radians(p1.getLongitude());
        double lat2 = (double)com.bbn.openmap.geo.Geo.radians(p2.getLatitude());
        double lon2 = (double)com.bbn.openmap.geo.Geo.radians(p2.getLongitude());
        
        double tc = 0.0;
        double dlon_W = mod(lon2-lon1,2*Math.PI);
        double dlon_E = mod(lon1-lon2, 2*Math.PI);
        double dphi = Math.log((1+Math.sin(lat2))/Math.cos(lat2)) -Math.log ((1+Math.sin(lat1))/Math.cos(lat1));
        if (dlon_W<dlon_E) {// West is the shortest
            tc = mod(Math.atan2(-dlon_W,dphi),2*Math.PI);
        } else {
            tc = mod(Math.atan2(dlon_E,dphi),2*Math.PI);
        }
        return (float)tc;
    }
    
}