星历表（DE406+Testeph.f）整理兼求教

最新推荐文章于 2024-05-13 16:02:52 发布

灵易联盟

最新推荐文章于 2024-05-13 16:02:52 发布

阅读量1.7k

点赞数

分类专栏： # 历法资料星历表

本文链接：https://blog.csdn.net/jiaguoxinzhi/article/details/89740704

版权

星历表同时被 2 个专栏收录

6 篇文章 1 订阅

订阅专栏

历法资料

3 篇文章 0 订阅

订阅专栏

之前自己也研究过，，就没研究懂。。只是作了翻译跟一些摘要，，另外没有转化成自己需要的数据。。

这次，，自己再次重拾Fortran+de406研究了下。。
现在卡在
J2000的坐标，，怎么转成实时的黄经钤度。。

求大神们指教。

整理资料如下。

C---------------------------------------------------------------
	!主程式开始
C---------------------------------------------------------------
       CHARACTER*6  NAMS(400)
       CHARACTER*3  ALF3

       DOUBLE PRECISION  DEL
       DOUBLE PRECISION  ET
       DOUBLE PRECISION  R(6)
       DOUBLE PRECISION  SS(3)
       DOUBLE PRECISION  VALS(400)
       DOUBLE PRECISION  XI
  
       INTEGER  I
       INTEGER  LINE
       data line/0/
       INTEGER  NVS,NTARG,NCTR,NCOORD
       INTEGER  NPT
       data npt/100/

	 INTEGER DFEBLK
!+++++++++++++++++++++++++++++++++++++++++
	!定义结束
!+++++++++++++++++++++++++++++++++++++++++

       CALL  Const (NAMS, VALS, SS, NVS)
       WRITE(*,*) ' JPL 测试-EPHEMERIS 程式. ' //
     .            ' 最后更新 2011-5.'
!       WRITE (*,'(/3F14.2/200(/2(A8,D24.16)))') SS,(NAMS(I)
!     . ,VALS(I),I=1,NVS)			!输出表头信息
       WRITE (*,*) SS,(NAMS(I)
     . ,VALS(I),I=1,NVS)			!输出表头信息
       WRITE (*,'(3F14.2)') SS			!输出表头信息


!   1  WRITE(*,'(a4)')'EOT?'
!	 READ(*,'(a3)')ALF3
!      IF(ALF3 .NE. 'EOT') GO TO 1	!寻找EOT标记

      WRITE(*,'(A91)') 
     . '  line -- jed --   t#   c#   x#   --- jpl value ---'//
     . '   --- user value --    -- difference --'
      WRITE(*,*)
   2  WRITE(*,'(a8)')'5个参数>'

!READ(*,'(15X,D10.1,3I3,F30.20)',END=9)  !跳过前缀、日期
	READ(*,'(15X,D10.1,3I3,F30.20)',END=9)  !跳过前缀、日期
     . ET,NTARG,NCTR,NCOORD,XI
!               2458602.12 11 12 0 0.0

!2458602.12000000               11          12 -1.713231114586650E-003  6.979541948858645E-003  2.990473552032142E-003 -8.349798584670847E-006  5.846299874843524E-007  4.893236050471883E-007

      IF(ET .LT. SS(1)) GO TO 902	!时间不在表内重读.
      IF(ET .GT. SS(2)) GO TO 902
