package com.base.core.operator;

import com.alibaba.fastjson.JSON;
import com.base.core.BizService;
import com.base.core.config.ImageLimitation;
import com.base.core.model.enity.*;
import com.base.core.source.SideAngleSurfaceGrid;
import com.base.core.utiles.SunElevationAngleCalculator;
import com.base.core.utiles.TimeFormatTool;
import com.base.core.utiles.TimeUtils;
import lombok.extern.slf4j.Slf4j;
import org.springframework.beans.factory.annotation.Value;
import org.springframework.stereotype.Service;
import org.springframework.util.CollectionUtils;

import javax.annotation.Resource;
import java.math.BigDecimal;
import java.util.*;

@Service
@Slf4j
public class AccessWindowService {

    @Resource
    private ImageLimitation rLimit;

    @Value("${maxWindowLimit:3600}")
    private Integer maxWindowLimit;

    @Resource
    BizService bizService;

    List<ImageAccessData> imageAccessDataList = new ArrayList<>();

    public List<ImageAccessData> radarCompute(String code,List<Eph> ephDatas, List<SphericalPoint> target) {
        List<ImageAccessData> windows_list = new ArrayList<>();
        for (SphericalPoint sphericalPoint_one : target) {
            List<ImageAccessData> imageAccessDatas = radarAccess_point(code,sphericalPoint_one, ephDatas);
            if (imageAccessDatas != null && imageAccessDatas.size() != 0) {
                windows_list.addAll(imageAccessDatas);
            }

        }
        return windows_list;
    }


    public List<ImageAccessData> radarAccess_point(String code,SphericalPoint sPoint, List<Eph> ephDatas) {

        bizService.getImageLimitation(code,rLimit);


        // 检查是否已初始化星历数据
        if (CollectionUtils.isEmpty(ephDatas) || ephDatas.size() < 2) {
            return null;
        }
        // 清空存储成像卫星访问计算结果的序列

        imageAccessDataList.clear();
        ImageAccessData accessData = new ImageAccessData();
        ImageAccessMessage accessMessage;
        int index = 0;

        // 定义pi值
        double pi = 3.14159265358979323846d;
        while (index < ephDatas.size() - 1) {

            Eph ephData = new Eph();
            String template=JSON.toJSONString(ephDatas.get(index));
            ephData = JSON.parseObject(template, Eph.class);

            // x为星下点到目标点球面距离，单位：弧度
            double x = Math.cos((90.0 - ephData.getLat()) * pi / 180.0)
                    * Math.cos((90.0 - sPoint.getLat()) * pi / 180.0)
                    + Math.sin((90.0 - ephData.getLat()) * pi / 180.0) * Math.sin((90.0 - sPoint.getLat()) * pi / 180.0)
                    * Math.cos(Math.abs(ephData.getLon() - sPoint.getLon()) * pi / 180.0);
            x = Math.acos(x);
            // 降轨成像：QWERT6A-1、QWERT6、QWERT9、QWERT10、QWERT14A-1
            if (rLimit.getImageSign() == 0) {
                if ((x <= rLimit.getMaxDistance()) && (ephData.getLat() > ephDatas.get(index + 1).getLat())) {
                    accessMessage = ImageComputeAll(sPoint, rLimit, ephData, index, ephDatas);

                    if ((accessMessage.getRollAngle() <= rLimit.getMaxRoll())
                            && (accessMessage.getRollAngle() >= rLimit.getMinRoll())
                            && (accessMessage.getSunElevation() >= rLimit.getMinSun())) {
                        accessData.setSatName(ephData.getCode());
                        accessData.setPass(Long.parseLong(ephData.getCirCle()));
                        accessData.setAccessTime(accessMessage.getAccessTime());
                        accessData.setAccessTimeend(accessMessage.getAccessTimeend());
                        accessData.setRollAngle(accessMessage.getRollAngle());
                        accessData.setSunElevation(accessMessage.getSunElevation());
                        imageAccessDataList.add(accessData);
                        // 如果计算得到第一次访问结果，按照重访周期直接跳到后面的星历开始计算，一般可见光卫星跳11小时
                        index = (int) (index + rLimit.getStepSpan());
                    } else {
                        // 如果计算没有得到访问结果，向后跳1圈左右的时间
                        index = (int) (index + (long) (1.3 * 60 * 60 / rLimit.getStepLength()));
                        // index++;
                    }
                }
            }
            // 成像与升降轨无关：QWERT5、QWERT7、QWERT11
            else if (rLimit.getImageSign() == 2) {
                if (x <= rLimit.getMaxDistance()) {
                    accessMessage = ImageComputeAll(sPoint, rLimit, ephData, index, ephDatas);
                    //20191123  侧摆角的绝对值做判断 gaoyue edit
                    double rollAngle2 = Math.abs(accessMessage.getRollAngle());
                    double rollAngle = accessMessage.getRollAngle();
                    if (rollAngle <= rLimit.getMaxRoll()
                            && (rollAngle >= rLimit.getMinRoll())
                            && (accessMessage.getSunElevation() >= rLimit.getMinSun())) {
                        accessData.setSatName(ephData.getCode());
                        accessData.setPass(Long.parseLong(ephData.getCirCle()));
                        accessData.setAccessTime(accessMessage.getAccessTime());
                        accessData.setAccessTimeend(accessMessage.getAccessTimeend());
                        accessData.setRollAngle(rollAngle2);
                        accessData.setSunElevation(accessMessage.getSunElevation());
                        imageAccessDataList.add(accessData);
                        index = (int) (index + rLimit.getStepSpan());
                    } else {
                        index = (int) (index + (long) (1.3 * 60 * 60 / rLimit.getStepLength()));
                        // index++;
                    }
                }
            }
            // 升轨成像：目前没有
            else if (rLimit.getImageSign() == 1) {
                if ((x <= rLimit.getMaxDistance())
                        && (ephDatas.get(index).getLat() < ephDatas.get(index + 1).getLat())) {
                    accessMessage = ImageComputeAll(sPoint, rLimit, ephData, index, ephDatas);

                    if ((accessMessage.getRollAngle() <= rLimit.getMaxRoll())
                            && (accessMessage.getRollAngle() >= rLimit.getMinRoll())
                            && (accessMessage.getSunElevation() >= rLimit.getMinSun())) {
                        accessData.setSatName(ephData.getCode());
                        accessData.setPass(Long.parseLong(ephData.getCirCle()));
                        accessData.setAccessTime(accessMessage.getAccessTime());
                        accessData.setAccessTimeend(accessMessage.getAccessTimeend());
                        accessData.setRollAngle(accessMessage.getRollAngle());
                        accessData.setSunElevation(accessMessage.getSunElevation());
                        imageAccessDataList.add(accessData);
                        index = (int) (index + rLimit.getStepSpan());
                    } else {
                        index = (int) (index + (long) (1.3 * 60 * 60 / rLimit.getStepLength()));
                    }
                }
            }

            index++;
        }
        return imageAccessDataList;
    }


    private double SphericalAngle(SphericalPoint source, SphericalPoint target) {
        // double pi=3.14159265358979323846;
        double back;
        SphericalPoint source1 = new SphericalPoint();
        SphericalPoint target1 = new SphericalPoint();
        source1.setLon(source.getLon() * Math.PI / 180.0);
        source1.setLat(source.getLat() * Math.PI / 180.0);
        target1.setLon(target.getLon() * Math.PI / 180.0);
        target1.setLat(target.getLat() * Math.PI / 180.0);

        back = Math.cos(Math.PI / 2 - source1.getLat()) * Math.cos(Math.PI / 2 - target1.getLat())
                + Math.sin(Math.PI / 2 - source1.getLat()) * Math.sin(Math.PI / 2 - target1.getLat())
                * Math.cos(Math.abs(source1.getLon() - target1.getLon()));

        return Math.acos(back) * 180.0 / Math.PI;
    }


    /**
     * 用于求各时刻太阳高度角
     *
     * @param sunTime
     * @param pSat
     * @return
     */
    public double sunAngle(Calendar sunTime, SphericalPoint pSat) {
        SphericalPoint sunPos = new SphericalPoint();
        // 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;
    }


