blob: b4fec40535c8fcc1f0a7742a0330888fca44cfc8 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
|
! "flmoon" function to compute the phases of the moon
!
! REFERENCE: "Numerical recipes in FORTRAN,
! The Art of Scientific Computing, Second Edition"
! Press, Teukolsky, Vetterling, Flannery
! Cambridge University Press
! ISBN 0-521-43064-X
! PAGES: 1-2 ("SUBROUTINE flmoon")
!
! INPUT: "n" the n-th such phase since January, 1900
! "nph" the phase desired
! 0: new moon
! 1: first quarter
! 2: full moon
! 3: last quarter
! OUTPUT: "jd" the Julian day number
! "frac" the fractional part of day to be added to it
INTEGER jd n nph ;
REAL frac ;
REAL RAD := $Pi_R 180. / ; ! NOTE: $Pi_R is a parametric constant
INTEGER i ;
REAL am as c t t2 xtra ;
:: >>n<< >>nph<< ;
EVALUATE c := n I_TO_R nph I_TO_R 4. / + ;
EVALUATE t := c 1236.85 / ;
EVALUATE t2 := t t * ;
EVALUATE as := 359.2242 29.105356 c * + ;
EVALUATE am := 306.0253 385.816918 c * 0.010730 t2 * + + ;
EVALUATE jd := 2415020 28 n * 7 nph * + + ;
EVALUATE xtra := 0.75933 1.53058868 c * +
1.178E-4 1.55E-7 t * - t2 * + ;
IF nph 0 = nph 2 = + THEN
EVALUATE xtra := xtra
0.1734 3.93E-4 t * - RAD as * SIN *
0.4068 RAD am * SIN * - + ;
ELSEIF nph 1 = nph 3 = + THEN
EVALUATE xtra := xtra
0.1721 4.E-4 t * - RAD as * SIN *
0.6280 RAD am * SIN * - + ;
ELSE
ECHO "*nph* is unknown in *flmoon*" ;
ENDIF ;
IF xtra 0. >= THEN
EVALUATE i := xtra R_TO_I ;
ELSE
EVALUATE i := xtra 1. - R_TO_I ;
ENDIF ;
EVALUATE jd := jd i + ;
EVALUATE frac := xtra i I_TO_R - ;
:: <<jd>> <<frac>> ;
QUIT " Routine *flmoon* XREF " .
|
