;+
; NAME:
;       iris_get_back()
;
; PURPOSE:
;       A function to estimate the FUV spectrograph stray light
;       background image.  A normalized background image is loaded and
;       scaled based on the telescope pointing, exposure time, and
;       binning.  The pointing-dependent background count rate is
;       estimated using a smoothed limb darkening function for a
;       wavelength of 500 nm.  The parameters for this function were
;       determined empirically for each filter (filter back-reflection
;       of visible-IR light onto slit prism seems to be the culprit). 
;
;       Some remaining questions:
;          Do level 1 data coordinates need to be updated?
;          Which index variable contains the proper roll angle?  
;          Do I have the sign right on roll angle?
;
; CALLING SEQUENCE:
;       image = iris_get_back(index, bgfile=bgfile)
;
; INPUTS:
;       index - An IRIS FUV spectrograph image index structure for
;               level 1 or 1.5 data.
;      bgfile - a filename pointing to the fits file to load
;
; OPTIONAL INPUTS:
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;       Returns an array containing the estimate of the stray light image. 
;
; OPTIONAL OUTPUTS:
;
; COMMON BLOCKS:
;       None
;
; PROCEDURES USED:
;       iris_back_func - function for calculating the smoothed limb
;                        darkening function
;
; COMMENTS:
;       Function must be given one index at a time, images should be
;       dark subtracted, but not yet flat fielded.
;
; EXAMPLES:
;          IDL> image = iris_get_back(index)
;
; MODIFICATION HISTORY:
;       Started 2013-Nov-14 by Sarah A. Jaeggli, Montana State University
;-


FUNCTION iris_get_back, index, bgfile=bgfile

  ;check if image is FUV
  if index.img_path ne 'FUV' then begin
     print, 'This is not FUV data, exiting...'
     return, 0
  endif

  ;relevant properties of the data which need to be read from the index are:
  ;pointing, roll angle
  xcen0 = index.xcen
  ycen0 = index.ycen
  roll  = index.sat_rot ;roll angle from solar north
  rsun  = index.rsun_obs

  ;rotate coordinates so they are w.r.t roll
  xcen = xcen0*cos(roll/180.*!pi) + ycen0*sin(roll/180.*!pi)
  ycen = xcen0*sin(roll/180.*!pi) + ycen0*cos(roll/180.*!pi)

  ;filter
  fw = index.ifwpos

  ;exposure time
  exp=index.exptime

  ;binning
  ybin=float(index.sumspat)
  xbin=float(index.sumsptrl)

  ;load normalized background image
  bg_img=readfits(bgfile, hdr)

  nx = (size(bg_img))[1]
  ny = (size(bg_img))[2]

  case fw of
     ;ld_pars =  x0[arcsec], y0[arcsec], gaussian width[arcsec], 
     ;             amplitude[DN/sec], offset[DN/sec]

     ;1330
     31  : ld_pars = [fxpar(hdr, 'LD31_0'), fxpar(hdr, 'LD31_1'), $
                      fxpar(hdr, 'LD31_2'), fxpar(hdr, 'LD31_3'), $
                      fxpar(hdr, 'LD31_4')]

     ;2796
     61  : ld_pars = [fxpar(hdr, 'LD61_0'), fxpar(hdr, 'LD61_1'), $
                      fxpar(hdr, 'LD61_2'), fxpar(hdr, 'LD61_3'), $
                      fxpar(hdr, 'LD61_4')]

     ;1400
     91  : ld_pars = [fxpar(hdr, 'LD91_0'), fxpar(hdr, 'LD91_1'), $
                      fxpar(hdr, 'LD91_2'), fxpar(hdr, 'LD91_3'), $
                      fxpar(hdr, 'LD91_4')]

     ;2832
     121 : ld_pars = [fxpar(hdr, 'LD121_0'), fxpar(hdr, 'LD121_1'), $
                      fxpar(hdr, 'LD121_2'), fxpar(hdr, 'LD121_3'), $
                      fxpar(hdr, 'LD121_4')]

     ;5000
     1   : begin
        print, 'No data available for the 5000A filter, exiting...'
        return, 0
     end

     ;1600
     151 : begin
        print, 'No data available for the 1600A filter, exiting...'
        return, 0
     end
  endcase

  ;calculate smoothed limb darkening function
  r_obs=sqrt((xcen-ld_pars[0])^2 + (ycen-ld_pars[1])^2)
  bg_level = iris_bgcal_func(r_obs, [ld_pars[2:*], rsun])

  ;scale intensity
  back = bg_img*bg_level*exp

  ;account for binning
  back = rebin(back, nx/xbin, ny/ybin, /SAMPLE)*xbin*ybin

return, back

end