package com.base.core.operator;

import com.base.api.model.SubPointDataRpcDTO;
import com.base.core.model.enity.*;
import com.base.core.source.LatLon;
import com.base.core.utiles.*;
import com.base.orbit.utils.DateUtils;
import org.springframework.stereotype.Component;

import java.util.ArrayList;
import java.util.Calendar;
import java.util.List;

/**
 * 条带操作常用函数实现
 * @author vrss2
 *
 */
@Component
public class StripRollCalc {
	
	/**
	 * 地球半径,值为6378.137,单位为km
	 */
	private static double	R					= 6378.137;						// 6378.137;
	
	/**
	 * 计算卫星物理上向观测区域西北点旋转对应的侧摆角是多少
	 *	通过方位角判断角度正负值
	 *	如果方位角为正,表示物理上向右旋转，如果卫星为降轨，则表示向左侧摆
	 *方位角为负,表示物理上向左旋转，如果卫星为降轨，则表示向右侧摆
	 * @param region
	 * 观测区域
	 * @param ephtrack
	 * 星下点轨迹
	 * @return
	 */
	public double calTurnedNWdegree(Region region, EphTrackCoreDTO ephtrack, double maxrollangle){
		double lat1, lat2, lon1, lon2, altitude;
		lat1 = region.getNwNode().getLat();
		lon1 = region.getNwNode().getLon();
		double dis = 1000;
		int count = 0;
		for (int i = 0; i < ephtrack.getSubPointlist().size(); i++) {
			double temp = Math.abs(ephtrack.getSubPointlist().get(i).getRsLat() - lat1);
			if (temp < dis) {
				dis = temp;
				count = i;
			}
		}

		lat2 = ephtrack.getSubPointlist().get(count).getRsLat();
		lon2 = ephtrack.getSubPointlist().get(count).getRsLon();
		altitude = ephtrack.getSubPointlist().get(count).getAlt();
		double distance = CalDistance.Distance(lat1, lon1, lat2, lon2);

		// 由最大侧摆转动的距离
		double turneddegree = CalDistance.CalRollAngle(distance, altitude);
		if(turneddegree == 0){
			return turneddegree;
		}
		LatLon l1 = LatLon.fromDegrees(lat1, lon1);
		LatLon l2 = LatLon.fromDegrees(lat2, lon2);
		//计算方位角
		double azimuth = LatLon.greatCircleAzimuth(l2, l1).getDegrees();
		//对侧摆值进行限制,如超过最大侧摆值限制,则以最大侧摆值限制为值
		double rollangle = 0;
		//如果卫星为下降轨，方位角为正表示物理向右侧摆，真实卫星向左侧摆
		if(ephtrack.getUpdown().equalsIgnoreCase("1")){
			if(azimuth > 0){
				rollangle = - turneddegree;
			}else{
				rollangle = turneddegree;
			}
		}else{
			//如果卫星为上升轨，方位角为正表示物理向右侧摆，真实卫星向右侧摆
			if(azimuth > 0){
				rollangle = turneddegree;
			}else{
				rollangle = -turneddegree;
			}
		}
		if(rollangle > 0 && rollangle > maxrollangle){
			rollangle = maxrollangle;
		}
		if(rollangle < 0 && rollangle < - maxrollangle){
			rollangle = - maxrollangle;
		}
		return rollangle;
		
	}

	/**
	 * 计算卫星物理上向观测区域东南点旋转对应的侧摆角是多少
	 *	通过方位角判断角度正负值
	 *	如果方位角为正,表示物理上向右旋转，如果卫星为降轨，则表示向左侧摆
	 *方位角为负,表示物理上向左旋转，如果卫星为降轨，则表示向右侧摆
	 * @param region
	 * 观测区域
	 * @param ephtrack
	 * 星下点轨迹
	 * @return
	 */
	
