; --------------------------------------------------------------------

FUNCTION n_lines,file

; This function opens the specified file, determines its length in
; bytes, and reads it in as a BYTE-valued array.  It then returns
; the number of occurrences of 10B in the file.  10B is the ASCII
; code for a line feed, which marks the end of a line of text (in
; UNIX, at least).

; Modified 10 January 1995 by D P Steele:
; If called under VMS, this function dodges the whole issue of
; counting lines and SPAWNs a call to WC.COM, which uses EVE (the
; Extensible Vax Editor) to count the number of lines and return
; the result.  Obviously, WC.COM must be accessible to the function.

IF !VERSION.OS EQ 'vms' THEN BEGIN
    SPAWN,'@USER$:[STEELE.IDL]WC '+file,n
    n=LONG(n(0))
    ENDIF ELSE BEGIN
    OPENR,unit,file,/GET_LUN
    f=FSTAT(unit)
    IF f.SIZE EQ 0 THEN n=0L ELSE BEGIN
        b=BYTARR(f.size)
        READU,unit,b
        n=LONG(TOTAL(b EQ 10B))
        ENDELSE
    CLOSE,unit
    FREE_LUN,unit
    ENDELSE

RETURN,n
END

;-------------------------------------------------------------

;+
; NAME:
;       QGETFILE
; PURPOSE:
;       Read a text file into a string array.
; CATEGORY:
; CALLING SEQUENCE:
;       s = qgetfile(f)
; INPUTS:
;       f = text file name.      in
; KEYWORD PARAMETERS:
;       Keywords:
;         ERROR=err  error flag: 0=ok, 1=file not opened.
;         /QUIET means give no error message.
; OUTPUTS:
;       s = string array.        out
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;       R. Sterner, 20 Mar, 1990
;	D. Steele, 18 Mar, 1994 - added line count, changed WHILE loop
;				  to FOR loop, and changed name from
;				  GETFILE.PRO to QGETFILE.PRO.
;
; Copyright (C) 1990, 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.
;-
;-------------------------------------------------------------

	function qgetfile, file, error=err, help=hlp, quiet=quiet

	if (n_params(0) lt 1) or keyword_set(hlp) then begin
	  print,' Read a text file into a string array.'
	  print,' s = getfile(f)'
	  print,'   f = text file name.      in'
	  print,'   s = string array.        out'
	  print,' Keywords:'
	  print,'   ERROR=err  error flag: 0=ok, 1=file not opened.'
	  print,'   /QUIET means give no error message.'
	  return, -1
	endif

	get_lun, lun
	on_ioerror, err

; The following lines were added by D P Steele, 94/03/18
	nl=n_lines(file)
	s=STRARR(nl)

	openr, lun, file

;	s = [' ']
	t = ''

;	while not eof(lun) do begin
	FOR i=0L,nl-1 DO BEGIN		; replace EOF search with FOR loop
	  readf, lun, t
;	  s = [s,t]
	  s(i)=t			; avoid slow string concatenation
;	endwhile
	  ENDFOR

	close, lun
	free_lun, lun
	err = 0
;	return, s(1:*)
	RETURN,s

err:	if !err eq -168 then begin
	  if not keyword_set(quiet) then print,' Non-standard text file format.'
	  free_lun, lun
;	  return, s(1:*)
	  RETURN,s
	endif
	if not keyword_set(quiet) then print,$
	  ' Error in getfile: File '+file+' not opened.'
	free_lun, lun
	err = 1
	return, -1

	end

; --------------------------------------------------------------------

;+
; NAME:
;       MPUTFILE
; PURPOSE:
;       Write a text file from a string array.
; CATEGORY:
; CALLING SEQUENCE:
;       mputfile, f, s
; INPUTS:
;       f = text file name.      in
;       s = string array.        in
; KEYWORD PARAMETERS:
;       Keywords:
;         ERROR=err  error flag: 0=ok, 1=invalid string array.
; OUTPUTS:
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;       R. Sterner, 20 Mar, 1990
;       R. Sterner,  4 Nov, 1992 --- allowed scalar strings.
;
; Copyright (C) 1990, 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 mputfile, file, s, error=err, help=hlp

	if (n_params(0) lt 1) or keyword_set(hlp) then begin
	  print,' Write a text file from a string array.'
	  print,' mputfile, f, s'
	  print,'   f = text file name.      in'
	  print,'   s = string array.        in'
	  print,' Keywords:'
	  print,'   ERROR=err  error flag: 0=ok, 1=invalid string array.'
	  return
	endif

	if datatype(s) ne 'STR' then begin
	  print,' Error in mputfile: argument must be a string array.'
	  err = 1
	  return
	endif

	get_lun, lun
	openw, lun, file

	for i = 0L, n_elements(s)-1 do begin
	  t = s(i)
	  if t eq '' then t = ' '
	  printf, lun, t
	endfor

	close, lun
	free_lun, lun
	err = 0
	return

	end

; --------------------------------------------------------------------

;+
; NAME:
;       MONTHNAMES
; PURPOSE:
;       Returns a string array of month names.
; CATEGORY:
; CALLING SEQUENCE:
;       mnam = monthnames()
; INPUTS:
; KEYWORD PARAMETERS:
; OUTPUTS:
;       mnam = string array of 13 items:     out
;         ['Error','January',...'December']
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;       R. Sterner, 18 Sep, 1989
;
; Copyright (C) 1989, 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.
;-
;-------------------------------------------------------------

	function monthnames, help=hlp

	if keyword_set(hlp) then begin
	  print,' Returns a string array of month names.'
	  print,' mnam = monthnames()'
	  print,'   mnam = string array of 13 items:     out'
	  print,"     ['Error','January',...'December']"
	  return, -1
	endif

	mn = ['Error','January','February','March','April','May',$
	      'June','July','August','September','October',$
	      'November','December']

	return, mn
	end

; --------------------------------------------------------------------

;+
; NAME:
;       YMD2JD
; PURPOSE:
;       From Year, Month, and Day compute Julian Day number.
; CATEGORY:
; CALLING SEQUENCE:
;       jd = ymd2jd(y,m,d)
; INPUTS:
;       y = Year (like 1987).                    in
;       m = month (like 7 for July).             in
;       d = month day (like 23).                 in
; KEYWORD PARAMETERS:
; OUTPUTS:
;       jd = Julian Day number (like 2447000).   out
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;       R. Sterner,  23 June, 1985 --- converted from FORTRAN.
;       Johns Hopkins University Applied Physics Laboratory.
;       RES 18 Sep, 1989 --- converted to SUN
;
; Copyright (C) 1985, 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.
;-
;-------------------------------------------------------------

	function ymd2jd, iy, im, id, help=hlp

	if (n_params(0) LT 3) or keyword_set(hlp) then begin
	  print,' From Year, Month, and Day compute Julian Day number.'
	  print,' jd = ymd2jd(y,m,d)'
	  print,'   y = Year (like 1987).                    in'
	  print,'   m = month (like 7 for July).             in'
	  print,'   d = month day (like 23).                 in'
	  print,'   jd = Julian Day number (like 2447000).   out'
	  return, -1
	endif

	y = long(iy)
	m = long(im)
	d = long(id)
	jd = 367*y-7*(y+(m+9)/12)/4-3*((y+(m-9)/7)/100+1)/4 $
             +275*m/9+d+1721029

	return, jd

	end

; --------------------------------------------------------------------

;+
; NAME:
;       WEEKDAY
; PURPOSE:
;       Compute weekday given year, month, day.
; CATEGORY:
; CALLING SEQUENCE:
;       wd = weekday(y,m,d,[nwd])
; INPUTS:
;       y, m, d = Year, month, day (Like 1988, 10, 31).      in
; KEYWORD PARAMETERS:
; OUTPUTS:
;       wd = Returned name of weekday.                       out
;       nwd = optional Weekday number.                       out
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;       R. Sterner. 31 Oct, 1988.
;       Johns Hopkins University Applied Physics Laboratory.
;       RES 18 Sep, 1989 --- converted to SUN
;
; Copyright (C) 1988, 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.
;-
;-------------------------------------------------------------

	FUNCTION WEEKDAY, Y, M, D, NWD, help=hlp

	IF (N_PARAMS(0) LT 3) or keyword_set(hlp) THEN BEGIN
	  PRINT,' Compute weekday given year, month, day.'
	  PRINT,' wd = weekday(y,m,d,[nwd])'
	  PRINT,'   y, m, d = Year, month, day (Like 1988, 10, 31).      in'
	  PRINT,'   wd = Returned name of weekday.                       out'
	  PRINT,'   nwd = optional Weekday number.                       out'
	  RETURN, -1
	ENDIF

	NAMES = ['','Sunday','Monday','Tuesday',$
		 'Wednesday','Thursday','Friday','Saturday']

	NWD = ((YMD2JD(Y,M,D) + 1) MOD 7) + 1
	RETURN, NAMES(NWD)

	END

