;+
; NAME
;   ml_evalplatemags
;
; PURPOSE:
;   Given a plateMags structure, evaluate it at a set of locations.
;
;   The input structure ismagnitudes with associated x,y
;   locations in units of arcsec, where +x=+RA, +y=+DEC
;
;   The input magnitudes are fiber2mags, and so are the outputs.
;
; CALLING SEQUENCE:
;   mags=ml_evalplatemags(platemags,xloc,yloc,band,[rlim=,rscale=])
;
; INPUTS:
;   platemags: plateMags structure for a particular IFU
;   xloc: Vector of x positions to get the mag at
;   yloc: Vector of y positions to get the mag at
;     (xloc and yloc must have same size)
;   band: g,r,i,z
;
; OPTIONAL INPUTS:
;   rlim: Cutoff radius outside which to ignore values in platemags grid
;     (default is 2.5 arcsec, the ideal fiber spacing)
;   rscale: Radial scalelength for exponential weighting function.
;     (default is 1.6/2.35 arcsec, 1-sigma width of 1.6 arcsec seeing profile)
;   /visual: Display an info plot showing fibers and magnitudes to X window
;
; OUTPUTS:
;   mags: Vector of magnitudes with same length as xloc and yloc
;
; PROCEDURES CALLED:
;   valid_num()
;   ml_drawcirc()
;
; BUGS:
;  
;
; REVISION HISTORY:
;   13-Nov-2013  David R. Law (Dunlap Institute; drlaw@di.utoronto.ca)
;       First written.
;   12-Mar-2014  Updated for real plateMags files.  Added /visual
;       keyword for debugging and verification, added ability to handle
;       nan values in input files. (D. Law)
;-

function ml_evalplatemags,platemags,xloc,yloc,band,rlim=rlim,rscale=rscale,visual=visual

; Default rlim is 2.6 arcsec (i.e., slightly bigger than fiber sep so
; that we can handle missing fibers)
if (keyword_set(rlim)) then rlim=rlim $
else rlim=2.6

; Default rscale is 1.6/2.35
; (i.e., 1-sigma width of gaussian assuming 1.6 arcsec seeing)
if (keyword_set(rscale)) then rscale=rscale $
else rscale=1.6/2.35

; How many points to evaluate at?
nloc=n_elements(xloc)

; Create a magnitude vector of that length
mags=fltarr(nloc)

; Check that xloc and yloc have the same size
if (n_elements(yloc) ne nloc) then begin
  splog,'WARNING- xloc and yloc have different sizes!'
  return,0
endif

; Define the band to evaluate magnitudes in
; If bad band specified return zeroes
if ((band ne 'g')and(band ne 'r')and(band ne 'i')and(band ne 'z')) then begin
  splog,'WARNING- bad band specified for platemags'
  return, mags 
endif else begin
  if (band eq 'g') then jband=0
  if (band eq 'r') then jband=1
  if (band eq 'i') then jband=2
  if (band eq 'z') then jband=3
endelse

; Loop over nloc evaluation locations
for i=0,nloc-1 do begin
  ; Radius of platemags defined locations from xloc,yloc
  r=sqrt((xloc[i]-platemags.raoff)^2+(yloc[i]-platemags.decoff)^2)
  ; Consider only grid values from platemags within a given
  ; cutoff radius
  index=where(r le rlim,nindex)

  ; If no grid values within cutoff radius, then no
  ; valid information for the extrapolation, simply return
  ; magnitude of nearest non-nan point
  if (nindex eq 0) then begin
    validnum=valid_num(platemags.fiber2mag[jband])
    wherevalid=where(validnum eq 1,nvalid)
    ; If no valid numbers, platemags is 0
    if (nvalid eq 0) then mags[i]=0. $
    ; Otherwise give magnitude of nearest valid point
    else begin
      rtemp=r[wherevalid]
      mtemp=platemags[wherevalid].fiber2mag[jband]
      minloc=where(rtemp eq min(rtemp))
      mags[i]=mtemp[minloc[0]]
    endelse

  ; If there was one or more grid values within the cutoff radius
  ; then combine them using a modified Shephards method
  ; with exponential radial weights
  endif else begin
    w=exp(-0.5*r[index]^2/rscale^2)
    ; Figure out where the magnitudes are valid numbers
    validnum=valid_num(platemags[index].fiber2mag[jband])
    wherevalid=where(validnum eq 1,nvalid)
    ; If no valid numbers, platemags is 0
    if (nvalid eq 0) then mags[i]=0. $
    ; Otherwise computed weighted mean of valid numbers
    else begin
      mags[i]=total(w[wherevalid]*(platemags[index].fiber2mag[jband])[wherevalid])/total(w[wherevalid])
    endelse
  endelse
endfor

; Show a plot to xterm
if (keyword_set(visual)) then begin
  set_plot,'x'
  device,decomposed=0
  loadct,1

  maxmag=max(platemags.fiber2mag[jband])
  minmag=min(platemags.fiber2mag[jband])
  platemags_col=250-(platemags.fiber2mag[jband]-minmag)*250./(maxmag-minmag)
  plot,[min(platemags.raoff),max(platemags.raoff)],[min(platemags.decoff),max(platemags.decoff)],/nodata
  ml_drawcirc,platemags.raoff,platemags.decoff,diam=replicate(2.0,n_elements(platemags)),color=platemags_col
  xyouts,platemags.raoff,platemags.decoff,platemags.fiber2mag[jband]

  mags_col=250-(mags[*]-minmag)*250./(maxmag-minmag)
  ml_drawcirc,xloc,yloc,diam=replicate(2.0,nloc),color=mags_col
endif

return,mags
end