!	GO TO 902

      CALL  Pleph ( ET, NTARG, NCTR, R )	!时间点,两星位置速度至R
	GO TO 903

      DEL=DABS(R(NCOORD)-XI)	!误差计算
      if(ntarg .eq. 15 .and. ncoord .eq. 3)then
        del=del/(0.23d0*(et-2451545.d0))
      endif
      if(ntarg .eq. 15 .and. ncoord .eq. 6)then
        del=del*0.01/(1.0d0+(et-2451545.d0)/365.25e2)
      endif

      LINE=LINE+1	!逢百输出
      IF ( MOD(LINE,npt) .EQ. 0 )  WRITE(*,200)
     . LINE,ET,NTARG,NCTR,NCOORD,XI,R(NCOORD),DEL
 200  FORMAT(I6,F10.1,3I5,2F25.13,1E13.5)

      IF (DEL .GE. 1.D-13) WRITE(*,201) !误差太大
     . LINE,ET,NTARG,NCTR,NCOORD,XI,R(NCOORD),DEL
 201  FORMAT(/ '*****  WARNING : next difference >= 1.D-13  *****'//
     . I6,F10.1,3I3,2F25.13,1E13.5/' ')

      GO TO 2		!下一条测试数据

 902	WRITE(*,*)ET,NTARG,NCTR,NCOORD,XI
	Go To 2 !下一条测试数据

 903	WRITE(*,*) "执行结果",ET, NTARG, NCTR, R
	Go To 2 !下一条测试数据

   9  STOP "OK"
	Go To 2 !下一条测试数据
      END
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	!主程式结束
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


C---------------------------------------------------------------
	!参数生成器（字量，项量，文件号，文件名）
C---------------------------------------------------------------
      SUBROUTINE FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL)

      CHARACTER*80  NAMFIL

      NAMFIL= 'JPLEPH'
      NRECL= 4	!or 1
      NRFILE=12
      IRECSZ=0

      INQUIRE(FILE=NAMFIL,RECL=IRECSZ)
	 !取得文件名为NAMFIL的数据块大小


      IF(IRECSZ .LE. 0) write(*,*)
     . ' 未查询到正在作业的文件' // NAMFIL

      KSIZE=IRECSZ/NRECL

      RETURN

      END
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



C---------------------------------------------------------------
	!参数生成器（字量，项量，文件号，文件名）
C---------------------------------------------------------------
      SUBROUTINE FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)

      SAVE

      CHARACTER*6 TTL(14,3),CNAM(400)
      CHARACTER*80 NAMFIL

      DIMENSION SS(3)	!动态数组

      INTEGER IPT(3,13)


      NRECL=4
      NRFILE=12
      NAMFIL= 'JPLEPH'


      MRECL=NRECL*1000

        OPEN(NRFILE,
     *       FILE=NAMFIL,
     *       ACCESS='DIRECT',
     *       FORM='UNFORMATTED',
     *       RECL=MRECL,
     *       STATUS='OLD')

      READ(NRFILE,REC=1)TTL,CNAM,SS,NCON,AU,EMRAT,
     * ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)

      CLOSE(NRFILE)




      KMX = 0
      KHI = 0

      DO I = 1,13
         IF (IPT(1,I) .GT. KMX) THEN
            KMX = IPT(1,I)
            KHI = I
         ENDIF
      ENDDO

      ND = 3
      IF (KHI .EQ. 12) ND=2

      KSIZE = 2*(IPT(1,KHI)+ND*IPT(2,KHI)*IPT(3,KHI)-1)

      RETURN

      END
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



C---------------------------------------------------------------
	!位置和速度(儒略日，目标行星，中心行星，返回位置速度（6）)
