import org.geotools.referencing.datum.DefaultEllipsoid; //导入依赖的package包/类
/**
 * Finds the shortest distance between two points on the surface of the Earth. This uses Vincenty's method
 * (assumes the Earth is an oblate ellipsoid). If this method fails to converge, it falls back to the Haversine
 * formula (assumes the Earth is a sphere).
 * @param point1 The first point.
 * @param point2 The second point.
 * @return The orthodromic distance, in kilometres.
 */
public static double findOrthodromicDistance(Point point1, Point point2) {
    // Ideally we would use GeodeticCalculator.getOrthodromicDistance() for this, but its internal call
    // to computeDirection() fails to converge for many real-life points. Instead we call the method on
    // DefaultEllipsoid directly. Its results largely match PostGIS's ST_DISTANCE function applied to geography
    // types.
    try {
        return DefaultEllipsoid.WGS84.orthodromicDistance(point1.getX(), point1.getY(),
                point2.getX(), point2.getY()) / METRES_IN_A_KILOMETRE;
    } catch (ArithmeticException e) {
        // Failed to converge, so fall back to using a sphere (like PostGIS does). Uses the Haversine formula.
        double latDistance = Math.toRadians(point2.getY() - point1.getY());
        double lonDistance = Math.toRadians(point2.getX() - point1.getX());
        double a = Math.sin(latDistance / 2) * Math.sin(latDistance / 2) +
                   Math.cos(Math.toRadians(point1.getY())) * Math.cos(Math.toRadians(point2.getY())) *
                   Math.sin(lonDistance / 2) * Math.sin(lonDistance / 2);
        double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
        return RADIUS_OF_THE_EARTH_IN_KILOMETRES * c;
    }
}
 

import org.geotools.referencing.datum.DefaultEllipsoid; //导入依赖的package包/类
public static void main(String[]args){
	Random random = new Random(123);
	ArrayList<LatLong> lls = new ArrayList<>();
	for(int i =0 ; i<100 ; i++){
		LatLongImpl ll = new LatLongImpl(-90 + 180*random.nextDouble(), -180 + 360*random.nextDouble());
		lls.add(ll);
	}
	
	//while(true){
		for(LatLong from:lls){
			for(LatLong to:lls){
				double a = greatCircle(from, to,false);
				
				DefaultEllipsoid elipsoid = DefaultEllipsoid.WGS84;
				double b = elipsoid.orthodromicDistance(from.getLongitude(), from.getLatitude(), to.getLongitude(), to.getLatitude());
				System.out.println("From " + from + " To "+ to + " a=" + a + " b=" + b);
			}
		}			
	//}
}
 

import org.geotools.referencing.datum.DefaultEllipsoid; //导入依赖的package包/类
/**
 * Returns great circle distance in metres
 * @param from
 * @param to
 * @param highAccuracy
 * @return
 */
public static double greatCircle(LatLong from, LatLong to, boolean highAccuracy) {
	double d=0;
	if(highAccuracy){
		DefaultEllipsoid elipsoid = DefaultEllipsoid.WGS84;
		d= elipsoid.orthodromicDistance(from.getLongitude(), from.getLatitude(), to.getLongitude(), to.getLatitude());

	}else{
		d = greatCircleApprox(from, to);
	}
	return d;
}
 

import org.geotools.referencing.datum.DefaultEllipsoid; //导入依赖的package包/类
@Override
public double measure(
		final Coordinate c1,
		final Coordinate c2 ) {
	try {
		return JTS.orthodromicDistance(
				c1,
				c2,
				getCRS());
	}
	catch (final TransformException e) {
		throw new RuntimeException(
				"Failed to transform coordinates to provided CRS",
				e);
	}
	catch (final java.lang.AssertionError ae) {
		// weird error with orthodromic distance..when distance is too close
		// (0.05 meter), it fails the tolerance test
		LOGGER.info(
				"when distance is too close(0.05 meter), it fails the tolerance test",
				ae);

		final GeodeticCalculator calc = new GeodeticCalculator(
				getCRS());
		calc.setStartingGeographicPoint(
				c1.x,
				c1.y);
		calc.setDestinationGeographicPoint(
				c2.x,
				c2.y);
		return ((DefaultEllipsoid) calc.getEllipsoid()).orthodromicDistance(
				calc.getStartingGeographicPoint().getX(),
				calc.getStartingGeographicPoint().getY(),
				calc.getDestinationGeographicPoint().getX(),
				calc.getDestinationGeographicPoint().getY());
	}
}
 

import org.geotools.referencing.datum.DefaultEllipsoid; //导入依赖的package包/类
/**
   * Construct a default geo viewpoint instance
   */
  protected BaseGeoViewpoint() {
      super("GeoViewpoint");

      vfFieldOfView = 0.785398f;
      vfJump = true;
      vfHeadlight = true;
      vfOrientation = new float[] { 0, 0, 1, 0 };
      vfPosition = new double[] { 0, 0, 100000 };
      vfSpeedFactor = 1;
      vfGeoSystem = new String[] {"GD","WE"};
      vfNavType = new String[] { NAV_TYPE_EXAMINE, NAV_TYPE_ANY };
      correctedAvatarSize = new float[] { 0.25f, 1.6f, 0.75f };
      correctedSpeed = 1;
      posVec = new Vector3d();

      qx = new Quat4d();
      qz = new Quat4d();
      qr = new Quat4d();
      tmpAxis = new AxisAngle4d();
      axis = new AxisAngle4f();

      local_quat = new Quat4d();
      rel_quat   = new Quat4d();
      comb_quat  = new Quat4d();

      hasChanged = new boolean[NUM_FIELDS];
      localPosition = new double[3];

      viewpointListeners = new ArrayList<ViewpointListener>(1);
      changeListener = new ArrayList<NavigationInfoChangeListener>();

gt = new GeocentricTransform( DefaultEllipsoid.WGS84, true );
double[] wgs84 = new double[3];
double[] ecef = new double[3];
  }
 

import org.geotools.referencing.datum.DefaultEllipsoid; //导入依赖的package包/类
/**
 * Returns the x value for the google mercator projection (EPSG:3857)
 *
 * @param x the Slippy Map xTile
 * @param z the Slippy Map zoom Level
 * @return a x coordinate in google mercator
 */
public static double tile2xMercator(int x, int z) {
    double n = Math.pow(2.0, z);
    return DefaultEllipsoid.WGS84.getSemiMajorAxis() * Math.PI * (2 * x / n - 1);
}
 

import org.geotools.referencing.datum.DefaultEllipsoid; //导入依赖的package包/类
/**
 * Returns the y value in mercator projection (EPSG:3857)
 *
 * @param y the Slippy Map y tile
 * @param z the Slippy Map zzom Level
 * @return y coordinate in google mercator
 */
public static double tile2yMercator(int y, int z) {
    double n = Math.pow(2.0, z);
    return DefaultEllipsoid.WGS84.getSemiMajorAxis() * Math.PI * (1 - (2*y/n));
}