;+
; Project     : SOHO - CDS     
;                   
; Name        : WAVE2PIX
;               
; Purpose     : Calculate the CDS detector pixel given a wavelength.
;               
; Explanation : Uses dummy transformations to translate a wavelength to a
;               pixel value for any of the CDS's six spectral regions.
;               
; Use         : IDL> pixel = wave2pix(spectrum,wavelength [,/limit])
;    
; Inputs      : spectrum  -  the spectrum identifier (string).  Only the first
;                            and last characters are used to identify the
;                            spectrum so inputs such as NIS1 and N1 are
;                            equivalent.  Valid entries are NIS1, NIS2, GIS1,
;                            GIS2, GIS3, GIS4 (and their abbreviations).
;               
;               wavelength - the wavelength to translate. Can be an array.  
;               
;
; Opt. Inputs : None
;
; Outputs     : The function value returned is the corresponding pixel.
;               A value of -1 is returned for an error condition.
;               
; Opt. Outputs: None
;               
; Keywords    : LIMIT  - if present an appropriate limiting pixel is returned
;                        for out of range wavelengths, else -1 is returned.
;
;		NOLIMIT- If input wavelength is out of range, then return a
;			 value anyway.
;
;               SLIT   - slit number to cater for slit dependence
;
;               REAL   - if given the output is returned as a floating point
;                        value
;
; Calls       : None
;               
; Restrictions: Only dummy transformations at present, to be updated when
;               calibrations are known.
;               
; Side effects: None
;               
; Category    : Calibration, GDS, VDS, Wavelength
;               
; Prev. Hist. : None
;
; Written     : C D Pike, RAL, 28-May-1993
;               
; Modified    : Change of constants in VDS wavelengths  CDP 16-Nov-93
;               Include Limit keyword, CDP, 2-Feb-94
;               Updated wavelength ranges.  CDP, 30-Jan-95
;               Replaced calls to NINT by ROUND.  CDP, 17-Jun-95
;               Incorporate first in-flight results by BJIB.  CDP, 27-Feb-96
;               Use common for communication of coefficients. CDP, 11-Mar-96
;               Make NIS quadratic.  23-Jul-96
;               Add REAL keyword.    3-Sep-96
;		Version 10, 13-Jan-1998, William Thompson, GSFC
;			Added keyword nolimit
;		Version 11, 19-Jan-2000, William Thompson, GSFC
;			Allow SPECTRUM to be an array
;
; Version     : Version 11, 19-Jan-2000
;-            

function wave2pix, spectrum, w, limit=limit, nolimit=nolimit, slit=slit, $
	real=real


;
;  common for communication of coefficients
;
common cds_wavecal, ncoff, gcoff

;
;  check parameters
;
if n_params() lt 2 then begin
   print,"Use: IDL> l = wave2pix(spectrum_id, wavelength [,limit]) "
   print,"eg   IDL> l = wave2pix('NIS1',[310.123,311.234])"
   return,-1
endif

;
;  If spectrum was passed as an array, and matches the dimensions of w,
;  then call wave2pix individually for each point.
;
if n_elements(spectrum) gt 1 then begin
    if n_elements(spectrum) eq n_elements(w) then begin
	result = 0.*w
	for i=0,n_elements(result)-1 do result(i) = wave2pix(spectrum(i), $
		w(i), limit=limit, nolimit=nolimit, slit=slit, real=real)
	return, result
    end else begin
	print, 'SPECTRUM_ID must be either a scalar'
	print, 'or have the same number of elements as PIXEL'
	return, 0
    endelse
endif

;
;  validity flag
;
valid = 1

;
;   parse the spectrum identifier
;
s = strtrim(spectrum,2)
sg = strupcase(strmid(s,0,1))
snum = fix(strmid(s,strlen(s)-1,1))

case sg of
  'N': begin
        if datatype(ncoff,1) ne 'Structure' then begin
           if load_wavecal('n') then begin
              print,'Loaded default NIS coefficients.'
           endif else begin
              print,'Error loading default NIS coefficients.'
              return, -1
           endelse
        endif

        c = ncoff.coeff
 
        case snum of
           1   : begin
                  cm = w - c(0,0,0)
                  b =      c(0,1,0)
                  a =      c(0,2,0)
                 end
;
           2   : begin
                  cm = w - c(1,0,0)
                  b =      c(1,1,0)
                  a =      c(1,2,0)
                 end
           else: begin
                  bell
                  print,'Normal incidence spectrum number not recognised.'
                  print,'Value was ',snum
                  return,-1
                 end
        endcase
       end
  'G': begin
        if datatype(gcoff,1) ne 'Structure' then begin
           if load_wavecal('g',gset_id=!def_gset_id) then begin
              print,'Loaded default (gset_id='+trim(!def_gset_id)+$
                     ') GIS coefficients.'
           endif else begin
              print,'Error loading default GIS coefficients.'
              return,-1
           endelse
        endif

        c = gcoff.coeff
 
        case snum of
           1   : begin
                  cm = w - c(0,0,0)
                  b =      c(0,1,0)
                  a =      c(0,2,0)
                 end
;
           2   : begin
                  cm = w - c(1,0,0)
                  b =      c(1,1,0)
                  a =      c(1,2,0)
                 end
;
           3   : begin
                  cm = w - c(2,0,0)
                  b =      c(2,1,0)
                  a =      c(2,2,0)
                 end
;
           4   : begin
                  cm = w - c(3,0,0)
                  b =      c(3,1,0)
                  a =      c(3,2,0)
                 end
;
           else: begin
                  bell
                  print,'Grazing incidence spectrum number not recognised.'
                  print,'Value was ',snum
                  return,-1
                 end
        endcase
;
       end
 else: begin
         bell
         print,'Spectrum identifier not recognised. Must be N or G.'
         print,'Character was ',sg
         return,-1
       end
endcase

;
;  invert the forward calibration
;
if a eq 0.0 then begin
   if not keyword_set(real) then begin
      pixel = fix(round(cm/b))
   endif else begin
      pixel = cm/b
   endelse
endif else begin
   if not keyword_set(real) then begin
      pixel = fix(round((-b + sqrt(b*b+4.*a*cm))/2./a))
   endif else begin
      pixel = (-b + sqrt(b*b+4.*a*cm))/2./a
   endelse
endelse
small = where(pixel lt 0)
big = where(pixel gt 2047)
if small(0) ge 0 or big(0) ge 0  then valid = 0

;
;  check limits on returned value
;
if not valid then begin
   if not keyword_set(limit) then begin
      if not keyword_set(nolimit) then begin
         bell,2
         print,'Invalid wavelength for spectrum.  Approximate limits are:'
         print,' '
         print,' NIS1  305 < wave < 379'
         print,' NIS2  513 < wave < 633'
         print,' '
         print,' GIS1  151 < wave < 221'
         print,' GIS2  256 < wave < 341'
         print,' GIS3  392 < wave < 492'
         print,' GIS4  657 < wave < 785'
         print,' '
         return, -1
      endif
   endif else begin
      if small(0) ge 0 then pixel(small) = 0
      if big(0) ge 0 then begin
         if sg eq 'N' then pixel(big) = 1023 else pixel(big) = 2047
      endif
   endelse
endif

;
;  successful return
;
if n_elements(pixel) eq 1 then pixel = pixel(0)
return, pixel

end