C---------------------------------------------------------------
      SUBROUTINE PLEPH ( ET, NTARG, NCENT, RRD )
	
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      DIMENSION RRD(6),ET2Z(2),ET2(2),PV(6,13)
      DIMENSION SS(3),CVAL(400),PVSUN(6),zips(2)
      data zips/2*0.d0/

      LOGICAL BSAVE,KM,BARY
      LOGICAL FIRST
      DATA FIRST/.TRUE./

      INTEGER LIST(12),IPT(39),DENUM


      COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT

      COMMON/STCOMX/KM,BARY,PVSUN

	WRITE(*,'(a8)')' PLEPH>'
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	!定义结束!初始化
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      ET2(1)=ET
      ET2(2)=0.D0
      GO TO 11

	!给定儒略日数组的新入口
      ENTRY DPLEPH(ET2Z,NTARG,NCENT,RRD)

      ET2(1)=ET2Z(1)
      ET2(2)=ET2Z(2)

  11  WRITE(*,'(a8)')' DPLEPH>'
	IF(FIRST) CALL STATE(zips,list,pv,pnut)
      FIRST=.FALSE.	!第一次取参数

  96  IF(NTARG .EQ. NCENT) RETURN	!目标同于中心退出不计算

      DO I=1,12
		LIST(I)=0
      ENDDO	!初始化LIST

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	!取两星参数
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	!第14星代码 章动
      IF(NTARG.NE.14) GO TO 97
        IF(IPT(35).GT.0) THEN
          LIST(11)=2
          CALL STATE(ET2,LIST,PV,RRD)
          RETURN
        ELSE
          do i=1,4
			rrd(i)=0.d0
          enddo
          WRITE(6,297)
  297     FORMAT(' *****  星历表中没有章动的信息  *****')
          STOP
        ENDIF


  97  do i=1,6
		rrd(i)=0.d0
      enddo

	!第15星代码 天平动
      IF(NTARG.NE.15) GO TO 98
        IF(IPT(38).GT.0) THEN
          LIST(12)=2
          CALL STATE(ET2,LIST,PV,RRD)
          DO I=1,6
			RRD(I)=PV(I,11)
          ENDDO
          RETURN
        ELSE
          WRITE(6,298)
  298     FORMAT(' *****  星历表中没有天平动的信息  *****')
          STOP
        ENDIF


	!第1-13星代码 第10星L3 3星L10 13星L3 其它10星不变
  98  BSAVE=BARY
      BARY=.TRUE.


      DO I=1,2
		K=NTARG
		IF(I .EQ. 2) K=NCENT
		IF(K .LE. 10) LIST(K)=2
		IF(K .EQ. 10) LIST(3)=2
		IF(K .EQ. 3) LIST(10)=2
		IF(K .EQ. 13) LIST(3)=2
      ENDDO


      CALL STATE(ET2,LIST,PV,RRD)

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	!各星数据整理
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      IF(NTARG .EQ. 11 .OR. NCENT .EQ. 11) THEN
		DO I=1,6
			PV(I,11)=PVSUN(I)
		ENDDO
      ENDIF	!11太阳参数

      IF(NTARG .EQ. 12 .OR. NCENT .EQ. 12) THEN
		DO I=1,6
			PV(I,12)=0.D0
		ENDDO
      ENDIF	!12星置0

      IF(NTARG .EQ. 13 .OR. NCENT .EQ. 13) THEN
		DO I=1,6
			PV(I,13)=PV(I,3)
		ENDDO
      ENDIF	!13星从L3取

      IF(NTARG*NCENT .EQ. 30 .AND. NTARG+NCENT .EQ. 13) THEN
		DO I=1,6
			PV(I,3)=0.D0
		ENDDO
		GO TO 99
      ENDIF	!地月系,3星地球置0

      IF(LIST(3) .EQ. 2) THEN
		DO I=1,6
			PV(I,3)=PV(I,3)-PV(I,10)/(1.D0+EMRAT)
		ENDDO
      ENDIF	!月球或地月置扣去地球参数

      IF(LIST(10) .EQ. 2) THEN
		DO I=1,6
			PV(I,10)=PV(I,3)+PV(I,10)
		ENDDO
      ENDIF	!地球置地球及月球之和
	Write(*,*) '修正结果',PV

  99  DO I=1,6
		RRD(I)=PV(I,NTARG)-PV(I,NCENT)
      ENDDO	!求出两星差值	
	Write(*,*) '最终结果',RRD

      BARY=BSAVE
      RETURN
      END

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



C---------------------------------------------------------------
	!插值（数据项13表1，偏量，13表2（即系数个数），2 or 3(即数据项数)，13表3行即区段，2，返回pv）
	!NCM*IFL=PV的长度NCF NCM L决定了EPH位置，即系数。
C---------------------------------------------------------------
      SUBROUTINE INTERP(BUF,T,NCF,NCM,NA,IFL,PV)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      SAVE
C
      DOUBLE PRECISION BUF(NCF,NCM,*),T(2),PV(NCM,*),PC(18),VC(18)

C
      DATA NP/2/
      DATA NV/3/
      DATA TWOT/0.D0/
      DATA PC(1),PC(2)/1.D0,0.D0/
      DATA VC(2)/1.D0/


	Write(*,*) "T",T,"NCF",NCF,"NCM",NCM,"NA",NA,"IFL",IFL
