; Compute the fiber input losses for the benchmark spectra cases
; Fiber diameter 1.5 arcsec
; Seeing is 1.1 arcsec FWHM with Moffat beta=3.5, and scale with
; wavelength as (6355./wave)^0.3
; Lateral offset is 0.2 arcsec from fiber center
; Cases are:
;   * Point source
;   * ELG with half-light radius 0.35 arcsec
;   * LRG with half-light radius 1.0 arcsec
pro fiberloss

   seeing = 1.1 ; FWHM in arcsec at 6000 Ang
   pixscale = 0.025 ; arcsec/pix
   npix = 501 ; number of pixels across image
   beta = 3.5 ; Moffat profile parameter
   fibrad = 0.75 ; fiber radius in arcsec
   offset = 0.20 ; lateral error in arcsec
   r_half = 0.45 ; half-light radius for ELGs
   r_e = 1.0 ; half-light radius for LRGs

   xcen = (npix-1)/2
   rr = shift(dist(npix,npix), xcen, xcen) * pixscale
   x1 = (npix-1)/2 - fix(0.2*npix)
   x2 = (npix-1)/2 + fix(0.2*npix)
   fibmask = shift(rr LT fibrad, offset/pixscale)

   for wave=9200., 10000, 400. do begin
      fwhm = seeing * (6355./wave)^0.3
      alpha = 0.5 * fwhm / sqrt(2.^(1./beta) -1)
      img_seeing = (beta - 1) * (!pi * alpha^2)^(-1) * (1 + rr^2/alpha^2)^(-beta)

      ;;; print, 'Wave= ', wave, ' Point source= ', $
      ;;;  total(img_seeing*fibmask)/total(img_seeing)

      ; Convolve w/out /edge_zero because it's too slow
      img_elg = exp(-1.678*rr/r_half)
      img_elg2 = convol(img_elg,img_seeing[x1:x2,x1:x2],/center)
      print, 'Wave= ', wave, ' ELG= ', $
      total(img_elg2*fibmask)/total(img_elg2)

      img_lrg = exp(-7.67*((rr/r_e)^0.25 -1))
      img_lrg2 = convol(img_lrg,img_seeing[x1:x2,x1:x2],/center)
      print, 'Wave= ', wave, ' LRG= ', $
       total(img_lrg2*fibmask)/total(img_lrg2)
       
      ; /edge_zero does matter for more extended objects like LRGs,
      ; giving a ~3% absolute difference in acceptance.  501 pixel window
      ; may be too small too.
      ;;; img_lrg3 = convol(img_lrg,img_seeing[x1:x2,x1:x2],/center,/edge_zero)
      ;;; print, 'Wave= ', wave, ' LRG3= ', $
      ;;;  total(img_lrg3*fibmask)/total(img_lrg3)
   endfor
end