	public double calTurnedSEdegree(Region region, EphTrackCoreDTO ephtrack, double maxrollangle) {
		double lat1, lat2, lon1, lon2, altitude;
		lat1 = region.getSeNode().getLat();
		lon1 = region.getSeNode().getLon();
		double dis = 1000;
		int count = 0;
		for (int i = 0; i < ephtrack.getSubPointlist().size(); i++) {
			double temp = Math.abs(ephtrack.getSubPointlist().get(i).getLsLat() - lat1);
			if (temp < dis) {
				dis = temp;
				count = i;
			}
		}

		lat2 = ephtrack.getSubPointlist().get(count).getLsLat();
		lon2 = ephtrack.getSubPointlist().get(count).getLsLon();
		altitude = ephtrack.getSubPointlist().get(count).getAlt();
		double distance = CalDistance.Distance(lat1, lon1, lat2, lon2);

		// 由最大侧摆转动的距离
		double turneddegree = CalDistance.CalRollAngle(distance, altitude);
		if(turneddegree == 0){
			return turneddegree;
		}
		LatLon l1 = LatLon.fromDegrees(lat1, lon1);
		LatLon l2 = LatLon.fromDegrees(lat2, lon2);
		double azimuth = LatLon.greatCircleAzimuth(l2, l1).getDegrees();
		// 对侧摆值进行限制,如超过最大侧摆值限制,则以最大侧摆值限制为值
		double rollangle = 0;
		// 如果卫星为下降轨，方位角为正表示物理向右侧摆，真实卫星向左侧摆
		if (ephtrack.getUpdown().equalsIgnoreCase("1")) {
			if (azimuth > 0) {
				rollangle = -turneddegree;
			} else {
				rollangle = turneddegree;
			}
		} else {
			// 如果卫星为上升轨，方位角为正表示物理向右侧摆，真实卫星向右侧摆
			if (azimuth > 0) {
				rollangle = turneddegree;
			} else {
				rollangle = -turneddegree;
			}
		}
		if (rollangle > 0 && rollangle > maxrollangle) {
			rollangle = maxrollangle;
		}
		if (rollangle < 0 && rollangle < - maxrollangle) {
			rollangle = - maxrollangle;
		}
		return rollangle;
	}
	
	/**
	 * 计算左右最大侧摆角值
	 * @param ephrelation
	 * @param region
	 * @param ephtrack
	 * @param maxrollangle
	 * @return
	 */
	public double[] calrollangle(String ephrelation, Region region, EphTrackCoreDTO ephtrack, double maxrollangle) throws Exception{
		double[] rollangle = {0,0};
		double leftrollangle = - maxrollangle;
		double rightrollangle = maxrollangle;
		double taskrollreq = 0;
		if(region.getCallRequireVisibility()==null){
			taskrollreq = Math.abs(maxrollangle);
		}else{
			taskrollreq = region.getCallRequireVisibility().getRollAngle().doubleValue();
		}		
		if(Math.abs(taskrollreq) > Math.abs(maxrollangle)){
			taskrollreq = Math.abs(maxrollangle);
		}
		leftrollangle = - taskrollreq;
		rightrollangle = taskrollreq;
/*		Log.getLogger().debug("debug", "leftrollangle:" + leftrollangle,region.getBizRequire().getCode());
		Log.getLogger().debug("debug", "rightrollangle:" + rightrollangle,region.getBizRequire().getCode());*/
		rollangle[0] = leftrollangle;
		rollangle[1] = rightrollangle;
		return rollangle;
	}
	
	

	/**
	 * 按角度Degree_piece平均拆分
	 * @param leftrollangle
	 * 最小侧摆角值（左侧摆角值）
	 * @param rightrollangle
	 * 最大侧摆角值（右侧摆角值）
	 * @param roll_lmt
	 * 侧摆角限制值
	 * @return
	 */
	public List<Double> calRoll(double leftrollangle, double rightrollangle, double roll_lmt,double degree_split) {
		if (leftrollangle > 0 && leftrollangle > roll_lmt) {
			leftrollangle = roll_lmt;
		}
		if (rightrollangle > 0 && rightrollangle > roll_lmt) {
			rightrollangle = roll_lmt;
		}
		if (leftrollangle < 0 && Math.abs(leftrollangle) > roll_lmt) {
			leftrollangle = -roll_lmt;
		}
		if (rightrollangle < 0 && Math.abs(rightrollangle) > roll_lmt) {
			rightrollangle = -roll_lmt;
		}
		double minroll = 0, maxroll = 0, rollangle = 0;
		if (leftrollangle > rightrollangle) {
			minroll = rightrollangle;
			maxroll = leftrollangle;

		} else {
			minroll = leftrollangle;
			maxroll = rightrollangle;
		}

		List<Double> roll_list = new ArrayList<Double>();
		if (maxroll - minroll != 0) {
			rollangle = minroll;
			while (rollangle < maxroll) {
				roll_list.add(rollangle);
				rollangle = rollangle + degree_split;
			}
			roll_list.add(maxroll);
		} else {
			roll_list.add(maxroll);
		}
		return roll_list;
	}
	
