package commons.gis;

/**
 * 星际旅行概论（坐标系偏移工具类）
 * 支持GC02、BD09和WGS84坐标系之间互相转化
 * @author marquis
 *
 */
public class CrsShiftUtils {
	
	/**
	 * π
	 */
	private final static double PI = 3.1415926535897931;
	
	/**
	 * 长半轴
	 */
	private final static double aa = 6378245.0;
	
	/**
	 * 偏心率平方
	 */
	private final static double ee = 0.00669342162296594323;
	private final static double ee1 = 0.0066934216229659433;
	
	/**
	 * 默认偏移中心点坐标
	 */
	private final static double centerLng = 105.0;
	private final static double centerLat = 35.0;
	private final static double centerLng1 = 100.0;
	private final static double centerLat1 = 30.0;
	
	/** 
	 * 偏移/反偏移方向
	 * @author marquis
	 *
	 */
	public enum Method {
		WGS84_2_GCJ02,
		GCJ02_2_WGS84,
		WGS84_2_BD09,
		BD09_2_WGS84,
		GCJ02_2_BD09,
		BD09_2_GCJ02
	}
	
	/** 
	 * 禁止实例化
	 */
	private CrsShiftUtils() {
	}
	
	/**
	 * 二维点数据的坐标系转换
	 * @param lng
	 * @param lat
	 * @param method
	 * @return
	 */
	public static double[] convert(double lng, double lat, Method method) {
		double ret[];
		switch (method) {
			case WGS84_2_GCJ02:
				ret = shift2Mars(lng, lat, Math.PI, aa, ee, lng, lat);
				ret[0] = lng + ret[0];
				ret[1] = lat + ret[1];
				return ret;
			case GCJ02_2_WGS84:
				ret = shift2Mars(lng, lat, Math.PI, aa, ee, centerLng, centerLat);
				ret[0] = lng - ret[0];
				ret[1] = lat - ret[1];
				return ret;
			case WGS84_2_BD09:
				ret = shift2Mars(lng, lat, Math.PI, aa, ee, centerLng, centerLat);
				ret[0] = lng + ret[0];
				ret[1] = lat + ret[1];
				return mars2Saturn(ret[0], ret[1], Math.PI);
			case BD09_2_WGS84:
				ret = saturn2Mars(lng, lat, Math.PI);
				ret = shift2Mars(ret[0], ret[1], Math.PI, aa, ee, centerLng, centerLat);
				ret[0] = lng - ret[0];
				ret[1] = lat - ret[1];
				return ret;
			case GCJ02_2_BD09:
				return mars2Saturn(lng, lat, Math.PI);
			case BD09_2_GCJ02:
				return saturn2Mars(lng, lat, Math.PI);
			default:
				ret = new double[2];
				ret[0] = lng;
				ret[1] = lat;
				return ret;
		}
	}
	
	/**
	 * 计算火星坐标系的偏移量，正向偏移时坐标加偏移量，反向偏移时减偏移量
	 * @param lng
	 * @param lat
	 * @param pi
	 * @param a
	 * @param e
	 * @param cLng
	 * @param cLat
	 * @return
	 */
	private static double[] shift2Mars(double lng, double lat, double pi, double a, double e, double cLng, double cLat) {
	    if (outOfChina(lng, lat)) {
			double[] ret = new double[2];
			ret[0] = lng;
			ret[1] = lat;
			return ret;
	    }
	    double dlat = transformlat(lng - cLng, lat - cLat, pi);
	    double dlng = transformlng(lng - cLng, lat - cLat, pi);
	    double radlat = lat / 180.0 * pi;
	    double magic = Math.sin(radlat);
	    magic = 1 - e * magic * magic;
	    double sqrtmagic = Math.sqrt(magic);
	    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi);
	    dlng = (dlng * 180.0) / (a / sqrtmagic * Math.cos(radlat) * pi);
	    return new double[]{dlng, dlat};
	}
	
	/**
	 * 计算火星坐标系到土星坐标系
	 * @param lng
	 * @param lat
	 * @param pi
	 * @return
	 */
	private static double[] mars2Saturn(double lng, double lat, double pi) {
		if (outOfChina(lng, lat)) {
			double[] ret = new double[2];
			ret[0] = lng;
			ret[1] = lat;
			return ret;
		}
		double magic = Math.sqrt(lng * lng + lat * lat) + 0.00002 * Math.sin(lat * pi);
		double theta = Math.atan2(lat, lng) + 0.000003 * Math.cos(lng * pi);
		double bdLng = magic * Math.cos(theta) + 0.0065;
		double bdLat = magic * Math.sin(theta) + 0.006;
		return new double[] { bdLng, bdLat };	
	}
	
	/**
	 * 计算土星坐标系到火星坐标系
	 * @param lng
	 * @param lat
	 * @param pi
	 * @return
	 */
	private static double[] saturn2Mars(double lng, double lat, double pi) {
		if (outOfChina(lng, lat)) {
			return new double[] { lng, lat };
		}
		double bdLng = lng - 0.0065, bdLat = lat - 0.006;
		double magic = Math.sqrt(bdLng * bdLng + bdLat * bdLat) - 0.00002 * Math.sin(bdLat * pi);
		double theta = Math.atan2(bdLat, bdLng) - 0.000003 * Math.cos(bdLng * pi);
		double gcjLng = magic * Math.cos(theta);
		double gcjLat = magic * Math.sin(theta);
		return new double[] { gcjLng, gcjLat };
	}

	/**
	 * 转换经度
	 * @param lng
	 * @param lat
	 * @param pi
	 * @return
	 */
	private static double transformlng(double lng, double lat, double pi) {
		double ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * Math.sqrt(Math.abs(lng));
		ret += (20.0 * Math.sin(6.0 * lng * pi) + 20.0 * Math.sin(2.0 * lng * pi)) * 2.0 / 3.0;
		ret += (20.0 * Math.sin(lng * pi) + 40.0 * Math.sin(lng / 3.0 * pi)) * 2.0 / 3.0;
		ret += (150.0 * Math.sin(lng / 12.0 * pi) + 300.0 * Math.sin(lng / 30.0 * pi)) * 2.0 / 3.0;
		return ret;
	}	
	
	/**
	 * 转换维度
	 * @param lng
	 * @param lat
	 * @param pi
	 * @return
	 */
	private static double transformlat(double lng, double lat, double pi) {
		double ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat
				+ 0.2 * Math.sqrt(Math.abs(lng));
		ret += (20.0 * Math.sin(6.0 * lng * pi) + 20.0 * Math.sin(2.0 * lng * pi)) * 2.0 / 3.0;
		ret += (20.0 * Math.sin(lat * pi) + 40.0 * Math.sin(lat / 3.0 * pi)) * 2.0 / 3.0;
		ret += (160.0 * Math.sin(lat / 12.0 * pi) + 320 * Math.sin(lat * pi / 30.0)) * 2.0 / 3.0;
		return ret;
	}
	
	/**
	 * 简单粗暴的判断是否在国内，不在国内不做偏移
	 * @param lng
	 * @param lat
	 * @return
	 */
	private static boolean outOfChina(double lng, double lat) {
		//return false;
		return !(lng > 73.66 && lng < 135.05 && lat > 3.86 && lat < 53.55);
	}

}