;-------------------------------------------------------------
;+
; NAME:
;       JD2YMD
; PURPOSE:
;       Find year, month, day from julian day number.
; CATEGORY:
; CALLING SEQUENCE:
;       jd2ymd, jd, y, m, d
; INPUTS:
;       jd = Julian day number (like 2447000).     in
; KEYWORD PARAMETERS:
; OUTPUTS:
;       y = year (like 1987).                      out
;       m = month number (like 7).                 out
;       d = day of month (like 23).                out
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;       R. Sterner.  21 Aug, 1986.
;       Johns Hopkins Applied Physics Lab.
;       RES 18 Sep, 1989 --- converted to SUN
;       R. Sterner, 30 Apr, 1993 --- cleaned up and allowed arrays.
;
; Copyright (C) 1986, Johns Hopkins University/Applied Physics Laboratory
; This software may be used, copied, or redistributed as long as it is not
; sold and this copyright notice is reproduced on each copy made.  This
; routine is provided as is without any express or implied warranties
; whatsoever.  Other limitations apply as described in the file disclaimer.txt.
;-
;-------------------------------------------------------------
 
	pro jd2ymd, jd, y, m, d, help=hlp
 
	if (n_params(0) lt 4) or keyword_set(hlp) then begin
	  print,' Find year, month, day from julian day number.'
	  print,' jd2ymd, jd, y, m, d'
	  print,'   jd = Julian day number (like 2447000).     in'
	  print,'   y = year (like 1987).                      out'
	  print,'   m = month number (like 7).                 out'
	  print,'   d = day of month (like 23).                out'
	  return
	endif
 
	y = fix((jd - 1721029)/365.25)		; Estimated year.
	jd0 = ymd2jd(y, 1, 0)			; JD for day 0.
	days = jd - jd0				; Day of year.
	w = where(days le 0, cnt)		; Find where year is wrong.
	if cnt gt 0 then begin
	  y(w) = y(w) - 1			; Year was off by 1.
	  jd0(w) = ymd2jd( y(w), 1, 0)		; New JD for day 0.
	  days(w) = jd(w) - jd0(w)		; New day of year.
	endif
 
	;---  Correct for leap-years.  -----
	ly = (((y mod 4) eq 0) and ((y mod 100) ne 0)) $
            or ((y mod 400) eq 0)
 
	;---  Days before start of each month.  -----
	ydays = [0,0,31,59,90,120,151,181,212,243,273,304,334,366]
	off   = [0,0, 0, 1, 1,  1,  1,  1,  1,  1,  1,  1,  1,  1]
 
	;----------------  Find which month.  --------------------------
	;     Algorithm: ydays has cumulative # days up to start of each month
	;     (13 elements so month number may be used as an index).
	;     This number needs 1 added for Mar to Dec if it is a leap year.
	;     This is done by adding the offset, off, times the leap year flag.
	;     The larger the day of year (days) the fewer elements of this
	;     corrected ydays array will be greater than days.  The 
	;     entries in the corrected ydays gt days are located by where and
	;     counted by n_elements.  Ydays has 13 elements so subtract result
	;     from 13 to get month number.
	;---------------------------------------------------------------
	njd = n_elements(jd)	; # of JDs to convert.
	m = intarr(njd)		; Set up storage for months.
	d = intarr(njd)		; Set up storage for day of month.
	for i = 0, njd-1 do begin	; Loop through each JD.
	  ydays2 = ydays+ly(i)*off	; Correct cumulative days for year.
	  dy = days(i)			; Days into year for i'th JD.
	  mn = 13-n_elements(where(dy le ydays2))  ; i'th month number.
	  m(i) = mn			; Store month.
	  d(i) = fix(dy - ydays2(mn))	; Find and store i'th day of month.
	endfor
 
	;---------  Make sure scalars are returned as scalars  -------
	if n_elements(m) eq 1 then begin
	  m = m(0)
	  d = d(0)
	endif
 
	return
	end