package com.base.core.utiles;


import com.base.core.model.enity.*;

import java.util.Calendar;

public class GeoUtils {

    private double[] convertLLAToECI(EphData ephData, double currentJulianDate){
        double[] lla_deg_m_sat = new double[]{ephData.getLat(),ephData.getLon(),ephData.getHigh()};
        double thetaDeg = Sidereal.Mean_Sidereal_Deg(currentJulianDate-AstroConst.JDminusMJD, lla_deg_m_sat[1]);
        return GeoFunctions.calculateECIposition(lla_deg_m_sat,thetaDeg);
    }


    //卫星对目标的瞬时视场角（访问侧摆角），参数sat：卫星空间位置测地坐标，target：目标点测地坐标
    public static double accessLLA_new(PositionLLA sat, SphericalPoint target) {

        //sat：卫星空间位置测地坐标，target：目标点测地坐标
        //Sa、S：星下点地心坐标，T：目标点地心坐标
        double pi=Math.PI;
        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 a=SphericalAngle(S,T);

        //OT：目标到地心距离，OS：星下点到地心距离
        double OT=EarthRadiusA(T.getLat());
        double OS=Sa.getR();
        double a0=pi/2-Math.asin(OT/OS);
        double f = Math.atan(Math.cos(a0)*Math.sin(a)/(1-Math.cos(a0)*Math.cos(a)));


        double ST=Math.sqrt( Math.pow(OT,2)+Math.pow(OS,2)-2*OT*OS*Math.cos(a*pi/180.0) ) ;
        System.out.println("\n\nlla f:"+Math.asin( OT/ST*Math.sin(a*pi/180.0) )*180.0/pi);
        System.out.println("lla_new f:"+f*180/pi);
        return f*180/pi;
    }

    //卫星对目标的瞬时视场角（访问侧摆角），参数sat：卫星空间位置测地坐标，target：目标点测地坐标
    public static double accessLLA_EZR(String timeString,PositionLLA sat, SphericalPoint target) {
        Time time1_julian = new Time();
        time1_julian.set(TimeUtils.try2ParseTime(timeString).getTimeInMillis());
        double julian_time1 = time1_julian.getJulianDate();
        //常规EPH为J2000 TEME
        double[] lla_pos = new double[]{sat.getLat(),sat.getLon(),sat.getHigh()};
        double[] eci_pos = GeoFunctions.convertLLAToECI(lla_pos,julian_time1);

        //sat.calculateTemePositionFromUT(time1);
        //计算地面站对于卫星在指定时刻的高度角
        PositionLLA positionLLA_target = new PositionLLA();
        positionLLA_target.setLat(target.getLat());
        positionLLA_target.setLon(target.getLon());
        positionLLA_target.setHigh(0);
        double h1 = 90-GeoFunctions.calculate_AER(julian_time1,positionLLA_target.getLla_deg_m(),eci_pos)[1];
        return h1;

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

        //sat：卫星空间位置测地坐标，target：目标点测地坐标
        //Sa、S：星下点地心坐标，T：目标点地心坐标
        double pi=Math.PI;
        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();
        //
        double ST=Math.sqrt( Math.pow(OT,2)+Math.pow(OS,2)-2*OT*OS*Math.cos(coreangle*pi/180.0) ) ;
//        System.out.println("\n\nlla f:"+Math.asin( OT/ST*Math.sin(coreangle*pi/180.0) )*180.0/pi);

        return Math.asin( OT/ST*Math.sin(coreangle*pi/180.0) )*180.0/pi;
    }

    //函数用测地纬度求地表上点的地心纬度，返回值地心纬度，参数lat测地纬度，见《球面天文学》
    public static  double GeodesicToCore(double lat) {
        //a为地球赤道半径，f为地球扁率
        double a= AstroConst.R_Earth,f=AstroConst.f_Earth;
        double pi=Math.PI;
        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;
    }

