pro avgmed,out,images,filename=filename,log=log, $
		noweight=noweight,nocoadd=nocoadd,average=average, $
		evenmed=evenmed,maxrej=maxrej,minrej=minrej, $
		silent=silent

;+
; PURPOSE:
; Median average a stack of FITS images, with optional median weighting.
;
; NOTES:
; All frames are weighted by their medians, unless /noweight is set.
; If /log, then writes weights to a file as well and the options chosen.
;
; 5/20/94 M. Liu (UCB)
;
; All frames read from file are normalized by their # of coadds 
; unless /nocoadd is set
;
; 6/27/94 MCL
;
; Added ability to read images from IDL array (in a very inelegant way)
;   if want to work with IDL array, just call avgmed like:
;	IDL> avgmed, out, array
;   to read stuff from a file, just define the 'filename' parameter:
;	IDL> avgmed, out, filename='file.list'
; 7/01/94
;
; significant modifications:
;
;   1. 	modified median averaging so images are scaled to the average
;    	of all the medians, instead of the median of the first image
;   	-> this bombs if an image has a median = 0
;
;   2. 	when median filtering if have <= 6 images, will take the mean
;	of the middle 2 images when there is an even (2,4,6) number of images
;	(only if /evenmed set since this can be quite S-L-O-W); this probably
;	make little difference in practice
;
;   3.  added 'minrej' and 'maxrej': will toss these number of values
;	out before median averaging (no effect for arithmatic avg)
;	note: this is done AFTER any weighting.
;
;   4.  major editing of program to make it nicer
;
; 12/13/94
;
; added minrej/maxrej option for arithmatic averaging as well
; and modified to toss more than one pixel
;
; 6/2/95 MCL
;
; don't make a copy of the original images, instead doing scaling
; when do median filtering - runs about 10% faster for
; a stack of 30 images with simple median averaging
;
; 7/2/95 MCL
;
; corrected error in calculating the proper image weights
; (was incorrectly scaling the weights to 1.0 by multiplying by the
;   average value instead of dividing by it!)
; 12/17/96 MCL
;
; the above bug fix was wrong - no idea why I thought it was right!
; 01/29/97
;
; new algorithm for /evenmed - may be faster, but untested
; 05/04/97 MCL
;
; added if/then for /noweight so avoids multiplying by array of 1's
;   clumsier code but maybe faster operation
; used temporary() function when getting medians (does this help?)
; 02/16/01 MCL
;
; cleaned up comments & tabs
; minor bug fix to get code running (odd, since I don't remember
; hacking it lately, but whatever ...)
; 05/08/01 MCL
;
; **--> should rewrite /evenmed using the new /even flag in MEDIAN!!
;-


if (n_params() lt 2) then begin
    print,'avgmed,out,images,[filename=],[log=],'
    print,'       [noweight],[nocoadd],[average],'
    print,'       [evenmed],[maxrej=],[minrej=],[silent]'
    return
endif


; beginning messages
if keyword_set(silent) eq 0 then begin
    
    if keyword_set(noweight) then  $
      message, 'no weighting for averaging', /info $
    else message,'median weighting for averaging', /info
    
    if keyword_set(filename) then begin
        if keyword_set(nocoadd) then $
          message, 'not dividing each frame by number of coadds', /info $
        else message,'dividing each frame by number of coadds', /info
    endif
    
endif



;--------------------------------------------------
; (1) figure out how many images
;--------------------------------------------------
; (1a) either from an input file, not counting blank lines 
if keyword_set(filename) eq 1 then begin
    
    print,format='($,"* reading images from file ",A,": ")',filename 
    openr,unit0,filename,/get_lun
    n = 0
    file = ' ' 
    while eof(unit0) ne 1 do begin
        readf,unit0,file
        if (strlen(file) ne 0) then n = n+1
    endwhile
    print,strc(n),' files *'
    free_lun,unit0

; (1b) or from IDL array
endif else begin 
    sz = size(images)
    xs = sz(1)
    ys = sz(2)
    n = sz(3)