    private ImageAccessMessage ImageComputeAll(SphericalPoint target, ImageLimitation limit, Eph ephData, int index,
                                               List<Eph> ephDatas) {
        ImageAccessMessage message = new ImageAccessMessage();
        message.setRollAngle(-180.0);

        int id = index;

        Eph ephdata1, ephdata2, ephdata3;

        SphericalPoint s1 = new SphericalPoint(), s2 = new SphericalPoint(), s3 = new SphericalPoint();

        PositionLLA pos1 = new PositionLLA(), pos2 = new PositionLLA(), pos3 = new PositionLLA();

        double dis1, dis2, dis3;

        double RollAngle;

        while (id < ephDatas.size() - 2) {
            ephdata1 = new Eph();
            ephdata2 = new Eph();

            String ephdata1str=JSON.toJSONString(ephDatas.get(id));
            ephdata1 = JSON.parseObject(ephdata1str, Eph.class);

            String ephdata2str=JSON.toJSONString(ephDatas.get(id + 1));
            ephdata2 = JSON.parseObject(ephdata2str, Eph.class);

            s1.setLon(ephdata1.getLon());
            s1.setLat(ephdata1.getLat());
            s2.setLon(ephdata2.getLon());
            s2.setLat(ephdata2.getLat());

            pos1.setLon(ephdata1.getLon());
            pos1.setLat(ephdata1.getLat());
            pos1.setHigh(ephdata1.getAlt());
            pos2.setLon(ephdata2.getLon());
            pos2.setLat(ephdata2.getLat());
            pos2.setHigh(ephdata2.getAlt());

            dis1 = SphericalAngle(s1, target);// 角度
            dis2 = SphericalAngle(s2, target);// 角度


            Calendar time2 = Calendar.getInstance();
            time2.setTime(ephdata2.getEphTime());

            ephdata3 = new Eph();
            String ephdata3str=JSON.toJSONString(ephDatas.get(id + 2));
            ephdata3 = JSON.parseObject(ephdata3str, Eph.class);

            Calendar time1 = Calendar.getInstance();
            time1.setTime(ephdata3.getEphTime());

            s3.setLon(ephdata3.getLon());
            s3.setLat(ephdata3.getLat());
            pos3.setLon(ephdata3.getLon());
            pos3.setLat(ephdata3.getLat());
            pos3.setHigh(ephdata3.getAlt());
            dis3 = SphericalAngle(s3, target);// 角度

            if ((dis2 <= dis1) && (dis2 <= dis3) && (dis2 <= 90.0)) {
               /* Calendar calendar1 = Calendar.getInstance();
                calendar1.setTimeInMillis(time1.getTimeInMillis());
                calendar1.add(Calendar.SECOND, -(int) limit.getXSpan());
                message.setAccessTimeend(calendar1.getTime());*/

                Calendar calendar = Calendar.getInstance();
                calendar.setTimeInMillis(time2.getTimeInMillis());
                calendar.add(Calendar.SECOND, -(int) limit.getXSpan());
                message.setAccessTime(calendar.getTime());

                // 太阳高度角
                message.setSunElevation(sunAngle(time2, s2));
                message.setSunElevation(calcSunAngleNew(time2, s2));

                // 侧摆角（滚动角）的绝对值
                message.setRollAngle(accessLLA(pos2, target));
                message.setRollAngle(Math.abs(message.getRollAngle()));

                // 轨道计算步长为1秒或更少时，精确计算侧视角度，提高精度
                if (limit.getStepLength() <= 1) {
                    OrbitElement sat;
                    PositionLLA bmin = pos2;
                    PositionVelocity PV0 = new PositionVelocity();
                    PV0.getR().setX(ephdata2.getX().doubleValue());
                    PV0.getR().setY(ephdata2.getY().doubleValue());
                    PV0.getR().setZ(ephdata2.getZ().doubleValue());
                    PV0.getV().setX(ephdata2.getVX().doubleValue());
                    PV0.getV().setY(ephdata2.getVY().doubleValue());
                    PV0.getV().setZ(ephdata2.getVZ().doubleValue());
                    sat = RvToElement(PV0);
                    double dismin = dis2;
                    double t = 0.01 / 806.8116341;
                    // double RollAngle;
                    // 以0.01秒为步进，只考虑J2、J3、J4的1阶、2阶长期项摄动把ephdata2点对应的轨道根数外推
                    for (int i = 1; i <= 100; i++) {
                        // 以0.1秒为步进，只考虑J2、J3、J4的1阶、2阶长期项摄动把ephdata2点对应的轨道根数外推
                        OrbitElement delta1 = Delta_1(sat);
                        OrbitElement delta2 = Delta_2(sat);
                        sat = sat.plus(delta1.plus(delta2).multiply(t));
                        // 轨道根数转换为准地固系测地坐标
                        PositionVelocity PV = ElementToRv(sat);
                        Position a = PV.getR();
                        // 历元平赤道地心惯性系坐标经岁差矩阵转换为瞬时平赤道地心系坐标，1毫秒
                        a = Precession(time2, a, 1.0 * i);
                        // 瞬时平赤道地心系坐标经章动矩阵转换为瞬时真赤道地心系坐标，1毫秒
                        a = Nutation(time2, a, 1.0 * i);
                        // 求赤经章动
                        double chijingnutation = LonNutation(time2, 1.0 * i);
                        // 求格林尼制真恒星时
                        double sg = SgTime(time2, chijingnutation, 1.0 * i);
                        // 瞬时真赤道地心系坐标经地球自转矩阵转换为准地固系坐标
                        a = Rotz(a, sg);
                        // 准地固系坐标（地心直角坐标）转换为准地固系坐标（测地坐标：测地纬度、测地经度、高度）
                        PositionLLA b = XyzTolla(a);
                        SphericalPoint smin = new SphericalPoint();
                        smin.setLat(b.getLat());
                        smin.setLon(b.getLon());
                        double distance = SphericalAngle(smin, target);
                        if (distance <= dismin) {
                            dismin = distance;
                            bmin = b;
                        }
                    }
                    // 侧摆角（滚动角）的绝对值
                    message.setRollAngle(accessLLA(bmin, target));
                    message.setRollAngle(Math.abs(message.getRollAngle()));
                }

                /*
                 * //另一种方法计算侧摆角，计算目标与卫星连线、星下点与卫星连线2条直线夹角，精度较差，待改进
                 * message.RollAngle=AccessXYZ(pos2,Target,ephdata2);
                 * message.RollAngle=fabs(message.RollAngle);
                 */

                // 判断侧视角度正负（精确判断方法见纸上）
                PositionLLA CBpdpos = new PositionLLA();
                double Const;
                SphericalPoint Targetnew;
                PositionLLA pos1new, pos2new, pos3new;
                Targetnew = target;
                pos1new = pos1;
                pos2new = pos2;
                pos3new = pos3;

                // 轨道跨正负180度经线情况,把-180至180经度转为0-360经度，用于判断侧摆正负
                double CbjIndex1 = pos3.getLon() * pos1.getLon();
                double CbjIndex2 = CbjIndex1 * target.getLon();
                if (Math.abs(pos1.getLon()) >= 155.0 || Math.abs(pos3.getLon()) >= 155.0
                        || Math.abs(target.getLon()) >= 155.0) {
                    if (pos1.getLon() < 0.0) {
                        pos1new.setLon(pos1.getLon() + 360.0);
                    }
                    if (pos2.getLon() < 0.0) {
                        pos2new.setLon(pos2.getLon() + 360.0);
                    }
                    if (pos3.getLon() < 0.0) {
                        pos3new.setLon(pos3.getLon() + 360.0);
                    }
                    if (target.getLon() < 0.0) {
                        Targetnew.setLon(target.getLon() + 360.0);
                    }
                }

                // 侧摆角正负精确判断，适用于太阳同步轨道，侧摆角正负判断最好还是按照QWERT14的方法用矢量判断
                CBpdpos.setLon(target.getLon());
                Const = (CBpdpos.getLon() - pos2.getLon()) / (pos1.getLon() - pos2.getLon());
                CBpdpos.setLat(pos2.getLat() + Const * (pos1.getLat() - pos2.getLat()));
                if (target.getLat() >= CBpdpos.getLat()) {
                    // SAR卫星
                    if ((limit.getSatName().equals("QWERT5-2")) || (limit.getSatName().equals("QWERT5-3"))
                            || (limit.getSatName().equals("QWERT7-2")) || (limit.getSatName().equals("QWERT7-3"))
                            || (limit.getSatName().equals("QWERT13A-1"))) {
                        message.setRollAngle(message.getRollAngle());
                    }
                    // 光学卫星
                    else {
                        message.setRollAngle(-1.0 * message.getRollAngle());
                    }
                } else if (target.getLat() < CBpdpos.getLat()) {
                    // SAR卫星
                    if ((limit.getSatName().equals("QWERT5-2")) || (limit.getSatName().equals("QWERT5-3"))
                            || (limit.getSatName().equals("QWERT7-2")) || (limit.getSatName().equals("QWERT7-3"))
                            || (limit.getSatName().equals("QWERT13A-1"))) {
                        message.setRollAngle(-1.0 * message.getRollAngle());
                    }
                    // 光学卫星
                    else {
                        message.setRollAngle(message.getRollAngle());
                    }
                }

                // if((message.RollAngle<=limit.MaxRoll)&&(message.RollAngle>=limit.MinRoll)&&(message.SunElevation>=limit.MinSun))
                // {
                return message;
                // }

                // return message;
            }

            // ephdata1=ephdata2;
            // ephdata2=ephdata3;
            // time1=time2;
            // time2=time3;
            // s1=s2;
            // s2=s3;
            // pos1=pos2;
            // pos2=pos3;
            // dis1=dis2;
            // dis2=dis3;

            id++;
        }

        return message;
    }

    public List<SphericalPoint> regionToListSphericalPoint(Region region) {
        List<SphericalPoint> point_inarea = new ArrayList<>();
        SphericalPoint sphericalPointNW = new SphericalPoint();
        sphericalPointNW.setLat(region.getNwNode().getLat());
        sphericalPointNW.setLon(region.getNwNode().getLon());
        point_inarea.add(sphericalPointNW);

        SphericalPoint sphericalPointNE = new SphericalPoint();
        sphericalPointNE.setLat(region.getNeNode().getLat());
        sphericalPointNE.setLon(region.getNeNode().getLon());
        point_inarea.add(sphericalPointNE);

        SphericalPoint sphericalPointSW = new SphericalPoint();
        sphericalPointSW.setLat(region.getSwNode().getLat());
        sphericalPointSW.setLon(region.getSwNode().getLon());
        point_inarea.add(sphericalPointSW);

        SphericalPoint sphericalPointSE = new SphericalPoint();
        sphericalPointSE.setLat(region.getSeNode().getLat());
        sphericalPointSE.setLon(region.getSeNode().getLon());
        point_inarea.add(sphericalPointSE);


        return point_inarea;
    }


    // 函数作用见本上说明，参数sat：卫星空间位置测地坐标，target：目标点测地坐标
    private double accessLLA(PositionLLA sat, SphericalPoint target) {

        // sat：卫星空间位置测地坐标，target：目标点测地坐标
        // Sa、S：星下点地心坐标，T：目标点地心坐标
        double pi = 3.14159265358979323846d;
        PositionLLR Sa;
        SphericalPoint S = new SphericalPoint(), T = new SphericalPoint();

        // 测地坐标转换为地心坐标
        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();
        //根据经纬度坐标，从文件中查询该点的高程值
        /*String srtm = Parameters.getInstance().getCoufigur().getPPAS_PARA().get("srtm")
                .toString();
        DEMReader demReader = SRTMReader.getSRTM_Reader(srtm);
        double lat = target.getLat() ;
        double lon = target.getLon();
        int dem1 = demReader.getdem(lat, lon);*/
        //暂时先把高程设置为0
        int dem1 = 0;
        OT = OT + 0.001 * dem1;


//		double lats = Sa.getLat() ;
//     	double lons = Sa.getLon();
//		int dem2 = demReader.getdem(lats, lons);
//		OS = OS + 0.001*dem2;//添加高程结束
        //
        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;
    }