	/**
	 * 拆分侧摆角，按分数Piece平均拆分Piece
	 * @param leftrollangle
	 * 最小侧摆角值（左侧摆角值）
	 * @param rightrollangle
	 * 最大侧摆角值（右侧摆角值）
	 * @param roll_lmt
	 * 侧摆角限制值
	 * @param piece
	 * 拆分份数
	 * @return
	 */
	public List<Double> calRollbyPiece(double leftrollangle, double rightrollangle, double roll_lmt,int piece) {
		if (leftrollangle > 0 && leftrollangle > roll_lmt) {
			leftrollangle = roll_lmt;
		}
		if (rightrollangle > 0 && rightrollangle > roll_lmt) {
			rightrollangle = roll_lmt;
		}
		if (leftrollangle < 0 && Math.abs(leftrollangle) > roll_lmt) {
			leftrollangle = -roll_lmt;
		}
		if (rightrollangle < 0 && Math.abs(rightrollangle) > roll_lmt) {
			rightrollangle = -roll_lmt;
		}
		double minroll = 0, maxroll = 0, rollangle = 0;
		if (leftrollangle > rightrollangle) {
			minroll = rightrollangle;
			maxroll = leftrollangle;

		} else {
			minroll = leftrollangle;
			maxroll = rightrollangle;
		}

		List<Double> roll_list = new ArrayList<Double>();
		if (maxroll - minroll != 0) {
			int count = 0;
			rollangle = minroll;
			double disroll = (maxroll - minroll) / (piece - 1);
			while (count != piece) {
				roll_list.add(rollangle);
				rollangle = rollangle + disroll;
				count++;
			}
		} else {
			roll_list.add(leftrollangle);
		}
		return roll_list;
	}
		
