public static double ellipsoidalDistance(LatLon p1,LatLon p2,double equatorialRadius,double polarRadius){
  final double F=(equatorialRadius - polarRadius) / equatorialRadius;
  final double R=1.0 - F;
  final double EPS=0.5E-13;
  if (p1 == null || p2 == null) {
    throw new IllegalArgumentException("Position Is Null");
  }
  double GLAT1=p1.getLatitude().radians;
  double GLAT2=p2.getLatitude().radians;
  double TU1=R * Math.sin(GLAT1) / Math.cos(GLAT1);
  double TU2=R * Math.sin(GLAT2) / Math.cos(GLAT2);
  double CU1=1. / Math.sqrt(TU1 * TU1 + 1.);
  double SU1=CU1 * TU1;
  double CU2=1. / Math.sqrt(TU2 * TU2 + 1.);
  double S=CU1 * CU2;
  double BAZ=S * TU2;
  double FAZ=BAZ * TU1;
  double GLON1=p1.getLongitude().radians;
  double GLON2=p2.getLongitude().radians;
  double X=GLON2 - GLON1;
  double D, SX, CX, SY, CY, Y, SA, C2A, CZ, E, C;
  do {
    SX=Math.sin(X);
    CX=Math.cos(X);
    TU1=CU2 * SX;
    TU2=BAZ - SU1 * CU2 * CX;
    SY=Math.sqrt(TU1 * TU1 + TU2 * TU2);
    CY=S * CX + FAZ;
    Y=Math.atan2(SY,CY);
    SA=S * SX / SY;
    C2A=-SA * SA + 1.;
    CZ=FAZ + FAZ;
    if (C2A > 0.) {
      CZ=-CZ / C2A + CY;
    }
    E=CZ * CZ * 2. - 1.;
    C=((-3. * C2A + 4.) * F + 4.) * C2A * F / 16.;
    D=X;
    X=((E * CY * C + CZ) * SY * C + Y) * SA;
    X=(1. - C) * X * F + GLON2 - GLON1;
  }
 while (Math.abs(D - X) > EPS);
  X=Math.sqrt((1. / R / R - 1.) * C2A + 1.) + 1.;
  X=(X - 2.) / X;
  C=1. - X;
  C=(X * X / 4. + 1.) / C;
  D=(0.375 * X * X - 1.) * X;
  X=E * CY;
  S=1. - E - E;
  S=((((SY * SY * 4. - 3.) * S * CZ* D / 6. - X) * D / 4. + CZ) * SY * D + Y) * C * equatorialRadius* R;
  return S;
}

public static double ellipsoidalDistance ( LatLon p1 , LatLon p2 , double equatorialRadius , double polarRadius ) { final double F = ( equatorialRadius - polarRadius ) / equatorialRadius ; final double R = 1.0 - F ; final double EPS = 0.5E-13 ; if ( p1 == null || p2 == null ) { throw new IllegalArgumentException ( "Position Is Null" ) ; } double GLAT1 = p1 . getLatitude ( ) . radians ; double GLAT2 = p2 . getLatitude ( ) . radians ; double TU1 = R * Math . sin ( GLAT1 ) / Math . cos ( GLAT1 ) ; double TU2 = R * Math . sin ( GLAT2 ) / Math . cos ( GLAT2 ) ; double CU1 = 1. / Math . sqrt ( TU1 * TU1 + 1. ) ; double SU1 = CU1 * TU1 ; double CU2 = 1. / Math . sqrt ( TU2 * TU2 + 1. ) ; double S = CU1 * CU2 ; double BAZ = S * TU2 ; double FAZ = BAZ * TU1 ; double GLON1 = p1 . getLongitude ( ) . radians ; double GLON2 = p2 . getLongitude ( ) . radians ; double X = GLON2 - GLON1 ; double D , SX , CX , SY , CY , Y , SA , C2A , CZ , E , C ; do { SX = Math . sin ( X ) ; CX = Math . cos ( X ) ; TU1 = CU2 * SX ; TU2 = BAZ - SU1 * CU2 * CX ; SY = Math . sqrt ( TU1 * TU1 + TU2 * TU2 ) ; CY = S * CX + FAZ ; Y = Math . atan2 ( SY , CY ) ; SA = S * SX / SY ; C2A = - SA * SA + 1. ; CZ = FAZ + FAZ ; if ( C2A > 0. ) { CZ = - CZ / C2A + CY ; } E = CZ * CZ * 2. - 1. ; C = ( ( - 3. * C2A + 4. ) * F + 4. ) * C2A * F / 16. ; D = X ; X = ( ( E * CY * C + CZ ) * SY * C + Y ) * SA ; X = ( 1. - C ) * X * F + GLON2 - GLON1 ; } while ( Math . abs ( D - X ) > EPS ) ; X = Math . sqrt ( ( 1. / R / R - 1. ) * C2A + 1. ) + 1. ; X = ( X - 2. ) / X ; C = 1. - X ; C = ( X * X / 4. + 1. ) / C ; D = ( 0.375 * X * X - 1. ) * X ; X = E * CY ; S = 1. - E - E ; S = ( ( ( ( SY * SY * 4. - 3. ) * S * CZ * D / 6. - X ) * D / 4. + CZ ) * SY * D + Y ) * C * equatorialRadius * R ; return S ; }
computes the distance between two points on an ellipsoid iteratively . < p / > note : this method was copied from the unidata netcdf java library . http : / / www . unidata . ucar . edu / software / netcdf - java / < p / > algorithm from u . s . national geodetic survey , fortran program " inverse , " subroutine " inver1 , " by l . pfeifer and john g . gergen . see http : / / www . ngs . noaa . gov / tools / inv_fwd / inv_fwd . html < p / > original documentation : solution of the geodetic inverse problem after t . vincenty modified rainsford ' s method with helmert ' s elliptical terms effective in any azimuth and at any distance short of antipodal standpoint / forepoint must not be the geographic pole < p / > requires close to 1 . 4 e - 5 seconds wall clock time per call on a 550 mhz pentium with linux 7 . 2 .