    // 此函数用于测地坐标向地心坐标的转换，参数sat为测地坐标，见《球面天文学》37页图2.7
    private PositionLLR LLAToLLR(PositionLLA sat) {
        double a = 6378137, f = 1.0 / 298.257, b;
        double pi = 3.14159265358979323846d;
        double Radius, CoreAngle;
        PositionLLR back = new PositionLLR();
        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.getHigh();
        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));
    }


    /* 20190615 YHL增加 */
    public List<SphericalPoint> convertAreaToPoints(List<SphericalPoint> polygonAreaPoints, double distance) {
        if (polygonAreaPoints == null)
            return null;
        if (polygonAreaPoints.size() <= 2)
            return polygonAreaPoints;
        SphericalPoint firstAreaPoints = null;
        SphericalPoint previousAreaPoints = null;
        List<SphericalPoint> resultPoints = new ArrayList<>();
        // 将多边形转为外接矩形
        List<Double> rec_XY = convertPolygonToRectangleXY(polygonAreaPoints);
        // 将外接矩形转为点阵
        List<SphericalPoint> rec_pointsArray = convertRectangleToPointsArray(rec_XY.get(0), rec_XY.get(1),
                rec_XY.get(2), rec_XY.get(3), distance);
        // 对点阵中的逐点进行判断，是否在多边形中,在多边形中的，放入结果列表中。
        List<SphericalPoint> polygon_pointsArray = new ArrayList();
        List<SphericalPoint> points_outside_Array = new ArrayList();

        for (SphericalPoint point : rec_pointsArray) {
            if (!containsPoint(polygonAreaPoints, point)) {
                points_outside_Array.add(point);
                // 如果没有与顶点重和，添加进入，否则不加入
                // boolean repeated = false;
                // for (SphericalPoint vertexPoint : polygonAreaPoints) {
                // if (point.getLat() == vertexPoint.getLat() && point.getLon()
                // == vertexPoint.getLon()) {
                // repeated = true;
                // continue;
                // }
                // }

                // if (!repeated)
                // polygon_pointsArray.add(point);
            }
        }
        rec_pointsArray.removeAll(points_outside_Array);
        polygon_pointsArray.addAll(rec_pointsArray);
        return polygon_pointsArray;

    }

    /* 20190615 YHL增加 */
    public List<SphericalPoint> convertRectangleToPointsArray(double x_min, double x_max, double y_min, double y_max,
                                                              double distance) {
        List<List<SphericalPoint>> points_array = convertRectangleToPointsArray_rows(x_min, x_max, y_min, y_max,
                distance);
        List<SphericalPoint> points_output = new ArrayList<>();
        if (points_array != null) {
            for (List<SphericalPoint> points_one_row : points_array) {
                if (points_one_row != null)
                    points_output.addAll(points_one_row);
            }
        }
        return points_output;
    }

    /* 20190615 YHL增加 */
    public List<List<SphericalPoint>> convertRectangleToPointsArray_rows(double x_min, double x_max, double y_min,
                                                                         double y_max, double distance) {
        // double
        // x_min=rectangle.get(0),x_max=rectangle.get(1),y_min=rectangle.get(2),y_max=rectangle.get(3);
        List<SphericalPoint> list = new ArrayList<>();

        List<List<SphericalPoint>> points_array = new ArrayList<>();
        // points_array.add(new SphericalPoint(x_min, y_min));
        // points_array.add(new SphericalPoint(x_max, y_min));
        // points_array.add(new SphericalPoint(x_min, y_max));
        // points_array.add(new SphericalPoint(x_max, y_max));
        boolean x_tail, y_tail;
        boolean cross180 = (180 - Math.abs(x_min) + (180 - Math.abs(x_max)) < (x_max - x_min)) ? true : false;

        for (double y = y_min; y < y_max + distance; y += distance) {
            List<SphericalPoint> points_row = new ArrayList<>();
            boolean flag_end = false;
            if (y >= y_max) {
                flag_end = true;
                y = y_max;
            }
            double coord_y;
            if (y > 90) {
                coord_y = 90 - (y - 90);
            } else {
                coord_y = y;
            }

            if (cross180) {
                double x;
                for (x = x_min; x > -180; x += -distance) {
                    double coord_x;

                    // -175至175
                    if (x < -180) {
                        coord_x = 180 + (x + 180);
                    } else if (x > 180) {
                        coord_x = -180 + (x - 180);
                    } else {
                        coord_x = x;
                    }

                    points_row.add(new SphericalPoint(coord_x, coord_y));

                }
                for (x = 180 + (180 + x); x > x_max; x += -distance) {
                    double coord_x;

                    // -175至175
                    if (x < -180) {
                        coord_x = 180 + (x + 180);
                    } else if (x > 180) {
                        coord_x = -180 + (x - 180);
                    } else {
                        coord_x = x;
                    }

                    points_row.add(new SphericalPoint(coord_x, coord_y));

                }
            } else {
                for (double x = x_min; x < x_max; x += distance) {
                    double coord_x;

                    // -175至175
                    if (x < -180) {
                        coord_x = 180 + (x + 180);
                    } else if (x > 180) {
                        coord_x = -180 + (x - 180);
                    } else {
                        coord_x = x;
                    }

                    points_row.add(new SphericalPoint(coord_x, coord_y));

                }
            }
            // 每行最后一点
            points_row.add(new SphericalPoint(x_max, y));

            points_array.add(points_row);
            if (flag_end) {
                break;
            }
        }
        // //添加最后一行
        // List<SphericalPoint> points_row_last = new ArrayList<>();
        // for (double x = x_min ; x < x_max; x += distance) {
        // double coord_x;
        // x_tail = false;

        // if (x > 180) {
        // coord_x = -180 + (x-180 );
        // }else{
        // coord_x=x;
        // }

        // points_row_last.add(new SphericalPoint(coord_x, y_max));

        // }
        // //添加最后一行，最后一个点
        // points_row_last.add(new SphericalPoint(x_max, y_max));

        return points_array;
    }


    public List<Double> convertPolygonToRectangleXY(List<SphericalPoint> polygonPoints) {
        double x_min = 0, x_max = 0, y_min = 0, y_max = 0;
        int count = 0;
        for (SphericalPoint point : polygonPoints) {
            if (count == 0) {
                x_min = point.getLon();
                y_min = point.getLat();
                x_max = point.getLon();
                y_max = point.getLat();
                count++;
                continue;
            } else {
                if (x_min > point.getLon())
                    x_min = point.getLon();
                if (x_max < point.getLon())
                    x_max = point.getLon();
                if (y_min > point.getLat())
                    y_min = point.getLat();
                if (y_max < point.getLat())
                    y_max = point.getLat();
            }
        }

        List<Double> rectangle_xyminmaxlist = new ArrayList<>();
        rectangle_xyminmaxlist.add(new Double(x_min));
        rectangle_xyminmaxlist.add(new Double(x_max));
        rectangle_xyminmaxlist.add(new Double(y_min));
        rectangle_xyminmaxlist.add(new Double(y_max));
        return rectangle_xyminmaxlist;
    }


    /* 20190615 YHL增加 */
    public boolean containsPoint(List<SphericalPoint> polygPoints, SphericalPoint p) {
        double[] point = new double[2];
        double[][] vertices = new double[polygPoints.size()][2];

        //gaoyue edit 20190625
        double x_min = 0, x_max = 0, y_min = 0, y_max = 0;
        int count = 0;
        for (SphericalPoint pp : polygPoints) {
            if (count == 0) {
                x_min = pp.getLon();
                y_min = pp.getLat();
                x_max = pp.getLon();
                y_max = pp.getLat();
                count++;
                continue;
            } else {
                if (x_min > pp.getLon())
                    x_min = pp.getLon();
                if (x_max < pp.getLon())
                    x_max = pp.getLon();
                if (y_min > pp.getLat())
                    y_min = pp.getLat();
                if (y_max < pp.getLat())
                    y_max = pp.getLat();
            }
        }
        if (Math.abs(x_min) + Math.abs(x_max) > 180) {
            point[0] = p.getLat() + 90;
            point[1] = p.getLon() < 0 ? p.getLon() + 360 : p.getLon();
            for (int i = 0; i < polygPoints.size(); i++) {
                vertices[i][0] = polygPoints.get(i).getLat() + 90;
                vertices[i][1] = polygPoints.get(i).getLon() < 0 ? polygPoints.get(i).getLon() + 360 : polygPoints.get(i).getLon();
            }
        } else {
            point[0] = p.getLat() + 90;
            point[1] = p.getLon();
            for (int i = 0; i < polygPoints.size(); i++) {
                vertices[i][0] = polygPoints.get(i).getLat() + 90;
                vertices[i][1] = polygPoints.get(i).getLon();
            }
        }
        // edit end 20190625
        /*
        point[0] = p.getLat() + 90;
        point[1] = p.getLon() +180 < 0 ? p.getLon() + 360 : p.getLon();

        for (int i = 0; i < polygPoints.size(); i++) {
            vertices[i][0] = polygPoints.get(i).getLat() + 90;
            vertices[i][1] = polygPoints.get(i).getLon() < 0 ? polygPoints.get(i).getLon() + 360: polygPoints.get(i).getLon();
        }
        */
        int verticesCount = vertices.length;
        int nCross = 0;
        for (int i = 0; i < verticesCount; ++i) {
            double[] p1 = vertices[i];
            double[] p2 = vertices[(i + 1) % verticesCount];

            // 求解 y=p.y 与 p1 p2 的交点
            if (p1[0] == p2[0]) { // p1p2 与 y=p0.y平行
                continue;
            }
            if (point[0] < Math.min(p1[0], p2[0])) { // 交点在p1p2延长线上
                continue;
            }
            if (point[0] >= Math.max(p1[0], p2[0])) { // 交点在p1p2延长线上
                continue;
            }
            // 求交点的 X 坐标
            double x = (point[0] - p1[0]) * (p2[1] - p1[1]) / (p2[0] - p1[0]) + p1[1];
            if (x > point[1]) { // 只统计单边交点
                nCross++;
            }
        }
        // 单边交点为偶数，点在多边形之外
        return (nCross % 2 == 1);
    }

    /**
     * 调用可见光的太阳高度角计算方法
     *
     * @param time2
     * @param s2
     * @return
     */
    private double calcSunAngleNew(Calendar time2, SphericalPoint s2) {
        double sunangle = 0;
        double Longitude = s2.getLon();
        double Latitude = s2.getLat();
        sunangle = SunElevationAngleCalculator.calculateSunElevationAngle(time2, Longitude, Latitude);
//		return sunangle;201812sunqy
        if (sunangle < 0) {
            return -90;
        } else {
            return sunangle;
        }
    }

    private OrbitElement RvToElement(PositionVelocity Pos) {
        double rx, ry, rz, vx, vy, vz, a, e, i, raan, aop, m;
        double r0, v0, rv0, me, gm;
        double a0, time, pi;
        double p3, q3, r1, r2, r3;

        rx = Pos.getR().getX();
        ry = Pos.getR().getY();
        rz = Pos.getR().getZ();
        vx = Pos.getV().getX();
        vy = Pos.getV().getY();
        vz = Pos.getV().getZ();

        // gm：地心引力常数
        gm = 398600.5e+9;
        pi = 3.14159265358979323846;
        // a0：地球参考椭球体赤道半径
        a0 = 6378137.0;
        // time：时间单位
        time = Math.sqrt(Math.pow(a0, 3) / gm);

        rx = rx / a0;
        ry = ry / a0;
        rz = rz / a0;
        vx = vx * time / a0;
        vy = vy * time / a0;
        vz = vz * time / a0;

        r0 = Math.sqrt(rx * rx + ry * ry + rz * rz);
        v0 = Math.sqrt(vx * vx + vy * vy + vz * vz);
        rv0 = rx * vx + ry * vy + rz * vz;

        a = 1.0 / (2.0 / r0 - Math.pow(v0, 2));
        me = Math.atan(rv0 / Math.sqrt(a) / (1.00 - r0 / a));
        if (rv0 / Math.sqrt(a) >= 0 && (1.00 - r0 / a) >= 0)
            me = me;
        if (rv0 / Math.sqrt(a) >= 0 && (1.00 - r0 / a) <= 0)
            me = me + pi;
        if (rv0 / Math.sqrt(a) < 0 && (1.00 - r0 / a) < 0)
            me = me + pi;
        if (rv0 / Math.sqrt(a) < 0 && (1.00 - r0 / a) > 0)
            me = me + 2 * pi;
        e = (1.00 - r0 / a) / Math.cos(me);
        m = me - e * Math.sin(me);

        p3 = Math.cos(me) / r0 * rz - Math.sqrt(a) * Math.sin(me) * vz;
        q3 = (Math.sin(me) / r0 * rz + Math.sqrt(a) * (Math.cos(me) - e) * vz) / Math.sqrt(1.00 - e * e);
        r1 = (ry * vz - rz * vy) / Math.sqrt(a * (1.00 - e * e));
        r2 = (rz * vx - rx * vz) / Math.sqrt(a * (1.00 - e * e));
        r3 = (rx * vy - ry * vx) / Math.sqrt(a * (1.00 - e * e));

        i = Math.acos(r3);

        aop = Math.atan(p3 / q3);
        if (p3 >= 0 && q3 >= 0) {
            aop = aop;
        }
        if (p3 >= 0 && q3 < 0) {
            aop = aop + pi;
        }
        if (p3 < 0 && q3 < 0) {
            aop = aop + pi;
        }
        if (p3 < 0 && q3 >= 0) {
            aop = aop + 2 * pi;
        }

        raan = Math.atan(-r1 / r2);
        if (r1 >= 0 && -r2 >= 0) {
            raan = raan;
        }
        if (r1 >= 0 && -r2 < 0) {
            raan = raan + pi;
        }
        if (r1 < 0 && -r2 < 0) {
            raan = raan + pi;
        }
        if (r1 < 0 && -r2 >= 0) {
            raan = raan + 2 * pi;
        }

        a = a * a0;
        i = i * 180.0 / pi;
        raan = raan * 180.0 / pi;
        aop = aop * 180.0 / pi;
        m = m * 180.0 / pi;

        OrbitElement back = new OrbitElement();
        back.setA(a);
        back.setE(e);
        back.setI(i);
        back.setAop(aop);
        back.setRaan(raan);
        back.setM(m);

        return back;
    }


    // Delta_1为J2、J3、J4摄动的0阶和1阶长期项摄动
    private OrbitElement Delta_1(OrbitElement sat) {
        double p, a2, j2 = 1.082636e-03, n;
        double a, e, i, aop, raan, M;
        double pi = 3.14159265358979323846d, a0 = 6378137.0d;
        OrbitElement delta = new OrbitElement();
        a2 = 1.5 * j2;
        a = sat.getA() / a0;
        e = sat.getE();
        i = sat.getI() * pi / 180.0;
        raan = sat.getRaan() * pi / 180.0;
        aop = sat.getAop() * pi / 180.0;
        M = sat.getM() * pi / 180.0;
        p = a * (1.0 - e * e);
        n = Math.pow(a, -1.5);

        delta.setA(0);
        delta.setE(0);
        delta.setI(0);

        delta.setRaan(-a2 / Math.pow(p, 2) * n * Math.cos(i));
        delta.setAop(a2 / Math.pow(p, 2) * n * (2.0 - 2.5 * Math.pow(Math.sin(i), 2)));
        delta.setM(n + a2 / Math.pow(p, 2) * n * (1.0 - 1.5 * Math.pow(Math.sin(i), 2)) * Math.sqrt(1.0 - e * e));

        delta.setRaan(delta.getRaan() * 180.0 / pi);
        delta.setAop(delta.getAop() * 180.0 / pi);
        delta.setM(delta.getM() * 180.0 / pi);

        return delta;
    }


    // Delta_2为J2、J3、J4摄动的2阶长期项摄动
    private OrbitElement Delta_2(OrbitElement sat) {
        double p, a2, a4, j2 = 1.082636e-03, j4 = -1.619e-06, n;
        double a, e, i, aop, raan, M;
        double pi = 3.14159265358979323846, a0 = 6378137.0;
        OrbitElement delta = new OrbitElement();
        a2 = 1.5 * j2;
        a4 = -35 / 8.0 * j4;
        a = sat.getA() / a0;
        e = sat.getE();
        i = sat.getI() * pi / 180.0;
        raan = sat.getRaan() * pi / 180.0;
        aop = sat.getAop() * pi / 180.0;
        M = sat.getM() * pi / 180.0;
        p = a * (1.0 - e * e);
        n = Math.pow(a, -1.5);

        delta.setA(0);
        delta.setE(0);
        delta.setI(0);

        delta.setRaan(-a2 * a2 / Math.pow(p, 4) * n * Math.cos(i)
                * ((1.5 + 1 / 6.0 * e * e + Math.sqrt(1.0 - e * e))
                - Math.pow(Math.sin(i), 2) * (5.0 / 3.0 - 5 / 24.0 * e * e + 1.5 * Math.sqrt(1.0 - e * e))
                + (a4 / Math.pow(a2, 2))
                * ((6 / 7.0 + 9 / 7.0 * e * e) - Math.pow(Math.sin(i), 2) * (1.5 + 9 / 4.0 * e * e))));

        delta.setAop(a2 * a2 / Math.pow(p, 4) * n * ((4 + 7 / 12.0 * e * e + 2 * Math.sqrt(1.0 - e * e))
                - Math.pow(Math.sin(i), 2) * (103 / 12.0 + 3 / 8.0 * e * e + 11 / 2.0 * Math.sqrt(1 - e * e))
                + Math.pow(Math.sin(i), 4) * (215 / 48.0 - 15 / 32.0 * e * e + 15 / 4.0 * Math.sqrt(1.0 - e * e))
                + a4 / Math.pow(a2, 2)
                * ((12 / 7.0 + 27 / 14.0 * e * e) - Math.pow(Math.sin(i), 2) * (93 / 14.0 + 27 / 4.0 * e * e)
                + Math.pow(Math.sin(i), 4) * (21 / 4.0 + 81 / 16.0 * e * e))));

        delta.setM(
                a2 * a2 / Math.pow(p, 4)
                        * n * Math
                        .sqrt(1.0
                                - e * e)
                        * (0.5 * (1.0 - 1.5 * Math.pow(Math.sin(i), 2)) * (1.0 - 1.5 * Math.pow(Math.sin(i), 2))
                        * Math.sqrt(1.0 - e * e)
                        + (2.5 + 10 / 3.0 * e * e) - Math
                        .pow(Math.sin(i),
                                2)
                        * (19 / 3.0
                        + 26 / 3.0 * e * e)
                        + Math.pow(Math.sin(i), 4) * (233 / 48.0 + 103 / 12.0 * e * e)
                        + Math.pow(e, 4) / (1.0 - e * e)
                        * (35 / 12.0 - 35 / 4.0 * Math.pow(Math.sin(i), 2)
                        + 315 / 32.0 * Math.pow(Math.sin(i), 4))
                        + a4 / Math.pow(a2, 2) * e * e * (9 / 14.0 - 45 / 14.0 * Math.pow(Math.sin(i), 2)
                        + 45 / 16.0 * Math.pow(Math.sin(i), 4))));

        delta.setRaan(delta.getRaan() * 180.0 / pi);
        delta.setAop(delta.getAop() * 180.0 / pi);
        delta.setM(delta.getM() * 180.0 / pi);

        return delta;
    }


    // ?公式未知，此函数用地心纬度求地表不同纬度圈地心距，参数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);
    }


    // 此函数把轨道根数转化为瞬时真赤道地心系（与J2000.0惯性系差岁差、章动）空间位置、速度
    private PositionVelocity ElementToRv(OrbitElement sat) {
        // SatPoint位置、速度
        PositionVelocity SatPoint = new PositionVelocity();
        double a0 = 6378137.0d, a, e, i, aop, raan, M;
        double pi = 3.14159265358979323846d;
        double p1, p2, p3, q1, q2, q3;
        double Me0 = 0, Me, f;
        // 半长轴以地球半径为单位
        a = sat.getA() / a0;
        e = sat.getE();
        i = sat.getI() * pi / 180.0;
        raan = sat.getRaan() * pi / 180.0;
        aop = sat.getAop() * pi / 180.0;
        M = sat.getM() * pi / 180.0;
        // 迭代法用平近点角M求偏近点角Me
        Me = M;
        while (Math.abs(Me - Me0) > 0.00000000001) {
            Me0 = Me;
            Me = M + e * Math.sin(Me);
        }

        while (Me < 0) {
            Me = Me + 2 * pi;
        }
        // while(Me>2*pi)
        // {
        Me = 2 * pi;
        // Me=Me-2*pi;
        // }
        // 求真近点角f
        f = Metof(e, Me * 180.0 / pi);

        p1 = Math.cos(raan) * Math.cos(aop) - Math.sin(raan) * Math.sin(aop) * Math.cos(i);
        p2 = Math.sin(raan) * Math.cos(aop) + Math.cos(raan) * Math.sin(aop) * Math.cos(i);
        p3 = Math.sin(aop) * Math.sin(i);

        q1 = -Math.cos(raan) * Math.sin(aop) - Math.sin(raan) * Math.cos(aop) * Math.cos(i);
        q2 = -Math.sin(raan) * Math.sin(aop) + Math.cos(raan) * Math.cos(aop) * Math.cos(i);
        q3 = Math.cos(aop) * Math.sin(i);

        SatPoint.getR().setX(a * a0 * (Math.cos(Me) - e) * p1 + a * a0 * Math.sqrt(1.0 - e * e) * Math.sin(Me) * q1);
        SatPoint.getR().setY(a * a0 * (Math.cos(Me) - e) * p2 + a * a0 * Math.sqrt(1.0 - e * e) * Math.sin(Me) * q2);
        SatPoint.getR().setZ(a * a0 * (Math.cos(Me) - e) * p3 + a * a0 * Math.sqrt(1.0 - e * e) * Math.sin(Me) * q3);

        double r = a * (1.0 - e * Math.cos(Me));

        SatPoint.getV().setX(Math.sqrt(a) / r * (-Math.sin(Me) * p1 + Math.sqrt(1 - e * e) * Math.cos(Me) * q1));
        SatPoint.getV().setY(Math.sqrt(a) / r * (-Math.sin(Me) * p2 + Math.sqrt(1 - e * e) * Math.cos(Me) * q2));
        SatPoint.getV().setZ(Math.sqrt(a) / r * (-Math.sin(Me) * p3 + Math.sqrt(1 - e * e) * Math.cos(Me) * q3));
        // ct是使引力常数G=1的时间单位
        double ct = 806.8110649227;

        SatPoint.getV().setX(SatPoint.getV().getX() * a0 / ct);
        SatPoint.getV().setY(SatPoint.getV().getY() * a0 / ct);
        SatPoint.getV().setZ(SatPoint.getV().getZ() * a0 / ct);

        return SatPoint;
    }

    private double Metof(double e, double me) {
        double f;
        // double pi=3.14159265358979323846;
        me = me * Math.PI / 180.0;
        f = 2 * Math.atan(Math.sqrt((1 + e) / (1 - e)) * Math.tan(me / 2));

        while (f > 2 * Math.PI) {
            f = f - 2 * Math.PI;
        }
        while (f < 0) {
            f = f + 2 * Math.PI;
        }

        return f * 180.0 / Math.PI;
    }


    // 函数把历元平赤道地心系坐标经岁差矩阵转换为瞬时平赤道地心系坐标
    // 参数sat：历元平赤道地心系坐标，返回值sat：瞬时平赤道地心系坐标
    private Position Precession(Calendar newtime, Position sat, double ms) {
        double d1, d2, d3;
        Calendar starttime = TimeUtils.getCalendar(2000, 01, 01, 20, 0, 0);
        long int_time = (newtime.getTimeInMillis() - starttime.getTimeInMillis()) / 1000;
        double rlxj = (double) int_time + ms;
        // rlxj为星历历元时间到J2000.0的儒略世纪数
        rlxj = rlxj / 86400.00 / 36525.0;
        // 岁差量计算，单位：度
        d1 = 2306.21810 / 3600.00 * rlxj + 0.301880 / 3600.00 * Math.pow(rlxj, 2)
                + 0.0179980 / 3600.0 * Math.pow(rlxj, 3);
        d2 = 2004.31090 / 3600.00 * rlxj - 0.426650 / 3600.00 * Math.pow(rlxj, 2)
                - 0.0418330 / 3600.0 * Math.pow(rlxj, 3);
        d3 = 2306.21810 / 3600.00 * rlxj + 1.094680 / 3600.00 * Math.pow(rlxj, 2)
                + 0.0182030 / 3600.0 * Math.pow(rlxj, 3);

        sat = Rotz(sat, -d1);
        sat = Roty(sat, d2);
        sat = Rotz(sat, -d3);
        return sat;
    }


    // 函数用测地纬度求地表上点的地心纬度，返回值地心纬度，参数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;
    }


    private Position Rotx(Position sat, double sita) {
        Position sat0 = new Position();
        double x0, y0, z0;
        ;
        // double pi=3.14159265358979323846;
        sita = sita * Math.PI / 180.0;

        x0 = sat.getX();
        y0 = sat.getY() * Math.cos(sita) + sat.getZ() * Math.sin(sita);
        z0 = -sat.getY() * Math.sin(sita) + sat.getZ() * Math.cos(sita);
        sat0.setX(x0);
        sat0.setY(y0);
        sat0.setZ(z0);
        return sat0;
    }


    private Position Roty(Position sat, double sita) {
        double x0, y0, z0;
        // double pi=3.14159265358979323846;
        sita = sita * Math.PI / 180.0;
        Position sat0 = new Position();

        x0 = sat.getX() * Math.cos(sita) - sat.getZ() * Math.sin(sita);
        y0 = sat.getY();
        z0 = sat.getX() * Math.sin(sita) + sat.getZ() * Math.cos(sita);
        sat0.setX(x0);
        sat0.setY(y0);
        sat0.setZ(z0);
        return sat0;
    }

    private Position Rotz(Position sat, double sita) {
        double x0, y0, z0;
        Position sat0 = new Position();
        // double pi=3.14159265358979323846;
        sita = sita * Math.PI / 180.0;
        x0 = sat.getX() * Math.cos(sita) + sat.getY() * Math.sin(sita);
        y0 = -sat.getX() * Math.sin(sita) + sat.getY() * Math.cos(sita);
        z0 = sat.getZ();
        sat0.setX(x0);
        sat0.setY(y0);
        sat0.setZ(z0);
        return sat0;
    }


    // 函数把瞬时平赤道地心系坐标经章动矩阵转换为瞬时真赤道地心系坐标
    // 参数sat：瞬时平赤道地心系坐标，返回值sat：瞬时真赤道地心系坐标
    private Position Nutation(Calendar newtime, Position sat, double ms) {
        // double pi=3.14159265358979323846;
        double rlxj, a1, a2, a3, a4, a5, delta_lon, delta_obliquity, hcjj;
        double[] xh = new double[107];
        double[] c_lon = new double[107];
        double[] c_obliquity = new double[107];
        for (int i = 0; i < 107; i++) {
            xh[i] = 0;
            c_lon[i] = 0;
            c_obliquity[i] = 0;
        }
        Calendar starttime = TimeUtils.getCalendar(2000, 01, 01, 20, 0, 0);
        long int_time = (newtime.getTimeInMillis() - starttime.getTimeInMillis()) / 1000;
        rlxj = (double) int_time + ms;
        // rlxj为星历历元时间到J2000.0的儒略世纪数
        rlxj = rlxj / 86400.00 / 36525.00;
        // hcjj为考虑岁差影响的平黄赤交角
        hcjj = 84381.4480 / 3600.0 - 46.81500 / 3600.0 * rlxj - 0.00059 / 3600.0 * Math.pow(rlxj, 2)
                + 0.0018130 / 3600.0 * Math.pow(rlxj, 3);

        a1 = 134.96298139 + (1325 * 360 + 198.8673981) * rlxj + 8.69722e-3 * Math.pow(rlxj, 2)
                + 1.78e-5 * Math.pow(rlxj, 3);
        a2 = 357.5277230 + (99 * 360 + 359.05034) * rlxj - 1.6e-4 * Math.pow(rlxj, 2) - 3.33e-6 * Math.pow(rlxj, 3);
        a3 = 93.27191028 + (1342 * 360 + 82.017538) * rlxj - 0.0036825 * Math.pow(rlxj, 2)
                + 3.06e-6 * Math.pow(rlxj, 3);
        a4 = 297.8503630 + (1236 * 360 + 307.11148) * rlxj - 1.914167e-3 * Math.pow(rlxj, 2)
                + 5.28e-6 * Math.pow(rlxj, 3);
        a5 = 125.0445220 - (5 * 360 + 134.13626) * rlxj + 2.0708e-3 * Math.pow(rlxj, 2) + 2.222e-6 * Math.pow(rlxj, 3);

        xh[1] = a5;
        xh[2] = 2 * a5;
        xh[3] = -2 * a1 + 2 * a3 + a5;
        xh[4] = 2 * a1 - 2 * a3;
        xh[5] = -2 * a1 + 2 * a3 + 2 * a5;
        xh[6] = a1 - a2 - a4;
        xh[7] = -2 * a2 + 2 * a3 - 2 * a4 + a5;
        xh[8] = 2 * a1 - 2 * a3 + a5;
        xh[9] = 2 * a3 - 2 * a4 + 2 * a5;
        xh[10] = a2;
        xh[11] = a2 + 2 * a3 - 2 * a4 + 2 * a5;
        xh[12] = -a2 + 2 * a3 - 2 * a4 + 2 * a5;
        xh[13] = 2 * a3 - 2 * a4 + a5;
        xh[14] = 2 * a1 - 2 * a4;
        xh[15] = 2 * a3 - 2 * a4;
        xh[16] = 2 * a2;
        xh[17] = a2 + a5;
        xh[18] = 2 * a2 + 2 * a3 - 2 * a4 + 2 * a5;
        xh[19] = -a2 + a5;
        xh[20] = -2 * a1 + 2 * a4 + a5;
        xh[21] = -a2 + 2 * a3 - 2 * a4 + a5;
        xh[22] = 2 * a1 - 2 * a4 + a5;
        xh[23] = a2 + 2 * a3 - 2 * a4 + a5;
        xh[24] = a1 - a4;
        xh[25] = 2 * a1 + a2 - 2 * a4;
        xh[26] = -2 * a3 + 2 * a4 + a5;
        xh[27] = a2 - 2 * a3 + 2 * a4;
        xh[28] = a2 + 2 * a5;
        xh[29] = -a1 + a5;
        xh[30] = a2 + 2 * a3 - 2 * a4;
        xh[31] = 2 * a3 + 2 * a5;
        xh[32] = a1;
        xh[33] = 2 * a3 + a5;
        xh[34] = a1 + 2 * a3 + 2 * a5;
        xh[35] = a1 - 2 * a4;
        xh[36] = -a1 + 2 * a3 + 2 * a5;
        xh[37] = 2 * a4;
        xh[38] = a1 + a5;
        xh[39] = -a1 + a5;
        xh[40] = -a1 + 2 * a3 + a5;
        xh[41] = a1 + 2 * a3 + a5;
        xh[42] = 2 * a3 + a5;
        xh[43] = 2 * a1 + a5;
        xh[44] = a1 + 2 * a3 + a5;
        xh[45] = 2 * a1 + 2 * a3 + 2 * a5;
        xh[46] = 2 * a3 + a5;
        xh[47] = -a1 + 2 * a3 + a5;
        xh[48] = -a1 + a5;
        xh[49] = a1 - 2 * a4 + a5;
        xh[50] = -a1 + 2 * a3 + 2 * a4 + a5;

        c_lon[1] = -171996 - 174.2 * rlxj;
        c_lon[2] = 2062 + 0.2 * rlxj;
        c_lon[3] = 46;
        c_lon[4] = 11;
        c_lon[5] = -3;
        c_lon[6] = -3;
        c_lon[7] = -2;
        c_lon[8] = 1;
        c_lon[9] = -13187 - 1.6 * rlxj;
        c_lon[10] = 1426 - 3.4 * rlxj;
        c_lon[11] = -517 + 1.2 * rlxj;
        c_lon[12] = 217 - 0.5 * rlxj;
        c_lon[13] = 129 + 0.1 * rlxj;
        c_lon[14] = 48;
        c_lon[15] = -22;
        c_lon[16] = 17 - 0.1 * rlxj;
        c_lon[17] = -15;
        c_lon[18] = -16 + 0.1 * rlxj;
        c_lon[19] = -12;
        c_lon[20] = -6;
        c_lon[21] = -5;
        c_lon[22] = 4;
        c_lon[23] = 4;
        c_lon[24] = -4;
        c_lon[25] = 1;
        c_lon[26] = 1;
        c_lon[27] = -1;
        c_lon[28] = 1;
        c_lon[29] = 1;
        c_lon[30] = -1;
        c_lon[31] = -2274 - 0.2 * rlxj;
        c_lon[32] = 712 + 0.1 * rlxj;
        c_lon[33] = -386 - 0.4 * rlxj;
        c_lon[34] = -301;
        c_lon[35] = -158;
        c_lon[36] = 123;
        c_lon[37] = 63;
        c_lon[38] = 63 + 0.1 * rlxj;
        c_lon[39] = -58 - 0.1 * rlxj;
        c_lon[40] = -59;
        c_lon[41] = -51;
        c_lon[42] = -38;
        c_lon[43] = 29;
        c_lon[44] = 29;
        c_lon[45] = -31;
        c_lon[46] = 26;
        c_lon[47] = 21;
        c_lon[48] = 16;
        c_lon[49] = -13;
        c_lon[50] = -10;

        c_obliquity[1] = 92025 + 8.9 * rlxj;
        c_obliquity[2] = -895 + 0.5 * rlxj;
        c_obliquity[3] = -24;
        c_obliquity[4] = 0;
        c_obliquity[5] = 1;
        c_obliquity[6] = 0;
        c_obliquity[7] = 1;
        c_obliquity[8] = 0;
        c_obliquity[9] = 5736.00 - 3.1 * rlxj;
        c_obliquity[10] = 54 - 0.1 * rlxj;
        c_obliquity[11] = 224 + 0.6 * rlxj;
        c_obliquity[12] = -95 + 0.3 * rlxj;
        c_obliquity[13] = -70;
        c_obliquity[14] = 1;
        c_obliquity[15] = 0;
        c_obliquity[16] = 0;
        c_obliquity[17] = 9;
        c_obliquity[18] = 7;
        c_obliquity[19] = 6;
        c_obliquity[20] = 3;
        c_obliquity[21] = 3;
        c_obliquity[22] = -2;
        c_obliquity[23] = -2;
        c_obliquity[24] = 0;
        c_obliquity[25] = 0;
        c_obliquity[26] = 0;
        c_obliquity[27] = 0;
        c_obliquity[28] = 0;
        c_obliquity[29] = 0;
        c_obliquity[30] = 0;
        c_obliquity[31] = 977 - 0.5 * rlxj;
        c_obliquity[32] = -7;
        c_obliquity[33] = 200;
        c_obliquity[34] = 129 - 0.1 * rlxj;
        c_obliquity[35] = -1;
        c_obliquity[36] = -53;
        c_obliquity[37] = -2;
        c_obliquity[38] = -33;
        c_obliquity[39] = 32;
        c_obliquity[40] = 26;

        c_obliquity[41] = 27;
        c_obliquity[42] = 16;
        c_obliquity[43] = -1;
        c_obliquity[44] = -12;
        c_obliquity[45] = 13;
        c_obliquity[46] = -1;
        c_obliquity[47] = -10;
        c_obliquity[48] = -8;
        c_obliquity[49] = 7;
        c_obliquity[50] = 5;

        delta_lon = 0;
        delta_obliquity = 0;

        for (int j = 1; j <= 106; j++) {
            delta_lon = delta_lon + c_lon[j] / 10000.0 / 3600.0 * Math.sin(xh[j] * Math.PI / 180.0);
            delta_obliquity = delta_obliquity + c_obliquity[j] / 10000.0 / 3600.0 * Math.cos(xh[j] * Math.PI / 180.0);
        }
        // delta_lon：黄经章动，delta_obliquity：交角章动

        double Lon_Nutation = delta_lon * Math.cos(hcjj * Math.PI / 180.0);

        sat = Rotz(sat, -delta_lon * Math.cos(hcjj * Math.PI / 180.0));
        sat = Roty(sat, delta_lon * Math.sin(hcjj * Math.PI / 180.0));
        sat = Rotx(sat, -delta_obliquity);
        return sat;
    }


    // PositionLLA是椭球地固坐标系－测地坐标
    // 函数把空间直角坐标xyz转换为测地坐标经、纬、高度
    private PositionLLA XyzTolla(Position sat) {
        // sat可能为空间位置xyz
        // 0.0066946=1-(b/a)*(b/a)=e*e，b、a为地球半短轴，半长轴，e为地球偏心率
        PositionLLA point = new PositionLLA();
        double x, y, z, lat, lon = 0, h = 0;
        double pi = 3.14159265358979323846d;
        double a = 6378137.0d, lat0 = -100, n, w;
        x = sat.getX();
        y = sat.getY();
        z = sat.getZ();

        lat = Math.atan(z / Math.sqrt(x * x + y * y));

        if ((x >= 0) && (y >= 0)) {
            lon = Math.atan(y / x);
        }
        if ((x < 0) && (y >= 0)) {
            lon = Math.atan(y / x) + pi;
        }
        if ((x < 0) && (y < 0)) {
            lon = Math.atan(y / x) + pi;
        }
        if ((x > 0) && (y < 0)) {
            lon = Math.atan(y / x) + 2 * pi;
        }

        if (lon < 0.0) {
            lon = lon + 2 * pi;
        }
        if (lon > pi) {
            lon = lon - 2 * pi;
        }

        // 用迭代法求椭球大地纬度
        while (Math.abs(lat - lat0) > 0.00000001) {
            lat0 = lat;
            n = a / Math.sqrt(1.0 - 0.0066946 * Math.sin(lat0) * Math.sin(lat0));
            w = Math.sqrt(1.0 - 0.0066946 * Math.sin(lat0) * Math.sin(lat0));
            h = Math.sqrt(x * x + y * y) / Math.cos(lat0) - n;
            lat = Math.atan(z / Math.sqrt(x * x + y * y) * (1.00 + a * 0.0066946 * Math.sin(lat0) / z / w));
        }
        lon = lon * 180.0 / pi;
        lat = lat * 180.0 / pi;
        point.setLon(lon);
        point.setLat(lat);
        point.setHigh(h);
        return point;
    }


    // 函数计算赤经章动，参数：计算时刻，返回值：赤经章动（度）
    private double LonNutation(Calendar newtime, double ms) {
        // double pi=3.14159265358979323846;
        double rlxj, a1, a2, a3, a4, a5, delta_lon, delta_obliquity, hcjj;
        double[] xh = new double[107], c_lon = new double[107], c_obliquity = new double[107];
        for (int i = 0; i < 107; i++) {
            xh[i] = 0;
            c_lon[i] = 0;
            c_obliquity[i] = 0;
        }

        Calendar starttime = TimeUtils.getCalendar(2000, 01, 01, 20, 0, 0);
        long int_time = (newtime.getTimeInMillis() - starttime.getTimeInMillis()) / 1000;
        rlxj = (double) int_time + ms;
        // rlxj从J2000.0起算的儒略世纪数
        rlxj = rlxj / 86400.00 / 36525.00;
        hcjj = 84381.4480 / 3600.0 - 46.81500 / 3600.0 * rlxj - 0.00059 / 3600.0 * Math.pow(rlxj, 2)
                + 0.0018130 / 3600.0 * Math.pow(rlxj, 3);

        a1 = 134.96298139 + (1325 * 360 + 198.8673981) * rlxj + 8.69722e-3 * Math.pow(rlxj, 2)
                + 1.78e-5 * Math.pow(rlxj, 3);
        a2 = 357.5277230 + (99 * 360 + 359.05034) * rlxj - 1.6e-4 * Math.pow(rlxj, 2) - 3.33e-6 * Math.pow(rlxj, 3);
        a3 = 93.27191028 + (1342 * 360 + 82.017538) * rlxj - 0.0036825 * Math.pow(rlxj, 2)
                + 3.06e-6 * Math.pow(rlxj, 3);
        a4 = 297.8503630 + (1236 * 360 + 307.11148) * rlxj - 1.914167e-3 * Math.pow(rlxj, 2)
                + 5.28e-6 * Math.pow(rlxj, 3);
        a5 = 125.0445220 - (5 * 360 + 134.13626) * rlxj + 2.0708e-3 * Math.pow(rlxj, 2) + 2.222e-6 * Math.pow(rlxj, 3);

        xh[1] = a5;
        xh[2] = 2 * a5;
        xh[3] = -2 * a1 + 2 * a3 + a5;
        xh[4] = 2 * a1 - 2 * a3;
        xh[5] = -2 * a1 + 2 * a3 + 2 * a5;
        xh[6] = a1 - a2 - a4;
        xh[7] = -2 * a2 + 2 * a3 - 2 * a4 + a5;
        xh[8] = 2 * a1 - 2 * a3 + a5;
        xh[9] = 2 * a3 - 2 * a4 + 2 * a5;
        xh[10] = a2;
        xh[11] = a2 + 2 * a3 - 2 * a4 + 2 * a5;
        xh[12] = -a2 + 2 * a3 - 2 * a4 + 2 * a5;
        xh[13] = 2 * a3 - 2 * a4 + a5;
        xh[14] = 2 * a1 - 2 * a4;
        xh[15] = 2 * a3 - 2 * a4;
        xh[16] = 2 * a2;
        xh[17] = a2 + a5;
        xh[18] = 2 * a2 + 2 * a3 - 2 * a4 + 2 * a5;
        xh[19] = -a2 + a5;
        xh[20] = -2 * a1 + 2 * a4 + a5;
        xh[21] = -a2 + 2 * a3 - 2 * a4 + a5;
        xh[22] = 2 * a1 - 2 * a4 + a5;
        xh[23] = a2 + 2 * a3 - 2 * a4 + a5;
        xh[24] = a1 - a4;
        xh[25] = 2 * a1 + a2 - 2 * a4;
        xh[26] = -2 * a3 + 2 * a4 + a5;
        xh[27] = a2 - 2 * a3 + 2 * a4;
        xh[28] = a2 + 2 * a5;
        xh[29] = -a1 + a5;
        xh[30] = a2 + 2 * a3 - 2 * a4;
        xh[31] = 2 * a3 + 2 * a5;
        xh[32] = a1;
        xh[33] = 2 * a3 + a5;
        xh[34] = a1 + 2 * a3 + 2 * a5;
        xh[35] = a1 - 2 * a4;
        xh[36] = -a1 + 2 * a3 + 2 * a5;
        xh[37] = 2 * a4;
        xh[38] = a1 + a5;
        xh[39] = -a1 + a5;
        xh[40] = -a1 + 2 * a3 + a5;
        xh[41] = a1 + 2 * a3 + a5;
        xh[42] = 2 * a3 + a5;
        xh[43] = 2 * a1 + a5;
        xh[44] = a1 + 2 * a3 + a5;
        xh[45] = 2 * a1 + 2 * a3 + 2 * a5;
        xh[46] = 2 * a3 + a5;
        xh[47] = -a1 + 2 * a3 + a5;
        xh[48] = -a1 + a5;
        xh[49] = a1 - 2 * a4 + a5;
        xh[50] = -a1 + 2 * a3 + 2 * a4 + a5;

        c_lon[1] = -171996 - 174.2 * rlxj;
        c_lon[2] = 2062 + 0.2 * rlxj;
        c_lon[3] = 46;
        c_lon[4] = 11;
        c_lon[5] = -3;
        c_lon[6] = -3;
        c_lon[7] = -2;
        c_lon[8] = 1;
        c_lon[9] = -13187 - 1.6 * rlxj;
        c_lon[10] = 1426 - 3.4 * rlxj;
        c_lon[11] = -517 + 1.2 * rlxj;
        c_lon[12] = 217 - 0.5 * rlxj;
        c_lon[13] = 129 + 0.1 * rlxj;
        c_lon[14] = 48;
        c_lon[15] = -22;
        c_lon[16] = 17 - 0.1 * rlxj;
        c_lon[17] = -15;
        c_lon[18] = -16 + 0.1 * rlxj;
        c_lon[19] = -12;
        c_lon[20] = -6;
        c_lon[21] = -5;
        c_lon[22] = 4;
        c_lon[23] = 4;
        c_lon[24] = -4;
        c_lon[25] = 1;
        c_lon[26] = 1;
        c_lon[27] = -1;
        c_lon[28] = 1;
        c_lon[29] = 1;
        c_lon[30] = -1;
        c_lon[31] = -2274 - 0.2 * rlxj;
        c_lon[32] = 712 + 0.1 * rlxj;
        c_lon[33] = -386 - 0.4 * rlxj;
        c_lon[34] = -301;
        c_lon[35] = -158;
        c_lon[36] = 123;
        c_lon[37] = 63;
        c_lon[38] = 63 + 0.1 * rlxj;
        c_lon[39] = -58 - 0.1 * rlxj;
        c_lon[40] = -59;
        c_lon[41] = -51;
        c_lon[42] = -38;
        c_lon[43] = 29;
        c_lon[44] = 29;
        c_lon[45] = -31;
        c_lon[46] = 26;
        c_lon[47] = 21;
        c_lon[48] = 16;
        c_lon[49] = -13;
        c_lon[50] = -10;

        c_obliquity[1] = 92025 + 8.9 * rlxj;
        c_obliquity[2] = -895 + 0.5 * rlxj;
        c_obliquity[3] = -24;
        c_obliquity[4] = 0;
        c_obliquity[5] = 1;
        c_obliquity[6] = 0;
        c_obliquity[7] = 1;
        c_obliquity[8] = 0;
        c_obliquity[9] = 5736.00 - 3.1 * rlxj;
        c_obliquity[10] = 54 - 0.1 * rlxj;
        c_obliquity[11] = 224 + 0.6 * rlxj;
        c_obliquity[12] = -95 + 0.3 * rlxj;
        c_obliquity[13] = -70;
        c_obliquity[14] = 1;
        c_obliquity[15] = 0;
        c_obliquity[16] = 0;
        c_obliquity[17] = 9;
        c_obliquity[18] = 7;
        c_obliquity[19] = 6;
        c_obliquity[20] = 3;
        c_obliquity[21] = 3;
        c_obliquity[22] = -2;
        c_obliquity[23] = -2;
        c_obliquity[24] = 0;
        c_obliquity[25] = 0;
        c_obliquity[26] = 0;
        c_obliquity[27] = 0;
        c_obliquity[28] = 0;
        c_obliquity[29] = 0;
        c_obliquity[30] = 0;
        c_obliquity[31] = 977 - 0.5 * rlxj;
        c_obliquity[32] = -7;
        c_obliquity[33] = 200;
        c_obliquity[34] = 129 - 0.1 * rlxj;
        c_obliquity[35] = -1;
        c_obliquity[36] = -53;
        c_obliquity[37] = -2;
        c_obliquity[38] = -33;
        c_obliquity[39] = 32;
        c_obliquity[40] = 26;
        c_obliquity[41] = 27;
        c_obliquity[42] = 16;
        c_obliquity[43] = -1;
        c_obliquity[44] = -12;
        c_obliquity[45] = 13;
        c_obliquity[46] = -1;
        c_obliquity[47] = -10;
        c_obliquity[48] = -8;
        c_obliquity[49] = 7;
        c_obliquity[50] = 5;

        delta_lon = 0;
        delta_obliquity = 0;

        for (int j = 1; j <= 106; j++) {
            delta_lon = delta_lon + c_lon[j] / 10000.0 / 3600.0 * Math.sin(xh[j] * Math.PI / 180.0);
            delta_obliquity = delta_obliquity + c_obliquity[j] / 10000.0 / 3600.0 * Math.cos(xh[j] * Math.PI / 180.0);
        }

        double Lon_Nutation = delta_lon * Math.cos(hcjj * Math.PI / 180.0);
        // 返回值：赤经章动,单位度
        return Lon_Nutation;
    }


    public List<StripDTO> accessWindowCombine(List<StripDTO> striplist) {
        List<StripDTO> working_accessDatas = new LinkedList<>();
        working_accessDatas.addAll(striplist);
        List<StripDTO> combined_result_List = new LinkedList<>();
        while (!working_accessDatas.isEmpty()) {
            StripDTO current = working_accessDatas.get(0);
            StripDTO target = null;
            boolean combineflag = false;
            for (int i = 1; i < working_accessDatas.size(); i++) {
                // 如果有相交，进行合并，TBD:或其他合并条件
                target = working_accessDatas.get(i);
                // 可以检查如果需要合并，合并后的时长是否超出了最大连续侦察时长。防止无限合并
                //  侧摆角相同的两个条带才进行合并
                Date currentStarttime = combineTwoRadarStrip(current, target).getObservationStartTime();
                Date currentEndtime = combineTwoRadarStrip(current, target).getObservationEndTime();
                long windowSize = currentEndtime.getTime() - currentStarttime.getTime();
                windowSize /= 1000;// to seconds
                if (windowSize <= Long.valueOf(maxWindowLimit) && Objects.equals(current.getRollAngle(), target.getRollAngle())) {
                    combineflag = true;
                    break;
                }
            }
            // 如果经过比较,发现不需要合并，说明该时间段已经为最长时间段，从工作列表中移至结果列表中。
            if (!combineflag) {
                working_accessDatas.remove(current);
                combined_result_List.add(current);

            } else {
                // 如果发现需要进行合并，将发现目标窗口与工作窗口合并，并将合并后的新窗口放在列表第一个，等待再次访问
                StripDTO combined_single = combineTwoRadarStrip(current, target);
                working_accessDatas.remove(current);
                working_accessDatas.remove(target);
                working_accessDatas.add(0, combined_single);

            }

        }
        return combined_result_List;
    }



    private StripDTO combineTwoRadarStrip(StripDTO current, StripDTO target) {
        Date ealiestBegin, latestEnd;
        ealiestBegin = current.getObservationStartTime().before(target.getObservationStartTime()) ? current.getObservationStartTime() : target.getObservationStartTime();
        latestEnd = current.getObservationEndTime().after(target.getObservationEndTime()) ? current.getObservationEndTime() : target.getObservationEndTime();
        current.setObservationStartTime(ealiestBegin);
        current.setObservationEndTime(latestEnd);
        return current;
    }


    /**
     * 计算条带在某侧摆值下星下点轨迹
     * @param rollangle
     * 侧摆值
     * @param ephtrack
     * 原轨迹
     * @param maxswathwidth
     * 幅宽
     * @param maxoverlapwidth
     * 重叠区域
     * @return
     * @throws CloneNotSupportedException
     */
    public List<SubPoint> calSubPoint(double rollangle, EphTrackCoreDTO ephtrack, double maxswathwidth,
                                      double maxoverlapwidth) throws Exception {
        List<SubPoint> temp = new ArrayList<>();
        // side为侧摆方向，side为0时，表示向右侧摆，side为1时表示向左侧摆
        int side;

        if (rollangle > 0) {
            side = 0;
        } else {
            side = 1;
        }
        double roll = Math.abs(rollangle);
        int duration_lmt = (int) ephtrack.getPayloadLMT().getMin_duration()+1;
        if(ephtrack.getSubPointlist().size() >= duration_lmt){
            for (int i = 0; i < ephtrack.getSubPointlist().size() - 1; i++) {
                com.base.core.source.LatLon latLon = com.base.core.source.LatLon.fromDegrees(ephtrack.getSubPointlist().get(i).getLat(),
                        ephtrack.getSubPointlist().get(i).getLon());
                com.base.core.source.Position pos = new com.base.core.source.Position(latLon, ephtrack.getSubPointlist().get(i).getAlt());
                com.base.core.source.LatLon latLon1 = com.base.core.source.LatLon.fromDegrees(ephtrack.getSubPointlist().get(i + 1).getLat(),
                        ephtrack.getSubPointlist().get(i + 1).getLon());// 下一点星下点数据
                com.base.core.source.Position nextPos = new com.base.core.source.Position(latLon1, ephtrack.getSubPointlist().get(i + 1).getAlt());
                SideAngleSurfaceGrid fake = new SideAngleSurfaceGrid(pos, nextPos, maxswathwidth/2, maxoverlapwidth, roll,
                        side);
                SubPoint tp=new SubPoint();
                /* tp = ephtrack.getSubPointlist().get(i);*/
                String template=JSON.toJSONString(ephtrack.getSubPointlist().get(i));
                tp = JSON.parseObject(template, SubPoint.class);
                if(ephtrack.getUpdown().equalsIgnoreCase("1")){
                    tp.setLsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setRsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointThr().getLongitude().getDegrees());
                }else{
                    tp.setRsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setLsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointThr().getLongitude().getDegrees());
                }

                temp.add(tp);
            }
            //是否满足观测需求的时长要求
        } else {
            int count = ephtrack.getSubPointlist().size() / 2;
            Date center_time = ephtrack.getSubPointlist().get(count).getDate();
            Date starttime = new Date(center_time.getTime() - duration_lmt * 500);
            Date enttime = new Date(center_time.getTime() + duration_lmt * 500);
            List<SubPoint> subPoints =new ArrayList<>();
            if(!ephtrack.getPayloadLMT().getMaxInfoBasisPayload().getPayloadType().equals("sar")){
                subPoints= bizService.getBizSubPointList(ephtrack.getSatallitCode(), starttime, enttime,"1");
            }else {
                subPoints = bizService.getBizSubPointList(ephtrack.getSatallitCode(), starttime, enttime);
            }
            /*List<SubPoint> subPoints = bizService.getBizSubPointList(ephtrack.getSatallitCode(),starttime,enttime);*/
            if (subPoints == null || subPoints.isEmpty()) {
                throw new Exception("no suitable eph data for missionstrip!");
            }
            for (int i = 0; i < subPoints.size() - 1; i++) {
                SubPoint tp = new SubPoint();
                SubPoint tpNext = new SubPoint();
                String tpstr=JSON.toJSONString(subPoints.get(i));
                tp = JSON.parseObject(tpstr, SubPoint.class);

                String tpNextstr=JSON.toJSONString(subPoints.get(i+1));
                tpNext = JSON.parseObject(tpNextstr, SubPoint.class);

                com.base.core.source.LatLon latLon =  com.base.core.source.LatLon.fromDegrees(tp.getLat().floatValue(), tp.getLon().floatValue());
                com.base.core.source.Position pos = new  com.base.core.source.Position(latLon, tp.getAlt().floatValue());
                com.base.core.source.LatLon latLon1 =  com.base.core.source.LatLon.fromDegrees(tpNext.getLat().floatValue(), tpNext.getLon().floatValue());// 下一点星下点数据
                com.base.core.source.Position nextPos = new  com.base.core.source.Position(latLon1, tpNext.getAlt().floatValue());
                SideAngleSurfaceGrid fake = new SideAngleSurfaceGrid(pos, nextPos, maxswathwidth/2, maxoverlapwidth, roll,
                        side);

                if(ephtrack.getUpdown().equalsIgnoreCase("1")){
                    tp.setLsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setRsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointThr().getLongitude().getDegrees());
                }else{
                    tp.setRsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setLsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointThr().getLongitude().getDegrees());
                }

                temp.add(tp);

            }

        }
        return temp;
    }


    /**
     * 计算网格弧段最接近的轨迹位置
     */
    public Map<Integer,Integer> calcSublistPosition(List<SubPoint> TrackPointlist,Region region){
        Map<Integer,Integer> regionvalue = new HashMap<>();

        int size = TrackPointlist.size();
        double max_lat = -90;
        double min_lat = 90;

        double max_lon = -180;
        double min_lon = 180;
        double temp = 0;
        String updown = TrackPointlist.get(0).getUpdown();
        String updown_end = TrackPointlist.get(size - 1).getUpdown();
//		float distance = 0;
        /**
         * 降轨时左上右下（左点幅宽纬度大于右点幅宽纬度）
         升轨时左下右上（左点幅宽纬度小于右点幅宽纬度）
         */
        if(updown.equals(updown_end)){
            if(updown.equalsIgnoreCase("1")){
                //降轨时LsLat最大
                max_lat = TrackPointlist.get(0).getLat();
                min_lat = TrackPointlist.get(size - 1).getLat();
                max_lon = TrackPointlist.get(0).getRsLon();
                min_lon = TrackPointlist.get(size - 1).getLsLon();
//				distance = TrackPointlist.get(0).getLat() - TrackPointlist.get(size - 1).getLat();
//				min_sublat = TrackPointlist.get(size - 1).getLat();
            }else{
                //升轨时RsLat最大
                max_lat = TrackPointlist.get(size - 1).getLat();
                min_lat = TrackPointlist.get(0).getLat();
                max_lon = TrackPointlist.get(size - 1).getRsLon();
                min_lon = TrackPointlist.get(0).getLsLon();
//				distance = TrackPointlist.get(size - 1).getLat() - TrackPointlist.get(0).getLat();
//				min_sublat = TrackPointlist.get(0).getLat();
            }
            for(int i = 0;i < region.getOriginArea().size();i++){
                float lb_lat = (float) region.getOriginArea().get(i).getGridlist().get(0).getLbLat();
                float rt_lat = (float) region.getOriginArea().get(i).getGridlist().get(0).getRtLat();
                float c_lat =  (float) region.getOriginArea().get(i).getGridlist().get(0).getCLat();
                double dis = 1000;
                int count = 0;

                if(min_lat < lb_lat && rt_lat < max_lat){
                    for (int j = 0; j < TrackPointlist.size(); j++) {
                        temp = Math.abs(TrackPointlist.get(j).getLat() - c_lat);
                        if (temp < dis) {
                            dis = temp;
                            count = j;
                        }
                    }
                    regionvalue.put(i, count);

                }
            }
        }
        return regionvalue;

    }

    /**
     * 计算Sar星点目标星下点
     * @param rollAngle
     * @param accessTime
     * @param ephTrack
     * @param maxSwathWidth
     * @param maxOverLapWidth
     * @param region
     * @return
     * @throws Exception
     */
    public List<SubPoint> calSubPointForSarPoint(double rollAngle,Date accessTime, EphTrackCoreDTO ephTrack, double maxSwathWidth,
                                                 double maxOverLapWidth, Region region ) throws Exception {
        List<SubPoint> temp = new ArrayList<SubPoint>();
        // side为侧摆方向，side为0时，表示向右侧摆，side为1时表示向左侧摆
        int side;
        /*
         * 降轨时，侧摆大于0表示地球坐标系向左摆
         */
        /*
         * if(ephtrack.getUpdown().equalsIgnoreCase("1")){ if (rollangle > 0) {
         * side = 1; } else { side = 0; } }else{ if (rollangle > 0) { side = 0;
         * } else { side = 1; } }
         */
        if (rollAngle > 0) {
            side = 0;
        } else {
            side = 1;
        }
        double roll = Math.abs(rollAngle);
        int duration_lmt = (int) ephTrack.getPayloadLMT().getMin_duration() + 1;

        // 实传模式（认为是侧摆角为0，不对再改
        if (ephTrack.getSubPointlist().size() >= duration_lmt) {
            for (int i = 0; i < ephTrack.getSubPointlist().size() - 1; i++) {
                com.base.core.source.LatLon latLon = com.base.core.source.LatLon.fromDegrees(ephTrack.getSubPointlist().get(i).getLat(),
                        ephTrack.getSubPointlist().get(i).getLon());
                com.base.core.source.Position pos = new com.base.core.source.Position(latLon, ephTrack.getSubPointlist().get(i).getAlt());
                com.base.core.source.LatLon latLon1 = com.base.core.source.LatLon.fromDegrees(ephTrack.getSubPointlist().get(i + 1).getLat(),
                        ephTrack.getSubPointlist().get(i + 1).getLon());// 下一点星下点数据
                com.base.core.source.Position nextPos = new com.base.core.source.Position(latLon1, ephTrack.getSubPointlist().get(i + 1).getAlt());
                SideAngleSurfaceGrid fake = new SideAngleSurfaceGrid(pos, nextPos, maxSwathWidth/2, maxOverLapWidth, roll,
                        side);
                SubPoint tp=new SubPoint();
                /* tp = ephtrack.getSubPointlist().get(i);*/
                String template=JSON.toJSONString(ephTrack.getSubPointlist().get(i));
                tp = JSON.parseObject(template, SubPoint.class);
                if(ephTrack.getUpdown().equalsIgnoreCase("1")){
                    tp.setLsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setRsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointThr().getLongitude().getDegrees());
                }else{
                    tp.setRsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setLsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointThr().getLongitude().getDegrees());
                }

                temp.add(tp);
            }

        } else {
            Date centerTime = accessTime ;
            Date starttime = new Date(centerTime.getTime() - duration_lmt * 500);
            Date enttime = new Date(centerTime.getTime() + duration_lmt * 500);
            List<SubPoint> subPoints =new ArrayList<>();
            if(!ephTrack.getPayloadLMT().getMaxInfoBasisPayload().getPayloadType().equals("sar")){
                subPoints= bizService.getBizSubPointList(ephTrack.getSatallitCode(), starttime, enttime,"1");
            }else {
                subPoints = bizService.getBizSubPointList(ephTrack.getSatallitCode(), starttime, enttime);
            }
            /*List<SubPoint> subPoints = bizService.getBizSubPointList(ephTrack.getSatallitCode(),starttime,enttime);*/
            if (subPoints == null || subPoints.isEmpty()) {
                throw new Exception("no suitable eph data for missionstrip!");
            }
            for (int i = 0; i < subPoints.size() - 1; i++) {
                SubPoint tp = new SubPoint();
                SubPoint tpNext = new SubPoint();
                String tpstr=JSON.toJSONString(subPoints.get(i));
                tp = JSON.parseObject(tpstr, SubPoint.class);

                String tpNextstr=JSON.toJSONString(subPoints.get(i+1));
                tpNext = JSON.parseObject(tpNextstr, SubPoint.class);

                com.base.core.source.LatLon latLon =  com.base.core.source.LatLon.fromDegrees(tp.getLat().floatValue(), tp.getLon().floatValue());
                com.base.core.source.Position pos = new  com.base.core.source.Position(latLon, tp.getAlt().floatValue());
                com.base.core.source.LatLon latLon1 =  com.base.core.source.LatLon.fromDegrees(tpNext.getLat().floatValue(), tpNext.getLon().floatValue());// 下一点星下点数据
                com.base.core.source.Position nextPos = new  com.base.core.source.Position(latLon1, tpNext.getAlt().floatValue());
                SideAngleSurfaceGrid fake = new SideAngleSurfaceGrid(pos, nextPos, maxSwathWidth/2, maxOverLapWidth, roll, side);

                if(ephTrack.getUpdown().equalsIgnoreCase("1")){
                    tp.setLsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setRsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointThr().getLongitude().getDegrees());
                }else{
                    tp.setRsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setLsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointThr().getLongitude().getDegrees());
                }
                temp.add(tp);
            }
        }
        return temp;
    }

    public List<SubPoint> calSubPointForCcdPoint(double rollangle, EphTrackCoreDTO ephTrackCoreDTO, float maxswathwidth, float maxoverlapwidth, Region region) {
        List<SubPoint> temp = new ArrayList<SubPoint>();
        // side为侧摆方向，side为0时，表示向右侧摆，side为1时表示向左侧摆
        int side;
        /*
         * 降轨时，侧摆大于0表示地球坐标系向左摆
         */
        /*
         * if(ephtrack.getUpdown().equalsIgnoreCase("1")){ if (rollangle > 0) {
         * side = 1; } else { side = 0; } }else{ if (rollangle > 0) { side = 0;
         * } else { side = 1; } }
         */
        if (rollangle > 0) {
            side = 0;
        } else {
            side = 1;
        }
        double roll = Math.abs(rollangle);
        int duration_lmt = (int) ephTrackCoreDTO.getPayloadLMT().getMin_duration() + 1;

        // 是否满足观测需求的时长要求 && 实传模式（认为是侧摆角为0，不对再改）---zhanghy
        /*BigDecimal isRealtimeTrans = region.getCallRequireVisibility().getIsRealtimeTrans();*/
        BigDecimal isRealtimeTrans = BigDecimal.valueOf(0);
        if ((ephTrackCoreDTO.getSubPointlist().size() >= duration_lmt)
                || (isRealtimeTrans !=null && isRealtimeTrans.intValue()==1)) {
            for (int i = 0; i < ephTrackCoreDTO.getSubPointlist().size() - 1; i++) {
                com.base.core.source.LatLon latLon = com.base.core.source.LatLon.fromDegrees(ephTrackCoreDTO.getSubPointlist().get(i).getLat(),
                        ephTrackCoreDTO.getSubPointlist().get(i).getLon());
                com.base.core.source.Position pos = new com.base.core.source.Position(latLon, ephTrackCoreDTO.getSubPointlist().get(i).getAlt());
                com.base.core.source.LatLon latLon1 = com.base.core.source.LatLon.fromDegrees(ephTrackCoreDTO.getSubPointlist().get(i + 1).getLat(),
                        ephTrackCoreDTO.getSubPointlist().get(i + 1).getLon());// 下一点星下点数据
                com.base.core.source.Position nextPos = new com.base.core.source.Position(latLon1, ephTrackCoreDTO.getSubPointlist().get(i + 1).getAlt());
                SideAngleSurfaceGrid fake = new SideAngleSurfaceGrid(pos, nextPos, maxswathwidth / 2, maxoverlapwidth,
                        roll, side);
                SubPoint tp = new SubPoint();
                tp = ephTrackCoreDTO.getSubPointlist().get(i);
                if (ephTrackCoreDTO.getUpdown().equalsIgnoreCase("1")) {
                    tp.setLsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setRsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointThr().getLongitude().getDegrees());
                } else {
                    tp.setRsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setLsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointThr().getLongitude().getDegrees());
                }

                temp.add(tp);
            }

        } else {

            int count = ephTrackCoreDTO.getSubPointlist().size() / 2;
            Date centerTime = ephTrackCoreDTO.getSubPointlist().get(count).getDate();
            Date starttime = new Date(centerTime.getTime() - duration_lmt * 500);
            Date enttime = new Date(centerTime.getTime() + duration_lmt * 500);
            String firsttime = TimeFormatTool.dateToStringDefault(starttime);
            String lasttime = TimeFormatTool.dateToStringDefault(enttime);
            List<SubPoint> dbdata = new ArrayList<>();
            if(!ephTrackCoreDTO.getPayloadLMT().getMaxInfoBasisPayload().getPayloadType().equals("sar")){
                dbdata= bizService.getBizSubPointList(ephTrackCoreDTO.getSatallitCode(), starttime, enttime,"1");
            }else {
                dbdata = bizService.getBizSubPointList(ephTrackCoreDTO.getSatallitCode(), starttime, enttime);
            }
            /*List<SubPoint> dbdata = bizService.getBizSubPointList(ephTrackCoreDTO.getSatallitCode(), starttime, enttime);*/

            for (int i = 0; i < dbdata.size() - 1; i++) {
                SubPoint tp = new SubPoint();
                tp.setId(String.valueOf(dbdata.get(i).getId()));
                tp.setUpdown(String.valueOf(dbdata.get(i).getUpdown()));
                tp.setOrbit(String.valueOf(dbdata.get(i).getOrbit()));
                tp.setSatelliteId(String.valueOf(dbdata.get(i).getSatelliteId()));
                tp.setDate(dbdata.get(i).getDate());

                Double alt = dbdata.get(i).getAlt();
                Double lat = dbdata.get(i).getLat();
                Double lon = dbdata.get(i).getLon();
                tp.setAlt(alt);
                tp.setLat(lat);
                tp.setLon(lon);
                Double nextalt = dbdata.get(i + 1).getAlt();
                Double nextlat = dbdata.get(i + 1).getLat();
                Double nextlon = dbdata.get(i + 1).getLon();

                com.base.core.source.LatLon latLon = com.base.core.source.LatLon.fromDegrees(lat.floatValue(), lon.floatValue());
                com.base.core.source.Position pos = new com.base.core.source.Position(latLon, alt.floatValue());
                com.base.core.source.LatLon latLon1 = com.base.core.source.LatLon.fromDegrees(nextlat.floatValue(), nextlon.floatValue());// 下一点星下点数据
                com.base.core.source.Position nextPos = new com.base.core.source.Position(latLon1, nextalt.floatValue());
                SideAngleSurfaceGrid fake = new SideAngleSurfaceGrid(pos, nextPos, maxswathwidth / 2, maxoverlapwidth,
                        roll, side);

                if (ephTrackCoreDTO.getUpdown().equalsIgnoreCase("1")) {
                    tp.setLsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setLsLon(fake.getPointSec().getLongitude().getDegrees());
                    tp.setRsLat(fake.getPointThr().getLatitude().getDegrees());
                    tp.setRsLon(fake.getPointThr().getLongitude().getDegrees());
                } else {
                    tp.setRsLat(fake.getPointSec().getLatitude().getDegrees());
                    tp.setRsLon( fake.getPointSec().getLongitude().getDegrees());
                    tp.setLsLat( fake.getPointThr().getLatitude().getDegrees());
                    tp.setLsLon( fake.getPointThr().getLongitude().getDegrees());
                }

                temp.add(tp);

            }

        }
        return temp;
    }
}