	public SubPoint getNearPoint(NodeCoreDTO center, EphTrackCoreDTO ephtrack){
		
		double lat1 = center.getLat();
		double lon1 = center.getLon();
		SubPoint nearpoint = new SubPoint();
//		double tempAngle = 0;
		double tempAngle=90.0;
		for(int i = 0;i < ephtrack.getSubPointlist().size();i++){
			//计算公式
			// 计算cos(alpha)值
//			double value = Math.sin(ephtrack.getSubPointlist().get(i).getLat() * Math.PI / 180) * Math.sin(lat1 * Math.PI / 180)
//					+ Math.cos(ephtrack.getSubPointlist().get(i).getLat() * Math.PI / 180) * Math.cos(lat1 * Math.PI / 180)
//					* (Math.cos((ephtrack.getSubPointlist().get(i).getLon() - lon1) * Math.PI / 180));
			double value=SphericalAngle(ephtrack.getSubPointlist().get(i), center);
			
			if (value >= 90) {
				System.out.println(value);
			}
			
//			if (value > tempAngle) {
			if (value < tempAngle) {
				tempAngle = value;
				nearpoint = ephtrack.getSubPointlist().get(i);	
				nearpoint.setNearPointNum(i);
			} 
			
		}
		
		return nearpoint;
	}
	private double SphericalAngle(SubPoint subPointIn, NodeCoreDTO center)
	{
		double pi=3.14159265358979323846;
		double subPointLon=subPointIn.getLon()*pi/180.0;
		double subPointLat=subPointIn.getLat()*pi/180.0;
		double targetLon=center.getLon()*pi/180.0;
		double targetLat=center.getLat()*pi/180.0;
		
		double back=Math.cos(pi/2-subPointLat)*Math.cos(pi/2-targetLat)+
				Math.sin(pi/2-subPointLat)*Math.sin(pi/2-targetLat)*Math.cos(Math.abs(subPointLon-targetLon));
		return Math.acos(back)*180.0/pi;
	}
	/**
	 * 	计算点目标的侧摆角值
	 * @param ephtrack
	 * @return
	 */
	public double calPointRollangle(SubPoint nearpoint,NodeCoreDTO center, EphTrackCoreDTO ephtrack,String payloadType){
		double lat1 = center.getLat();
		double lon1 = center.getLon();

		double rollangle=90;
		if(nearpoint.getLat() == lat1 && nearpoint.getLon() == lon1){
			rollangle = 0 ;
		}else{
//			//计算公式
//			double cosAlpha = Math.sin(nearpoint.getLat() * Math.PI / 180) * Math.sin(lat1 * Math.PI / 180)
//					+ Math.cos(nearpoint.getLat() * Math.PI / 180) * Math.cos(lat1 * Math.PI / 180)
//					* (Math.cos((nearpoint.getLon() - lon1) * Math.PI / 180));
//			double sinAlpha = Math.sqrt(1 - Math.pow(cosAlpha, 2));		
//			double H =Double.valueOf(nearpoint.getHigh()/1000.0);
//			double sinRow = R / (R + H);
//			double turneddegree = Math.atan(sinRow * sinAlpha / (1 - sinRow * cosAlpha)) * (180 / Math.PI);
				
			//计算侧摆角
			double turneddegree=accessLLA(nearpoint,center);
			
//			LatLon l1 = LatLon.fromDegrees(lat1, lon1);
//			LatLon l2 = LatLon.fromDegrees(nearpoint.getLat(), nearpoint.getLon());
//			//计算方位角
//			double azimuth = LatLon.greatCircleAzimuth(l2, l1).getDegrees();
//			
//			//如果卫星为下降轨,方位角为正表示物理向右侧摆,真实卫星向左侧摆
//			if(ephtrack.getUpdown().equalsIgnoreCase("1")){
//				if(azimuth < 0){
//					rollangle =  -turneddegree;
//				}else{
//					rollangle =  turneddegree;
//				}
//			}else{
//				//如果卫星为上升轨,方位角为正表示物理向右侧摆,真实卫星向右侧摆
//				if(azimuth < 0){
//					rollangle = turneddegree;
//				}else{
//					rollangle = -turneddegree;
//				}
//			}
			
			
			int nearnum = nearpoint.getNearPointNum();
			
			if(nearnum==0)//第一个点最近  高越改 
			{
				if(nearnum == (ephtrack.getSubPointlist().size()-1))
				{
					//只有一个点
				}
				else
				{
					SubPoint nearpointafter=ephtrack.getSubPointlist().get(nearnum+1);
					
					Node CBpdpos=new Node();
				
					//侧摆角正负精确判断，适用于太阳同步轨道，侧摆角正负判断最好还是按照QWERT14的方法用矢量判断
					CBpdpos.setLon(center.getLon());
					
					//轨道跨正负180度经线情况,把-180至180经度转为0-360经度，用于判断侧摆正负		
					if(Math.abs(CBpdpos.getLon())>=155.0||Math.abs(nearpoint.getLon())>=155.0||Math.abs(nearpointafter.getLon())>=155.0)
					{
						if(nearpointafter.getLon()<0.0)
						{
							nearpointafter.setLon(nearpointafter.getLon()+360.0);
						}
						if(nearpoint.getLon()<0.0)
						{
							nearpoint.setLon(nearpoint.getLon()+360.0);
						}
						if(CBpdpos.getLon()<0.0)
						{
							CBpdpos.setLon(CBpdpos.getLon()+360.0);
						}				
					}
					
					double Const=(CBpdpos.getLon()-nearpoint.getLon())/(nearpoint.getLon()-nearpointafter.getLon());
					CBpdpos.setLat(nearpoint.getLat()+Const*(nearpoint.getLat()-nearpointafter.getLat()));
					if(center.getLat()>=CBpdpos.getLat())
					{
						
						// 载荷特性 0光学/1红外/2电子/3雷达
						//SAR卫星				
						if (payloadType.contains("3"))// 3雷达
						{
							rollangle = turneddegree;
						}
						//光学卫星
						else
						{
							rollangle = -turneddegree;
						}
					}
					else
					{
						//SAR卫星
						if (payloadType.contains("3"))// 3雷达
						{
							rollangle = -turneddegree;
						}
						//光学卫星
						else
						{
							rollangle = turneddegree;
						}
					}
				}
			}
			else//不是第一个点最近，且星下点大于1
			{
				SubPoint nearpointbefore=ephtrack.getSubPointlist().get(nearnum-1);
				
				Node CBpdpos=new Node();
			
				//侧摆角正负精确判断，适用于太阳同步轨道，侧摆角正负判断最好还是按照QWERT14的方法用矢量判断
				CBpdpos.setLon(center.getLon());
				
				//轨道跨正负180度经线情况,把-180至180经度转为0-360经度，用于判断侧摆正负		
				if(Math.abs(CBpdpos.getLon())>=155.0||Math.abs(nearpoint.getLon())>=155.0||Math.abs(nearpointbefore.getLon())>=155.0)
				{
					if(nearpointbefore.getLon()<0.0)
					{
						nearpointbefore.setLon(nearpointbefore.getLon()+360.0);
					}
					if(nearpoint.getLon()<0.0)
					{
						nearpoint.setLon(nearpoint.getLon()+360.0);
					}
					if(CBpdpos.getLon()<0.0)
					{
						CBpdpos.setLon(CBpdpos.getLon()+360.0);
					}				
				}
				
				double Const=(CBpdpos.getLon()-nearpoint.getLon())/(nearpointbefore.getLon()-nearpoint.getLon());
				CBpdpos.setLat(nearpoint.getLat()+Const*(nearpointbefore.getLat()-nearpoint.getLat()));
				if(center.getLat()>=CBpdpos.getLat())
				{
					
					// 载荷特性 0光学/1红外/2电子/3雷达
					//SAR卫星				
					if (payloadType.contains("3"))// 3雷达
					{
						rollangle = turneddegree;
					}
					//光学卫星
					else
					{
						rollangle = -turneddegree;
					}
				}
				else
				{
					//SAR卫星
					if (payloadType.contains("3"))// 3雷达
					{
						rollangle = -turneddegree;
					}
					//光学卫星
					else
					{
						rollangle = turneddegree;
					}
				}
			}
			
			
			
		}
		
		return NumberUtils.getRoundHalfUp(rollangle,4);//4舍5入,保留4位小数
	}
	
	
		//函数作用见本上说明，参数sat：卫星空间位置测地坐标，target：目标点测地坐标
		private double accessLLA(SubPoint sat, NodeCoreDTO target) {

			//sat：卫星空间位置测地坐标，target：目标点测地坐标
			//Sa、S：星下点地心坐标，T：目标点地心坐标
			double pi=3.14159265358979323846d;
			SubPointDataRpcDTO Sa;
			SubPoint S = new SubPoint();
			NodeCoreDTO T = new NodeCoreDTO();
			
			//测地坐标转换为地心坐标
			Sa=LLAToLLR(sat); 
			T.setLat(GeodesicToCore(target.getLat()));
			T.setLon(target.getLon());
			S.setLat(Sa.getLat());
			S.setLon(Sa.getLon());
			//调用函数SphericalAngle求星下点与目标球面距离coreangle：弧度	  
			double coreangle=SphericalAngle(S,T);
			//OT：目标到地心距离，OS：星下点到地心距离
			double OT=EarthRadiusA(T.getLat());
			double OS=Sa.getR();
			//	  
			double ST=Math.sqrt( Math.pow(OT,2)+Math.pow(OS,2)-2*OT*OS*Math.cos(coreangle*pi/180.0) ) ;
			return Math.asin( OT/ST*Math.sin(coreangle*pi/180.0) )*180.0/pi;
		}
		
