src/com/wherehoo/WHGeo.java
/**
     *Calculates the azimuth of a line from <tt>origin</tt> to <tt>end</tt><br>. 
     *@param origin the starting point
     *@param end the ending point
     *@return Angle between a meridian that passes through <tt>origin</tt> and a line from 
     *<tt>origin</tt> to <tt>end</tt>. There are two special cases:<br>
     *If <tt>origin == end</tt> then return 0<br>  
     *If <tt>origin</tt> is a pole then return -1.
     */
Calculates the azimuth of a line from origin to end.
@param origin the starting point
@param end the ending point
@return Angle between a meridian that passes through origin and a line from
origin to end. There are two special cases:
If origin == end then return 0
If origin is a pole then return -1.
Calculates the azimuth of a line from origin to end . @param origin the starting point @param end the ending point @return Angle between a meridian that passes through origin and a line from origin to end . There are two special cases : If origin == end then return 0 If origin is a pole then return - 1 .
public static double heading(Point2D.Double origin, Point2D.Double end){
	double lat1;
	double lat2;
	double lon1;
	double lon2;
	double dlat;
	double dlon;
	double rlat1;
	double rlat2;
	double rlon1;
	double rlon2;
	double drlat;
	double drlon;
	double drlon_gc;
	double pre_theta1=0;
	double pre_theta2=0;
	boolean zero_lon_crossing;
	int w_e;
	int s_n;
	lat1=origin.getX();
	lat2=end.getX();
	lon1=origin.getY();
	lon2=end.getY();
	rlat1=Math.toRadians(lat1); 
	rlat2=Math.toRadians(lat2);
	rlon1=Math.toRadians(lon1);
	rlon2=Math.toRadians(lon2);
	if (lat1==90 || rlat1==-90) return -1;
        zero_lon_crossing=(Math.abs(lon1-lon2)>(360-Math.abs(lon1-lon2)));
	dlat=Math.min(Math.abs(lat2-lat1),360-Math.abs(lat2-lat1));
	dlon=Math.min(Math.abs(lon2-lon1),360-Math.abs(lon2-lon1));
	drlat = Math.toRadians(dlat);
	drlon = Math.toRadians(dlon);
	if (drlat==0 && drlon==0)
	    return 0;
	if (rlat2>rlat1){
	    s_n=1;
	} else {
	    if (rlat2==rlat1){
		s_n=0;
	    } else {
		s_n=-1;
	    }
	}
	if (rlon2>rlon1){
	    w_e=1;
	} else {
	    if (rlon2==rlon1){
		w_e=0;
	    } else {
		w_e=-1;
	    }
	}
	if (zero_lon_crossing) {
	    w_e=w_e*-1;
	}
	try{
	    /*
	      double H=Math.asin(Math.sin(rlat1)*Math.sin(rlat2)+{Math.cos(rlat1)*Math.cos(rlat2)*Math.cos(drlon)});
	      //System.out.println("H :"+H);
	      double Z=Math.abs(Math.asin(Math.cos(rlat2)*Math.sin(drlon)/Math.cos(H)));
	      pre_theta1=Math.toDegrees(Z);
	    */
	    BigDecimal pre_h_1= (new BigDecimal(Math.sin(rlat1))).multiply(new BigDecimal(Math.sin(rlat2)));
	    BigDecimal pre_h_2= (new BigDecimal(Math.cos(rlat1))).multiply(new BigDecimal(Math.cos(rlat2))).multiply(new BigDecimal(Math.cos(drlon)));
	    double H=Math.asin((pre_h_1.add(pre_h_2)).doubleValue());
	    double pre_z=((new BigDecimal(Math.cos(rlat2))).multiply(new BigDecimal(Math.sin(drlon))).divide( new BigDecimal(Math.cos(H)),50,BigDecimal.ROUND_HALF_UP)).doubleValue();
	    double Z=Math.abs(Math.asin(pre_z));
	    pre_theta1=Math.toDegrees(Z);
	} catch (Exception e){
	    System.out.println(e);
	    pre_theta1=90;
	}
	if (!(pre_theta1>=0 && pre_theta1<=360)){
	    pre_theta1=90;
	}
	/*adjust according to mathematical quadrant of the end point
	 *cartesian ccordinates with "origin" point at the origin
	 *and y axis is parallel to a meridian passing through "origin" point
	 */
	if(s_n==1 && w_e==0) {
	    pre_theta2=pre_theta1;
	}
	if(s_n==1 && w_e==1){ 
	    pre_theta2=pre_theta1;
	}
	if(s_n==0 && w_e==1) {
	    pre_theta2=pre_theta1;
	}
	if(s_n==-1 && w_e==1) {
	    pre_theta2=180-pre_theta1;
	}
	if(s_n==-1 && w_e==0) {
	    pre_theta2=180;
	}
	if(s_n==-1 && w_e==-1) {
	    pre_theta2=180+pre_theta1;
	}
	if(s_n==0 && w_e==-1) {
	    pre_theta2=180+pre_theta1;
	}
	if(s_n==1 && w_e==-1) {
	    pre_theta2=360-pre_theta1;
	}
	/**
	 *calculate the precision of the result
	 *the precision is such that for a given set of coordinates
	 *if a client travels towards end point with resulting heading 
	 *returned with given precision, the client will end up no more than 1m off the end point
	 *it should really be calculated by asin(1/distance(origin,end))and finding order of first non zero digit
	 *but I found out that (log base 10 of distance)-1 is very close approximation
	 */
	double dist=Math.max(WHGeo.distance(origin,end),1);
	int precision = Math.max(((int) Math.rint(Math.log(dist)/Math.log(10))),0);
	return ((new BigDecimal(pre_theta2)).divide(new BigDecimal(1),precision,BigDecimal.ROUND_HALF_UP)).doubleValue();
    }