; --------------------------------------------------------------------

;+
; NAME:
;       STREP
; PURPOSE:
;       Edit a string by position. Precede, Follow, Replace, Delete.
; CATEGORY:
; CALLING SEQUENCE:
;       newstring = strep(string,cmd,p,ss,[iflg])
; INPUTS:
;       string = string to edit.                               in
;       cmd = edit command:                                    in
;         'P' = precede position p with substring ss.
;         'F' = follow position p with substring ss.
;         'R' = replace text starting at position p
;               with text from substring ss.
;         'D' = delete N characters starting at
;               position p.  The calling sequence for
;               this command is slightly different:
;               IFLG = STREP(string,'D',p,n,[iflg])
;               Where n = number of characters to delete.
;       p = character position to use.                         in
;           0 = first char.  Any number larger
;           than the string length = last char.
;       ss = substring to use.  For 'D' command                in
;            n is used instead of ss.
; KEYWORD PARAMETERS:
; OUTPUTS:
;       iflg = 0 for a successful edit,                        out
;       iflg = -1 for an error and no change to string.
;       newstring = edited string.                             out
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;       Written by R. Sterner, 27 Dec, 1984.
;       Converted to SUN 13 Aug, 1989 --- R. Sterner.
;       Johns Hopkins University Applied Physics Laboratory.
;
; Copyright (C) 1984, 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.
;-
;-------------------------------------------------------------

	FUNCTION strep,S,CMD,IP,SS,IFLG, help = h

	IF (N_PARAMS(0) LT 4) or keyword_set(h) THEN BEGIN
	  PRINT,' Edit a string by position. Precede, Follow, Replace, Delete.'
	  print,' newstring = strep(string,cmd,p,ss,[iflg])
	  print,'   string = string to edit.                               in'
	  print,'   cmd = edit command:                                    in'
	  print,"     'P' = precede position p with substring ss.
	  print,"     'F' = follow position p with substring ss.
	  print,"     'R' = replace text starting at position p
	  print,'           with text from substring ss.
	  print,"     'D' = delete N characters starting at
	  print,'           position p.  The calling sequence for
	  print,'           this command is slightly different:
	  print,"           IFLG = STREP(string,'D',p,n,[iflg])
	  print,'           Where n = number of characters to delete.
	  print,'   p = character position to use.                         in'
	  print,'       0 = first char.  Any number larger
	  print,'       than the string length = last char.
	  print,"   ss = substring to use.  For 'D' command                in"
	  print,'        n is used instead of ss.
	  print,'   iflg = 0 for a successful edit,                        out'
	  print,'   iflg = -1 for an error and no change to string.
	  print,'   newstring = edited string.                             out'
	  RETURN, -1
	ENDIF

	N = STRLEN(S) - 1
	IF N LT 0 THEN $
	  IF CMD EQ 'D' THEN RETURN, '' ELSE RETURN, SS
	P = IP>0<N
	BS = BYTE(S)
	IF (CMD NE 'D') THEN  BSS = BYTE(SS)
	JFLG = 0

	CASE CMD OF
  'P':	BEGIN
	IF P EQ 0 THEN T = [BSS,BS] ELSE $
		T = [BS(0:P-1),BSS,BS(P:*)]
	END
  'F':	BEGIN
	IF P EQ N THEN BEGIN
	  T = [BS,BSS]
	ENDIF ELSE BEGIN
	  IF IP GE 0 THEN BEGIN
	    T = [BS(0:P),BSS,BS(P+1:*)]
	  ENDIF ELSE BEGIN
	    T = [BSS, BS]
	  ENDELSE
	ENDELSE
	END
  'R':  BEGIN
	NSS = STRLEN(SS)
	IF NSS GT N - P + 1 THEN $
		T = [BS(0:P-1),BSS] $
	ELSE BEGIN $
		T = BS
		T(P) = BSS
	ENDELSE
	END
  'D':	BEGIN
	ND = SS < (N + 1)
	IF P EQ 0 THEN $
	  IF ND EQ N+1 THEN T = '' ELSE $
	    T = [BS(ND:*)]
	IF P GT 0 THEN $
	  IF ND+P-1 GE N THEN $
	    T = [BS(0:P-1)] ELSE $
	    T = [BS(0:P-1),BS(ND+P:*)]
	END
  ELSE: JFLG = -1
  ENDCASE

	IF JFLG NE 0 THEN T = S
	IF JFLG EQ 0 THEN T = STRING(T)
	IF N_PARAMS(0) EQ 5 THEN IFLG = JFLG


	RETURN,T

	END

; --------------------------------------------------------------------

;+
; NAME:
;       STRESS
; PURPOSE:
;       String edit by sub-string. Precede, Follow, Delete, Replace.
; CATEGORY:
; CALLING SEQUENCE:
;       new = stress(old,cmd,n,oldss,newss,ned)
; INPUTS:
;       old = string to edit.                               in
;       cmd = edit command:                                 in
;         'P' = precede.
;         'F' = follow.
;         'D' = delete.
;         'R' = replace.
;       n = occurrence number to process (0 = all).         in
;       oldss = reference substring.                        in
;       oldss may have any of the following forms:
;         1. s	  a single substring.
;         2. s...    start at substring s, end at end of string.
;         3. ...e    from start of string to substring e.
;         4. s...e   from subs s to subs e.
;         5. ...     entire string.
;       newss = substring to add. Not needed for 'D'        in
; KEYWORD PARAMETERS:
; OUTPUTS:
;       ned = number of occurrences actually changed.       out
;       new = resulting string after editing.               out
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;       Written by R. Sterner, 6 Jan, 1985.
;       RES --- 23 May, 1988 fixed a bug in SSTYP = 2.
;       Converted to SUN 13 Aug, 1989 --- R. Sterner. (FOR loop change).
;       Johns Hopkins University Applied Physics Laboratory.
;
; Copyright (C) 1985, 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.
;-
;-------------------------------------------------------------

	FUNCTION STRESS,STRNG,CMDX,N,OLD,NEW,NED, help = h

	if (n_params(0) lt 3) or keyword_set(h) then begin
	  print,' String edit by sub-string. Precede, Follow, Delete, Replace.'
	  print,' new = stress(old,cmd,n,oldss,newss,ned)
	  print,'   old = string to edit.                               in'
	  print,'   cmd = edit command:                                 in'
	  print,"     'P' = precede.
	  print,"     'F' = follow.
	  print,"     'D' = delete.
	  print,"     'R' = replace.
	  print,'   n = occurrence number to process (0 = all).         in'
	  print,'   oldss = reference substring.                        in'
	  print,'   oldss may have any of the following forms:
	  print,'     1. s	  a single substring.
	  print,'     2. s...    start at substring s, end at end of string.
	  print,'     3. ...e    from start of string to substring e.
	  print,'     4. s...e   from subs s to subs e.
	  print,'     5. ...     entire string.
	  print,"   newss = substring to add. Not needed for 'D'        in"
	  print,'   ned = number of occurrences actually changed.       out'
	  print,'   new = resulting string after editing.               out'
	  return, -1
	endif

	CMD = STRUPCASE(CMDX)
	PDOT = STRPOS(OLD,'...')
	SSL = STRLEN(OLD)
	SSTYP = 0
	POS1 = -1
	POS2 = -1
	RSTR = STRNG
	IF (PDOT EQ -1) THEN SSTYP = 1
;	IF ((PDOT>0) EQ SSL-3) THEN SSTYP = 2
	IF (PDOT GT 0) AND (PDOT EQ SSL-3) THEN SSTYP = 2
	IF (PDOT EQ 0) AND (SSL GT 3) THEN SSTYP = 3
	IF (PDOT GT 0) AND (PDOT LT SSL-3) THEN SSTYP = 4
	IF (PDOT EQ 0) AND (SSL EQ 3) THEN SSTYP = 5
	NED = 0		; number of occurrences actually changed.


	CASE SSTYP OF
1:	  BEGIN
	    S = OLD
	    E = ''
	  END
2:	  BEGIN
	    S = STRSUB(OLD,0,SSL-4)
	    E = ''
    	  END
3:  	  BEGIN
	    S = ''
	    E = STRSUB(OLD,3,SSL-1)
	  END
4:  	  BEGIN
	    S = STRSUB(OLD,0,PDOT-1)
	    E = STRSUB(OLD,PDOT+3,SSL-1)
	  END
5:  	  BEGIN
	    S = ''
	    E = ''
	  END
ELSE: 	  PRINT, 'ERROR IN SSTYP'
	ENDCASE


;---------------  Find substring # N start  ---------------
	POS = -1
	nfor = n>1