!+++++++++++++++++++++++++++++++++++++++++
	!定义结束!计算插项基础
!+++++++++++++++++++++++++++++++++++++++++
      DNA=DBLE(NA) !段数 DBLE转双精
      DT1=DINT(T(1))		!记录内偏量整数部
      TEMP=DNA*T(1)		!时间点偏量
      L=IDINT(TEMP-DT1)+1	!时间点位置

      TC=2.D0*(DMOD(TEMP,1.D0)+DT1)-1.D0

      IF(TC.NE.PC(2)) THEN
        NP=2
        NV=3
        PC(2)=TC
        TWOT=TC+TC
      ENDIF

!+++++++++++++++++++++++++++++++++++++++++
	!计算第一组属性!位置属性
!+++++++++++++++++++++++++++++++++++++++++
      IF(NP.LT.NCF) THEN !数据项数
        DO 1 I=NP+1,NCF
		  PC(I)=TWOT*PC(I-1)-PC(I-2)
    1   CONTINUE
        NP=NCF
      ENDIF
	Write(*,*) 'TC',TC,'第一切比',PC

      DO 2 I=1,NCM	!数据类数NCM*2
		PV(I,1)=0.D0
		DO 3 J=NCF,1,-1
			PV(I,1)=PV(I,1)+PC(J)*BUF(J,I,L)
			!Write(*,*) '库内系数',J,BUF(J,I,L)
    3		CONTINUE
	Write(*,*) '第一组属性',I,PV(I,1)
    2 CONTINUE

!+++++++++++++++++++++++++++++++++++++++++
	!计算第二组属性!速度属性
!+++++++++++++++++++++++++++++++++++++++++
      IF(IFL.LE.1) RETURN	!1组值
      VFAC=(DNA+DNA)/T(2)
      VC(3)=TWOT+TWOT

      IF(NV.LT.NCF) THEN
        DO 4 I=NV+1,NCF
		  VC(I)=TWOT*VC(I-1)+PC(I-1)+PC(I-1)-VC(I-2)
    4   CONTINUE
        NV=NCF
      ENDIF
	Write(*,*) 'VFAC',VFAC,'第二切比',VC

      DO 5 I=1,NCM
		PV(I,2)=0.D0
		DO 6 J=NCF,2,-1
			PV(I,2)=PV(I,2)+VC(J)*BUF(J,I,L)
			!Write(*,*) '库内系数',J,BUF(J,I,L)
    6     CONTINUE
		PV(I,2)=PV(I,2)*VFAC
		Write(*,*) '第二组属性',I,PV(I,2)
    5 CONTINUE
      RETURN
      END
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	!程序结束
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


C---------------------------------------------------------------
	!双精分隔(双精,结果)3.5 分隔成2 1.5 ;-0.5 --> -1 1.5
C---------------------------------------------------------------
      SUBROUTINE SPLIT(TT,FR)
		IMPLICIT DOUBLE PRECISION (A-H,O-Z)

		DIMENSION FR(2)

		FR(1)=DINT(TT)	!双精截成整数
		FR(2)=TT-FR(1)	!双精小数部分
		IF(TT.GE.0.D0 .OR. FR(2).EQ.0.D0) RETURN !小数部为零退出

		FR(1)=FR(1)-1.D0	!整数部-1到
		FR(2)=FR(2)+1.D0	!小数部+1
		RETURN
      END
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


C---------------------------------------------------------------
	!（儒略日，，返回值，返回值）