public static double heading ( Point2D . Double origin , Point2D . Double end ) { double lat1 ; double lat2 ; double lon1 ; double lon2 ; double dlat ; double dlon ; double rlat1 ; double rlat2 ; double rlon1 ; double rlon2 ; double drlat ; double drlon ; double drlon_gc ; double pre_theta1 = 0 ; double pre_theta2 = 0 ; boolean zero_lon_crossing ; int w_e ; int s_n ; lat1 = origin . getX ( ) ; lat2 = end . getX ( ) ; lon1 = origin . getY ( ) ; lon2 = end . getY ( ) ; rlat1 = Math . toRadians ( lat1 ) ; rlat2 = Math . toRadians ( lat2 ) ; rlon1 = Math . toRadians ( lon1 ) ; rlon2 = Math . toRadians ( lon2 ) ; if ( lat1 == 90 || rlat1 == - 90 ) return - 1 ; zero_lon_crossing = ( Math . abs ( lon1 - lon2 ) > ( 360 - Math . abs ( lon1 - lon2 ) ) ) ; dlat = Math . min ( Math . abs ( lat2 - lat1 ) , 360 - Math . abs ( lat2 - lat1 ) ) ; dlon = Math . min ( Math . abs ( lon2 - lon1 ) , 360 - Math . abs ( lon2 - lon1 ) ) ; drlat = Math . toRadians ( dlat ) ; drlon = Math . toRadians ( dlon ) ; if ( drlat == 0 && drlon == 0 ) return 0 ; if ( rlat2 > rlat1 ) { s_n = 1 ; } else { if ( rlat2 == rlat1 ) { s_n = 0 ; } else { s_n = - 1 ; } } if ( rlon2 > rlon1 ) { w_e = 1 ; } else { if ( rlon2 == rlon1 ) { w_e = 0 ; } else { w_e = - 1 ; } } if ( zero_lon_crossing ) { w_e = w_e * - 1 ; } try { /*
	      double H=Math.asin(Math.sin(rlat1)*Math.sin(rlat2)+{Math.cos(rlat1)*Math.cos(rlat2)*Math.cos(drlon)});
	      //System.out.println("H :"+H);
	      double Z=Math.abs(Math.asin(Math.cos(rlat2)*Math.sin(drlon)/Math.cos(H)));
	      //System.out.println("Z :"+Z);
	      pre_theta1=Math.toDegrees(Z);
	    */ BigDecimal pre_h_1 = ( new BigDecimal ( Math . sin ( rlat1 ) ) ) . multiply ( new BigDecimal ( Math . sin ( rlat2 ) ) ) ; BigDecimal pre_h_2 = ( new BigDecimal ( Math . cos ( rlat1 ) ) ) . multiply ( new BigDecimal ( Math . cos ( rlat2 ) ) ) . multiply ( new BigDecimal ( Math . cos ( drlon ) ) ) ; double H = Math . asin ( ( pre_h_1 . add ( pre_h_2 ) ) . doubleValue ( ) ) ; double pre_z = ( ( new BigDecimal ( Math . cos ( rlat2 ) ) ) . multiply ( new BigDecimal ( Math . sin ( drlon ) ) ) . divide ( new BigDecimal ( Math . cos ( H ) ) , 50 , BigDecimal . ROUND_HALF_UP ) ) . doubleValue ( ) ; double Z = Math . abs ( Math . asin ( pre_z ) ) ; pre_theta1 = Math . toDegrees ( Z ) ; } catch ( Exception e ) { System . out . println ( e ) ; pre_theta1 = 90 ; } if ( ! ( pre_theta1 >= 0 && pre_theta1 <= 360 ) ) { pre_theta1 = 90 ; } /*adjust according to mathematical quadrant of the end point
	 *cartesian ccordinates with "origin" point at the origin
	 *and y axis is parallel to a meridian passing through "origin" point
	 */ if ( s_n == 1 && w_e == 0 ) { pre_theta2 = pre_theta1 ; } if ( s_n == 1 && w_e == 1 ) { pre_theta2 = pre_theta1 ; } if ( s_n == 0 && w_e == 1 ) { pre_theta2 = pre_theta1 ; } if ( s_n == - 1 && w_e == 1 ) { pre_theta2 = 180 - pre_theta1 ; } if ( s_n == - 1 && w_e == 0 ) { pre_theta2 = 180 ; } if ( s_n == - 1 && w_e == - 1 ) { pre_theta2 = 180 + pre_theta1 ; } if ( s_n == 0 && w_e == - 1 ) { pre_theta2 = 180 + pre_theta1 ; } if ( s_n == 1 && w_e == - 1 ) { pre_theta2 = 360 - pre_theta1 ; } /**
	 *calculate the precision of the result
	 *the precision is such that for a given set of coordinates
	 *if a client travels towards end point with resulting heading 
	 *returned with given precision, the client will end up no more than 1m off the end point
	 *it should really be calculated by asin(1/distance(origin,end))and finding order of first non zero digit
	 *but I found out that (log base 10 of distance)-1 is very close approximation
	 */ double dist = Math . max ( WHGeo . distance ( origin , end ) , 1 ) ; int precision = Math . max ( ( ( int ) Math . rint ( Math . log ( dist ) / Math . log ( 10 ) ) ) , 0 ) ; return ( ( new BigDecimal ( pre_theta2 ) ) . divide ( new BigDecimal ( 1 ) , precision , BigDecimal . ROUND_HALF_UP ) ) . doubleValue ( ) ; }
