; program to add write velocities and line widths to eis simulated
; data;
; 30-Apr-2009 A.Gardini Fixed 1-pixel shift error
;
; read velocities into variable v:
restore,'~/solarb/eis/idl/eis_sim/velocity.dat'
;
; read eis fits file with 5 lines:
;
cc=2.997929e8
hdr=obj_new('fits_hdr')
data=obj_new('eis_data',datasource='fits')
data-> read,'~/solarb/eis/data/science/fits/eis_sim0.fits'
hdr-> read,'~/solarb/eis/data/science/fits/eis_sim0.fits'

line_id=hdr->getline_id()
line_wvl=hdr->getline_wvl()
nwin=data->getnwin()

; hardcoded temperatures of formation for the 5 lines as well as
; atomic weight for the lines:
;
te = [1.2e6, 1.3e6, 5.e4, 1.e6, 1.5e6]  ; Fe XII, Fe XIII, He II, Si X, Fe XV
awgt = [9.27426e-23, 9.27426e-23, 6.64226e-24, 4.66453e-23, 9.27426e-23]

; get wavelength scale
lam_str=data->getlambda()

for i=0,nwin-1 do begin
  data->getwin,i,wd,pos
  if pos[0] gt 2147 then begin
    pos[0] = pos[0] - 2148
    det_a=1
    lambda=lam_str.scale_a
  endif else begin
    lambda=lam_str.scale_b
  endelse 
  lstart=pos[0]   ; line start pixel
  lstop=pos[0]+pos[1]-1
  llam=lambda[lstart:lstop]
  ; find central wavelength of line from max of mean line profile:
  mline = mean_spec(wd)
  lmax = max(mline, mpix)
  lambda0 = llam[mpix]
  ny0=double(cc/lambda0*1.e10)       ; angstrom to meter
  ny=double(cc/llam*1.e10)           ; angstrom to meter
  dny = ny-shift(ny, 1)
  dny[0] = dny[1]
  nny = (size(ny))[1]
  ; calculate phiny with velocities and non-thermal widths:
  chi = 20.    ; non-thermal line width
  u = v
  phiny, ny0, ny, u, te[i], awgt[i], chi, phi
  ; compute total line intensity in old line profile:
  wtot = total(wd(mpix-3:mpix+3), 1)
  wnew = wd
  for j = 0, nny-1 do $
    wnew[j, *, *] = phi[j, *, *]*wtot*abs(dny[j]) > wd(j, *, *)
  ; normalize:
  t1 = total(wd)
  t2 = total(wnew)
  wnew = wnew*t1/t2
  data-> setvar, wnew, i
  print, 'Processing line window ' +strtrim(i, 2)
endfor

; write new data back to fits file:
  eis_modfits, data, hdr, fitsfile = '~/solarb/eis/data/science/fits/eis_1.fits'

end
