;+
; NAME:
;       iris_back_func()
;
; PURPOSE:
;       A function to calculate a limb-darkening function which is
;       arbitrarily offset from solar center, smoothed, scaled, and
;       off-set in intensity.  For use in iris_get_back.
;
; CALLING SEQUENCE:
;       inten = iris_get_back(x, y, params)
;
; INPUTS:
;       x - solar X coordinate in arcsec 
;       y - solar Y coordinate in arcsec
;       params - an array of the function parameters containing:
;                [x0, y0, dw, amp, off, rsun]
;                x0   - X center of limb darkening function in arcsec
;                y0   - Y center of limb darkening function in arcsec
;                dw   - Gaussian width of smoothing funtion in arcsec
;                amp  - max intensity of limb darkening function from
;                       pedestal
;                off  - intensity pedestal for limb darkening function
;                rsun - solar radius for time of observation in arcsec
;
; OPTIONAL INPUTS:
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;       Returns intensity array of dimensions input for X and Y
;
; OPTIONAL OUTPUTS:
;       None
;
; COMMON BLOCKS:
;       None
;
; PROCEDURES USED:
;
; COMMENTS:
;
; EXAMPLES:
;          IDL> inten = iris_get_back(x, y, params)
;
; MODIFICATION HISTORY:
;       Started 2013-Nov-14 by Sarah A. Jaeggli, Montana State University
;-

function ld_func, r, a
  t = asin(r/1.) ;mu angle

  ;limb darkening function at 0.5 micron
  u=0.97 & v=-0.22

  I_ld = 1.0 - u - v + u*cos(t) + v*(cos(t))^2

  ;set off-limb values to 0
  good=where(finite(I_ld) eq 1, complement=bad)
  I_ld[bad]=0.

  return, I_ld
end

function iris_back_func, x, y, a

  x0=a[0]
  y0=a[1]
  dw=a[2]

  amp=a[3]
  off=a[4]

  rsun=a[5] ;solar radius is a fixed parameter

  ;calculate limb darkening function over a larger number of points so it
  ;can be accurately smoothed with a gaussian
  nr=1500
  n=nr*2+1

  rarray=findgen(n)-nr ;the radius in arcsec, 1 arcsec bins

  out=ld_func(rarray/rsun)

  psf=exp(-0.5*((rarray-nr)/dw)^2) + exp(-0.5*((rarray+nr)/dw)^2)

  out1=real_part(fft(fft(out)*fft(psf), /inverse))
  out1=out1/max(out1)

  ;use x,y coordinates as indices in function
  ridx=sqrt((x-x0)^2+(y-y0)^2)+nr
  z=out1[ridx]

  z=z*amp + off

  return, z
end