Calculates the azimuth of a line from <tt>origin</tt> to <tt>end</tt><br>.
Calculates the azimuth of a line from <tt > origin< / tt > to <tt > end< / tt > <br > .
//longitudinal difference translated to angular difference along great circle
 //these are helper variables that will help with quadrant calculations, will take -1,0,1 values;
 // measures come in as degrees, the trig functions want radians
 //check if we have pole origin, if so, return -1
 //zero crossings?
 //calculate differences
 //System.out.println("dlat dlon :" + dlat +" "+dlon);
 //check if origin==end if so, return 0
 //System.out.println("Zero lon crossing "+zero_lon_crossing); 
 //calculate helper vars for determining quadrant (mathematical, has nothing to do with Earth coordinates)
 //s_n
 //w_e
 //System.out.println("pre_h_2 :"+pre_h_2.doubleValue());
 //System.out.println("pre_z: "+pre_z);
 //System.out.println("Z :"+Z);
 //System.out.println("pre_theta1  :"+pre_theta1);    
 //System.out.println("Pre_theta1: "+pre_theta1);
 //0 degrees
 //System.out.println("zero degrees");
 //first quadrant
 //System.out.println("first quadrant");
 //90 degrees
 //System.out.println("90 degrees");
 //second quadrant
 //System.out.println("second quadrant");
 //180 degrees
 //System.out.println("180 degrees");
 //third quadrant
 //System.out.println("third quadrant");
 //270 degrees
 //System.out.println("270 degrees");
 //fourth quadrant
 //System.out.println("fourth quadrant");
 //System.out.println("Precision: "+precision);
 //round off the result
