之前自己也研究过,,就没研究懂。。只是作了翻译跟一些摘要,,另外没有转化成自己需要的数据。。
这次,,自己再次重拾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那个,,是否需要章动数据?