PRO fg_reg_wave, iindex, idata, oindex, odata, wave_ref=wave_ref, $
  interp=interp, cubic=cubic, p_arr=p_arr, q_arr=q_arr, $
  display=display, quiet=quiet

;+
; NAME:
;       fg_reg_wave
; PURPOSE:
;       Register an FG image or image array to a common wavelength reference
; SAMPLE CALLING SEQUENCE:
;	fg_reg_wave, fg_reg_wave, iindex, idata, oindex, odata
; INPUTS:
;	iindex -   Index record(s) corresponding to input data
;	idata -	   Input data cube, unregistered
; OPIONAL KEYWORDS:
;	wave_ref
;	interp
;	cubic
; OUTPUTS
;	oindex -  Output index records
;	odata -   Output data array, scaled and shifted to the
;               reference wavelength. Default is 450.5nm Blue Continuum. 
; OPTIONAL OUTPUTS
;	P_ARR - Polynomial coefficients for x coordinate transformation via poly_2d
;	Q_ARR - Polynomial coefficients for y coordinate
;               transformation via poly_2d
;
; COMMON BLOCKS: none
;
; RESTRICTIONS:
;       Requires running under SolarSoft. Requires
;       FG_GET_REG_COEFF.PRO.
;
;MODIFICATION HISTORY
;v0.9  Written by G. Slater, LMSAL, Oct. 2007
;v1.0  Modifications to sum & bin checking, T. Berger, LMSAL
;      18-Oct-2007.
;v1.1  Added coefficient file bomb-out. TEB, LMSAL 19-Oct-2007.
;-
;-----------------------------------------------------------------------------
;;Name, version, and timing information
prognam = 'FG_REG_WAVE.PRO'
progver = 'V1.1'

; Set defaults section:
; ---------------------
if N_ELEMENTS(wave_ref) eq 0 then wave_ref = 'blue' 
if N_ELEMENTS(interp) eq 0 then interp = 2
if N_ELEMENTS(cubic) eq 0 then cubic = -0.5
; ---------------------

wave_arr    = ['gband', 'blue', 'red', 'green', 'ca', 'cn', 'mag']
wave_string = ['G band 4305', 'blue cont 4504', 'red cont 6684', $
               'green cont 5550', 'Ca II H line', 'CN bandhead 3883', 'TF FE I 6302']

if get_logenv('SOT_REG_COEFF') eq '' then begin
   MESSAGE,'Coefficient file not found. Returning...'
   RETURN
end

nimg = N_ELEMENTS(idata[0,0,*])
;the following assumes that summing and binning are mutually exclusive
;operations. I.e. 2x2 binning will never occur on a 2x2 summed image.
sum = 1.
if iindex[0].camssum ne 1 then begin
   if iindex[0].campsum ne iindex[0].camssum then begin
      MESSAGE,'Serial and parallel summming differs. Returning...'
      RETURN
   end 
   sum = FLOAT(iindex[0].camssum)
end
if iindex[0].fgbinx ne 1 then begin
   if iindex[0].fgbiny ne iindex[0].fgbinx then begin
      MESSAGE,'X and Y binning differ. Returning...'
      RETURN
   end
   sum = FLOAT(iindex[0].fgbinx)
end

;Currently not sure if this routine works with sub-field
;read-outs. Print warning...
if iindex[0].fgccdix0 ne 0 or iindex[0].fgccdix1 ne 4095 $
   or iindex[0].fgccdiy0 ne 0 or iindex[0].fgccdiy1 ne 2047 $
then begin 
   MESSAGE,/INFO,'Warning: FG_REG_WAVE is currently tested for full FOV images only. Use with caution on sub-field readouts...'
end

oindex = iindex
odata = idata

for i=0,nimg-1 do begin

  index0 = iindex(i)
  wave_string0 = index0.wave
  ss_match0 = WHERE(wave_string eq wave_string0, n_match0)
  if n_match0 eq 1 then wave0 = (wave_arr(ss_match0))[0] else $
    stop, ' No match to wavelength found.  Stopping.'

  if wave0 eq wave_ref then begin
    odata(0,0,i) = idata[*,*,i]
  endif else begin
    t0 = anytim(index0.date_obs,/int)
    pq_cube = fg_get_reg_coeff(t0, wave0=wave0, wave_ref=wave_ref, filnam_geny=filnam_geny, $
      p_arr=p0, q_arr=q0)
    p0[0,0] /= sum
    q0[0,0] /= sum
    p0[0,1] = 1. / (1. + (1./p0[0,1] - 1.)/sum)
    q0[1,0] = 1. / (1. + (1./q0[1,0] - 1.)/sum)

    if N_ELEMENTS(p_arr) eq 0 then p_arr = p0 else p_arr = [[[p_arr]],[[p0]]]
    if N_ELEMENTS(q_arr) eq 0 then q_arr = q0 else q_arr = [[[q_arr]],[[q0]]]

    odata(0,0,i) = $
      POLY_2D(REFORM(idata[*,*,i]), p0, q0, interp, cubic=cubic, missing=missing)

  endelse

endfor

RETURN
END