pro mphase,jd, k
;+
; NAME:
; MPHASE
; PURPOSE:
; Return the illuminated fraction of the Moon at given Julian date(s)
;
; CALLING SEQUENCE:
; MPHASE, jd, k
; INPUT:
; JD - Julian date, scalar or vector, double precision recommended
; OUTPUT:
; k - illuminated fraction of Moon's disk (0.0 < k < 1.0), same number
; of elements as jd. k = 0 indicates a new moon, while k = 1 for
; a full moon.
; EXAMPLE:
; Plot the illuminated fraction of the moon for every day in July
; 1996 at 0 TD (~Greenwich noon).
;
; IDL> jdcnv, 1996, 7, 1, 0, jd ;Get Julian date of July 1
; IDL> mphase, jd+dindgen(31), k ;Moon phase for all 31 days
; IDL> plot, indgen(31),k ;Plot phase vs. July day number
;
; METHOD:
; Algorithm from Chapter 46 of "Astronomical Algorithms" by Jean Meeus
; (Willmann-Bell, Richmond) 1991. SUNPOS and MOONPOS are used to get
; positions of the Sun and the Moon (and the Moon distance). The
; selenocentric elongation of the Earth from the Sun (phase angle)
; is then computed, and used to determine the illuminated fraction.
; PROCEDURES CALLED:
; MOONPOS, SUNPOS
; REVISION HISTORY:
; Written W. Landsman Hughes STX June 1996
; Converted to IDL V5.0 W. Landsman September 1997
; Use /RADIAN keywords to MOONPOS, SUNPOS internally W. Landsman Aug 2000
;-
On_error,2
if N_params() LT 2 then begin
print,'Syntax - MPHASE, jd, k'
return
endif
diss = 1.49598e8 ;Earth-Sun distance (1 AU)
moonpos, jd, ram, decm, dism, /RADIAN
sunpos, jd, ras, decs, /RADIAN
; phi - geocentric elongation of the Moon from the Sun
; inc - selenocentric (Moon centered) elongation of the Earth from the Sun
phi = acos( sin(decs)*sin(decm) + cos(decs)*cos(decm)*cos(ras-ram) )
inc = atan( diss * sin(phi), dism - diss*cos(phi) )
k = (1 + cos(inc))/2.
return
end