LOOP:
	FOR I = 1, nfor DO BEGIN
	  POS = POS + 1
	  CASE SSTYP OF
    1:      POS = STRPOS(RSTR,S,POS)
    2:      POS = STRPOS(RSTR,S,POS)
    4:      POS = STRPOS(RSTR,S,POS)
    3:      POS = STRPOS(RSTR,E,POS)
    5:      POS = 0
	  ENDCASE
  	  IF POS LT 0 THEN GOTO, DONE
	ENDFOR

;----------  Find substring # N END  ----------------
    	CASE SSTYP OF
1:  	  BEGIN
	    POS1 = POS
	    POS2 = POS + STRLEN(S) - 1
	  END
2:  	  BEGIN
	    POS1 = POS
	    POS2 = STRLEN(RSTR) - 1
	  END
3:  	  BEGIN
	    POS1 = 0
	    POS2 = POS + STRLEN(E) - 1
	  END
4:  	  BEGIN
	    POS1 = POS
	    POS2 = STRPOS(RSTR,E,POS+1)
	    IF (POS2 LT 0) THEN GOTO, DONE
	    POS2 = POS2 + STRLEN(E) - 1
	  END
5:  	  BEGIN
	    POS1 = 0
	    POS2 = STRLEN(RSTR) - 1
	  END
	ENDCASE

;------------  edit string  --------------
    	CASE CMD OF
'P':  	  BEGIN
	    RSTR = STREP(RSTR,CMD,POS1,NEW)
	    POS = POS + STRLEN(NEW)
	  END
'F':  	  BEGIN
	    RSTR = STREP(RSTR,CMD,POS2,NEW)
	    POS = POS + STRLEN(NEW)
	  END
'R':  	  BEGIN
	    RSTR = STREP(RSTR,'D',POS1,POS2-POS1+1)
	    IF (POS1 GT 0) THEN $
	      RSTR = STREP(RSTR,'F',POS1-1,NEW)
	    IF (POS1 EQ 0) THEN $
	      RSTR = STREP(RSTR,'P',0,NEW)
	    POS = POS + STRLEN(NEW) - 1
	  END
'D':  	  BEGIN
	    RSTR = STREP(RSTR,CMD,POS1,POS2-POS1+1)
	    POS = POS - 1
	  END
ELSE: 	  BEGIN
	    PRINT, 'Error in cmd'
	    RETURN,RSTR
	  END
ENDCASE

	NED = NED + 1
	IF SSTYP EQ 5 THEN RETURN,RSTR
	IF N EQ 0 THEN GOTO, LOOP

DONE:
	RETURN, RSTR
	END

; --------------------------------------------------------------------

;+
; NAME:
;       YMD2DATE
; PURPOSE:
;       Convert from year, month, day numbers to date string.
; CATEGORY:
; CALLING SEQUENCE:
;       date = ymd2date(Y,M,D)
; INPUTS:
;       y = year number (like 1986).                         in
;       m = month number (1 - 12).                           in
;       d = day of month number (1 - 31).                    in
; KEYWORD PARAMETERS:
;       Keywords:
;         FORMAT = format string.  Allows output date to be customized.
;            The following substitutions take place in the format string:
;         Y$ = 4 digit year.
;         y$ = 2 digit year.
;         N$ = full month name.
;         n$ = 3 letter month name.
;         d$ = day of month number.
;         W$ = full weekday name.
;         w$ = 3 letter week day name.
; OUTPUTS:
;       date = returned date string (like 24-May-1986).      out
; COMMON BLOCKS:
; NOTES:
;       Notes:
;         The default format string is 'd$-n$-Y$' giving 24-Sep-1989
;         Example: FORMAT='w$ N$ d$, Y$' would give 'Mon
; MODIFICATION HISTORY:
;       R. Sterner.  16 Jul, 1986.
;       RES 18 Sep, 1989 --- converted to SUN
;       R. Sterner, 28 Feb, 1991 --- modified format.
;       R. Sterner, 16 Dec, 1991 --- added space to 1 digit day.
;       Johns Hopkins University Applied Physics Laboratory.
;
; 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.
;-
;-------------------------------------------------------------

	function ymd2date, y, m, d, help=hlp, format=frmt

	if (n_params(0) lt 3) or keyword_set(hlp) then begin
	  print,' Convert from year, month, day numbers to date string.'
	  print,' date = ymd2date(Y,M,D)'
	  print,'   y = year number (like 1986).                         in'
	  print,'   m = month number (1 - 12).                           in'
	  print,'   d = day of month number (1 - 31).                    in'
	  print,'   date = returned date string (like 24-May-1986).      out'
	  print,' Keywords:'
	  print,'   FORMAT = format string.  Allows output date to be '+$
	    'customized.'
	  print,'      The following substitutions take place in the '+$
	    'format string:'
	  print,'   Y$ = 4 digit year.'
	  print,'   y$ = 2 digit year.'
	  print,'   N$ = full month name.'
	  print,'   n$ = 3 letter month name.'
	  print,'   d$ = day of month number.'
	  print,'   W$ = full weekday name.'
	  print,'   w$ = 3 letter week day name.'
	  print,' Notes:'
	  print,"   The default format string is 'd$-n$-Y$' giving 24-Sep-1989"
	  print,"   Example: FORMAT='w$ N$ d$, Y$' would give 'Mon "+$
	    "September 18, 1989'"
	  return, -1
	endif

	;---- error check -----
	w = where(y lt 0, cnt)
	if cnt gt 0 then begin
	  print,'Error in ymd2date: invalid year: ',y(w)
	  return, -1
	endif
	;-----  Handle 2 digit years  ------
	w = where(y lt 50, cnt)			; Y < 50 assumed 2xxx.
	if cnt gt 0 then y(w) = y(w) + 2000
	w = where(y lt 100, cnt)		; 50 < Y < 100 assumed 19xx.
	if cnt gt 0 then y(w) = y(w) + 1900
	w = where((m lt 1) or (m gt 12),cnt)
	if cnt gt 0 then begin
	  print,'Error in ymd2date: invalid month: ',m(w)
	  return, -1
	endif
	w = where((d lt 1) or (d gt 31), cnt)
	if cnt gt 0 then begin
	  print,'Error in ymd2date: invalid month day: ', d(w)
	  return, -1
	endif

	;-----  format string  ------
	fmt = 'd$-n$-Y$'
	if keyword_set(frmt) then fmt = frmt

	n_dates = n_elements(y)
	date_arr = strarr(n_dates)

	;-----  Loop through all dates  ---------
	for i = 0, n_dates-1 do begin

	  ;-----  Get all the allowed parts  -----
	  yu = strtrim(y(i),2)
	  yl = strtrim(fix(y(i)-100*fix(y(i)/100)),2)
	  mnames = monthnames()
	  mu = mnames(m(i))
	  ml = strmid(mu,0,3)
	  dl = strtrim(d(i),2)
	  if strlen(dl) eq 1 then dl = ' '+dl	; Add space to 1 digit day.
	  wu = weekday(y(i),m(i),d(i))
	  wl = strmid(wu,0,3)

	  ;----  Do substitutions  ------
	  date = fmt
	  date = stress(date, 'R', 0, 'Y$', yu)	; 4 digit year.
	  date = stress(date, 'R', 0, 'y$', yl)	; 2 digit year.
	  date = stress(date, 'R', 0, 'd$', dl)	; day of month.
	  date = stress(date, 'R', 0, 'N$', mu)	; Full month name.
	  date = stress(date, 'R', 0, 'n$', ml)	; 3 letter month name.
	  date = stress(date, 'R', 0, 'W$', wu)	; Full weekday name.
	  date = stress(date, 'R', 0, 'w$', wl)	; 3 letter weekday name.

	  date_arr(i) = date
	endfor

	if n_dates eq 1 then date_arr=date_arr(0)  ; Return scalars as scalars.

	return, date_arr

	end

; --------------------------------------------------------------------

;+
; NAME:
;       SETMINUS
; PURPOSE:
;       Eliminate elements from a set that are also in another set.
; CATEGORY:
; CALLING SEQUENCE:
;       c = setminus(a,b)
; INPUTS:
;       a = set (array) to process.                in
;       b = set (array) of elements to eliminate.  in
; KEYWORD PARAMETERS:
; OUTPUTS:
;       c = resulting set (array).                 out
; COMMON BLOCKS:
; NOTES:
;       Notes: works for any type arrays.
;         Will not remove the element of a one element array.
;         This is because IDL does not have empty arrays
;         (corresponding to empty sets).
; MODIFICATION HISTORY:
;       R. Sterner, 7 Aug, 1991
;
; Copyright (C) 1991, 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.
;-
;-------------------------------------------------------------

	function setminus, a, b, help=hlp

	if (n_params(0) lt 2) or keyword_set(hlp) then begin
	  print,' Eliminate elements from a set that are also in another set.'
	  print,' c = setminus(a,b)'
	  print,'   a = set (array) to process.                in'
	  print,'   b = set (array) of elements to eliminate.  in'
	  print,'   c = resulting set (array).                 out'
	  print,' Notes: works for any type arrays.'
	  print,'   Will not remove the element of a one element array.'
	  print,'   This is because IDL does not have empty arrays'
	  print,'   (corresponding to empty sets).'
	  return, -1
	endif

	c = a
	for i = 0, n_elements(b)-1 do begin
	  w = where(c ne b(i), count)
	  if count gt 0 then c = c(w)
	endfor
	return, c
	end

