;+
; Project     :	STEREO
;
; Name        :	WCS_FIT_GRISM_FUNCT()
;
; Purpose     :	Function used by WCS_FIT_GRISM
;
; Category    :	FITS, Coordinates, WCS
;
; Explanation :	This function is used by the WCS_FIT_GRISM procedure to fit
;               grism parameters to spectral coordinates.  See WCS_FIT_GRISM
;               for more information.
;
; Syntax      :	Result = WCS_FIT_GRISM( PIXREL, PARAM )
;
; Examples    :	If PIXEL contains the array of pixel locations, and WAVE
;               contains the array of measured wavelengths, in Angstroms, and
;               CRPIX is the desired value of the FITS keyword, then
;
;                       PIXREL = PIXEL - (CRPIX - 1)
;                       WCS_FIT_GRISM, PIXREL, WAVE, PARAM, COORD='WAVE', $
;                               CUNIT='Angstrom', /INIT, MAX_ITER=1000
;
;               will do an initial fit of up to 1000 iterations, and
;
;                       OPLOT, PIXEL, WAVE-WCS_FIT_GRISM_FUNCT(PIXREL, PARAM)
;
;               will display the residuals.
;
; Inputs      :	PIXREL  = Array of pixel locations, relative to the reference
;                         pixel.  The FITS keyword CRPIXia must be consistent
;                         with the reference pixel used in calculating PIXREL.
;                         Note that FITS pixels start with 1, while IDL pixels
;                         start with 0--that distinction must be taken into
;                         account.
;
;               PARAM   = The parameter array, with the following elements:
;
;                               CRVAL           Reference value
;                               CDELT           Pixel spacing
;                               GME             G * m / cos(epsilon)
;                               NRA             nr * sin(alpha)
;                               NRAPRIME        nrprime * sin(alpha)
;                               THETA           theta
;
;                         See WCS_FIT_GRISM for more information.
;
; Opt. Inputs :	None.
;
; Outputs     :	The result of the function is the spectral coordinate, as a
;               function of pixel position.
;
; Opt. Outputs:	None.
;
; Keywords    :	None.
;
; Calls       :	Only native IDL statements.
;
; Common      :	The common block wcs_fit_grism is used to pass information from
;               the calling routine WCS_FIT_GRISM.
;
; Restrictions:	Not all combinations of the grism parameters are valid.
;
; Side effects:	None.
;
; Prev. Hist. :	None.
;
; History     :	Version 1, 28 June 2005, William Thompson, GSFC 
;
; Contact     :	WTHOMPSON
;-
;
function wcs_fit_grism_funct, pixrel, param
;
on_error, 2
common wcs_fit_grism, base_units, factor, coord_type, restfrq, restwav, gra
;
;  Define various constants, including those for calculating the refractive
;  index of air.
;
c = 2.99792458d8        ;Speed of light
h = 6.6260693d-34       ;Planck constant
apar = [287.6155d0,  1.62887d-12,  0.01360d-24]
pder = [287.6155d0, -1.62887d-12, -0.04080d-24]
;
;  Parse the parameter array
;
crval    = param[0]
cdelt    = param[1]
gme      = param[2]
nra      = param[3]
nraprime = param[4]
theta    = param[5] * !dpi / 180.
;
case coord_type of
    'FREQ': begin
        nu0 = factor * crval
        lambda0 = c / nu0
        deriv = -c / nu0^2
    endcase
    'ENER': begin
        nu0 = factor * crval / h
        lambda0 = c / nu0
        deriv = -c * h / nu0^2
    endcase
    'WAVN': begin
        kappa0 = factor * crval
        lambda0 = 1.d0 / kappa0
        deriv = -1.d0 / (c * kappa0)^2
    endcase
    'VRAD': begin
        v0 = factor * crval
        nu0 = restfrq * (1.d0 - v0 / c)
        lambda0 = c / nu0
        deriv = restfrq / nu0^2
    endcase
    'WAVE': begin
        lambda0 = factor * crval
        deriv = 1.d0
    endcase
    'VOPT': begin
        z0 = factor * crval
        lambda0 = restwav * (1.d0 + z0 / c)
        deriv = restwav / c
    endcase
    'ZOPT': begin
        z0 = factor * crval
        lambda0 = restwav * (1.d0 + z0)
        deriv = restwav
    endcase
    'AWAV': begin
        x = factor * crval
        if keyword_set(gra) then begin
            lambda0 = x
            deriv = 1.d0
        end else begin
            lambda0 = x*(1.d0 + 1.d-6*(apar[0] + apar[1]/x^2 + apar[2]/x^4))
            deriv = (1.d0 + 1.d-6 * (pder[0] + pder[1]/x^2 + pder[2]/x^4))
        endelse
    endcase
    'VELO': begin
        v0 = factor * crval
        lambda0 = restwav * (c + v0) / sqrt(c^2 - v0^2)
        deriv = c * restwav / ((c - v0) * sqrt(c^2 - v0^2))
    end
    'BETA': begin
        v0 = factor * crval * c
        lambda0 = restwav * (c + v0) / sqrt(c^2 - v0^2)
        deriv = c^2 * restwav / ((c - v0) * sqrt(c^2 - v0^2))
    endcase