C---------------------------------------------------------------
      SUBROUTINE STATE(ET2,LIST,PV,PNUT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      SAVE

      DIMENSION ET2(2),PV(6,12),PNUT(4),T(2),PJD(4),BUF(1500),
     . SS(3),CVAL(400),PVSUN(6)

      INTEGER LIST(12),IPT(4,13)

      LOGICAL FIRST
      DATA FIRST/.TRUE./

      CHARACTER*6 TTL(14,3),CNAM(400)
      CHARACTER*80 NAMFIL

      LOGICAL KM,BARY
      DOUBLE PRECISION  ET,ET0

      COMMON/EPHHDR/CVAL,SS,AU,EMRAT,NUMDE,NCON,IPT
      COMMON/CHRHDR/CNAM,TTL
      COMMON/STCOMX/KM,BARY,PVSUN	!公共区STCOMX有KM,BARY,PVSUN
!+++++++++++++++++++++++++++++++++++++++++
	!定义结束！取得参数
!+++++++++++++++++++++++++++++++++++++++++

      IF(FIRST) THEN
        FIRST=.FALSE.
C	!取得参数
C        CALL FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL)
C        CALL FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)
        CALL FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)
		IF(NRECL .EQ. 0) WRITE(*,*)' 单元个数不能为0 '

		IRECSZ=NRECL*KSIZE	!记录项大小
		Call	FixRECL(IRECSZ)
		NCOEFFS=KSIZE/2		!记录子项数

        OPEN(NRFILE,
     *       FILE=NAMFIL,
     *       ACCESS='DIRECT',
     *       FORM='UNFORMATTED',
     *       RECL=IRECSZ,
     *       STATUS='OLD')

!       WRITE(12,REC=1,IOSTAT=OUT) TTL,CNAM,SS,NCON,AU,EMRAT,IPT,NUMDE,LPT
		READ(NRFILE,REC=1)TTL,CNAM,SS,NCON,AU,EMRAT,
     . ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)
		READ(NRFILE,REC=2) CVAL
       !WRITE(*,*) TTL,CNAM,SS,NCON,AU,EMRAT,IPT,NUMDE,LPT

		NRL=0

      ENDIF
!+++++++++++++++++++++++++++++++++++++++++
	!取参数结束!取PVSUN
!+++++++++++++++++++++++++++++++++++++++++



      IF(ET2(1) .EQ. 0.D0) RETURN
      ET=ET2(1)+ET2(2)	!得到儒略日
      IF(ET.LT.SS(1) .OR. ET.GT.SS(2)) GO TO 98
	!在星历表范围之外跳到98 ET2(1)+ET2(2)=儒略日

	ET0=(ET-SS(1))			!
	NR=IDINT(ET0/SS(3))+3	!第三条起才是数据
      IF(ET.EQ.SS(2)) NR=NR-1	!最后一刻回上条

	  tmp1 = DBLE(NR-3)*SS(3) + SS(1)!记录开始天数
        tmp2 = ET - tmp1	!得到记录内天数偏移
        T(1) = tmp2/SS(3)	!指向记录内偏量

      IF(NR.NE.NRL) THEN
        NRL=NR
        READ(NRFILE,REC=NR,ERR=99)(BUF(K),K=1,NCOEFFS)
      ENDIF

      IF(KM) THEN
		T(2)=SS(3)*86400.D0	!记录总秒数
		AUFAC=1.D0
      ELSE
		T(2)=SS(3)			!记录总天数
		AUFAC=1.D0/AU		!
      ENDIF


      CALL INTERP(BUF(IPT(1,11)),T,IPT(2,11),3,IPT(3,11),2,PVSUN)
      DO I=1,6
		PVSUN(I)=PVSUN(I)*AUFAC
      ENDDO
	Write(*,*)"太阳",AUFAC,T,IPT(2,11),3,IPT(3,11),2,PVSUN


!+++++++++++++++++++++++++++++++++++++++++
	!取LIST中指定的PV
	!第10星L3 3星L10 13星L3  其它10星不变
	!第11星这里没计算是因为上边已经计算了太阳。
!+++++++++++++++++++++++++++++++++++++++++
      DO 4 I=1,10
      IF(LIST(I).EQ.0) GO TO 4

      CALL INTERP(BUF(IPT(1,I)),T,IPT(2,I),3,IPT(3,I),
     & LIST(I),PV(1,I))

      DO J=1,6	!9以内BARY=FALSE时需扣去太阳参数
       IF(I.LE.9 .AND. .NOT.BARY) THEN
		 PV(J,I)=PV(J,I)*AUFAC-PVSUN(J)
       ELSE		!10即地球不变
		 PV(J,I)=PV(J,I)*AUFAC
       ENDIF
      ENDDO
	Write(*,*)"行星",I,"BARY",BARY,AUFAC,T,IPT(2,I),3,IPT(3,I),2,PV

   4  CONTINUE