; --------------------------------------------------------------------

FUNCTION PoCasite,year,month,day

yr=year MOD 100
ns=N_ELEMENTS(yr)
PCsite=REPLICATE(1,ns)			; default location is Calgary

; Decide where Polar Camera was located
rabb=WHERE((yr EQ 93) AND ((month EQ 1) OR (month EQ 2)),nrabb)
IF nrabb GT 0 THEN PCsite(rabb)=2	; Rabbit Lake in Jan, Feb 1993
eur=WHERE(((yr EQ 93) AND (month GE 10)) OR (yr GE 94),neur)
IF neur GT 0 THEN PCsite(eur)=3		; Eureka since October 1993

IF ns EQ 1 THEN PCsite=PCsite(0)	; scalar inputs so scalar output

RETURN,PCsite
END

; --------------------------------------------------------------------

FUNCTION ymdh2jdf,year,month,day,hour

; This function is basically a wrapper for ymd2jd.pro from the
; JHU/APL IDL library.  It includes the specification of a (possibly
; fractional) hour of the day, and returns a (possibly fractional)
; Julian date.

ijd=ymd2jd(year,month,day)
RETURN,DOUBLE(ijd)+(hour-12.)/24.
END

; --------------------------------------------------------------------

PRO MOONPOS, dte, ra, dec, x, y, z
;+
; NAME:
;	MOONPOS
; PURPOSE:
;	To compute the RA and Dec of the Moon at a given date.
;
; CALLING SEQUENCE:
;	MOONPOS, Dte, Ra, Dec, [X, Y, Z]
;
; INPUTS:
;	Dte - The Julian date of the day (and time) in question.
;
; OUTPUTS:
;	Ra  - The right ascension of the moon at that date in radians.
;	Dec - The declination of the moon at that date in radians.
;	X, Y, Z - the geocentric equatorial inertial coordinates of the
;	          center of the Moon w.r.t. the center of the Earth (Re).
;
; SIDE EFFECTS:
;	The user should be aware that the equations used are not highly
;	accurate.    The error in right ascension will rarely exceed 0.3
;	degrees; the error in declination will rarely exceed 0.2 degrees.
;
; PROCEDURE:
;     The lunar position is computed from the low precision formulae provided
;     in the Astronomical Almanac.
;
; EXAMPLE:
;	IDL> juldate,[1982,4,6],jd    ;Get reduced julian date
;	IDL> jd = jd + 2400000.       ;Convert to full julian date
;	IDL> moonpos, jd, ra ,dec     ;Get RA and Dec of moon
;	IDL> print,adstring(ra*!RADEG,dec*!RADEG)
;		==> 11 17  6.1  +09 17 56.0
;
;	This is about 5' from the position given in the Astronomical Almanac
; MODIFICATION HISTORY:
;	Written by Michael R. Greason, STX, 31 October 1988.
;       Modified by David P. Steele, CNSR, 17 February 1994:
;       Added code to calculate and return the GEI cartesian
;       coordinates of the Moon, using formulae from the Astronomical
;	Almanac.
;-
 On_error,2
 if N_params() EQ 0 then begin
     print,'Syntax - MOONPOS, dte, ra, dec'
     print,'dte = Julian Date, Output Ra and Dec are in radians'
     return
 endif
 t = (double(dte) - 2451545.0D0) / 36525.D0
 d2r = !DPI / 180.D0
 sz = size(dte)
 scalar = sz(0) eq 0               ;Scalar input?
 if scalar then t = dblarr(1) + t
;
 theta1 = (134.9D0 + (477198.85D0 * t)) * d2r
 lamda = 218.32D0 + (481267.883D0 * t) + (6.29D0 * sin(theta1))
 theta1 = (259.2D0 - (413335.38D0 * t)) * d2r
 theta2 = (235.7D0 + (890534.23D0 * t)) * d2r
 lamda = lamda - (1.27D0 * sin(theta1)) + (0.66D0 * sin(theta2))
 theta1 = (269.9D0 + (954397.70D0 * t)) * d2r
 theta2 = (357.5D0 + (35999.05D0 * t)) * d2r
 lamda = lamda + (0.21D0 * sin(theta1)) - (0.19D0 * sin(theta2))
 theta1 = (186.6D0 + (966404.05D0 * t)) * d2r
 lamda = lamda - (0.11D0 * sin(theta1))
 neg = where( lamda lt 0.D0, Nneg)
 if ( Nneg GT 0 ) then lamda(neg) = 360.D0 + ( lamda(neg) mod 360)
 lamda = (lamda mod 360.D0) * d2r
;
 theta1 = (93.3D0 + (483202.03D0 * t)) * d2r
 theta2 = (228.2D0 + (960400.87D0 * t)) * d2r
 beta = (5.13D0 * sin(theta1)) + (0.28D0 * sin(theta2))
 theta1 = (318.3D0 + (6003.18D0 * t)) * d2r
 theta2 = (217.6D0 - (407332.20D0 * t)) * d2r
 beta = beta - (0.28D0 * sin(theta1)) - (0.17D0 * sin(theta2))
 neg = where( beta LT 0.D0, Nneg )
 if ( Nneg GT 0 ) then beta(neg) = 360.D0 + ( beta(neg) mod 360)
 beta = (beta mod 360.D0) * d2r

 theta1 = (134.9 + 477198.85D0 * t)*d2r
 theta2 = (259.2 - 413335.38D0 * t)*d2r
 pi = (0.9508 + 0.0518*COS(theta1) + 0.0095*COS(theta2))
 theta1 = (235.7 + 890534.23D0 * t)*d2r
 theta2 = (269.9 + 954397.70D0 * t)*d2r
 pi = pi + 0.0078*COS(theta1) + 0.0028*COS(theta2)

 r = 1.0D0/SIN(pi*d2r)

 l = cos(beta) * cos(lamda)
 m = (0.9175D0 * cos(beta) * sin(lamda)) - (0.3978D0 * sin(beta))
 n = (0.3978D0 * cos(beta) * sin(lamda)) + (0.9175D0 * sin(beta))

 x = r*l
 y = r*m
 z = r*n

 ra = atan(m, l)
 neg = where( ra LT 0.0D, Nneg )
 if ( Nneg GT 0 ) then ra(neg) = 2.D0*!DPI + ra(neg)
 dec = asin(n)
 if scalar then begin
    ra = ra(0) & dec = dec(0) & x = x(0) & y = y(0) & z = z(0)
 endif

 return
 end

; --------------------------------------------------------------------

;+
; NAME:
;       LMST
; PURPOSE:
;       Give local mean sidereal time.
; CATEGORY:
; CALLING SEQUENCE:
;       st = lmst( jd, ut, long)
; INPUTS:
;       jd = Julian Day (starting at noon).         in
;       ut = Universal time as fraction of day.     in
;       long = observer longitude (deg, East is +). in
; KEYWORD PARAMETERS:
; OUTPUTS:
;       st = sidereal time as fraction of day.      out
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;       R. Sterner.  16 Sep, 1986.
;       R. Sterner, 15 Jan, 1991 --- converted to V2.
;       Johns Hopkins University Applied Physics Laboratory.
;
; 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.
;-
;-------------------------------------------------------------

	FUNCTION LMST, JD, UT, LNG, help=hlp

	if (n_params(0) lt 3) or keyword_set(hlp) then begin
	  print,' Give local mean sidereal time.'
	  print,' st = lmst( jd, ut, long)'
	  print,'   jd = Julian Day (starting at noon).         in'
	  print,'   ut = Universal time as fraction of day.     in'
	  print,'   long = observer longitude (deg, East is +). in'
	  print,'   st = sidereal time as fraction of day.      out'
	  return, -1
	endif

	SOL2SID = 1.0027379093D0	; solar to sidereal rate.
	T0 = 6D0/24. + 38D0/1440. + 45.836D0/86400.	; 0 pt sidereal time.

	DLONG = -LNG/360.D0	; time diff from Greenwich.
	TDAYS = JD - 2415020.5	; days since noon Jan 0, 1986.
	T = T0 + (SOL2SID - 1.0)*TDAYS + SOL2SID*UT - DLONG
	RETURN, T - FLOOR(T)	; only want fraction.
	END

