;+
; NAME:
;       IRESTORE_WMAP()
;
; PURPOSE:
;       Restore an IDL save set two dimensional wavelength map
;       generated by IARCFIT.
;
; CALLING SEQUENCE:
;       wavemap = irestore_wmap(wmapname,datapath=,minwave_in=,$
;          maxwave_in=,dwave=,minwave_out=,maxwave_out=,wheader=,$
;          wset=,/silent)
;
; INPUTS:
;       wmapname - wavelength map to restore [NWMAP]
;
; OPTIONAL INPUTS:
;       datapath    - path name to WMAPNAME
;       minwave_in  - input minimum wavelength
;       maxwave_in  - input maximum wavelength
;       minwave_out - output minimum wavelength
;       maxwave_out - output maximum wavelength
;
; KEYWORD PARAMETERS:
;       silent - suppress messages to STDOUT
;
; OUTPUTS:
;       wavemap - wavelength map corresponding to WMAPNAME [NWMAP] 
;
; OPTIONAL OUTPUTS:
;       dwave   - output wavelength dispersion
;       wheader - simple wavelength map header
;       wset    - trace set corresponding to WAVEMAP 
;
; PROCEDURES USED:
;       CMRESTORE, TRACESET2XY, MKHDR, SXPAR(), SXADDPAR(),
;       DJS_LAXISGEN(), ARM_DOUBLE()
;
; COMMENTS:
;
; EXAMPLES:
;
; MODIFICATION HISTORY:
;       J. Moustakas, 2003 December 7, U of A - excised from
;          ICALIBRATE, documented, and developed substantially  
;       jm05apr19uofa - documentation improved
;
; Copyright (C) 2003, 2005, John Moustakas
; 
; This program is free software; you can redistribute it and/or modify 
; it under the terms of the GNU General Public License as published by 
; the Free Software Foundation; either version 2 of the License, or
; (at your option) any later version. 
; 
; This program is distributed in the hope that it will be useful, but 
; WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
; General Public License for more details. 
;-

function irestore_wmap, wmapname, datapath=datapath, minwave_in=minwave_in, $
  maxwave_in=maxwave_in, dwave=dwave, minwave_out=minwave_out, $
  maxwave_out=maxwave_out, wheader=wheader, wset=wset, silent=silent

    nwmap = n_elements(wmapname)
    if nwmap eq 0L then begin
       print, 'Syntax - wavemap = irestore_wmap(wmapname,datapath=,$'
       print, '   minwave_in=,maxwave_in=,dwave=,minwave_out=,$'
       print, '   maxwave_out=,wheader=,wset=,/silent)'
       return, -1L
    endif

    if n_elements(datapath) eq 0L then datapath = cwd()

; call this routine recursively

    if nwmap gt 1L then begin

       for iw = 0L, nwmap-1L do begin
          wavemap1 = irestore_wmap(wmapname[iw],minwave_in=minwave_in,$
            maxwave_in=maxwave_in,dwave=dwave1,minwave_out=minwave1,$
            maxwave_out=maxwave1,wset=wset1,wheader=wheader1)
          if (iw eq 0L) then begin
             wavemap = wavemap1
             dwave = dwave1
             minwave_out = minwave1
             maxwave_out = maxwave1
             wheader = wheader1
             wset = wset1
          endif else begin
             wavemap = [ [ [wavemap] ], [ [wavemap1] ] ]
             dwave = [dwave,dwave1]
             minwave_out = [minwave_out,minwave1]
             maxwave_out = [maxwave_out,maxwave1]
             wheader = [ [wheader] ,[wheader1] ]
             wset = [ [wset], [wset1] ]
          endelse
          
       endfor
       wset = reform(wset)
       
       return, wavemap
       
    endif

    if file_test(datapath+wmapname,/regular) eq 0L then begin
       splog, 'Wavelength map '+datapath+wmapname+' not found.'
       return, -1L
    endif

    if not keyword_set(silent) then splog, 'Restoring '+datapath+wmapname+'.'
    cmrestore, wmapname, names=['wset','header'], /quiet

; restore the wavelength map in double precision

    nx = long(wset.xmax-wset.xmin+1L)
    ntrace = (size(wset.coeff,/dimension))[1]
    xtemp = djs_laxisgen([nx,ntrace],iaxis=0L) + arm_double(wset.xmin)
    
    traceset2xy, wset, xtemp, wavemap
    wsize = size(wavemap,/dimension)
    ncols = wsize[0] & nrows = wsize[1]
    
; default linearized wavelength scale.  if MINWAVE is not specified
; then choose the minimum starting wavelength.  DWAVE is computed as
; the average dispersion of the middle row down to two decimal places 
    
    dwave = fix(100.0*(max(wavemap[*,nrows/2L])-min(wavemap[*,nrows/2L]))/(ncols-1))/100.0 ; [Angstrom/pixel]

    if n_elements(minwave_in) eq 0L then minwave_in = float(ceil(max(wavemap[0,*])))
    minwave_out = minwave_in > float(ceil(max(wavemap[0,*]))) ; [Angstrom]

    if n_elements(maxwave_in) eq 0L then maxwave_in = minwave_out + dwave*(ncols-1) 
    maxwave_out = maxwave_in < floor(min(wavemap[ncols-1,*]))   ; [Angstrom]

; create a wavelength map header

    mkhdr, wheader, wavemap
    jd = sxpar(header,'JD',count=jdcount)
    ra = sxpar(header,'RA',count=racount)
    de = sxpar(header,'DEC',count=decount)
    if jdcount eq 0L then jd = systime(/julian)
    if racount eq 0L then ra = '00:00:00'
    if decount eq 0L then de = '00:00:00'
    
    sxaddpar, wheader, 'JD', jd, ' Julian date', before='DATE'
    sxaddpar, wheader, 'RA', ra, ' right ascension', before='DATE'
    sxaddpar, wheader, 'DEC', de, ' declination', before='DATE'
    
return, wavemap
end