!+++++++++++++++++++++++++++++++++++++++++
	!章动L11
!+++++++++++++++++++++++++++++++++++++++++


      IF(LIST(11).GT.0 .AND. IPT(2,12).GT.0)
     * CALL INTERP(BUF(IPT(1,12)),T,IPT(2,12),2,IPT(3,12),
     * LIST(11),PNUT)

!+++++++++++++++++++++++++++++++++++++++++
	!天平动L12
!+++++++++++++++++++++++++++++++++++++++++

      IF(LIST(12).GT.0 .AND. IPT(2,13).GT.0)
     * CALL INTERP(BUF(IPT(1,13)),T,IPT(2,13),3,IPT(3,13),
     * LIST(12),PV(1,11))

      RETURN
!+++++++++++++++++++++++++++++++++++++++++
	!正常结束
!+++++++++++++++++++++++++++++++++++++++++

  98  WRITE(*,198)ET2(1)+ET2(2),SS(1),SS(2)
 198  format(' ***  查询的儒略日,',f12.2,
     * ' 不在星表范围,',2f12.2,'之内 ***')
      stop	"超出范围"

   99 WRITE(*,'(2F12.2,A80)') NRL,'记录读取失败'
      STOP	"星历表文件错误"

      END
C---------------------------------------------------------------
	!非正常结束
C---------------------------------------------------------------


C---------------------------------------------------------------
	!返回（名称表，值表，范围表，返回名表项数）
C---------------------------------------------------------------
      SUBROUTINE CONST(NAM,VAL,SSS,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      SAVE	!保存变量

      CHARACTER*6 NAM(*),TTL(14,3),CNAM(400)

      DOUBLE PRECISION VAL(*),SSS(3),SS(3),CVAL(400),zips(2)
      DOUBLE PRECISION xx(99)
      data zips/2*0.d0/

      INTEGER IPT(3,13),DENUM,list(11)
      logical first
      data first/.true./

      COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT
      COMMON/CHRHDR/CNAM,TTL

!+++++++++++++++++++++++++++++++++++++++++
	!定义结束
!+++++++++++++++++++++++++++++++++++++++++

      IF(FIRST) CALL STATE(zips,list,xx,xx)
      first=.false.

      N=NCON	!名称项数

      DO I=1,3
		SSS(I)=SS(I)
      ENDDO

      DO I=1,N
		NAM(I)=CNAM(I)
		VAL(I)=CVAL(I)
      ENDDO

      RETURN

      END
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


C---------------------------------------------------------------
	!参数生成器（字量，项量，文件号，文件名）
C---------------------------------------------------------------
      SUBROUTINE FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)

      SAVE	!保存变量表，有静态变量的味道

      CHARACTER*80 NAMFIL

      NRECL=4
      NRFILE=12
      NAMFIL= 'JPLEPH.406'
      KSIZE = 1456

      RETURN
      END
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



C---------------------------------------------------------------
	!记录块大小校正（字节数）返回实际大小
C---------------------------------------------------------------
	SUBROUTINE FixRECL(IRECSZ)
	 INTEGER OUT,IN,DFEBLK,IRECSZ
	 CHARACTER*64	BLANK
	 BLANK='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'

	 !打开一个临时文件
       OPEN ( UNIT   = 10001,
     .        ACCESS = 'DIRECT',
     .        FORM   = 'UNFORMATTED',
     .        RECL   = 1,
     .        STATUS = 'SCRATCH' )
	!NEW OLD REPLACE UNKNOW SCRATCH ,SEQUENTIAL DIRECT,UNFORMATTED FORMATTED
	 !写入一个数据块
	 WRITE (10001,REC=1,IOSTAT=OUT) BLANK

	 !读取一个数据块
	 BLANK=''
	 read(10001,REC=1,IOSTAT=IN) BLANK

	 DFEBLK=len_trim(BLANK)	!求出实际长度

	!关闭文件
       CLOSE (10001)
	 IRECSZ=IRECSZ / DFEBLK

       RETURN	!返回值
      END
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++






