地球上两点距离与航向

网上的文章，计算过程写的太过反人类，写个简单的。
就是地球表面的弧，再变成等腰梯形，用等腰梯形近似计算两点长度、角度就OK了

---

贴一段觉得网上写的还算不错的部分
Ea 表示赤道半径，
Eb表示极半径，
Ec的作用就是修正因为纬度不断变化的球半径长度。
地球是一个近似球体，Ea与Eb稍微有点差距。如果在GLAT=0，即在赤道上的时候，Ec=Eb+(Ea-Eb)(90-0)/90=Ea，那Ec就刚好是赤道半径Ea；如果在极点GLAT=90，Ec=Eb+(Ea-Eb)(90-90)/90=Eb，那Ec 就刚好是极半径Eb。
Ed是GLAT所在纬度的纬度圈的半径
Ed = ec * Math.Cos(GLAT * Math.PI / 180)

---

Ed就是一个切面的半径。

double dx=(B.m_RadLo-A.m_RadLo)*A.Ed;
double dy=(B.m_RadLa-A.m_RadLa)*A.Ec;//我感觉这是一个近似计算
这是一种算法

还有用一种计算方式是用梯形上边、下边、高来确定梯形，计算起来会准确些。

附上一段我用的java代码

package com.landsea.webgis.tools;

import java.math.BigDecimal;
import java.text.DecimalFormat;

/**
 * 　　* @Description: GisUtil 经纬度坐标计算
 * 　　* @throws
 * 　　* @author gaow
 * 　　* @date 2018/9/3 13:16
 */
public class GisUtil {

	private static double EARTH_RADIUS = 6378.137;

	static DecimalFormat df = new DecimalFormat("######0.00");

	private static double rad(double d) {
		return d * Math.PI / 180.0;
	}


	public static final double leftTopLongitude = -180;
	public static final double leftTopLatitude = 85;
	public static final double rightBottomLongitude = 180;
	public static final double rightBottomLatitude = -85;

	/**
	 * 获取船速(节)
	 *
	 * @param s 距离（千米）
	 * @param t 时间段(小时)
	 * @return
	 */
	public static Double getShipSog(double s, double t) {
		Double v = s / t;
		Double params = 1.852;//一节=一海里/小时=1.852公里/小时
		return v / params;
	}