endelse


; set up output array and make sure things are ok
out = fltarr(xs,ys)
if keyword_set(maxrej) or keyword_set(minrej) then begin
    if (maxrej+minrej) gt n then message,'(minrej+maxrej) too large!'
    if minrej gt n then message,'minrej too large!' 
    if maxrej gt n then message,'maxrej too large!'
endif else begin
    if keyword_set(minrej) eq 0 then minrej = 0 
    if keyword_set(maxrej) eq 0 then maxrej = 0 
endelse
	


;--------------------------------------------------
; (2) read in images 
;--------------------------------------------------
; (2a) from a file
if keyword_set(filename) eq 1 then begin

   openr,unit0,filename,/get_lun
   for i =0,(n-1) do begin

	readf,unit0,file
        
        ; first image defines array sizes
        if (i eq 0) then begin
            out = readfits(file,head,/silent)
            sz = size(out)
            xs = sz(1)
            ys = sz(2)
            imgs = fltarr(xs,ys,n)
            imgs(*,*,0) = out
            print,'image size = ',strc(sz(1)),' x ',strc(sz(2))
        endif else begin
            ;imgs(*,*,i) = readfits(file,head,/silent)
            images(*,*,i) = readfits(file,head,/silent)
	endelse

        ; normalize by # of coadds if desired
        if (keyword_set(nocoadd) eq 0) then begin 
            extract_head,head,'M56COAD0',coadd
            imgs(*,*,i) = imgs(*,*,i) / float(coadd)
        endif
        
   endfor
   free_lun,unit0

; (2b) or images already in a 3d array, just get the medians
endif else begin

    ;;;imgs = images			; don't change the originals!

endelse



;-------------------------------------------------------
; (3) median weighting so medians for all images
;     are equal to the average of all the image medians
;-------------------------------------------------------
wtlist = fltarr(n)+1.0
if keyword_set(noweight) eq 0 then begin
    
    ; get medians
    medlist = fltarr(n,/nozero)		
    for k=0,n-1 do begin
        medlist(k) = median(images(*,*,k))
    endfor
    
    ; scale the images
    mm = total(medlist) / n
    for k = 0,n-1 do $
      wtlist(k) = mm / medlist(k)   
    ;wtlist(k) = 1./medlist(k) / mm  ; -> 12/17/96 edit is wrong!

    ; fancy two column printing of weights
    if keyword_set(silent) eq 0 then begin
        for k = 0,n/2-1 do begin
            print,'image ',strc(k),' weight = ',strc(wtlist(k)),$
              '      image ',strc(k+n/2),' weight = ', $
              strc(wtlist(k+n/2))
        endfor
        if (n-1)/2 eq n/2 then $
          print,'image ',strc(n-1),' weight = ',strc(wtlist(n-1))
    endif
endif