; --------------------------------------------------------------------

;+
; NAME:
;	MOONAZEL
; PURPOSE:
;	Calculation of the direction to the Moon at a given time and
;	place.
; CATEGORY:
;
; CALLING SEQUENCE:
;	MOONAZEL, YEAR, MONTH, DAY, HOUR, MINUTE, SECOND $
;		, AZIMUTH, ELEVATION [, MRA, MDEC]
; INPUTS:
;       YEAR:   year, in either 4-digit or 2-digit representation
;       MONTH:  month, a number from 1 to 12
;       DAY:    day of month
;       HOUR:   hour, in Universal Time
;       MINUTE: minute, in UT
;       SECOND: second, in UT
; OPTIONAL INPUTS:
;	None.
; OUTPUTS:
;       AZIMUTH: The angle East of geographic North of the lunar
;                meridian
;       ELEVATION: The angle up from the horizon to the Moon in the
;                  meridian specified by AZIMUTH (ELEVATION < 90 deg)
; OPTIONAL OUTPUTS:
;       SRA:    The lunar right ascension
;       SDEC:   The lunar declination
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	None.
; RESTRICTIONS:
;       The algorithm giving the lunar position is valid for an
;       uncertain period of time.  The accuracy of the position is
;       typically better than 0.3 deg in RA, and 0.2 deg in
;       declination.  The place for which the lunar direction is
;       calculated is assumed to be the location of the Polar
;	Camera, and this is inferred from the specified date.
; PROCEDURE:
;       The location of the Polar Camera is established from the
;       date specified.  The direction to the Moon in inertial
;       coordinates is calculated using a routine coded by Michael
;       R. Greason, Hughes STX Corporation, October 1988.  The
;	rest of the calculation is based on equations in the
;	Observer's Handbook.
; EXAMPLE:
;       moonazel,94,10,12,0,0,0,az,el,/ref      ; calculate lunar
;                                               ; az and el for
;                                               ; October 12, 1994,
;                                               ; 00:00 UT (Eureka).
; SEE ALSO:
;	moonpos.pro
; MODIFICATION HISTORY:
; 	Written by:	D P Steele, ISR, U of Calgary, March 1994.
;       Documented:     September 1994.
;-
PRO moonazel,year,month,day,hour,minute,second,azimuth,elevation $
	    ,mra,mdec

; This routine computes the APPROXIMATE lunar azimuth and elevation
; for a given observing location.  If enough arguments are specified,
; the lunar right ascension and declination returned by moonpos.pro
; are also returned.

fhour=hour+(second/60.+minute)/60.

yr=year
l1900=WHERE(yr LT 1900.,nl1900)
IF nl1900 GT 0 THEN BEGIN
    l50=WHERE(yr(l1900) LT 50.,nl50)
    IF nl50 GT 0 THEN yr(l1900(l50))=yr(l1900(l50))+2000.
    IF nl50 LT nl1900 THEN BEGIN
	g50=setminus(INDGEN(N_ELEMENTS(l1900)),l50)
        yr(l1900(g50))=yr(l1900(g50))+1900.
	ENDIF
    ENDIF

PCsite=PoCasite(year,month,day)

ns=N_ELEMENTS(PCsite)
lat=REPLICATE(51.08D0,ns)	; default to Calgary
lon=REPLICATE(245.8681D0,ns)

rabb=WHERE(PCsite EQ 2,nrabb)	; Rabbit Lake coordinates
IF nrabb GT 0 THEN BEGIN
    lat(rabb)=58.23D0
    lon(rabb)=256.333D0
    ENDIF
eur=WHERE(PCsite EQ 3,neur)	; Eureka coordinates
IF neur GT 0 THEN BEGIN
    lat(eur)=80.0533D0
    lon(eur)=273.5839D0
    ENDIF

other=WHERE(ABS(PCsite-2) GT 1,nother)
IF nother GT 0 THEN BEGIN
    PRINT,'An unrecognized Polar Camera site is indicated!'
    PRINT,'Enter the N lat and E long (deg)'
    READ,olat,olon
    lat(other)=olat
    lon(other)=olon
    END

IF ns EQ 1 THEN BEGIN		; scalar inputs, scalar outputs
    lat=lat(0)
    lon=lon(0)
    ENDIF

jdf=ymdh2jdf(yr,month,day,fhour)
moonpos,jdf,mra,mdec,mx,my,mz	; mra, mdec in radians; mx, my, mz in Re

t2000=(DOUBLE(jdf)-2451545.0D0)/36525.0D0

locmeanst=lmst(ROUND(jdf),fhour/24.0D0,lon)*360.D0	; lmst in degrees

xp=mx-COS(lat*!DTOR)*COS(locmeanst*!DTOR)		; topocentric
yp=my-COS(lat*!DTOR)*SIN(locmeanst*!DTOR)		; EI coordinates
zp=mz-SIN(lat*!DTOR)					; wrt lat, lon

rp=SQRT(xp*xp+yp*yp+zp*zp)
alphap=ATAN(yp,xp)					; apparent RA
wneg=WHERE(alphap LT 0,nwneg)
IF nwneg GT 0 THEN alphap(wneg)=alphap(wneg)+2.0*!DPI
deltap=ASIN(zp/rp)					; apparent dec

mra=alphap
mdec=deltap

hourangle=(locmeanst-!RADEG*mra) MOD 360.0D0		; ditto for hourangle
whan=WHERE(hourangle LT 0.0D0,nwhan)
IF nwhan GT 0 THEN hourangle(whan)=hourangle(whan)+360.0D0

elevation=!RADEG*ASIN(SIN(mdec)*SIN(!DTOR*lat) $
	             +COS(!DTOR*hourangle)*COS(mdec)*COS(!DTOR*lat))

azimuth=!RADEG*ASIN(-COS(mdec)*SIN(!DTOR*hourangle)/COS(!DTOR*elevation))

ca=WHERE((((0.0D0 LE azimuth) AND (azimuth LT 90.0D0)) AND $
          ((270.0D0 LE hourangle) AND (hourangle LT 360.0D0))) OR $
         (((-90.0D0 LE azimuth) AND (azimuth LT 0.0D0)) AND $
          ((0.0D0 LE hourangle) AND (hourangle LT 90.0D0))),nca)
IF nca GT 0 THEN azimuth(ca)=180.0D0-azimuth(ca)

lz=WHERE(azimuth LT 0.0D0,nlz)
IF nlz GT 0 THEN azimuth(lz)=azimuth(lz)+360.0D0

RETURN
END

; --------------------------------------------------------------------

;+
; Name:
;	legend
; Purpose:
;	This procedure makes a legend for a plot.  The legend can contain
;	a mixture of symbols, linestyles, and filled polygons (usersym).
; Examples:
;	The call:
;		legend,['diamond','asterisk','square'],psym=[4,2,6]
;	produces:
;		-----------------
;		|		|
;		|  <> diamond	|
;		|  *  asterisk	|
;		|  [] square	|
;		|		|
;		-----------------
;	Each symbol is drawn with a plots command, so they look OK.
;	Other examples are given in usage.
; Usage:
;	legend,items,linestyle=linestyle	; vertical legend at upper left
;	legend,items,psym=psym			; ditto except using symbols
;	legend,items,psym=psym,/horizontal	; horizontal format
;	legend,items,psym=psym,box=0		; sans border
;	legend,items,psym=psym,delimiter='='	; embed an '=' betw psym & text
;	legend,items,psym=psym,margin=2		; 2-character margin
;	legend,items,psym=psym,position=pos	; position of legend
;	legend,items,psym=psym,number=2		; plot two symbols, not one
;	legend,items,/fill,psym=[8,8,8],colors=[10,20,30]; 3 filled squares
; Inputs:
;	items = text for the items in the legend, a string array
; Optional Inputs:
;	linestyle = array of linestyle numbers
;	psym = array of plot symbol numbers.  If psym(i) is negative, then a
;		line connects pts for ith item.  If psym(i) = 8, then the
;		procedure usersym is called with vertices define in the
;		keyword usersym.
;	N. B.: Choose either linestyle, psym or both.  If both linestyle and
;		psym parameters are present, normal plot behaviour occurs.
;		By default, if psym is positive, you get one point so there is
;		no connecting line.
; Optional Keywords:
;	/help = flag to print header
;	/horizontal = flag to make the legend horizontal
;	/vertical = flag to make the legend vertical (D=vertical)
;	box = flag to include/omit box around the legend (D=include)
;	delimiter = embedded character(s) between symbol and text (D=none)
;	colors = array of colors for plot symbols/lines (D=!color)
;	textcolors = array of colors for text (D=!color)
;	margin = margin around text measured in characters and lines
;	spacing = line spacing (D=bit more than character height)
;	pspacing = psym spacing (D=2 characters)
;	charsize = just like !p.charsize for plot labels
;	position = normalized coordinates of the upper left of the legend
;	number = number of plot symbols to plot or length of line (D=1)
;	usersym = 2-D array of vertices, cf. usersym in IDL manual. (D=square)
;	/fill = flag to fill the usersym
; Outputs:
;	legend to current plot device
; Common blocks:
;	none
; Procedure:
;	If keyword help is set, call doc_library to print header.
; Restrictions:
;	Here are some things that aren't implemented.
;	It would be nice to allow data and device coords as well.
;	An orientation keyword would allow lines at angles in the legend.
;	An array of usersyms would be nice---simple change.
;	An order option to interchange symbols and text might be nice.
;	Somebody might like double boxes, e.g., with box = 2.
; Side Effects:
; Modification history:
;	write, 24-25 Aug 92, F K Knight (knight@ll.mit.edu)
;-
pro legend,help=help,items,linestyle=linestyle,psym=psym $
  ,horizontal=horizontal,vertical=vertical,box=box,margin=margin $
  ,delimiter=delimiter,spacing=spacing,charsize=charsize,pspacing=pspacing $
  ,position=position,number=number,colors=colors,textcolors=textcolors $
  ,fill=fill,usersym=usersym
