星历表ASC2EPH修改且汉化

最新推荐文章于 2022-11-04 19:40:13 发布

灵易联盟

最新推荐文章于 2022-11-04 19:40:13 发布

阅读量2.1k

点赞数

分类专栏：星历表文章标签： header character integer subroutine db2 access

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

版权

星历表专栏收录该内容

6 篇文章 1 订阅

订阅专栏

C **********************************************************************************
C       PARAMETER ( NRECL = 4 )	!基本是用重定义输入
C       PARAMETER ( NRECL = 1 )	!所以最好直接设定或改HEAER
C **********************************************************************************
C	Linyee	2011-05		修改
C	增加NRECL参数 增加了文件名参数 增加了数据转换自动适应能力。增加了文件自动查找转换功能。
C	与官方最大的不同是，可以直接放到DE405目录下执行 ASC2EPH < header.405
C	可以直接生成HEAER中指定的文件名JPLEPH.EPH，唯一需要的是将header.405第一行修改成下边的形式
C---------------------------------------------------------------
CKSIZE=  2036    NCOEFF=  1018     NRECL=    4    JPLEPH.EPH
C    20  ascp1600.405    ascp2200.405
C    00  ascp1600.405    ascp2200.405
C---------------------------------------------------------------
CKSIZE=  1456    NCOEFF=   728     NRECL=    4    JPLEPH.EPH
C  -100  ascm3000.406    ascm0100.406
C   100  ascp0000.406    ascp2900.406
C---------------------------------------------------------------
C	转换以后只原本1/3不到的数据大小。第二行是第一步循回。第三行是第二次循回。不用时00。
C	要用时就根据实际设置比如DE406设为-100，即公元前及公元后。两组。
C	单个文件小点操作比较快，所以本法比原法用COPY生成新.405再转换会快非常多。
C **********************************************************************************
C---------------------------------------------------------------
	!主程式开始
C---------------------------------------------------------------

	!定义部分+++++++++++++++++

        IMPLICIT NONE


       CHARACTER*6  CNAM (400)
       CHARACTER*6  TTL  (14,3)
       CHARACTER*12 HEADER
       CHARACTER*20 DFEFLE	!文件名
	 CHARACTER*4  ASCNME(6), ASCNM2(6)
	 CHARACTER*12 ASCFNM(2)	!起始终止文件名

       DOUBLE PRECISION  AU
       DOUBLE PRECISION  CVAL (400)
       DOUBLE PRECISION  DB(3000)
       DOUBLE PRECISION  DB2Z
       DOUBLE PRECISION  EMRAT
       DOUBLE PRECISION  SS   (3)
       DOUBLE PRECISION  T1, T2

	 INTEGER  ASCLEN, ASCSTT, ASCEND, ASCLE2, ASCST2, ASCEN2
	 INTEGER  DFENCO, DFEBLK
       INTEGER  I, IPT (3,12), IRECSZ, IN
       INTEGER  J
       INTEGER  K,KK,KP2,KSIZE
       INTEGER  LPT (3)
       INTEGER  N, NCON, NROUT, NCOEFF, NRW, NUMDE
       INTEGER  NRECL
       INTEGER  OUT


       LOGICAL  FIRST, IWHILE,FTRUE
       DATA     FIRST / .TRUE. /
       DATA     FTRUE / .TRUE. /

C       PARAMETER ( NRECL = 4 )	!基本是用重定义输入
C       PARAMETER ( NRECL = 1 )	!所以最好直接设定或修HEAER

       db2z=0.d0			!起始儒略日
       T1 = 0.D0
       T2 = 9999999.d0	!儒略日最大限制
	 

 	!代码部分+++++++++++++++++

       write(*,*) ' JPL ASCII-TO-DIRECT-I/O 程式. ' //
     &            ' 最后修改 2011-05.'
     	write(*,*)	'Linyee 修改兼汉化，NRECL参数改为HEADER文件中获取'


 	!读取参数+++++++++++++++++
C      读取main ephemeris格式大小的记录.(从 header.xxx)

	 !跳过6个字符取6位数整数并输入!NRECL=4
       READ  (*,'(6X,I6,2(11X,I6),A20)') KSIZE,DFENCO,NRECL,DFEFLE	
       WRITE (*,100)              KSIZE,DFENCO,NRECL,DFEFLE