endcase
;
;  For grisms-in-air, unless the coordinate type is AWAV, correct for the
;  distinction between vacuum and air wavelengths.
;
if keyword_set(gra) and (coord_type ne 'AWAV') then begin
    lambda0 = lambda0 / $
      (1.d0 + 1.d-6 * (apar[0] + apar[1]/lambda0^2 + apar[2]/lambda0^4))
    deriv = deriv * $
      (1.d0 + 1.d-6 * (pder[0] + pder[1]/lambda0^2 + pder[2]/lambda0^4))
endif
;
;  Calculate the reference angle gamma_r, and the derivative dGamma/dw.
;
gamma0 = gme * lambda0 - nra
if abs(gamma0) le 1 then gamma0 = asin(gamma0) else message, $
  'Incompatible grism parameters'
denom = cos(gamma0) * cos(theta)^2
if denom eq 0 then message, 'Incompatible grism parameters'
dgdw = deriv * (gme - nraprime) / denom
;
;  Calculate the grism parameter, and use it to calculate the wavelength.
;
intermediate = cdelt*pixrel
gamma = -tan(theta) + dgdw * factor * intermediate
gamma = atan(gamma) + gamma0 + theta
denom = gme - nraprime
if denom eq 0 then message, 'Incompatible grism parameters'
lambda = (nra - nraprime * lambda0 + sin(gamma)) / denom
;
;  Unless the coordinate type is AWAV, correct for the distinction between
;  vacuum and air wavelengths.
;
if keyword_set(gra) and (coord_type ne 'AWAV') then lambda = lambda * $
  (1.d0 + 1.d-6 * (apar[0] + apar[1]/lambda^2 + apar[2]/lambda^4))
;
;  Convert into the final spectral coordinate variable.
;
case coord_type of
    'FREQ': s = c / lambda
    'ENER': s = h * c / lambda
    'WAVN': s = 1.d0 / lambda
    'VRAD': begin
        freq = c / lambda
        s = c * (restfrq - freq) / restfrq
    endcase
    'WAVE': s = lambda
    'VOPT': s = c * (lambda - restwav) / restwav
    'ZOPT': s = (lambda - restwav) / restwav
    'AWAV': if keyword_set(gra) then s = lambda else s = lambda * $
      (1.d0 + 1.d-6 * (apar[0] + apar[1]/lambda^2 + apar[2]/lambda^4))
    'VELO': s = c * (lambda^2 - restwav^2) / (lambda^2 + restwav^2)
    'BETA': s = (lambda^2 - restwav^2) / (lambda^2 + restwav^2)
endcase
;
;  Apply the appropriate units conversion, and return.
;
return, s / factor
end