;   MODIFICATION HISTORY
;          18-Mar-2008 A. Gardini - Changed the column's order of moments.
;           4-Jul-2008 A. Gardini - Added the init keyword.
;          11-Jul-2008 A. Gardini - Added group_leader in xmoment's call.
;          19-Sep-2008 A. Gardini - Added the error of intensity
;          23-Sep-2008 A. Gardini - Added the error of velocity
;          12-Mar-2009 V. Hansteen - Added possibility of setting aeffrelerr.
;          18-Jun-2011 V. Hansteen - Changed orbital variation and
;                                    tilt correction to use Kamio and
;                                    Fredviks method.
; $Id: eis_moment__moment.pro 13809 2022-04-03 13:26:04Z steinhh $
pro eis_moment::moment, data, indx, init=init
  if n_elements(init) eq 1 then return

cc = 2.99792e5  ; light speed
missing=-100.

; read line window data
  data-> getwin, indx, wd, pos
  err=data->geterr(indx)
  badpix=err*0.
  badpix[where(err eq missing)]=-100.
  nwave = (data-> getxw())[indx]
  nspat = (data-> getyw())[indx]
  nexp = data-> getnexp()
  nraster = data-> getnraster()

;; get wavelength scales
;  lstr = data-> getlambda(lamtype=self.lam_type) ; gives scale_a and scale_b
;
;; get linear dispersion:
;  cal = obj_new('eis_cal')
;  dispstr = cal-> getdisp()
;  obj_destroy,cal
;
;  if pos[0] gt 2147 then begin
;    lambda = lstr.scale_a
;    disp = dispstr.a
;    pos[0] = pos[0]-2148
;  endif else begin
;    lambda = lstr.scale_b
;    disp = dispstr.b
;  endelse
;  lambda = lambda[pos[0]:pos[0]+pos[1]-1]
  filename = data -> getfilename()
  header = data -> getwinheader(filename,indx)
  if pos[0] gt 2147 then begin
    pos[0] = pos[0]-2148-50
    ccd=0
  endif else begin
    pos[0] = pos[0]-50
    ccd=1
  endelse
  lambda = eis_pixel_to_wavelength(pos[0]+indgen(pos[1]),ccd,header.slit_ind)

; Launch widget tool to mark line and continuum pixels
  if self.line_def then begin
    file='line_def_'+strtrim(string(indx), 2)+'.tmp'
    on_ioerror, nofile
    restore, file
  endif else begin
    xmoment, data, wd, indx, lambda, groupl = self.group_leader
    if not(obj_valid(data)) then return ; if the xcontrol window closes
    line_px = (data-> getwd_def())[indx].line_px
    cont_px = (data-> getwd_def())[indx].cont_px
    save, line_px, cont_px, file='line_def_'+strtrim(string(indx), 2)+'.tmp'
  endelse
  lpx = indgen(line_px[1]-line_px[0]+1)+line_px[0]
  cpx = indgen(cont_px[1]-cont_px[0]+1)+cont_px[0]
  nlpx= float(n_elements(lpx))
  ncpx= float(n_elements(cpx))
  sum = fltarr(nwave)
  sumc= fltarr(nwave)
  sum[lpx] = 1.
  sumc[cpx] = 1.

; set up of moment array completed, now calculate
; set "laboratory wavelength"
; self.p0 = (line_px[0]+line_px[1])/2
; dvpx = (lpx-self.p0)*disp*cc/lambda[self.p0]
; dv = fltarr(nwave)
; dv[lpx] = dvpx
; totdv = total(dv)
  pix = fltarr(nwave) 
  pix[lpx] = lpx 
  pix2= fltarr(nwave) 
  pix2[lpx] = lpx*lpx 
; lambda0 = (lambda[line_px[0]]+lambda[line_px[1]])/2.
  lam = fltarr(nwave) 
  lam[lpx] = float(lambda[lpx])
  lambda1 = lambda[line_px[0]]
  disp=lambda[line_px[0]+1]-lambda[line_px[0]]
   
  hdr=data->gethdr()
  line_id=strtrim((hdr->getline_id())[indx],2)
  refwvl_id=(strsplit(line_id,' ',/extract))[2]
  refwvl=float(refwvl_id)