	/**
	 * 通过经纬度获取距离(单位：千米)
	 *
	 * @param lat1 第一点的纬度
	 * @param lng1 第一点的精度
	 * @param lat2 第二点的纬度
	 * @param lng2 第二点的精度
	 * @return
	 */
	public static Double getDistance(double lat1, double lng1, double lat2, double lng2) {
		double radLat1 = rad(lat1);
		double radLat2 = rad(lat2);
		double a = radLat1 - radLat2;
		double b = rad(lng1) - rad(lng2);
		double s = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(a / 2), 2)
				+ Math.cos(radLat1) * Math.cos(radLat2)
				* Math.pow(Math.sin(b / 2), 2)));
		s = s * EARTH_RADIUS;
		s = Math.round(s * 10000d) / 10000d;
		//s = s*1000;
		return s;
	}

	/**
	 * 基于余弦定理求两经纬度距离
	 *
	 * @param lat1 第一点的纬度
	 * @param lng1 第一点的精度
	 * @param lat2 第二点的纬度
	 * @param lng2 第二点的精度
	 * @return 返回的距离，单位km
	 */
	public static double lantitudeLongitudeDist(double lat1, double lng1, double lat2, double lng2) {
		double EARTH_RADIUS = 6378.137;
		double dist = EARTH_RADIUS * Math.acos(Math.sin(Math.toRadians(lat1)) * Math.sin(Math.toRadians(lat2)) + Math.cos(Math.toRadians(lat1)) * Math.cos(Math.toRadians(lat2)) * Math.cos(Math.toRadians(lng1 - lng2)));
		String distStr = df.format(dist);
		return Double.parseDouble(distStr);
	}

	/**
	 * 判断是否穿越经度边界
	 *
	 * @param oldLng 旧纬度
	 * @param newLng 新纬度
	 * @return true:穿越；false 未穿越
	 */
	public static Boolean crossEdge(double oldLng, double newLng) {
		return (Math.abs(oldLng - newLng) > 180);
	}

	/**
	 * 　　* @Description: 验证经纬度
	 *
	 * @param lat 纬度
	 * @param lng 经度
	 * @return true:有效；false 无效
	 */
	public static Boolean verifyLatLng(Object lat, Object lng) {
		Boolean res = false;
		if (lat instanceof BigDecimal && lng instanceof BigDecimal) {
			if (lat != null && lng != null) {
				res = verifyLatLng(((BigDecimal) lat).doubleValue(), ((BigDecimal) lng).doubleValue());
			}
		} else if (lat instanceof Double && lng instanceof Double) {
			if (Math.abs((Double) lat) > 0 && Math.abs((Double) lat) < 85 && Math.abs((Double) lng) > 0 && Math.abs((Double) lng) <= 180) {
				res = true;
			}
		}
		return res;
	}

	/**
	 * 获取角度
	 *
	 * @param lat1 第一点的纬度
	 * @param lng1 第一点的精度
	 * @param lat2 第二点的纬度
	 * @param lng2 第二点的精度
	 * @return
	 */
	public static double getAngle(double lat1, double lng1, double lat2, double lng2) {

		double y = Math.sin(lng2 - lng1) * Math.cos(lat2);
		double x = Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1) * Math.cos(lat2) * Math.cos(lng2 - lng1);
		double brng = Math.atan2(y, x);

		brng = Math.toDegrees(brng);
		if (brng < 0)
			brng = brng + 360;
		return brng;

	}
	/**
	 * 获取AB连线与正北方向的角度
	 * @param A  A点的经纬度
	 * @param B  B点的经纬度
	 * @return  AB连线与正北方向的角度（0~360）
	 */
	public  static double getAngle(MyLatLng A,MyLatLng B){
		double dx=(B.m_RadLo-A.m_RadLo)*A.Ed;
		double dy=(B.m_RadLa-A.m_RadLa)*A.Ec;
		double angle=0.0;
		if(dy!=0)
		{
            angle=Math.atan(Math.abs(dx/dy))*180./Math.PI;
        }
		double dLo=B.m_Longitude-A.m_Longitude;
		double dLa=B.m_Latitude-A.m_Latitude;
		if(dLo>0&&dLa<=0){
			angle=(90.-angle)+90;
		}
		else if(dLo<=0&&dLa<0){
			angle=angle+180.;
		}else if(dLo<0&&dLa>=0){
			angle= (90.-angle)+270;
		}
		return angle;
	}
	static class MyLatLng {
		final static double Rc=6378137;
		final static double Rj=6356725;
		double m_LoDeg,m_LoMin,m_LoSec;
		double m_LaDeg,m_LaMin,m_LaSec;
		double m_Longitude,m_Latitude;
		double m_RadLo,m_RadLa;
		double Ec;
		double Ed;
		public MyLatLng(double latitude,double longitude){
			m_LoDeg=(int)longitude;
			m_LoMin=(int)((longitude-m_LoDeg)*60);
			m_LoSec=(longitude-m_LoDeg-m_LoMin/60.)*3600;

			m_LaDeg=(int)latitude;
			m_LaMin=(int)((latitude-m_LaDeg)*60);
			m_LaSec=(latitude-m_LaDeg-m_LaMin/60.)*3600;

			m_Longitude=longitude;
			m_Latitude=latitude;
			m_RadLo=longitude*Math.PI/180.;
			m_RadLa=latitude*Math.PI/180.;
//			Ec=Rj+(Rc-Rj)*(90.-m_Latitude)/90.;
			Ec=Rc;
			Ed=Rc*Math.cos(m_RadLa);
//			Ed = Rj;
		}
	}

	// Test
	public static void main(String[] args) {
		/*double lat1 = 38.9158809764;
		double lng1 = 121.6274209139;
		double lat2 = 43.8282528952;
		double lng2 = 126.5457496619;

		double lat1 = 16.8046207000;
		double lng1 = 117.0072556000;
		double lat2 = 16.8153400000;
		double lng2 = 117.0159225000;
		System.out.println(getAngle(lat1, lng1, lat2, lng2));


		System.out.println(getDistance(lat1, lng1, lat2, lng2));
		System.out.println(lantitudeLongitudeDist(lat1, lng1, lat2, lng2));*/
		
		double lat1 = -31.2768879000;
		double lng1 = 50.3630295000;
		double lat2 = -31.2733192000;
		double lng2= 50.3776588000;


		MyLatLng A=new MyLatLng(lat1,lng1);
		MyLatLng B=new MyLatLng(lat2,lng2);
		System.out.println(getAngle(A,B));
	}
	public static double getAngleNew(double lat1, double lng1, double lat2, double lng2) {
		MyLatLng A=new MyLatLng(lat1,lng1);
		MyLatLng B=new MyLatLng(lat2,lng2);
		return  getAngle(A,B);
	}
}