;
;	=====>> HELP
;
if keyword_set(help) then begin & doc_library,'legend' & return & endif
;
;	=====>> SET DEFAULTS
;
n = n_elements(items)
szt = size(items)
if (szt(szt(0)+1) ne 7) then message,'First parameter must be a string array.  For help, type legend,/help'
doline = (n_elements(psym) ne n)
if (n_elements(linestyle) ne n) and doline then message,'You must have either linestyle or psym keyword with the same number of elements as items.  For help, type legend,/help'
if n_elements(linestyle) eq 0 then linestyle = intarr(n)	; D=SOLID
if n_elements(psym) eq 0 then psym = intarr(n)			; D=SOLID

if n_elements(horizontal) eq 0 then begin 			; D=VERTICAL
  if n_elements(vertical) eq 0 then vertical = 1
endif else begin
  if n_elements(vertical) eq 0 then vertical = not horizontal
endelse

if n_elements(box) eq 0 then box = 1
if n_elements(margin) eq 0 then margin = 1
if n_elements(delimiter) eq 0 then delimiter = ''
if n_elements(charsize) eq 0 then charsize = 1
if n_elements(spacing) eq 0 then spacing = 1.2
if n_elements(pspacing) eq 0 then pspacing = 3
xspacing = !d.x_ch_size/float(!d.x_size) * (spacing > charsize)
yspacing = !d.y_ch_size/float(!d.y_size) * (spacing > charsize)
if !x.window(0) eq !x.window(1) then begin
  plot,/nodata,xstyle=4,ystyle=4,[0],/noerase
endif
;	next line takes care of weirdness with small windows
pos = [min(!x.window),min(!y.window),max(!x.window),max(!y.window)]
if n_elements(position) eq 0 then position = [pos(0),pos(3)] + [xspacing,-yspacing]
if n_elements(number) eq 0 then number = 1
if n_elements(colors) eq 0 then colors = !color + intarr(n)
if n_elements(textcolors) eq 0 then textcolors = !color + intarr(n)
fill = keyword_set(fill)
if n_elements(usersym) eq 0 then usersym = [[0,0],[0,1],[1,1],[1,0]]-0.5
;
;	=====>> INITIALIZE POSITIONS
;
yoff = 0.2*yspacing
maxwidth = 0					; SAVED WIDTH FOR DRAWING BOX
x0 = position(0) + (margin)*xspacing		; INITIAL X & Y POSITIONS
y0 = position(1) - (margin-0.5)*yspacing
y = y0						; STARTING X & Y POSITIONS
x = x0
if vertical then begin				; CALC OFFSET FOR TEXT START
  xt = 0					; DEFAULT X VALUE
  for i = 0,n-1 do begin
    if psym(i) eq 0 then num = (number + 1) > 3 else num = number
    if psym(i) lt 0 then num = number > 2	; TO SHOW CONNECTING LINE
    if psym(i) eq 0 then expand = 1 else expand = 2
    xt = (expand*pspacing*(num-1)*xspacing) > xt
    endfor
endif	; NOW xt IS AN X OFFSET TO ALIGN ALL TEXT ENTRIES
;
;	=====>> OUTPUT TEXT FOR LEGEND, ITEM BY ITEM
;	=====>> FOR EACH ITEM, PLACE SYM/LINE, THEN DELIMITER,
;	=====>> THEN TEXT---UPDATING X & Y POSITIONS EACH TIME.
;
for i = 0,n-1 do begin
  if vertical then x = x0 else y = y0		; RESET EITHER X OR Y
  x = x + xspacing
  y = y - yspacing
  if psym(i) eq 0 then num = (number + 1) > 3 else num = number
  if psym(i) lt 0 then num = number > 2		; TO SHOW CONNECTING LINE
  if psym(i) eq 0 then expand = 1 else expand = 2
  xp = x + expand*pspacing*indgen(num)*xspacing
  yp = y + yoff + intarr(num)
  if psym(i) eq 0 then begin
    xp = [min(xp),max(xp)]			; TO EXPOSE LINESTYLES
    yp = [min(yp),max(yp)]			; DITTO
    endif
  if psym(i) eq 8 then usersym,usersym*charsize,fill=fill,color=colors(i)
  plots,xp,yp,color=colors(i),/normal $
    ,linestyle=linestyle(i),psym=psym(i),symsize=charsize
  if vertical then x = x + xt else x = max(xp)
  x = x + xspacing
  xyouts,x,y,delimiter,width=width,/norm,color=textcolors(i),size=charsize
  x = x + width	+ xspacing
  xyouts,x,y,items(i),width=width,/norm,color=textcolors(i),size=charsize
  x = x + width
  if not vertical and (i lt (n-1)) then x = x+2*xspacing; ADD INTER-ITEM SPACE
  maxwidth = (x + xspacing*margin) > maxwidth	; UPDATE MAXIMUM X
  endfor
;
;	=====>> OUTPUT BORDER
;
if box then begin
  x = position(0)
  y = position(1)
  if vertical then bottom = n else bottom = 1
  ywidth = - (2*margin+bottom-0.5)*yspacing
  plots,[x,maxwidth,maxwidth,x,x],y + [0,0,ywidth,ywidth,0],/norm
  endif
return
end

; --------------------------------------------------------------------

;+
; NAME:
;     PSOPEN
; PURPOSE:
;     Opens a file for PostScript graphics output
; CATEGORY:
;     Utility
; CALLING SEQUENCE:
;     PSOPEN,pfile,font=font,portrait=portrait,landscape=landscape $
;                 ,color=color,longpage=longpage
; INPUTS:
;     pfile = string specifying name of PostScript output file to create
; KEYWORD PARAMETERS:
;     font:      specifies PostScript font other than default (e.g., 'times')
;     portrait:  specifies portrait orientation on page (default)
;     landscape: specifies landscape orientation on page
;     color:     specifies color PostScript output, 8 bits per pixel
;     longpage:  specifies full length of page available for plotting
;
;     Available fonts include AVANTGARDE, BKMAN, COURIER, HELVETICA,
;        PALATINO, SCHOOLBOOK, TIMES, ZAPFCHANCERY, ZAPFDINGBATS
;
; OUTPUTS:
;     None
; OPTIONAL OUTPUT PARAMETERS:
;     None
; COMMON BLOCKS:
;     PS_STUFF,Pfont,Dname,Ps_filename
;        used by CLOSE_PS to reset font, graphics device type, and to enable
;        it to print a ^D to the file for compatibility with QMS printer
; SIDE EFFECTS:
;     Opens a file for plotting; sets plot type to PostScript;
;     sets !P.FONT to 0 (if a font change requested)
; RESTRICTIONS:
;     None known.
; PROCEDURE:
;
; MODIFICATION HISTORY:
;     T.J.Hughes, HIA/NRCC , November 30, 1990
;     D. P. Steele, CNSR/UofC, April 21 1993, modified to accept command
;	line parameters; added color output option
;-

PRO PSOPEN,pfile,font=font,portrait=portrait,landscape=landscape $
                ,color=color,longpage=longpage