		//函数用测地纬度求地表上点的地心纬度，返回值地心纬度，参数lat测地纬度，见《球面天文学》
		private double GeodesicToCore(double lat) {
			//a为地球赤道半径，f为地球扁率
			double a=6378137,f=1.0/298.257;
			double pi=3.14159265358979323846d;
			double C,S;
		    lat=lat*pi/180.0;
			C=a/( Math.sqrt(Math.pow(Math.cos(lat),2)+(1-f)*(1-f)*Math.pow(Math.sin(lat),2)));
			S=C*(1-f)*(1-f);
			//lat为测地纬度
		    return Math.atan(S/C*Math.tan(lat))*180.0/pi; 
		}

		//?公式未知，此函数用地心纬度求地表不同纬度圈地心距，参数lat：地心纬度，函数返回值地心距
		private double EarthRadiusA(double lat) {
			double a=6378137,f=1.0/298.257,b;
			double pi=3.14159265358979323846d;
			double back;
			lat=lat*pi/180.0;
			b=a*(1-f);
		    back=(1+Math.pow(Math.tan(lat),2))*a*a*b*b/(b*b+a*a*Math.pow(Math.tan(lat),2));
			return Math.sqrt(back);
		}

		//此函数用于测地坐标向地心坐标的转换，参数sat为测地坐标，见《球面天文学》37页图2.7
		private SubPointDataRpcDTO LLAToLLR(SubPoint sat) {
			double a=6378137,f=1.0/298.257,b;
			double pi=3.14159265358979323846d;
			double Radius,CoreAngle;
			SubPointDataRpcDTO back = new SubPointDataRpcDTO();
			b=a*(1-f);
		    Radius=EarthRadiusB(sat.getLat()); 
		    CoreAngle=Math.atan( Math.tan(sat.getLat()*pi/180.0)*b*b/a/a ) ;
			double X=Radius*Math.cos(CoreAngle);
			double Y=Radius*Math.sin(CoreAngle);
			//ST即图中NM，OS即图中ON
		    double ST=Y/Math.sin(sat.getLat()*pi/180.0)+sat.getAlt();
			double OS=X-Y/Math.tan(sat.getLat()*pi/180.0);			
			//R为卫星到地心距离
			back.setR(Math.sqrt( OS*OS+ST*ST-2*OS*ST*Math.cos((180.0-sat.getLat())/180.0*pi) )); 
			//back.lon为星下点地心经度
		    back.setLon(sat.getLon());
			//back.lat为星下点地心纬度
			back.setLat(Math.asin ( ST/back.getR()*Math.sin(sat.getLat()*pi/180.0) )*180.0/pi);
		    return back;
		}

