;+
; NAME:
;       EIS_GETMOMENTS
;
; PURPOSE:
;       EIS_GETMOMENTS is a script to calculate moments from one,
;       several or all spectral line windows from a (fits or CCSDS) data file. 
;
; CATEGORY:
;       Hansteen/Wikstol EIS data analysis SW
;
; CALLING SEQUENCE:
;       eis_getdata, file, iwin, int, vel, wid, (iwin=iwin, all=all,
;                    /gauss, /moment, /line_def, /lamtype, /ccsds)
;
; INPUTS:
;       file:  Name of file (string) containing EIS data (and
;              headers).
;
; KEYWORD PARAMETERS:
;       iwin: The line number indexes (starting with 0, max 24.)
;       all: process all line windows
;       gauss: set this to use gaussian fits to calculate int, vel, wid.
;       moment: (Default) Use moments to calculate int, vel, wid
;       line_def: Use pre-stored pixel definitions for line and
;                 continuum. These are stored in text files, as pixels defining
;                 where on the CCD the spectral line and it's (to be subtracted)
;                 continuum is. If the line ID is 'Line 0' then these pixel
;                 definitions should be stores in a file called
;                 'line_def_0.tmp'. If no line definitions are available, the
;                 user will be presented with a widget where an average line
;                 profile is plotted, and the pixel defintions are set
;                 interactively. 
;       lamptype: Can take two values: lamtype='pix' (DEFAULT) for
;                 pixel wavelength scale or lamtype='wav' for wavelength scale
;                 in Ångstrøm.  
;       ccsds: set this keyword if telemetry file. Standard is fits.
;
; OUTPUTS:
;       int: intensity (0th order moment). Normally int[x,y, nwin]
;       vel: velocity  (1st order momemt). Normally vel[x,y, nwin]
;       wid: width     (2nd order moment). Normally wid[x,y, nwin]
;
; CALLS:
;       EIS_DATA__READ
;       XCFIT (OPTIONAL)
;
; COMMON BLOCKS:
;
;
; PROCEDURE:
;	example 1: Calculate int, vel wid, from Line window 0 
;                  from CCSDS file1= 'file1.dat'. Use pre-defined 
;                  pixel values:
;                  IDL> eis_getmoments, file1, int,vel,wid, /ccsds,$
;                       /line_def
;       example 2: Get line windows 2 and 3 from fitsfile
;                  file2='file2.fits'. No predined pixel defintions
;                  exist, and use wavelength scale in stead of
;                  pixels. Use gaussian fits:
;                  IDL> eis_getmoments, file2,int,vel, wid,
;                  iwin=[2,3], lamtype='wav', /gauss
; RESTRICTIONS:
;
;       With the /gauss keyword, the program calls 'XCFIT' which is a
;       part of the CDS package distributed with solar soft. 
;
; MODIFICATION HISTORY:
;       01-Mar-2004: Oivind Wikstol.
;       15-Dec-2004: Oivind Wikstol.
;-
pro eis_getmoments, file, int, vel, wid, iwin =iwin, all = all, gauss = gauss, $
                    moment = moment, ccsds = ccsds, lamtype = lamtype, $
                    line_def = line_def

  if n_params() lt 1 then begin
    print, 'eis_getmoments, file,[int, vel, wid], (iwin = iwin, /all'
    print, '                ,/ccsds, /gauss, /moment, /line_def'
    print, '                , lamtype = lamtype'
    return
endif

; create data object:
  data = obj_new('eis_data', file, datasource = datasource)
  ; check how many line windows there are in this data set:
  nwin = data-> getnwin()

  if keyword_set(all) then iwin = indgen(nwin)

  ; if both iwin and all is NOT set: 
  if n_elements(iwin) eq 0 then begin
    print, 'eis_getmoments: No lines specified. Reading line 0'
    iwin=0
 endif

  if keyword_set(ccsds) then datasource = 'ccsds' else datasource = 'fits'
  if keyword_set(gauss) then fit = 'gauss' else fit = 'moment'
  
  if max(iwin) gt nwin then begin
    print, 'eis_getmoments: only '+strtrim(nwin, 2)+' line(s) in this file'
    return
  endif

  if min(iwin) lt 0 then begin
    print, 'eis_getmoments: minimum value for iwin is 0'
    return
  endif

  
  ; extract data 
  xnwin = n_elements(iwin)
  nslit = (data-> getyw())[0]
  nraster = data->getnraster()
  int = fltarr(xnwin, nslit, nraster)
  vel = int
  wid = int

  for i = 0, xnwin -1 do begin
  ; calculate moments
    if keyword_set(line_def) then begin
      filedef = 'line_def_'+strtrim(iwin[i],2)+'.tmp
      finfo = file_info(filedef)
      if(finfo.exists) then begin
        moment = obj_new('eis_moment', data, iwin = iwin[i], fit = fit, $
                          lamtype = lamtype, line_def = filedef)
    endif else begin
      print,'No pixel definitions. Define manually'
      moment = obj_new('eis_moment', data, iwin = iwin[i], fit = fit, $
                            lamtype = lamtype)
    endelse
  endif else begin
      print,'No pixel definitions. Define manually'
      moment = obj_new('eis_moment', data, iwin = iwin[i], fit = fit, $
                            lamtype = lamtype)
  endelse

;   extract moments:
    wd = moment-> getvar(iwin[i])
    int[i, *, *] = reform(wd[*, *, 0])
    vel[i, *, *] = reform(wd[*, *, 1])
    wid[i, *, *] = reform(wd[*, *, 2])
  endfor

  ; take away nwin dimension if only 1 line window:
  int = reform(int)
  vel = reform(vel)
  wid = reform(wid)

  return
end