COMMON ps_stuff,pfont,dname,ps_filename
nfont=''
IF KEYWORD_SET(font) THEN nfont=',/'+font
orient=',/PORTRAIT'
IF KEYWORD_SET(landscape) THEN orient=',/LANDSCAPE'
ystuff=''
IF KEYWORD_SET(longpage) THEN ystuff=',YSIZE=22.6,YOFFSET=3.46'
ps_filename=pfile
dvfile="'"+pfile+"'"
dname=!D.NAME
pfont=!P.FONT
SET_PLOT,'ps'
dev_string='DEVICE,FILE='+dvfile+orient+nfont+ystuff
IF KEYWORD_SET(color) THEN dev_string=dev_string+',/color,bits_per_pixel=8'
istat=EXECUTE(dev_string)
RETURN
END

; --------------------------------------------------------------------

PRO PSCLOSE, NORMPS3=NORMPS3
;+
; NAME:
;     PSCLOSE
; PURPOSE:
;     Closes a PostScript output file (opened by OPEN_PS). It then reopens
;     the file for Update, and appends a <CTRL>D to the file for compatability
;     with the QMS PostSCript printer
; CATEGORY:
;     Utility
; CALLING SEQUENCE:
;     PSCLOSE
; INPUTS:
;     None
; OPTIONAL INPUT PARAMETERS:
;     None
; KEYWORD PARAMETERS:
;       /NORMPS3: If set, prevents the output PostScript file being
;               reopened to remove the spurious PS-Adobe-3.0
;               characters in the first line that mess up the
;		PostScript printers.
; OUTPUTS:
;     None
; OPTIONAL OUTPUT PARAMETERS:
;     None
; COMMON BLOCKS:
;     PS_STUFF,pfont,dname,ps_filename
; SIDE EFFECTS:
;     Closes the file PS_FILENAME
;     Appends a <CTRL>D to that file
; RESTRICTIONS:
;     None known
; PROCEDURE:
;
; MODIFICATION HISTORY:
;     T.J.Hughes, HIA/NRCC , November 30, 1990
;     D.P.Steele, CNSR/UofC, April 21, 1993, modified from CLOSE_PS.PRO
;-
COMMON ps_stuff,pfont,dname,ps_filename

;      +-------------------------------------------------------------------+
;      |Close PS_FILENAME so that all the graphics buffers are cleared.    |
;      |EMPTY didn't seem to do it, and we can't use FLUSH,UNIT because    |
;      |we don't know what unit has been used.                             |
;      +-------------------------------------------------------------------+
DEVICE,/CLOSE
ctrl_d=string(4b)

;      +-------------------------------------------------------------------+
;      |Open the same file again, for updating, and append a <CTRL>D to    |
;      |it                                                                 |
;      +-------------------------------------------------------------------+
GET_LUN,unit
OPENU,unit,ps_filename,/APPEND
PRINTF,unit,ctrl_d
CLOSE,unit
FREE_LUN,unit

;      +-------------------------------------------------------------------+
;      |Reset the output graphics type to what it was before OPEN_PS was   |
;      |called, and reset !P.FONT similarly                                |
;      +-------------------------------------------------------------------+
SET_PLOT,dname
!P.FONT=pfont

;      +-------------------------------------------------------------------+
;      |Edit the first line of the file to remove the PS-Adobe-3.0 string  |
;      |that messes up (at least) the b/w PostScript printer.              |
;      +-------------------------------------------------------------------+
IF KEYWORD_SET(normps3) THEN GOTO,done

pstext=qgetfile(ps_filename)
pstext(0)=STRMID(pstext(0),0,2)
mputfile,ps_filename,pstext

done:

RETURN
END

; --------------------------------------------------------------------
;+
; NAME:
;	QLOOK
; PURPOSE:
;	Display of zenith brightnesses from QLOOK.## quick-look files
;	created by the Polar Camera.
; CATEGORY:
;	Data display
; CALLING SEQUENCE:
;	QLOOK, FILE [, /DEBUG ] [, /PS ]
; INPUTS:
;	FILE:	The full pathname of the file to be processed.
; OPTIONAL INPUTS:
;	None.
; KEYWORD PARAMETERS:
;	/DEBUG:	Displays debugging output from the procedure.
;	/PS:	Sends the output plot to a PostScript file, instead of the
;		default plot window.
; OUTPUTS:
;	None.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	Reads the specified file; otherwise, this baby is self-contained.
; RESTRICTIONS:
;	The specified file must exist (duhhh!).
; PROCEDURE:
;	Later.
; EXAMPLE:
;	Later.
; SEE ALSO:
;	Later.
; MODIFICATION HISTORY:
; 	Written by:	D P Steele, ISR, January 1995.
;-

PRO qlook,year,month,day,binfac,debug=debug,ps=ps

; First convert HH:MM:SS times to UT seconds.
; Then separate data by camera.
; Then determine cycles and identify dark counts.
; Save dark counts and remove them from the other data.
; Subtract interpolated dark counts from all data.
; Plot data by camera/filter.

;	Print usage message

IF N_PARAMS() LT 4 THEN BEGIN
    MESSAGE,"Usage: qlook,year,month,day,binfac,debug=debug,ps=ps"
    RETURN
    ENDIF

;	Set up OS-dependent parameters explicitly

dos=(STRLOWCASE(STRMID(!VERSION.OS,0,3)) EQ 'win')

IF dos THEN BEGIN
    dd='\'
    iroot='c:\images'
    gzipcmd='c:\datacomp\gzip386 -f '
    psroot='c:\ps'
    ENDIF ELSE BEGIN
    dd='/'
    iroot='/jasper/cnsr3_data1'
    gzipcmd='/usr/local/bin/gzip -f '
    psroot=root+dd+'ps'
    ENDELSE
logroot=iroot+dd+'log'

;	Load the quick-look data file

dayst=STRING(day,FORMAT='(I2.2)')
file='qlook.'+dayst
IF month LT 10 THEN hexmo=STRING(month,FORMAT='(I1)') $
               ELSE hexmo=STRING(BYTE(month-10)+BYTE('a'))
hexdate=hexmo+dayst
path=logroot+dd+STRING(year MOD 100,FORMAT='(I2.2)') $
            +dd+hexdate+dd+file

f=rfindfile(path+'*',count=nf)
IF nf EQ 0 THEN BEGIN
    MESSAGE,path+' not found'
    READ,'Enter alternate path to qlook.## file:',path
    path=path+dd+file
    f=rfindfile(path,count=nf)
    IF nf EQ 1 THEN f=f(0) ELSE BEGIN
        MESSAGE,file+' not found; aborting' $
               ,/INFORMATIONAL
        RETURN
        ENDELSE
    ENDIF ELSE BEGIN
    f=f(0)
    IF STRPOS(f,'z') NE -1 THEN BEGIN	; gzipped!
        zipped=1B
        SPAWN,gzipcmd+'-d '+f
        f=rfindfile(path,count=nf)
        IF nf EQ 1 THEN f=f(0) ELSE BEGIN
            MESSAGE,'Can''t find gunzipped file!' $
                   ,/INFORMATIONAL
            RETURN
            ENDELSE
        ENDIF ELSE zipped=0B
    ENDELSE

q=qgetfile(f)
nq=N_ELEMENTS(q)
IF nq LT 2 THEN STOP
IF zipped THEN SPAWN,gzipcmd+f
IF KEYWORD_SET(debug) THEN help,q
; READ,'Enter the CCD pixel binning factor for image pixels (1/2/4): ',binfac

;	Extract numeric data and convert HH:MM:SS times to UT seconds

cam=INTARR(nq)
filt=cam
expos=cam
data=LONARR(nq)
utsec=data
FOR i=0,nq-1 DO BEGIN
    qi=q(i)
    cam(i)=FIX(STRMID(qi,0,1))
    filt(i)=FIX(STRMID(qi,2,1))
    expos(i)=FIX(STRMID(qi,5,3))
    hut=LONG(STRMID(qi,9,2))
    mut=LONG(STRMID(qi,12,2))
    sut=LONG(STRMID(qi,15,2))
    utsec(i)=(hut*60+mut)*60+sut
    data(i)=LONG(STRMID(qi,18,STRLEN(qi)-18))
    ENDFOR

;	Separate data by camera.

c0=WHERE((cam EQ 0) AND (expos NE 0),nc0)
c1=WHERE((cam EQ 1) AND (expos NE 0),nc1)

;	Merge data from two cameras and fill gaps

nc2=MAX([nc0,nc1])+10	; allow a good margin of error

filt2=BYTARR(2,nc2)	; arrays to hold merged data
expos2=INTARR(2,nc2)
utsec2=LONARR(2,nc2)
data2=utsec2

p0=0			; pointers to current data
p1=p0
p2=p0

REPEAT BEGIN		; process all data entries
    pc0=c0(p0<(nc0-1))
    pc1=c1(p1<(nc1-1))
    CASE 1 OF

;	simultaneous entries from both cameras

