Fortran计算月食初亏——以2017年第一次月食为例

本程序示例是计算2017年第一次月食初亏时刻，程序从2017年8月7日0点开始迭代。运行环境：Fortran95，VS2010，WIN764bit

原理示意图如下。所有矢量都以地球为中心。S为太阳，E为地球，M为月球，O为地球本影影锥。当角ANGLE_EO+角ANGLE_MO=角ANGLE_EOM时，开始发生月食。

流程图如下：

源码包含文件：

1、testeph.f：JPL提供的官方程序，包括获取目标天体位置速度等子程序

2、compute.f：JPL其他官方程序中的子程序组合而成

3、Source1.f90：主程序

4、JPLEPH：JPL提供的DE421历表数据文件，原版为ASCII码版本，已转为二进制形式

子程序说明（只解释主程序中用到的子程序）：

1、testeph.f

• SUBROUTINE PLEPH ( ET, NTARG, NCENT, RRD )

功能：给定中心天体和目标天体，得到目标天体相对于中心天体的位置和速度

输入：ET代表时刻，儒略日，double型

   NTARG代表目标天体编号，NCENT代表中心天体编号，整型

输出： RRD代表返回参数，6元素数组，double型，前三个为位置分量（单位AU），后三个为速度分量（单位AU/s）

  天体编号：

 1-水星     2-金星       3-地球4-火星5-木星6-土星7-天王星8-海王星 9-冥王星10-月球11-太阳12-太阳系质心

13-地月系质心 14-章动 15-天平动

2、compute.f

• SUBROUTINE iau_CAL2JD ( IY, IM, ID, DJM0, DJM, J )

功能：从格里历到儒略历

输入：IY、IM、ID，年月日，整型

输出：DJM0，DJM，double型。其中DJM0固定为2400000.5，DJM0+DJM=儒略日

J为输出状态，0为正常，-1为错误

•  SUBROUTINE iau_JD2CAL ( DJ1, DJ2, IY, IM, ID, FD, J )

功能：儒略日到格里历

输入：DJ1，DJ2分别为儒略日的两部分，double型

输出：IY、IM、ID，年月日，整型

    FD日的小数部分，double型

   J为输出状态，0为正常，-1为错误

•  SUBROUTINE iau_PDP ( A, B, ADB )

功能：两个向量做点乘

输入：A，B，两个向量，3元素数组，double型

输出：ADB点乘结果，double型

•  SUBROUTINE sla_DD2TF (NDP, DAYS, SIGN, IHMSF)

功能：将小数天转为时分秒

输入：NDP，秒的小数点后位数

   DAYS，小数形式的天，double型

  SIGN：+或-号，char型

  IHMSF，4元素数组，整型，分别为时、分、秒、秒的小数部分

• subroutine selconQ(nams,vals)

功能：选择常数

输入：nams，常数编号，6元素char型

输出：vals，常数值，double型，单位Km

常数编号：AU-天文单位 ASUN-太阳半径 RE-地球半径 AM-月球半径 CLIGHT-光速

【附】源码地址：https://github.com/5663015/compute-eclipse-of-the-moon

    DE421历表数据文件：http://pan.baidu.com/s/1dEVTOqt 密码：xypi