public double calculateDistance() throws IllegalStateException,
			ArithmeticException {
		double U1, U2, w, sigma, ssig, csig, salp, c2alp, c2sigm;
		double u, A, B, C, dsig, b;
		double lambda, dlam = 0.0, lastdlam = 0.0;
		double p, q;
		int nIters = 0, nMaxIters = 10;

		double a = ellipsoid.getSemiMajorAxis();
		double f = ellipsoid.getFlattening();

		if (!havePositionA)
			throw new IllegalStateException("positionA is not defined");

		if (!havePositionB)
			throw new IllegalStateException("positionB is not defined");

		b = a * (1.0 - f);

		U1 = Math.atan((1.0 - f) * Math.tan(latitudeA));
		U2 = Math.atan((1.0 - f) * Math.tan(latitudeB));

		w = longitudeB - longitudeA;
		if (w < -Math.PI)
			w += TwoPi;
		if (w > Math.PI)
			w -= TwoPi;

		lambda = w;

		do {
			lastdlam = dlam;

			p = Math.cos(U2) * Math.sin(lambda);
			q = Math.cos(U1) * Math.sin(U2) - Math.sin(U1) * Math.cos(U2)
					* Math.cos(lambda);

			ssig = Math.sqrt(p * p + q * q);

			csig = Math.sin(U1) * Math.sin(U2) + Math.cos(U1) * Math.cos(U2)
					* Math.cos(lambda);

			sigma = Math.atan2(ssig, csig);

			salp = Math.cos(U1) * Math.cos(U2) * Math.sin(lambda) / ssig;

			c2alp = 1.0 - salp * salp;

			c2sigm = csig - 2.0 * Math.sin(U1) * Math.sin(U2) / c2alp;

			C = (f / 16.0) * c2alp * (4.0 + f * (4.0 - 3.0 * c2alp));

			dlam = (1.0 - C)
					* f
					* salp
					* (sigma + C
							* ssig
							* (c2sigm + C * csig
									* (2.0 * c2sigm * c2sigm - 1.0)));

			lambda = w + dlam;

			nIters += 1;
		} while ((Math.abs(dlam - lastdlam) > EPSILON) && (nIters <= nMaxIters));

		if (nIters > nMaxIters)
			throw new ArithmeticException("Exceeded iteration limit");

		u = (1.0 - salp * salp) * (a * a - b * b) / (b * b);

		A = 1.0 + (u / 16384.0)
				* (4096.0 + u * (-768.0 + u * (320.0 - 175.0 * u)));

		B = (u / 1024.0) * (256.0 + u * (-128.0 + u * (74.0 - 47.0 * u)));

		dsig = B
				* ssig
				* (c2sigm + (B / 4.0)
						* (csig * (2.0 * c2sigm * c2sigm - 1.0) - (B / 6.0)
								* c2sigm * (4.0 * ssig * ssig - 3.0)
								* (4.0 * c2sigm * c2sigm - 3.0)));

		p = Math.cos(U2) * Math.sin(lambda);
		q = Math.cos(U1) * Math.sin(U2) - Math.sin(U1) * Math.cos(U2)
				* Math.cos(lambda);

		azimuthA = Math.atan2(p, q);
		haveAzimuthA = true;

		p = -Math.cos(U1) * Math.sin(lambda);
		q = Math.cos(U2) * Math.sin(U1) - Math.sin(U2) * Math.cos(U1)
				* Math.cos(lambda);

		azimuthB = Math.atan2(p, q);
		haveAzimuthB = true;

		geodesicDistance = b * A * (sigma - dsig);
		haveDistance = true;

		return geodesicDistance;
	}