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资源上。需要可以去找下。