; moments and related quantities
  intc= fltarr(nspat, nexp)
  int = fltarr(nspat, nexp)
  bad = fltarr(nspat, nexp)
  vel = fltarr(nspat, nexp)
  wid = fltarr(nspat, nexp)
  wdl = fltarr(nwave, nspat, nexp) 
  mu  = fltarr(nspat, nexp) 
  interr= fltarr(nspat, nexp)
  intcerr= fltarr(nspat, nexp)
  err2= fltarr(nwave, nspat, nexp) 
  err2=err*err 
  err2[where(err eq missing)]=0.
  velb = fltarr(nspat, nexp)
  lambdac = fltarr(nspat, nexp)
  velerrfac = fltarr(nwave, nspat, nexp)
  velerr= fltarr(nspat, nexp)
  for i = 0, nexp-1 do begin
    intc[*, i]= reform(sumc#wd[*,*,i])/ncpx ; specif. intensity of the cont.
    intcerr[*, i] = reform(sumc#err2[*, *, i])
    for k = 0, nwave-1 do $
; specif. intens. continuum subtracted
      wdl[k,*,i]= wd[k,*,i]-intc[*,i] > 0.
; 0th order moment
;;   int[*, i] = (reform(sum#wd[*,*,i])-intc[*,i]*nlpx) > 1
;   int[*, i] = reform(sum#wd[*, *, i])-intc[*, i]*nlpx
    int[*, i] = reform(sum#wdl[*, *, i])
    interr[*, i] = reform(sum#err2[*, *, i])
; 1st order moment
;;   vel[*, i] = (dv#reform(wd[*, *, i])-intc[*,i]*totdv)/int[*, i]
    mu[*, i] = reform(pix#wdl[*, *, i])/int[*, i] 
    vel[*,i] = (lambda1+(mu[*,i]-line_px[0])*disp - refwvl)/refwvl*cc ;AG
    lambdac[*, i] = reform(lam#wdl[*, *, i])/int[*, i] 
    velb[*,i] = (lambdac[*,i] - refwvl)/refwvl*cc ;AG
; 2nd order moment
;; wid[*, i] = ((dv*dv)#reform(wd[*, *, i])-intc[*,i]*total(dv*dv))/int[*, i]
    wid[*, i]= reform(pix2#wdl[*,*,i])/int[*,i]-mu[*,i]*mu[*,i]
    bad[*, i] = total(badpix[lpx,*,i],1)+total(badpix[cpx,*,i],1)
  endfor
; Computation of the velocity errors
  for i = 0, nexp-1 do $
    for k = lpx[0], lpx[nlpx-1] do $
      for j = 0, nspat-1 do $
        velerrfac[k,j,i]= lam[k]/int[j,i]*(1.-lambdac[j,i]/lam[k])
  velerrfac=velerrfac*err
  velerrfac=velerrfac^2
  for i = 0, nexp-1 do $
    velerr[*, i] = reform(sum#velerrfac[*,*,i])
  velerr = sqrt(velerr)
  velerr = velerr/refwvl*cc ;AG
; Computation of the intensity errors (square root of the sum of the squares)
  gain=*self.cal->getgain()
  gainerr=*self.cal->getgainerr()
  aeffrelerr=*self.cal->getaeffrelerr()
  interr=sqrt(interr+intcerr*(float(nlpx)/float(ncpx))^2)
  interr=int*sqrt((interr/int)^2+(gainerr/gain)^2+aeffrelerr^2)
; Final scalings:
  int=int*disp
  interr=interr*disp

; Final computation of the line width: the factor 2.35 is introduced to get 
; the FWHM ( = 2*sqrt(2ln2) * sigma for a Gaussian distribution).
  wid = sqrt(wid)*disp*1000.*2.35
; Now put everything in place in the w array
  (*self.w[indx])[0,*,*] = int
  (*self.w[indx])[1,*,*] = vel 
  (*self.w[indx])[2,*,*] = wid 
  (*self.w[indx])[3,*,*] = interr
  (*self.w[indx])[4,*,*] = velerr
;
; correct for orbital variations and tilt
  corr=data->gethkwavecorr(refwvl)
  (*self.w[indx])[1,*,*] = vel-rotate(corr.im/refwvl*cc,3)
; 
; OLD method:
;  self->tilt_correction,data,indx,refwvl
;  self->corrvel,indx,bad 
; set fitslev to 2 (Level 2 fits data)
  fitslev = 2
  self-> setfitslev, fitslev
  return

  nofile:
  ok = dialog_message('No pixel information file')
  return
end