C---------------------------------------------------------------
	!默认记录单元字节数（字节数）返回默认记录单元字节数
C---------------------------------------------------------------
	Function DFERECL()
		 INTEGER OUT,IN,DFEBLK,DFERECL
		 CHARACTER*64	BLANK
		 BLANK='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'

		 !打开一个临时文件
		 OPEN ( UNIT   = 10001,
     .        ACCESS = 'DIRECT',
     .        FORM   = 'UNFORMATTED',
     .        RECL   = 1,
     .        STATUS = 'SCRATCH' )
		!NEW OLD REPLACE UNKNOW SCRATCH ,SEQUENTIAL DIRECT,UNFORMATTED FORMATTED
		 !写入一个数据块
		 WRITE (10001,REC=1,IOSTAT=OUT) BLANK

		 !读取一个数据块
		 BLANK=''
		 read(10001,REC=1,IOSTAT=IN) BLANK

		 DFEBLK=len_trim(BLANK)	!求出实际长度

		!关闭文件
		 CLOSE (10001)
		 DFERECL=DFEBLK
		 RETURN	!返回值
      END
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

测试数据

注意事项（有个别.eph文件的double是内存低位是数值的高位即大数序，，我测试的都是内存低位是数值的低位即小数序）

               2458602.98 11 12 0 0.0
 执行结果   2458602.98000000               11          12
 -1.720412519355159E-003  6.980040683812857E-003  2.990892612037133E-003
 -8.351124054440001E-006  5.752328444190759E-007  4.852403668307290E-007

测试数据2

其中赤经"23h59m2.9s"已与其它软件比较过了是正确的，但这个只是J2000的值。

怎么转换为Of Date (黄经钤度呢)。。

===== 查询数据 ===== JDEPOC=2440400.5 1969/6/28 星期六 00:00:00
求值行星编号11 中文[太阳] 英文[SUN]
参照行星编号3 中文[地球] 英文[EARTH]
查询时间：儒略日：2458563.41550926      TZ(08:00:00)2019/3/21 星期四 05:58:20   UTC:2019/3/20 星期三 21:58:20]
日期范围：[2458544.5-2458576.5] [2019/3/2 星期六 00:00:00-2019/4/3 星期三 00:00:00] 块内偏量：18.9155092593282/32(块天数长度)
=== 不是章动，也不是天平动，继续... ===
段数：2 当前所在段号1   段起始日：2458560.5     段内偏量：2.91550925932825/16(每段天数) t=-0.635561342583969
段数：2 当前所在段号1   段起始日：2458560.5     段内偏量：2.91550925932825/16(每段天数) t=-0.635561342583969
段数：8 当前所在段号4   段起始日：2458560.5     段内偏量：2.91550925932825/4(每段天数) t=0.457754629664123
最终结果(J2000)：[[0.99586250598589188,-0.004133541887724009,-0.0017924704779106767,0.00035769200667771237,0.015837173965159212,0.0068657454664686054]]
转换球坐标：[-0.0041506916115086976,1.572596226972848,0.99587269766142916,0.00035769200667771237,0.015837173965159212,0.0068657454664686054]
转换为角度：[-0.23781711139980269,90.103126683741465,0.99587269766142916,0.00035769200667771237,0.015837173965159212,0.0068657454664686054]
转换为赤经：["23h59m2.92389326404731s","90.1031266837415","0.995872697661429","0.000357692006677712","0.0158371739651592","0.00686574546646861"]
转换为罗盘度：[315.23781711139986,90.103126683741465,0.99587269766142916,0.00035769200667771237,0.015837173965159212,0.0068657454664686054]
转换为角秒：["315°14'16.1416010395078''","90.1031266837415","0.995872697661429","0.000357692006677712","0.0158371739651592","0.00686574546646861"]
Ok [315.23781711139986,90.103126683741465,0.99587269766142916,0.00035769200667771237,0.015837173965159212,0.0068657454664686054]

下图 RA/DE(J2000) 就是现在计算出来的值。。想转为of date那个，，是否需要章动数据？