;--------------------------------------------------
; (4) median averaging
;--------------------------------------------------
if keyword_set(average) eq 0 then begin

  message, 'median averaging frames', /info

  ; (4a) evenmed: if number of images is even,
  ;               instead of using IDL's median (the higher value), 
  ;               take average of the two values straddling the center
  if (n/2 eq n/2.0) and keyword_set(evenmed) then begin
      print,'even # of frames: using modified median'
      if (n eq 2) then begin
          out = (images(*,*,0) + images(*,*,1))/2.0
      endif else begin
          for j = 0,(ys-1) do $
            for i = 0,(xs-1) do begin
              
              cut = images(i,j,*)*wtlist
              
              ; this method is more elegant, but may actually be slower,
              ; at least for only 4 images (maybe faster for more images)
              mcut = median(cut)
              ww = where(cut ne mcut, nw)
              if (nw gt 0) then $
                out(i, j) =  (mcut+median(cut(ww)))/2.0 $
              else $
                out(i, j) = mcut
              
              ;old method:         
              ;cut = cut(sort(cut))
              ;out(i,j) = (cut((n-1)/2) + cut(n/2)) / 2.0
              
              if i eq 40 and j eq 40 then $
                print,'compare ',out(i,j),' with ',median(images(i,j,*)*wtlist)
              
          endfor
      endelse


  ; (4b) minrej & maxrej: toss the highest and/or lowest pixels before
  ;                       median averaging (this only matters if maxrej != minrej)
  endif else if keyword_set(maxrej) and keyword_set(minrej) then begin
      aa = 0
      zz = n-1
      if keyword_set(maxrej) then begin
          message,'* maxrej: tossing highest pixel *', /info
          zz = n-2
      endif
      if keyword_set(minrej) then begin
          message,'* minrej: tossing lowest pixel *', /info
          aa = 1
      endif
      if keyword_set(minrej) and keyword_set(maxrej) then begin
          if minrej eq maxrej then  $
            message, '* minrej=maxrej w/median averaging does nothing! *', /info
      endif
      for j = 0,(ys-1) do begin
          for i = 0,(xs-1) do begin
              cut = images(i,j,*)*wtlist
              cut = cut(sort(cut))
              out(i,j) = median(cut(aa:zz))
          endfor
      endfor
      
  ; (4c) just median average each pixel
  endif else begin
      if keyword_set(noweight) then begin
          for j = 0,(ys-1) do $
            for i = 0,(xs-1) do begin
              out(i,j) = median(images(i,j,*))
          endfor
      endif else begin
          for j = 0,(ys-1) do $
            for i = 0,(xs-1) do begin
              out(i,j) = median(images(i,j,*)*wtlist)
          endfor
      endelse

  endelse       

; --------------------------------------------------
; (5) arithmatic averaging
; --------------------------------------------------
endif else begin
  message, 'arithmatic averaging frames', /info


  ; (5a) with minrej and/or maxrej before averaging
  if keyword_set(maxrej) or keyword_set(minrej) then begin
      aa = 0
      zz = n-1
      nn = n
      if keyword_set(maxrej) then begin
          print,'* maxrej: tossing highest ',strc(maxrej),' pixels *'
          zz = zz - maxrej
          nn = nn - maxrej
      endif
      if keyword_set(minrej) then begin
          print,'* minrej: tossing lowest ',strc(minrej),' pixels *'
          aa = minrej
          nn = nn - minrej
      endif
      for j = 0,(ys-1) do begin
          for i = 0,(xs-1) do begin
              cut = images(i,j,*)*wtlist
              cut = cut(sort(cut))
              out(i,j) = total(cut(aa:zz))/nn
          endfor
      endfor
      
  ; (5b) plain arithmatic averaging
  endif else begin
      for j = 0,(ys-1) do $
        for i = 0,(xs-1) do $
        out(i,j) = total(images(i,j,*)*wtlist)/n
  endelse
  
endelse


; --------------------------------------------------
; (6) write to log file
; --------------------------------------------------
if keyword_set(log) then begin
    
    print,'writing to log file: ',log
    openw,unitlog,log,/get_lun
    printf,'** log file from avgmed.pro **'
    
    if keyword_set(filename) eq 1 then begin
        openr,unit0,filename,/get_lun
        file = ' ' 
        while eof(unit0) ne 1 do begin
            readf,unit0,file
            if (strlen(file) ne 0) then n = n+1
        endwhile
        printf,unitlog,file,'	',strc(wtlist(i))  
        free_lun,unit0
    endif else for i=0,(n-1) do printf,unitlog,'image ',strc(i),$
      '	',strc(wtlist(i))  
    
    if keyword_set(nocoadd) then $
      printf,unitlog,'* not normalized to counts/coadd *' $
    else printf,unitlog,'* normalized to counts/coadd *'
    
    if keyword_set(noweight) then $
      printf,unitlog,'* no weighting for frames *' $
    else printf,unitlog,'* median weighting for frames *'
    
    if keyword_set(average) eq 0 then $
      printf,unitlog,'* median averaged frames *' $
    else printf,unitlog,'* arithmatic averaged frames *'
    
    close,unitlog
endif

end