public double surfaceDistance(final GeoPoint pt1,final GeoPoint pt2){
  final double L=pt2.getLongitude() - pt1.getLongitude();
  final double U1=Math.atan((1.0 - flattening) * Math.tan(pt1.getLatitude()));
  final double U2=Math.atan((1.0 - flattening) * Math.tan(pt2.getLatitude()));
  final double sinU1=Math.sin(U1);
  final double cosU1=Math.cos(U1);
  final double sinU2=Math.sin(U2);
  final double cosU2=Math.cos(U2);
  final double dCosU1CosU2=cosU1 * cosU2;
  final double dCosU1SinU2=cosU1 * sinU2;
  final double dSinU1SinU2=sinU1 * sinU2;
  final double dSinU1CosU2=sinU1 * cosU2;
  double lambda=L;
  double lambdaP=Math.PI * 2.0;
  int iterLimit=0;
  double cosSqAlpha;
  double sinSigma;
  double cos2SigmaM;
  double cosSigma;
  double sigma;
  double sinAlpha;
  double C;
  double sinLambda, cosLambda;
  do {
    sinLambda=Math.sin(lambda);
    cosLambda=Math.cos(lambda);
    sinSigma=Math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) + (dCosU1SinU2 - dSinU1CosU2 * cosLambda) * (dCosU1SinU2 - dSinU1CosU2 * cosLambda));
    if (sinSigma == 0.0) {
      return 0.0;
    }
    cosSigma=dSinU1SinU2 + dCosU1CosU2 * cosLambda;
    sigma=Math.atan2(sinSigma,cosSigma);
    sinAlpha=dCosU1CosU2 * sinLambda / sinSigma;
    cosSqAlpha=1.0 - sinAlpha * sinAlpha;
    cos2SigmaM=cosSigma - 2.0 * dSinU1SinU2 / cosSqAlpha;
    if (Double.isNaN(cos2SigmaM))     cos2SigmaM=0.0;
    C=flattening / 16.0 * cosSqAlpha * (4.0 + flattening * (4.0 - 3.0 * cosSqAlpha));
    lambdaP=lambda;
    lambda=L + (1.0 - C) * flattening * sinAlpha* (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1.0 + 2.0 * cos2SigmaM * cos2SigmaM)));
  }
 while (Math.abs(lambda - lambdaP) > Vector.MINIMUM_RESOLUTION && ++iterLimit < 40);
  final double uSq=cosSqAlpha * this.squareRatio;
  final double A=1.0 + uSq / 16384.0 * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq)));
  final double B=uSq / 1024.0 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq)));
  final double deltaSigma=B * sinSigma * (cos2SigmaM + B / 4.0 * (cosSigma * (-1.0 + 2.0 * cos2SigmaM * cos2SigmaM) - B / 6.0 * cos2SigmaM * (-3.0 + 4.0 * sinSigma * sinSigma) * (-3.0 + 4.0 * cos2SigmaM * cos2SigmaM)));
  return c * A * (sigma - deltaSigma);
}

public double surfaceDistance ( final GeoPoint pt1 , final GeoPoint pt2 ) { final double L = pt2 . getLongitude ( ) - pt1 . getLongitude ( ) ; final double U1 = Math . atan ( ( 1.0 - flattening ) * Math . tan ( pt1 . getLatitude ( ) ) ) ; final double U2 = Math . atan ( ( 1.0 - flattening ) * Math . tan ( pt2 . getLatitude ( ) ) ) ; final double sinU1 = Math . sin ( U1 ) ; final double cosU1 = Math . cos ( U1 ) ; final double sinU2 = Math . sin ( U2 ) ; final double cosU2 = Math . cos ( U2 ) ; final double dCosU1CosU2 = cosU1 * cosU2 ; final double dCosU1SinU2 = cosU1 * sinU2 ; final double dSinU1SinU2 = sinU1 * sinU2 ; final double dSinU1CosU2 = sinU1 * cosU2 ; double lambda = L ; double lambdaP = Math . PI * 2.0 ; int iterLimit = 0 ; double cosSqAlpha ; double sinSigma ; double cos2SigmaM ; double cosSigma ; double sigma ; double sinAlpha ; double C ; double sinLambda , cosLambda ; do { sinLambda = Math . sin ( lambda ) ; cosLambda = Math . cos ( lambda ) ; sinSigma = Math . sqrt ( ( cosU2 * sinLambda ) * ( cosU2 * sinLambda ) + ( dCosU1SinU2 - dSinU1CosU2 * cosLambda ) * ( dCosU1SinU2 - dSinU1CosU2 * cosLambda ) ) ; if ( sinSigma == 0.0 ) { return 0.0 ; } cosSigma = dSinU1SinU2 + dCosU1CosU2 * cosLambda ; sigma = Math . atan2 ( sinSigma , cosSigma ) ; sinAlpha = dCosU1CosU2 * sinLambda / sinSigma ; cosSqAlpha = 1.0 - sinAlpha * sinAlpha ; cos2SigmaM = cosSigma - 2.0 * dSinU1SinU2 / cosSqAlpha ; if ( Double . isNaN ( cos2SigmaM ) ) cos2SigmaM = 0.0 ; C = flattening / 16.0 * cosSqAlpha * ( 4.0 + flattening * ( 4.0 - 3.0 * cosSqAlpha ) ) ; lambdaP = lambda ; lambda = L + ( 1.0 - C ) * flattening * sinAlpha * ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( - 1.0 + 2.0 * cos2SigmaM * cos2SigmaM ) ) ) ; } while ( Math . abs ( lambda - lambdaP ) > Vector . MINIMUM_RESOLUTION && ++ iterLimit < 40 ) ; final double uSq = cosSqAlpha * this . squareRatio ; final double A = 1.0 + uSq / 16384.0 * ( 4096.0 + uSq * ( - 768.0 + uSq * ( 320.0 - 175.0 * uSq ) ) ) ; final double B = uSq / 1024.0 * ( 256.0 + uSq * ( - 128.0 + uSq * ( 74.0 - 47.0 * uSq ) ) ) ; final double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4.0 * ( cosSigma * ( - 1.0 + 2.0 * cos2SigmaM * cos2SigmaM ) - B / 6.0 * cos2SigmaM * ( - 3.0 + 4.0 * sinSigma * sinSigma ) * ( - 3.0 + 4.0 * cos2SigmaM * cos2SigmaM ) ) ) ; return c * A * ( sigma - deltaSigma ) ; }
compute surface distance between two points .
