The purpose of the program is to compute the initial eclipse time of the first lunar eclipse in 2019, when the moon just entered the umbra, based on the DE421 and SOFA.
The coder wrote the main function in C language and Fortran language respectively. they call the same library functions but have different makefile configuration files. Please chone or download the version your need as follows:
- Fortran：1，3，4，5，7，8
- C：1，2，4，5，6，8
| NO. | FILENAME | DESCRIPTION |
|---|---|---|
| 1 | JPLEPH | the binary JPL ephemeris file. |
| 2 | MOON_ECLIPSE.c | the main function file in C language. |
| 3 | MOON_ECLIPSE.f95 | the main function file in Fortran language. |
| 4 | README.md | readme file. |
| 5 | SOFALIB.f | SOFA subroutines used in this project. |
| 6 | makefile_C | C version makefile (the file directs make on how to compile and link a program). |
| 7 | makefile_F | Fortran version makefile (the file directs make on how to compile and link a program). |
| 8 | selcon.f | the select constant file in the DE ephemeris. |
Please check that gfortran and gcc are installed before running. If not, refer to 'here' please.
- Fortran Version
$ make -f makefile_F
$ ./MOON_F.exe
- C Version
$ make -f makefile_C
$ ./MOON_C.exe
- The figure below shows the spatial relationship between the earth, the moon and the sun before the eclipse,
Figure 1. Schematic diagram of the spatial relationship between the earth, the moon and the sun at the time before the eclipse
Where S, E and M are the sun, the earth and the moon; O is the earth's umbra cone; ∠EOM is the Angle between the vectors EO and MO; ∠OE and ∠OM are the angular radii of the earth and moon relative to the shadow cone, respectively. When an eclipse occurs, ∠EOM is equal to ∠OE plus ∠OM.
- Since it takes time for the sun's rays to leave the sun and pass through the earth to form a shadow cone, it is necessary to calculate the light travel time. The figure below is a schematic diagram for calculating the light line.
Figure 2. Iteratively solve for light travel time
Steps as follows:
(1). Calculate the distance L0 from the sun to the earth according to the JPEPH of given time T0, then compute the light trivel time DT0.
(2). Calculate the new site of T0+DT0 based on the velocity of the object being measured,re-calculate the distance L1 from the sun to the earth,compute the light trivel time DT1.
(3). Calculate the new site of T0+DT1 based on the velocity of the object being measured,re-calculate the distance L2 from the sun to the earth,compute the light trivel time DT2.
(4). When the |DTn - DTn-1| approximately equal zero, We get the precise light trivial time DTn (n > = 1), and then we can calculate the measured body's location accurately.
- Strategies used to speed up the computation process
Figure 3. Locate the most likely day
-- Since the eclipse occurs in a straight line from the date to the earth and the moon, the Angle between the sun and the earth and the moon is close to 180 degrees, the first step is to take a day as the step length, filter the time when the Angle is too large (the time when the eclipse is impossible), and locate the day when the eclipse occurs.
--The distance between the shadow cone and the earth is equal to about 10000 earth's radius, the moon to the earth's distance is equal to about 60 radius of the earth, When an eclipse occurs the Angle between moon-earth and earth-cone less than 1 °
-- After locating the most likely time, You can continue to approach the eclipse at smaller intervals, such as half a day.
-
The procedure flow is as follows:
Figure 4. The flow chart
- Read astronomical units "AU" light speed "CLIGHT" and other constants from the binary almanac file JPLEPH.
- Converte the Gregorian calendar to the Julian calendar.
- Iterate to the day when an eclipse is likely to occur
- Iterate over the light travel time.
- Calculate the vector STE from the sun to the earth.
- Calculate the vector ETO from the earth to the shadow cone according to the similar triangle principle.
- Calculate the MTE(vector from the moon to the earth).
- Calculate the MTO according to the vector addition rule.
- Calculate the Angle between vector ETO and MTO according to cosine theorem.
- Calculate the angular radius of the shadow cone relative to the earth and the moon according to the principle of inverse trigonometric function.
- Compute ERR, ERR = ∠EOM - ∠OM - ∠ OE.
- Judge ERR < = 0? , if yes, jump to 13; if not, time +1 then jump to 4.
- Converte the Julian calendar to the Gregorian calendar.
- Converte the decimal part to UTC time.
- Output the results.
-
The calculation results of the program need to be referred to. Timeanddata of stavanger, Norway gives the time of the first lunar eclipse in 2019 as shown in the figure below:
Source of the image:https://www.timeanddate.com/eclipse/lunar/2019-january-21
Figure 5. Reference time
- The following is the result of running the program,it is in good agreement with Figure 5.
Figure 6. The result of the program
The subroutines invoked by the program are listed below:
SUBROUTINE iau_CAL2JD ( IY, IM, ID, DJM0, DJM, J )- 所在文件：cal2jd.for
- 功能说明：格里高历转为儒略历
- 参数说明：
Given:
IY,IM,ID i year, month, day in Gregorian calendar (Note 1)
Returned:
DJM0 d MJD zero-point: always 2400000.5
DJM d Modified Julian Date for 0 hrs
J i status:
0 = OK
-1 = bad year (Note 3: JD not computed)
-2 = bad month (JD not computed)
-3 = bad day (JD computed)
SUBROUTINE iau_JD2CAL ( DJ1, DJ2, IY, IM, ID, FD, J )- 所在文件：jd2cal.for
- 功能说明：儒略历转为格里高历
- 参数说明：
Given:
DJ1,DJ2 d Julian Date (Notes 1, 2)
Returned:
IY i year
IM i month
ID i day
FD d fraction of day
J i status:
0 = OK
-1 = unacceptable date (Note 1)
SUBROUTINE iau_PDP ( A, B, ADB )- 所在文件：pdp.for
- 功能说明：向量做点积
- 参数说明：
Given:
A d(3) first p-vector
B d(3) second p-vector
Returned:
ADB d A . B
SUBROUTINE sla_DD2TF (NDP, DAYS, SIGN, IHMSF)- 所在文件：selcon.f
- 功能说明：天到时分秒
- 参数说明：
Given:
NDP int number of decimal places of seconds
DAYS dp interval in days
Returned:
SIGN char '+' or '-'
IHMSF int(4) hours, minutes, seconds, fraction
subroutine selconQ(nams,vals)- 所在文件：selcon.f
- 功能说明：常量提取
- 参数说明：
Input
nams : names (character*6)
Output
vals : values corresponding to the input names
QI-ZHAOXIANG@20130402
SUBROUTINE PLEPH ( ET, NTARG, NCENT, RRD )- 所在文件：selcon.f
- 功能说明：返回给定时间中心天体到目标天体的矢量
- 参数说明：
THIS SUBROUTINE READS THE JPL PLANETARY EPHEMERIS
AND GIVES THE POSITION AND VELOCITY OF THE POINT 'NTARG'
WITH RESPECT TO 'NCENT'.
CALLING SEQUENCE PARAMETERS:
ET = D.P. JULIAN EPHEMERIS DATE AT WHICH INTERPOLATION
IS WANTED.
** NOTE THE ENTRY DPLEPH FOR A DOUBLY-DIMENSIONED TIME **
THE REASON FOR THIS OPTION IS DISCUSSED IN THE
SUBROUTINE STATE
NTARG = INTEGER NUMBER OF 'TARGET' POINT.
NCENT = INTEGER NUMBER OF CENTER POINT.
THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS:
1 = MERCURY 8 = NEPTUNE
2 = VENUS 9 = PLUTO
3 = EARTH 10 = MOON
4 = MARS 11 = SUN
5 = JUPITER 12 = SOLAR-SYSTEM BARYCENTER
6 = SATURN 13 = EARTH-MOON BARYCENTER
7 = URANUS 14 = NUTATIONS (LONGITUDE AND OBLIQ)
15 = LIBRATIONS, IF ON EPH FILE
(IF NUTATIONS ARE WANTED, SET NTARG = 14. FOR LIBRATIONS,
SET NTARG = 15. SET NCENT=0.)
RRD = OUTPUT 6-WORD D.P. ARRAY CONTAINING POSITION AND VELOCITY
OF POINT 'NTARG' RELATIVE TO 'NCENT'. THE UNITS ARE AU AND
AU/DAY. FOR LIBRATIONS THE UNITS ARE RADIANS AND RADIANS
PER DAY. IN THE CASE OF NUTATIONS THE FIRST FOUR WORDS OF
RRD WILL BE SET TO NUTATIONS AND RATES, HAVING UNITS OF
RADIANS AND RADIANS/DAY.
The option is available to have the units in km and km/sec.
For this, set km=.true. in the STCOMX common block.
this project：https://github.com/Sardingfish/MOON_ECLIPSE
SOFA：http://www.iausofa.org/index.html
DE/LE EPHEMERIS：ftp://ssd.jpl.nasa.gov
This project is licensed under the terms of the MIT license.