utsec(pc0) EQ utsec(pc1): BEGIN
            filt2(0,p2)=filt(pc0)
            filt2(1,p2)=filt(pc1)
            expos2(0,p2)=expos(pc0)
            expos2(1,p2)=expos(pc1)
            utsec2(0,p2)=utsec(pc0)
            utsec2(1,p2)=utsec(pc1)
            data2(0,p2)=data(pc0)
            data2(1,p2)=data(pc1)
            p0=(p0+1)<nc0
            p1=(p1+1)<nc1
            p2=p2+1
            END

;	good cam 0 entry and missing cam 1 entry

utsec(pc0) LT utsec(pc1): BEGIN
            filt2(0,p2)=filt(pc0)
            filt2(1,p2)=9		; 9 => bad filter!
            expos2(0,p2)=expos(pc0)
            expos2(1,p2)=expos(pc0)	; assume same exposure
            utsec2(0,p2)=utsec(pc0)
            utsec2(1,p2)=utsec(pc0)	; assume same time
            data2(0,p2)=data(pc0)
            data2(1,p2)=0		; 0 => no data
            p0=(p0+1)<nc0
            p2=p2+1
            END

;	good cam 1 entry and missing cam 0 entry

utsec(pc0) GT utsec(pc1): BEGIN
            filt2(0,p2)=9		; 9 => bad filter!
            filt2(1,p2)=filt(pc1)
            expos2(0,p2)=expos(pc1)	; assume same exposure
            expos2(1,p2)=expos(pc1)
            utsec2(0,p2)=utsec(pc1)	; assume same time
            utsec2(1,p2)=utsec(pc1)
            data2(0,p2)=0		; 0 => no data
            data2(1,p2)=data(pc1)
            p1=(p1+1)<nc1
            p2=p2+1
            END

        ENDCASE

    ENDREP UNTIL (p0 EQ nc0) AND (p1 EQ nc1)

npairs=p2
utsec=REFORM(utsec2(0,0:npairs-1))		; one utsec array will do

;	Determine cycles and identify dark counts

dutsec=(utsec-SHIFT(utsec,1))
newday=WHERE(dutsec LT -23*3600.,nnewday)	; find day wraparound

IF (nnewday EQ 1) AND (newday(0) NE 0) THEN BEGIN
    utsec(newday(0):*)=utsec(newday(0):*)+24*3600.
    dutsec(newday(0))=dutsec(newday(0))+24*3600.
    ENDIF

dutsec=dutsec>0
cycstart=[0,WHERE(dutsec GT 68),npairs]
ncycle=N_ELEMENTS(cycstart)

nincycle=(SHIFT(cycstart,-1)-cycstart)>0
nincycle=nincycle(0:ncycle-2)
cychist=HISTOGRAM(nincycle,MIN=0,BIN=1)
nincnorm=MAX(cychist,normlen)
longer=WHERE(nincycle GT normlen)

wdark=cycstart(longer+1)-1		; indices of dark counts
dkut=utsec(wdark)			; times of dark counts

;	Save dark counts and remove them from the other data

dark=LONARR(2,N_ELEMENTS(wdark))
dark(0,0)=data2(0,wdark)		; values of dark counts
dark(1,0)=data2(1,wdark)
; PRINT,[TRANSPOSE(dkut/3600.),dark]

;	Retain only open-shutter data and parameters

allcts=REPLICATE(1,npairs)
allcts(wdark)=0				; 1 if shutter open
whimg=WHERE(allcts,nimg)
utsec=utsec(whimg)
expos2=expos2(*,whimg)
filt2=filt2(*,whimg)
data2=data2(*,whimg)			; delete dark counts

;	Smooth dark counts to remove hot-pixel noise

FOR i=0,1 DO BEGIN
    tmp=dark(i,*)
    tmp=REFORM(MEDIAN(REFORM(tmp),3),1,N_ELEMENTS(tmp))
    dark(i,*)=tmp
    ENDFOR

;	Subtract interpolated dark counts from all data

idark=LONARR(2,nimg)
idark(0,0)=REFORM(LONG(interpol(dark(0,*),dkut,utsec)),1,nimg)
idark(1,0)=REFORM(LONG(interpol(dark(1,*),dkut,utsec)),1,nimg)
data2(0,*)=(data2(0,*)-idark(0,*))>0
data2(1,*)=(data2(1,*)-idark(1,*))>0
IF KEYWORD_SET(debug) THEN hlp,data2,idark

;	Plot data by camera/filter

;	Calibration values (rayleigh-seconds per DN for 2x2 binning)

cal=FLTARR(2,5)
cal(0,*)=[129.4, 62.3, 72.4, 62.7,  0.0]*(2./binfac)^2
cal(1,*)=[113.3,103.8, 88.0, 69.2, 49.1]*(2./binfac)^2

;	Wavelengths and features

fname=STRARR(2,5)
fname(0,*)=['!17829 nm   OH','!17608 nm   bg','!17630.0 nm OI' $
           ,'!17668 nm   N!i2!n','!17520.0 nm NI']
fname(1,*)=['!17835 nm   OH','!17820 nm   bg','!17785 nm   N!i2!e+!n' $
           ,'!17557.7 nm OI','!17714 nm   bg']

;	Prepare plot device

IF KEYWORD_SET(ps) THEN BEGIN
    yymm=STRING(year MOD 100,month,FORMAT='(2I2.2)')
    psfile=yymm+dayst+'ql.ps'
    psopen,psroot+dd+psfile,/long,font='times'
    ENDIF ELSE WINDOW,0,XSIZE=530,YSIZE=740

;	Set up plot space

!P.MULTI=[0,1,2,0,0]
date=ymd2date(year,month,day,format='d$ n$ y$')
!X.TITLE='!17Universal Time (hours, '+date+')'
!X.RANGE=[0,24]
!X.STYLE=1
!X.TICKS=4
!X.TICKV=INDGEN(5)*6
!X.MINOR=6
!Y.TITLE='!17Zenith Intensity (R or R/nm)'
!Y.RANGE=[0,1000]
!Y.STYLE=1

IF KEYWORD_SET(debug) THEN BEGIN
    STOP
    resp=''
    READ,'Enter an alternate full-scale brightness? ',resp
    IF resp NE '' THEN !Y.RANGE(1)=FLOAT(resp)
    ENDIF

;	Plot data for each camera, starting with 0

FOR ic=0,1 DO BEGIN
    ptitle='!17Quick Look Plot for Camera '+STRTRIM(ic,2)+', '+date
    PLOT,[0,1],/NODATA,TITLE=ptitle
    found=REPLICATE(0,5)
    FOR jf=0,4 DO BEGIN
        wjf=WHERE(filt2(ic,*) EQ jf,nwjf)
        IF nwjf GT 0 THEN BEGIN
	    found(jf)=nwjf
	    tf=utsec(wjf)
	    df=data2(ic,wjf)*cal(ic,jf)/expos2(ic,wjf)
            CASE ic OF
;	        R => R/nm for 608 nm bg channel
             0: IF jf EQ 1 THEN df=df/4.4
             1: BEGIN
;	            R => R/nm for 820 nm bg channel
	            IF jf EQ 1 THEN df=df/8.2
;	            R => R/nm for 714 nm bg channel
	            IF jf EQ 4 THEN df=df/6.9
                    END
                ENDCASE
	    IF KEYWORD_SET(debug) THEN STOP
	    OPLOT,tf/3600.,df,LINESTYLE=jf
	    ENDIF
        ENDFOR

;	If any data were plotted, create a legend.

    IF TOTAL(found) GT 0 THEN BEGIN
        wfound=WHERE(found)
        labels=REFORM(fname(ic,wfound))
        lines=wfound
        syms=REPLICATE(0,N_ELEMENTS(wfound))

;	If Moon up, show its elevation

    IF ic EQ 0 THEN BEGIN
        y=REPLICATE(year,25)
        mo=REPLICATE(month,25)
        d=REPLICATE(day,25)
        h=FINDGEN(25)
        mi=FLTARR(25)
        s=mi
        moonazel,y,mo,d,h,mi,s,maz,mel
        moonup=(MAX(mel) GT 0.)
        IF moonup THEN BEGIN
	    OPLOT,h,mel*10.,PSYM=-6
	    labels=[labels,'!17Moon el.*10']
	    lines=[lines,0]
	    syms=[syms,-6]
	    ENDIF
        ENDIF

;	Determine the normal coordinates of the upper left plot corner
;	and plot the legend

        xdvpos=!P.CLIP(0)
        ydvpos=!P.CLIP(3)
        npos=CONVERT_COORD([xdvpos,ydvpos],/DEVICE,/TO_NORMAL)
        npos=npos(0:1)
        legend,labels,linestyle=lines,psym=syms,position=npos
        ENDIF
    timestmp
    ENDFOR

;	Close the PS device, if necessary.

IF KEYWORD_SET(ps) THEN psclose

;	We're done

RETURN
END