    //此函数用于测地坐标向地心坐标的转换，参数sat为测地坐标，见《球面天文学》37页图2.7
    public static  PositionLLR LLAToLLR(PositionLLA sat) {
        double a= AstroConst.R_Earth,f=AstroConst.f_Earth,b;
        double pi=Math.PI;
        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：测地纬度，函数返回值地心距
    public static  double EarthRadiusB(double lat) {
        //a为地球赤道半径，f为地球扁率
        double a= AstroConst.R_Earth,f=AstroConst.f_Earth;
        double pi= Math.PI;
        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);
        return Math.sqrt( Math.pow(C*Math.cos(lat),2)+Math.pow(S*Math.sin(lat),2) );
    }
    public static  double SphericalAngle(SphericalPoint source_arg, SphericalPoint destination_arg) {
        double  back;

        SphericalPoint source = new SphericalPoint();
        SphericalPoint destination = new SphericalPoint();
        source.setLon(source_arg.getLon() * Math.PI / 180.0);
        source.setLat(source_arg.getLat() * Math.PI / 180.0);
        destination.setLon(destination_arg.getLon() * Math.PI / 180.0);
        destination.setLat(destination_arg.getLat() * Math.PI / 180.0);

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

        return Math.acos(back) * 180.0 / Math.PI;
    }
    //?公式未知，此函数用地心纬度求地表不同纬度圈地心距，参数lat：地心纬度，函数返回值地心距
    public static  double EarthRadiusA(double lat) {
        double a= AstroConst.R_Earth,f=AstroConst.f_Earth,b;
        double pi=Math.PI;
        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);
    }


    //PositionLLA是椭球地固坐标系－测地坐标
    //函数把空间直角坐标xyz转换为测地坐标经、纬、高度
    public static 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;
    }

    //函数计算赤经章动，参数：计算时刻，返回值：赤经章动（度）
    public static 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]*pi/180.0);
            delta_obliquity=delta_obliquity + c_obliquity[j]/10000.0/3600.0 * Math.cos(xh[j]*pi/180.0);
        }

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

    //函数把瞬时平赤道地心系坐标经章动矩阵转换为瞬时真赤道地心系坐标
    //参数sat：瞬时平赤道地心系坐标，返回值sat：瞬时真赤道地心系坐标
    public static 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]*pi/180.0);
            delta_obliquity=delta_obliquity + c_obliquity[j]/10000.0/3600.0 * Math.cos(xh[j]*pi/180.0);
        }
        //delta_lon：黄经章动，delta_obliquity：交角章动

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

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

    //函数把历元平赤道地心系坐标经岁差矩阵转换为瞬时平赤道地心系坐标
    //参数sat：历元平赤道地心系坐标，返回值sat：瞬时平赤道地心系坐标
    public static 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;
    }

    //此函数把轨道根数转化为瞬时真赤道地心系（与J2000.0惯性系差岁差、章动）空间位置、速度
    public static 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;
    }

    //Delta_2为J2、J3、J4摄动的2阶长期项摄动
    public static 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;
    }

    //Delta_1为J2、J3、J4摄动的0阶和1阶长期项摄动
    public static 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;
    }

    public static Position Con_Precession(Calendar ephTime, Position sat) {
        double d1,d2,d3;
        Calendar startTime = TimeUtils.getCalendar((int)2000, (int)1, (int)1, (int)20, (int)0, (int)0);
        long int_time=(ephTime.getTimeInMillis() - startTime.getTimeInMillis()) / 1000;
        //rlxj为星历历元时间到J2000.0的儒略世纪数
        double rlxj=int_time/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,d3);
        sat=Roty(sat,-1.0*d2);
        sat=Rotz(sat,d1);
        //	sat=Rotz(sat,-d1);
        //	sat=Roty(sat,d2);
        //	sat=Rotz(sat,-d3);
        return sat;
    }

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

        while(f>2*pi) {f=f-2*pi;}
        while(f<0)    {f=f+2*pi;}

        return f*180.0/pi;
    }

    public static double Mtoe(double e, double M) {
        double  me,me0;
        double pi=3.14159265358979323846d;
        M=M*pi/180.0;
        me0=0d;
        me=M;
//		while(
//				Math.abs(me-me0) > 0.00000000040
//			)
        while(
                Math.abs(me-me0) > 0.00000000001
                )
        {
            me0=me;
            me=M + e*Math.sin(me);
        }
        while(me>2*pi) { me=me-2*pi;}
        while(me<0)    {me=me+2*pi;}

        return me*180.0/pi;
    }

    public static Position Con_Nutation(Calendar ephTime, Position sat) {
        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]=0d;c_lon[i]=0d;c_obliquity[i]=0d;}
        Calendar startTime = TimeUtils.getCalendar((int)2000, (int)1, (int)1, (int)20, (int)0, (int)0);
        long int_time=(ephTime.getTimeInMillis() - startTime.getTimeInMillis()) / 1000;
        //rlxj为星历历元时间到J2000.0的儒略世纪数
        rlxj=int_time/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]=	46d;
        c_lon[4]= 11d;
        c_lon[5]= -3d;
        c_lon[6]= -3d;
        c_lon[7]= -2d;
        c_lon[8]= 1d;
        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]= 48d;
        c_lon[15]= -22d;
        c_lon[16]= 17 - 0.1*rlxj;
        c_lon[17]= -15d;
        c_lon[18]= -16 +0.1*rlxj;
        c_lon[19]= -12d;
        c_lon[20]= -6d;
        c_lon[21]=-5d;
        c_lon[22]=4d;
        c_lon[23]=4d;
        c_lon[24]=-4d;
        c_lon[25]=1d;
        c_lon[26]=1d;
        c_lon[27]=-1d;
        c_lon[28]=1d;
        c_lon[29]=1d;
        c_lon[30]=-1d;
        c_lon[31]=-2274 -0.2*rlxj;
        c_lon[32]=712 +0.1*rlxj;
        c_lon[33]=-386 -0.4*rlxj;
        c_lon[34]=-301d;
        c_lon[35]=-158d;
        c_lon[36]=123d;
        c_lon[37]=63d;
        c_lon[38]=63 +0.1*rlxj;
        c_lon[39]=-58 -0.1*rlxj;
        c_lon[40]=-59d;
        c_lon[41]=-51d;
        c_lon[42]=-38d;
        c_lon[43]=29d;
        c_lon[44]=29d;
        c_lon[45]=-31d;
        c_lon[46]=26d;
        c_lon[47]=21d;
        c_lon[48]=16d;
        c_lon[49]=-13d;
        c_lon[50]=-10d;



        c_obliquity[1]= 92025 + 8.9*rlxj;
        c_obliquity[2]= -895 + 0.5*rlxj;
        c_obliquity[3]= -24d;
        c_obliquity[4]= 0d;
        c_obliquity[5]= 1d;
        c_obliquity[6]= 0d;
        c_obliquity[7]= 1d;
        c_obliquity[8]= 0d;
        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]= -70d;
        c_obliquity[14]= 	1d;
        c_obliquity[15]= 	0d;
        c_obliquity[16]= 0d;
        c_obliquity[17]= 9d;
        c_obliquity[18]= 7d;
        c_obliquity[19]= 6d;
        c_obliquity[20]=3d;
        c_obliquity[21]= 3d;
        c_obliquity[22]= -2d;
        c_obliquity[23]= -2d;
        c_obliquity[24]= 0d;
        c_obliquity[25]= 0d;
        c_obliquity[26]= 0d;
        c_obliquity[27]= 0d;
        c_obliquity[28]= 0d;
        c_obliquity[29]= 0d;
        c_obliquity[30]= 0d;
        c_obliquity[31]= 977 -0.5*rlxj;
        c_obliquity[32]= -7d;
        c_obliquity[33]= 200d;
        c_obliquity[34]= 129 -0.1*rlxj;
        c_obliquity[35]= -1d;
        c_obliquity[36]= -53d;
        c_obliquity[37]= 	-2d;
        c_obliquity[38]= 	-33d;
        c_obliquity[39]= 	32d;
        c_obliquity[40]= 	26d;

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

        delta_lon=0;
        delta_obliquity=0;


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

        //delta_lon：黄经章动，delta_obliquity：交角章动，hcjj黄赤交角
        //Lon_Nutation=delta_lon*Math.cos(hcjj*pi/180.0);

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

    private static Position Roty(Position sat, double sita) {
        double x0,y0,z0;
        double pi=3.14159265358979323846;
        sita=sita*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 ;
    }

    public static Position Rotx(Position sat, double sita) {
        Position sat0 = new Position();
        double x0,y0,z0;;
        double pi=3.14159265358979323846;
        sita=sita*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;
    }

    public static Position Rotz(Position sat, double sita) {
        double x0,y0,z0;
        Position sat0 = new Position();
        double pi=3.14159265358979323846;
        sita=sita*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;
    }

    public static double LonNutation(Calendar ephTime) {
        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]=0d;c_lon[i]=0d;c_obliquity[i]=0d;}
        Calendar startTime = TimeUtils.getCalendar(2000, 1, 1, 20, 0, 0);
        long int_time=(ephTime.getTimeInMillis() - startTime.getTimeInMillis()) / 1000;
        //rlxj从J2000.0起算的儒略世纪数
        rlxj=int_time/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]=	46d;
        c_lon[4]= 11d;
        c_lon[5]= -3d;
        c_lon[6]= -3d;
        c_lon[7]= -2d;
        c_lon[8]= 1d;
        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]= 48d;
        c_lon[15]= -22d;
        c_lon[16]= 17 - 0.1*rlxj;
        c_lon[17]= -15d ;
        c_lon[18]= -16 +0.1*rlxj;
        c_lon[19]= -12d;
        c_lon[20]= -6d;
        c_lon[21]=-5d;
        c_lon[22]=4d;
        c_lon[23]=4d;
        c_lon[24]=-4d;
        c_lon[25]=1d;
        c_lon[26]=1d;
        c_lon[27]=-1d;
        c_lon[28]=1d;
        c_lon[29]=1d;
        c_lon[30]=-1d;
        c_lon[31]=-2274 -0.2*rlxj;
        c_lon[32]=712 +0.1*rlxj;
        c_lon[33]=-386 -0.4*rlxj;
        c_lon[34]=-301d;
        c_lon[35]=-158d;
        c_lon[36]=123d;
        c_lon[37]=63d;
        c_lon[38]=63 +0.1*rlxj;
        c_lon[39]=-58 -0.1*rlxj;
        c_lon[40]=-59d;
        c_lon[41]=-51d;
        c_lon[42]=-38d;
        c_lon[43]=29d;
        c_lon[44]=29d;
        c_lon[45]=-31d;
        c_lon[46]=26d;
        c_lon[47]=21d;
        c_lon[48]=16d;
        c_lon[49]=-13d;
        c_lon[50]=-10d;



        c_obliquity[1]= 92025 + 8.9*rlxj;
        c_obliquity[2]= -895 + 0.5*rlxj;
        c_obliquity[3]= -24d;
        c_obliquity[4]= 0d;
        c_obliquity[5]= 1d;
        c_obliquity[6]= 0d;
        c_obliquity[7]= 1d;
        c_obliquity[8]= 0d;
        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]= -70d;
        c_obliquity[14]= 	1d;
        c_obliquity[15]= 	0d;
        c_obliquity[16]= 0d;
        c_obliquity[17]= 9d;
        c_obliquity[18]= 7d;
        c_obliquity[19]= 6d;
        c_obliquity[20]=3d;
        c_obliquity[21]= 3d;
        c_obliquity[22]= -2d;
        c_obliquity[23]= -2d;
        c_obliquity[24]= 0d;
        c_obliquity[25]= 0d;
        c_obliquity[26]= 0d;
        c_obliquity[27]= 0d;
        c_obliquity[28]= 0d;
        c_obliquity[29]= 0d;
        c_obliquity[30]= 0d;
        c_obliquity[31]= 977 -0.5*rlxj;
        c_obliquity[32]= -7d;
        c_obliquity[33]= 200d;
        c_obliquity[34]= 129 -0.1*rlxj;
        c_obliquity[35]= -1d;
        c_obliquity[36]= -53d;
        c_obliquity[37]= 	-2d;
        c_obliquity[38]= 	-33d;
        c_obliquity[39]= 	32d;
        c_obliquity[40]= 	26d;

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

        delta_lon=0;
        delta_obliquity=0;


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

        return delta_lon*Math.cos(hcjj*pi/180.0);
        //返回值：赤经章动
    }

    public static Position LLAToXyz(PositionLLA pos) {
        Position Xyz = new Position();
        double pi=3.14159265358979323846;
        double Re=6378140;
        //f参考椭球体扁率
        double f=0.003352813178;
        double N=Re/Math.sqrt(Math.pow(Math.cos(pos.getLat()*pi/180.0),2)+Math.pow(1.0-f,2)*Math.pow(Math.sin(pos.getLat()*pi/180.0),2));
        Xyz.setX(N*Math.cos(pos.getLat()*pi/180.0)*Math.cos(pos.getLon()*pi/180.0));
        Xyz.setY(N*Math.cos(pos.getLat()*pi/180.0)*Math.sin(pos.getLon()*pi/180.0));
        Xyz.setZ(N*Math.pow(1.0-f,2)*Math.sin(pos.getLat()*pi/180.0));
        return Xyz;
    }

    public static 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;
    }

    /**
     * 用于求各时刻太阳高度角
     *
     * @param sunTime
     * @param pSat
     * @return
     */
    public static double sunAngle(Calendar sunTime, SphericalPoint pSat) {
        SphericalPoint sunPos = new SphericalPoint();
        double pi=3.14159265358979323846;
        //UTC时间 2000年1月1日 11:58:55.816
        Calendar utcJ2000Time = TimeUtils.getCalendar(2000, 1, 1, 11, 58, 55);
        long count = (sunTime.getTimeInMillis() - utcJ2000Time.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*pi/180)*Math.cos(hcAngle*pi/180))*180/pi);}
        if(l>90 && l<=180)  { sunPos.setLon(180 +Math.atan(Math.tan(l*pi/180)*Math.cos(hcAngle*pi/180))*180/pi);}
        if(l>180 && l<=270) { sunPos.setLon(180 +Math.atan(Math.tan(l*pi/180)*Math.cos(hcAngle*pi/180))*180/pi);}
        if(l>270 && l<=360) { sunPos.setLon(360 +Math.atan(Math.tan(l*pi/180)*Math.cos(hcAngle*pi/180))*180/pi);}
        //通过太阳黄经L、黄赤交角hcangle求太阳赤纬llat
        sunPos.setLat(Math.asin(Math.sin(l*pi/180)*Math.sin(hcAngle*pi/180))*180/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;}
    }


    public static double SgTime(Calendar sunTime, double lonNutation, double delta) {
        //  J2000.0标准历元的起始时间
        Calendar startTime = TimeUtils.getCalendar(2000, 1, 1, 11, 58, 55);
        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;
    }

  //meter
    public static double computeDistanceRadBySatFovAngle(double altitude,double angle){

        double R = AstroConst.R_Earth;
        double F = 90-Math.toDegrees(Math.acos(R/(R+altitude)));//bigest Hfov angle
        if(angle>F)
            angle = F-0.0001;
        double n = 1+altitude/R;
        double a_rad = -Math.toRadians(angle)+Math.asin(n*Math.sin(Math.toRadians(angle)));
        return a_rad;
    }

    //altitude in meter
    public static double computeMaxCoverageDistanceKMBySatAndEarth(double altitude){
        double R = AstroConst.R_Earth;

        double theta = Math.acos(R/(R+altitude));
        double distance_km_to_stationView = theta*R;
        return distance_km_to_stationView/1000;
    }
    


    
    public static void main(String[] arg){
        double angle = GeoUtils.sunAngle(TimeUtils.try2ParseTime("2018-03-03 2:17:00"),new SphericalPoint(113,38));
        System.out.println("angle:"+angle);
    }

}