		//此函数用测地纬度求地表不同纬度圈地心距，参数lat：测地纬度，函数返回值地心距
		private double EarthRadiusB(double lat) {
			//a为地球赤道半径，f为地球扁率
			double a=6378137,f=1.0/298.257;
			// double pi=3.14159265358979323846;
			double C,S;
		    lat=lat*Math.PI/180.0;
			C=a/( Math.sqrt(Math.pow(Math.cos(lat),2)+(1-f)*(1-f)*Math.pow(Math.sin(lat),2)));
			S=C*(1-f)*(1-f);
		    return Math.sqrt( Math.pow(C*Math.cos(lat),2)+Math.pow(S*Math.sin(lat),2) ); 
		}

		
		/**
		 * 用于求各时刻太阳高度角
		 * 
		 * @param sunTime
		 * @param pSat
		 * @return
		 */
		public double sunAngle(Calendar sunTime, NodeCoreDTO pSat) {
			SubPoint sunPos = new SubPoint();
			// double pi=3.14159265358979323846;
			Calendar cstTime = TimeUtils.getCalendar(2000, 1, 1, 20, 0, 0);
			long count = (sunTime.getTimeInMillis() - cstTime.getTimeInMillis()) / 1000;
			// T: 从J2000.0起算的儒略世纪数
			double t = (double)count / 3600.0 / 24.0 / 36525.0;
			// L: 纽康的太阳几何平黄线（以消除光行差影响），见《球面天文学》 20页
			double l = 279.0+41.0/60+24.68/3600+129602768.13/3600*(1+t)+1.089*(1+t)*(1+t)/3600;
			while (l > 360) l = l - 360;
			while (l < 0) l = l + 360;
			// 平黄赤交角
			double hcAngle = 84381.4480/3600.0-46.81500/3600.0*t-0.00059/3600.0*Math.pow(t,2)+0.0018130/3600.0*Math.pow(t,3);
			//通过太阳黄经L、黄赤交角hcangle求太阳赤经lon
		    if(l<=90)           { sunPos.setLon(Math.atan(Math.tan(l*Math.PI/180)*Math.cos(hcAngle*Math.PI/180))*180/Math.PI);}
			if(l>90 && l<=180)  { sunPos.setLon(180 +Math.atan(Math.tan(l*Math.PI/180)*Math.cos(hcAngle*Math.PI/180))*180/Math.PI);}
			if(l>180 && l<=270) { sunPos.setLon(180 +Math.atan(Math.tan(l*Math.PI/180)*Math.cos(hcAngle*Math.PI/180))*180/Math.PI);}
			if(l>270 && l<=360) { sunPos.setLon(360 +Math.atan(Math.tan(l*Math.PI/180)*Math.cos(hcAngle*Math.PI/180))*180/Math.PI);}
			//通过太阳黄经L、黄赤交角hcangle求太阳赤纬llat
			sunPos.setLat(Math.asin(Math.sin(l*Math.PI/180)*Math.sin(hcAngle*Math.PI/180))*180/Math.PI);
			//赤经从春分点起算，所以减去格林尼治恒星时转化为太阳经度
		    sunPos.setLon(sunPos.getLon()-SgTime(sunTime, 0, 0));
			while(sunPos.getLon()<0) { sunPos.setLon(sunPos.getLon()+360);}
			while(sunPos.getLon()>360) { sunPos.setLon(sunPos.getLon()-360);}
			//计算太阳与目标的球面距离：弧度
			double angle=SphericalAngle(sunPos,pSat);
		    if(angle>90) { return -90;}
			//返回值：太阳高度角
		    else if (angle > 0) {return 90-angle;}
		    else {return angle;}
		}
		
