!  "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 " .