405
基于JPL DE405和IAU SOFA计算太阳系行星运动状态的Fortran平台搭建与测试
一些铺垫
我们为什么要得到行星运动状态
行星的运动状态是指行星的位置和速度，这两个信息对我们非常重要，比如想要为人造航天器设计一条从地球去往火星的轨道，那么至少需要知道一段时间内地球以及火星相对于太阳的位置和速度，才能找到合适的时机（包括时间与位置）脱离地球以及进入火星。
JPL DE405是什么
一种由美国Jet Propulsion Laboratory （JPL）提出的精确星历表，包括行星的章动和天平动。
IAU SOFA是什么
Standard of fundamental astronomy （SOFA）是由国际组织International Astronomical Union（IAU）提出的一套基础天文学标准，可以通过其提供的软件来调用许多重要的功能。
Fortran语言的优势
一种能够高效执行科学计算的编程语言。
Fortran平台的构成
在windows 10操作系统上，采用了最新的Microsoft Visual Studio Enterprise (VS) 2017+Intel parallel Studio XE (PSXE) 2018的方案，VS提供界面与项目管理，PSXE提供编译支持，这个方案虽然新但也不失稳定性。
2. Visual Studio Enterprise 2017安装
Visual Studio Enterprise 2017是Microsoft公司开发工具包系列产品，我们想要使用它的集成开发环境(IDE)。
安装的VS版本为Visual Studio Enterprise 2017，需要激活（在网上容易搜到产品秘钥进行激活），
下载地址为：https://visualstudio.microsoft.com/zh-hans/downloads/
安装过程中，遇到以下选择工作负载的类似界面时，选择“Windows（3）”中的“使用C++的桌面开发”，确认后继续接下来的步骤，完成安装即可
安装结束可以在windows菜单栏找到VS
3. Intel Parallel Studio XE 2018安装
Parallel Studio XE是Intel公司推出的一款IDE软件，该程序具有强大的多线程编程功能，能够提供高效的Fortran编译支持。
安装的版本为：parallel_studio_xe_2018_update1_cluster_edition。PSXE可以通过高校邮箱申请许可的序列码，具体步骤按照Intel官网的步骤操作即可，
地址为：https://software.intel.com/en-us/parallel-studio-xe/choose-download
开始安装PSXE，下面的步骤与帖子（感谢Springer_）基本一致：
将申请的序列码填入：
安装结束后，在windows菜单栏可以找到PSXE，如下图
可以看到，PSXE提供了两个DOS窗口：for IA-32 VS以及for Intel 64 VS，使用起来没有发现区别，我们下面测试JPL提供的DE405例程就从该窗口进入。
4. JPL DE405的下载与测试
JPL DE405是我们对JPL的行星星历工具包的称呼，它包括一定时间跨度的行星运动状态数据，以及调用的程序，其中最为核心的星历数据保存在DExxx.xxx（比如，比较常用的DE405.405，405是数据版本的编号）文件中。为了方便，我们用星历文件名来代表工具包。
JPL DE405工具包下载的官网为：https://ssd.jpl.nasa.gov/?planet_eph_export
星历DE405的下载地址为：ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/
DE405转换为二进制存储以及读取DE405的程序的下载地址为：ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran
程序的说明文件下载的地址为：ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/userguide.txt
至此，我们得到了所有需要的文件：
文件名ascYXXXX.405表示公元后(p)XXXX年的数据，一个文件存储了20年内的数据，读者根据自己的需要下载相应时间跨度的文件；asc2eph.f是ascii转二进制的程序文件；testeph.f是读取二进制文件并且相关计算的程序文件；userguide.txt是说明文件。
在windows开始菜单，找到Intel Parallel Studio XE 2018，打开compiler 18.0 for Intel 64 VS的DOS窗口对JPL提供的Fortran原码进行测试：
首先操作进入存放JPL DE405工具包的文件夹内，合并文件：
copy header.405+ascp2000.405+ascp2020.405+ascp2040.405 infile.405
此命令会产生一个infile.405的文件，继续在命令行界面输入：
ifort asc2eph.f
ifort是PSXE对应的表示编译的指令，通过asc2eph.f程序将infile读取（<）到控制台：
asc2eph < infile.405
此命令会产生一个名为JPLEPH的文件，并有执行结果：
打开testeph.f文件，给所有的NRECL赋值为1（对应windows操作系统），FSIZER3中的KSIZE赋值为2036，并且取消子程序STATE中的调用FSIZER3的注释符，编译testeph.f文件：
ifort testeph.f
通过testeph.f程序将testpo.405文件读取到控制台：
testeph < testpo.405
正常情况下，控制台会输出一系列的常用参数，以及和标准数据的比较结果，至此，测试成功。
5. IAU SOFA的下载
因为JPL提供的查询行星位置和速度的程序输入为儒略日（Julian Date），这与我们习惯使用的记时法格里高利历（Gregorian calendar）不同，所以我们需要将GC转化为JD，幸运的是，IAU SOFA提供了Fortran的程序：cal2jd.f
程序的下载地址为：http://www.iausofa.org/index.html
值得注意的是，cal2jd.f输出为简约儒略日（MJD），但是可以通过简单的修改（+2400000.5）得到JD。在上述的网站，IAU SOFA还提供了许多有用的程序，比如epv00.f，它可以直接返回地球相对于太阳的速度与位置（稍后我们就采用了这个程序与JPL程序进行对比）。
6. 基于VS项目管理的Fortran77程序
在VS下新建Fortran工程，程序的结构为：
其中：
- JPLEPH：存储星历的二进制文件
- cal2jd.for：格里高利历转化为儒略历
- const.for：查找星历文件里所有的常数
- epv00.for：计算地球相对于太阳的位置和速度（对比程序）
- FSIZER3.for：初始化
- Get_ephemeris.for：主程序
- INTERP.for：微分并插值Chebyshev系数
- PLEPH.for：计算目标星相对于参考星的位置和速度
- SPLIT.for：将双精度的实数拆分为整数与小数部分
- STATE.for：读取并插值星历文件
除了主程序Get_ephemeris.for，cal2jd.for与epv00.for来源于IAU SOFA原文件，其余的子程序均来自testeph.f文件。
主程序的代码为：
program Get_ephemeris
implicit none
C++++++++++++++++++++++++++++
C For routine iau_cal2jd(year,month,day,djm0,djm,djm_status)
C year,month,day integer year, month, day in Gregorian calendar
C djm0 d MJD zero-point: always 2400000.5
C djm d Julian Date for 0 hrs
C djm_status i status:
INTEGER(4) year,month,day,djm_status
DOUBLE PRECISION djm0,djm
C++++++++++++++++++++++++++++
C For routine iau_EPV00(djm1,djm2,pvh,pvb,jstat)
C djm1 d TDB date part A (Note 1)
C djm2 d TDB date part B (Note 1)
C pvh d(3,2) heliocentric Earth position/velocity (au,au/day)
C pvb d(3,2) barycentric Earth position/velocity (au,au/day)
C jstat i status: 0 = OK
C +1 = warning: date outside 1900-2100 AD
INTEGER jstat
DOUBLE PRECISION djm1,djm2
DIMENSION pvh(3,2),pvb(3,2)
DOUBLE PRECISION pvh,pvb
C++++++++++++++++++++++++++++
C For routine PLEPH(djm3,ntarg,ncent,rrd)
C ntarg = INTEGER NUMBER OF 'TARGET' POINT.
C
C ncent = INTEGER NUMBER OF CENTER POINT.
C
C THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS:
C
C 1 = MERCURY 8 = NEPTUNE
C 2 = VENUS 9 = PLUTO
C 3 = EARTH 10 = MOON
C 4 = MARS 11 = SUN
C 5 = JUPITER 12 = SOLAR-SYSTEM BARYCENTER
C 6 = SATURN 13 = EARTH-MOON BARYCENTER
C 7 = URANUS 14 = NUTATIONS (LONGITUDE AND OBLIQ)
C 15 = LIBRATIONS, IF ON EPH FILE
C rrd = OUTPUT 6-word D.P. ARRAY CONTaiNING POSITION AND VELOCITY
C OF POINT 'NTARG' RELATIVE TO 'NCENT'.
DOUBLE PRECISION djm3
DIMENSION rrd(6)
DOUBLE PRECISION rrd
year = 2011
month = 5
day = 1
call iau_cal2jd(year,month,day,djm0,djm,djm_status)
djm1 = djm
djm2 = 0D0
call iau_EPV00(djm1,djm2,pvh,pvb,jstat)
djm3 = djm
call PLEPH(djm3,3,11,rrd)
write(*,*) "查询地球相对于太阳位置速度的时间为（年/月/日）："
write(*,*) year,month,day
write(*,100)
write(*,*) "地球相对于太阳的位置和速度为(au,au/y) (by IAU EPV00):"
write(*,*) pvh(1,1),pvh(2,1),pvh(3,1),pvh(1,2),pvh(2,2),pvh(3,2)
write(*,100)
write(*,*) "地球相对于太阳的位置和速度为(au,au/y) (by JPL PLEPH):"
write(*,*) rrd(1),rrd(2),rrd(3),rrd(4),rrd(5),rrd(6)
100 format(/)
stop
end程序的运行结果为：
可以看到，通过IAU epv00与JPL PLEPH计算的结果基本一致，因此结果是可信的。
相关阅读
在研究用户需求上没有什么捷径可以走，不要以为自己可以想当然地猜测用户习惯产品研发中心最容易犯的一个错误是：研发者往往对自己挖
你手里有5万、50万、500万，创业都不一定成功。现在只有一千块，怎么创业?看看他们都是怎么回答的。彻底否决型来自山东的创业者罗文
短短7天，相互保已经吸引了超过1000万的用户加入，可以肯定的是，相互保的火热绝对会给整个行业带来巨大的改变。那么，本文就和大家一起
世界瞬息万变，有目的的试验比有目的的实践更重要。故事要从爱迪生谈起，绝大多数人都会认为是爱迪生发明了电灯泡。其实，这是一个天大
看起来是打包了一堆权益，其实最值钱的商品是会员自己。如果你是个淘宝剁手党，最近可能在手淘注意到 88 VIP 会员的促销 banner。这