100    FORMAT(/'KSIZE =',I6,' DFENCOEFF=',I6,'     NRECL=',I6,a20)

       IRECSZ = NRECL * KSIZE	!单元模块大小 文本1字节 二进制VFT默认4字节一单元 
	 !stop 'debug'
 	!打开文件+++++++++++++++++
	!文件号12,JPLEPH文件名,任意读写,二进制方式,模块大小

       OPEN ( UNIT   = 12,
     .        FILE   = DFEFLE,
     .        ACCESS = 'DIRECT',
     .        FORM   = 'UNFORMATTED',
     .        RECL   = 1,
     .        STATUS = 'REPLACE' )
	!NEW OLD REPLACE UNKNOW SCRATCH ,SEQUENTIAL DIRECT,UNFORMATTED FORMATTED
	 WRITE (12,REC=1,IOSTAT=OUT) (DB(K),K=1,1)
	 HEADER=''
	 read(12,REC=1,IOSTAT=IN) HEADER
	 DFEBLK=len_trim(HEADER)
	 write (*,*) 	"默认单元大小：" ,DFEBLK

	!重新打开文件
        CLOSE (12)
	IRECSZ=IRECSZ / DFEBLK
       OPEN ( UNIT   = 12,
     .        FILE   = DFEFLE,
     .        ACCESS = 'DIRECT',
     .        FORM   = 'UNFORMATTED',
     .        RECL   = IRECSZ,
     .        STATUS = 'REPLACE' )
	 WRITE (12,REC=1,IOSTAT=OUT) (DB(K),K=1,1)
	 !stop 'debug'

	!文件名自适应功能
       READ  (*,'(I8,3a4,4x,3a4)') ASCLEN,ASCNME	!每行一读
	 read(ASCNME(2),*) ASCSTT	!当前文件号
	 read(ASCNME(5),*) ASCEND	!终止文件号

	Write(ASCFNM(1),'(a4,i4.4,a4)') ASCNME(1),ASCSTT,ASCNME(3)
	 ASCFNM(2)=ASCNME(4) // ASCNME(5) // ASCNME(6)

       WRITE (*,'(3a15)') '起始-终止文件名:',ASCFNM(1),ASCFNM(2)
     	 WRITE (*,*)ASCLEN


       READ  (*,'(I8,3a4,4x,3a4)') ASCLE2,ASCNM2	!每行一读
	 ASCFNM(1)=ASCNM2(1) // ASCNM2(2) // ASCNM2(3)
	 ASCFNM(2)=ASCNM2(4) // ASCNM2(5) // ASCNM2(6)

	 read(ASCNM2(2),*) ASCST2	!当前文件号
	 read(ASCNM2(5),*) ASCEN2	!终止文件号
       WRITE (*,'(3a15)') '起始-终止文件名:',ASCFNM(1),ASCFNM(2)
     	 WRITE (*,*)ASCLE2

	!stop "debug"


	!检测参数
       if(nrecl .ne. 1 .and. nrecl .ne. 4) then
         write(*,*)'*** 错误: 未设置或使用了未定义的 NRECL ***'
         stop
       endif

	!读(GROUP 1010)组.并输出

       CALL  NXTGRP ( HEADER )
       IF (HEADER .NE. 'GROUP   1010')  THEN
          CALL  ERRPRT ( 1010, '未找到组号' )
       ENDIF

       READ (*,'(14A6)')    TTL
       WRITE(*,'(/(14A6))') TTL

	!读(GROUP 1030)组.三个数据开始儒略日，结束日，数据块长度

       CALL  NXTGRP ( HEADER )
       IF ( HEADER .NE. 'GROUP   1030' ) THEN
          CALL ERRPRT ( 1030, '未找到组号' )
       ENDIF

       READ (*,'(3D12.0)')  SS
       WRITE(*,*) SS

	!读(GROUP 1040)组。名称

       CALL  NXTGRP ( HEADER )
       IF ( HEADER .NE. 'GROUP   1040 ') THEN
          CALL ERRPRT ( 1040, '未找到组号' )
       ENDIF

	!隐式循环(I),I=3,11,4即3 7 11三次。
	!读取8字符只能存6字符
       READ (*,'(I6)')    N
       READ (*,'(10A8)') (CNAM(I),I=1,N)	
       !WRITE(*,*) CNAM
       NCON = N

	!读取数据忽略空格。值

       CALL  NXTGRP ( HEADER )
       IF ( HEADER .NE. 'GROUP   1041' ) THEN
          CALL ERRPRT ( 1041, '未找到组号' )
       ENDIF

       READ (*,'(I6)')       N
       READ (*,*)  (CVAL(I),I=1,N)

       DO  I = 1, N
           IF ( CNAM(I) .EQ. 'AU    ' )  AU    = CVAL(I)
           IF ( CNAM(I) .EQ. 'EMRAT ' )  EMRAT = CVAL(I)
           IF ( CNAM(I) .EQ. 'DENUM ' )  NUMDE = CVAL(I)
       END DO

	!格式输出最多500行两列。每列前8为名后为双精
      WRITE (*,'(500(/2(A8,D24.16)))') (CNAM(I),CVAL(I),I=1,N)

	!读(GROUP 1050)组.13列

       CALL  NXTGRP ( HEADER )
       IF ( HEADER .NE. 'GROUP   1050' ) THEN
          CALL ERRPRT ( 1050, '未找到组号' )
       ENDIF

	!每行13列，三行。前12列到IPT，第13列到LPT。
       READ (*,'(13I6)')   ((IPT(I,J),J=1,12),LPT(I),I=1,3)
	!Linyee 增加属性个数，因为除了章动2其它都3，所以又去掉
	!do i=1 , 11
	!	IPT(4,i)=(IPT(1,i+1)-IPT(1,i))/IPT(2,i)/IPT(3,i)
	!end do
	!	IPT(4,12)=(LPT(1)-IPT(1,12))/IPT(2,12)/IPT(3,12)
	!	LPT(4)=(DFENCO+1-LPT(1))/LPT(2)/LPT(3)

      WRITE(*,'(/(3I5))') IPT, LPT	!输完第一数再输第二数.

 	!读取数据+++++++++++++++++
	!读(GROUP 1070)组.数据组

      CALL  NXTGRP ( HEADER )
      IF ( HEADER .NE. 'GROUP   1070' ) CALL ERRPRT(1070,'未找到组号')