		private double SgTime(Calendar sunTime, double lonNutation, double delta) {
			//  J2000.0标准历元的起始时间
			Calendar startTime = TimeUtils.getCalendar(2000, 1, 1, 20, 0, 0);
			long span = (sunTime.getTimeInMillis() - startTime.getTimeInMillis()) / 1000;
			//  GetTotalSecond返回总秒数
			//? 此公式有问题	
			//  double sg=(double(18.69737460) + (int_time+delta_t) * double(1.002737909350795) /3600.0) *15.0;
		    double tt=(double)span/(24.0*3600.0)/36525;
			//格林尼治平恒星时
			double sg=(18.6973746+tt*879000.0513367+0.093104/3600.0*tt*tt-6.2/1000000.0/3600.0*tt*tt*tt)*15.0;
			//加上赤经章动的修正，格林尼治真恒星时
			sg=sg+lonNutation;
			while (sg>=360)
			{
				sg=sg-360;
			}
			/*
		    while(sg>=360.0)
			{
			sg=fmod(sg,360.0);
			}
			*/
			return sg;
		}

		
		
		/**
		 * 計算太陽高度角
		 * @param sublist
		 * @return
		 */
		public double calcSunAngle(List<SubPoint> sublist){
			double sunangle = 0;
			int center_count = sublist.size() / 2;
			String centerTime = DateUtils.DateToString(sublist.get(center_count).getDate());
			double centerLongitude = sublist.get(center_count).getLon();
			double centerLatitude = sublist.get(center_count).getLat();
			sunangle = SunElevationAngleCalculator.calculateSunElevationAngle
					(centerTime, centerLongitude, centerLatitude);
			return sunangle;
		}
	
		/**
		 * 計算目标点太陽高度角
		 * @return
		 */
		public double calcSunAngle(String nearPointTime,double Longitude,double Latitude){
//			double sunangle = 0;
//			int center_count = sublist.size() / 2;
//			String centerTime = sublist.get(center_count).getTime();			
//			sunangle = SunElevationAngleCalculator.calculateSunElevationAngle
//					(centerTime, Longitude, Latitude);
			
			TimePattern dbTimePattern=CalcPara.INNER_TIME_PATTERN;//数据库的时区，北京时
			TimePattern innerTimePattern=CalcPara.DB_TIME_PATTERN_BJ;//轨道计算算法的时区，UTC	
			//TimePattern dbTimePattern=CalcPara.INNER_TIME_PATTERN;//UTC, gaoyue edit 20190725
			//TimePattern innerTimePattern=CalcPara.INNER_TIME_PATTERN;//UTC, gaoyue edit 20190725
								
			String sunAngleTime=DateFormatter.srcTimePattern2DstTimePattern(nearPointTime,dbTimePattern,innerTimePattern);
			
			Calendar sunTime=TimeUtils.try2ParseTime(sunAngleTime);
			NodeCoreDTO pSat=new NodeCoreDTO();
			pSat.setLat(Latitude);
			pSat.setLon(Longitude);
			double sunangle=sunAngle(sunTime, pSat);
			
			return NumberUtils.getRoundHalfUp(sunangle,4);//4舍5入，保留4位小数;
		}
		
		
		public double calcSunAngleNew(){
			double sunangle = 0;
			
			return sunangle;
		}
}