C-----------------------------------------------------------------
	!循环读数据！ASCLEN, ASCSTT, ASCEND
C-----------------------------------------------------------------
      NROUT  = 0
      IN     = 0
      OUT    = 0

	!第一循回
77	if (ASCLEN .NE. 0) then
		IWHILE =.TRUE.
		go to 55
	endif
		go to 66
	stop	'第一循回失败转第二循回'


55	do while(IWHILE)

	Write(ASCFNM(1),'(a4,i4.4,a4)') ASCNME(1),ASCSTT,ASCNME(3)
	!打开文件
       OPEN ( UNIT   = 22,
     .        FILE   = ASCFNM(1),
     .        ACCESS = 'SEQUENTIAL',
     .        FORM   = 'FORMATTED',
     .        IOSTAT   = IN,
C     .        RECL   = IRECSZ,
     .        STATUS = 'OLD' )

       if (IN .ne. 0) then
		call  ERRPRT(ASCSTT,ASCFNM(1) // ' 未找到后续数据文件')
	 endif
C---------------------------------------------------------------------------
	!文件记录遍历
C---------------------------------------------------------------------------
	!读数据。记录号，记录项。
  11  read(22,'(2i6)',IOSTAT =IN)nrw,ncoeff
          if(IN .NE. 0)then
             CALL ERRPRT(1070,' 未找到数据块')
          endif
      if(nrw .eq. 0) go to 11

	!读共1018-1个数据，每行3个。
      do k=1,ncoeff,3
        kp2=min(k+2,ncoeff)	!每行结束号
        read(22,*,IOSTAT = IN)(DB(kk),kk=k,kp2)
          if(IN .NE. 0)then
             CALL ERRPRT(k,' 此条数据读取错误或文件已结束')
          endif
      enddo


      DO WHILE (       ( IN    .EQ. 0 )
     .           .AND. ( DB(2) .LT. T2) )

          IF ( 2*NCOEFF .NE. KSIZE ) THEN
             CALL ERRPRT(NCOEFF,' 记录项数，
     .	   不符合规定大小')
          ENDIF

          IF  ( (DB(2) .GE. T1) .AND. (DB(1) .GE. DB2Z) ) THEN
             IF ( FIRST ) THEN
                DB2Z  = DB(1)
                FIRST = .FALSE.
             ENDIF

             IF (DB(1) .NE. DB2Z ) THEN
                CALL ERRPRT (NRW, '记录初始日期错误')
             ENDIF

             DB2Z  = DB(2)	!置下一次初始日期
             NROUT = NROUT + 1
C	保存到文件
             WRITE (12,REC=NROUT+2,IOSTAT=OUT) (DB(K),K=1,NCOEFF)

             IF ( OUT .NE. 0 ) THEN
                CALL ERRPRT (NROUT, 
     &                     '记录保存错误')
             ENDIF
C	开始日期，结束日期，日期长度
             IF (NROUT .EQ. 1) THEN
                SS(1) = DB(1)
                SS(3) = DB(2) - DB(1)
             ENDIF
             SS(2) = DB(2)

C	IF ( NROUT .EQ. 1 ) THEN
C		WRITE (*,271) NROUT, DB(2)
C	ENDIF
271          FORMAT (I6,
     &              ' 条EPH记录.  结束JED = ',
     &              F12.2," :From: ",A12)

C	输出个位为1的记录信息

			
             IF ( MOD(NROUT,200) .EQ. 1 ) THEN
                IF ( DB(1) .GE. T1 ) THEN
                   WRITE (*,271) NROUT, DB(2),ASCFNM(1)
                ELSE
                   WRITE (*,*) ' 查找第一记录中……'
                ENDIF
             ENDIF

          ENDIF
	!查找下一条记录
          READ (22,'(2I6)',IOSTAT =IN) NRW, NCOEFF

          if(IN .eq. 0)then
            do k=1,ncoeff,3
              kp2=min(k+2,ncoeff)
              read(22,*,IOSTAT = IN)(DB(kk),kk=k,kp2)
			if(IN .NE. 0)then
			   CALL ERRPRT(k,' 此条数据读取错误或文件已结束')
			endif
            enddo
          endif
      END DO
 
	!下一文件完全判别
       CLOSE (22)
	ASCSTT=ASCSTT+ASCLEN
	IWHILE= .FALSE.
	if ((ASCSTT .GE. ASCEND) .and. (ASCLEN .LT. 0)) then
	 IWHILE=.true.
          endif
	if ((ASCSTT .LE. ASCEND) .and. (ASCLEN .GT. 0)) then
	 IWHILE=.true.
          endif
      END DO

C---------------------------------------------------------------------------
	!第二循回
C---------------------------------------------------------------------------
66	if (ASCLE2 .NE. 0) then
		ASCSTT=ASCST2
		ASCLEN=ASCLE2
		ASCEND=ASCEN2
		ASCNME=ASCNM2

		ASCLE2=0
		 IWHILE=.true.
		go to 55
	endif

C---------------------------------------------------------------------------
	!循回结束.
C---------------------------------------------------------------------------


C	输出最后的记录信息
       WRITE (*,271) NROUT, DB(2),ASCFNM(1)

C      输出头信息

       NROUT = 1
       WRITE(12,REC=1,IOSTAT=OUT) TTL,CNAM,SS,NCON,AU,
     &                            EMRAT,IPT,NUMDE,LPT
       IF ( OUT .NE. 0 ) THEN
          CALL ERRPRT ( NROUT, '记录保存错误')
       ENDIF

C      输出头信息中的值表

       NROUT = 2
       WRITE(12,REC=2,IOSTAT=OUT)CVAL
       IF ( OUT .NE. 0 ) THEN
          CALL ERRPRT ( NROUT, '记录保存错误')
       ENDIF


       CLOSE (12)
       STOP ' OK'
       END
C---------------------------------------------------------------
	!主程式结束
C---------------------------------------------------------------



C---------------------------------------------------------------
	!输出错误信息并中止程式
C---------------------------------------------------------------
       SUBROUTINE  ERRPRT (I, MSG)

       CHARACTER*(*)  MSG
       INTEGER        I

       WRITE (*,200)  I, MSG
200    FORMAT('ERROR #',I8,2X,A50)

       STOP ' ERROR '
       END
C---------------------------------------------------------------
C---------------------------------------------------------------


C---------------------------------------------------------------
	!进入下一组数据区
	!'GROUP   1070'行作为结束行.
C---------------------------------------------------------------
      SUBROUTINE  NXTGRP ( HEADER )

      CHARACTER*(*)  HEADER
      CHARACTER*12   BLANK

 44     HEADER = ' '
      DO WHILE ( HEADER .EQ. ' ' )	!跳过空格行
          READ (*,'(A)') HEADER
      ENDDO

	READ(HEADER,'(A5)') BLANK		!不是组标识再找过
	if (trim(BLANK) .NE. "GROUP") goto 44

      IF ( HEADER .NE. 'GROUP   1070' )READ (*, '(A)') BLANK	!跳过空格

      RETURN
      END
C---------------------------------------------------------------
C---------------------------------------------------------------

WINDOWS下直接用 VisualFortran6.5 编译 DF ASC2eph.f就会生成 ASC2eph.exe

然后将程序复制到405 406下，将头文件第一行改成

KSIZE= 2036 NCOEFF= 1018 NRECL= 4 JPLEPH.405 20 ascp1600.405 ascp2200.405 00 ascp1600.405 ascp2200.405

406改成

KSIZE= 1456 NCOEFF= 728 NRECL= 4 JPLEPH.406 -100 ascm3000.406 ascm0100.406 100 ascp0000.406 ascp2900.406

直接运行 ASC2eph <heaer.405

就会根据HEADER前三行声明生成JPLEPH.405的文件。每个记录项大小为 Ksize*4即 NCOEFF*8。双精度单项大小。

源码跟程式都上传到了CSDN资源上。需要可以去找下。

灵易联盟

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
打赏
2
评论
星历表ASC2EPH修改且汉化

C **********************************************************************************C PARAMETER ( NRECL = 4 ) !基本是用重定义输入C PARAMETER ( NRECL = 1 ) !所以最好直接设定或改HEAERC ********************...
复制链接

扫一扫