pro fiss_conf_wvcalib_v2, FILE, wvband, refx, dldw, wwc, wvcen

;+
;  Name : FISS_CONF_WVCALIB_v2
;  Purpose :
;       Confirm wavelengths and calibrate wavelengths using the earth Atm line
;
;  Syntax : fiss_conf_wvcalib_v2, FILE, wvband, refx, dldw, wwc, wvcen
;
;  Input :
;      wvband : '6562', Ha line or '8542' Ca II 8542 line
;      file : the spectral profile to be used for the calibration.
;      refx : reference frame
;
;  Output :
;      dldw : angstrom per pixel
;      wwc : pixel position of line center (6562 or 8542)
;      wvcen : angstrom of line center (6562 or 8542)
;
;  Keywords :
;
;  Required routines : fiss_wv_calib
;
;  History
;       ref. fiss_wv_calib.pro made by J. Chae
;       2010 Oct, Dong-uk Song
;
;-

; Input file *************************************************************
;DATA_DIR='/data/FISS_DATA/2010/07/29/comp/fil/'
;FILE=DATA_DIR+'FISS_20100729_192517_B1_c.fts'
;wvband = '8542'
;*************************************************************************

ref=fiss_read_frame(FILE,refx,h)

nx=fxpar(h,'NAXIS3')
ny=fxpar(h,'NAXIS2')
nw=fxpar(h,'NAXIS1')

profile=fiss_read_frame(FILE,refx,h);fiss_read_frame(FILE,i,h)

;*************************************************************************

;; Determine the wavelength center using the average of the spectral profile

  sprofile=total(profile, 2)    ; Average of the spectral profile
  sprofile=sprofile/n_elements(profile[0,*])
  fiss_wv_calib, wvband, sprofile, wvpar
;;**************************************************************************

; Output of the fiss_wv_calib.pro******************************************
  wc=wvpar[0]
;**************************************************************************



case strmid(wvband,0,4) of
  '6562': begin
    wave=[6559.5880d0, 6560.555d0, 6560.68d0, 6561.097d0, 6562.817d0, 6563.521d0, 6565.545]
    wvname=['Ti II','Atm H2O','Si I','Atm H2O','H I','Atm H2O','Atm H2O']

    lines=[wvname[4],wvname[3], wvname[6]]
    lambda_line=[wave[4], wave[3], wave[6]]
    lamc=lambda_line[0]
    dlambda_line=lambda_line-lamc
    dldw=0.019086
    end

  '8542': begin
    wave=[8536.165d0, 8538.0152d0, 8539.888d0, 8540.817d0, 8542.089d0, 8546.222d0]
    wvname=['Si I','Fe I','Atm H2O','Atm H2O','Ca II','Atm H2O']

    lines=[wvname[4], wvname[3],wvname[5]];, wvname[5]];[wvname[4],wvname[3], wvname[5]]
    lambda_line=[wave[4], wave[3], wave[5]];, wave[5]][[wave[4], wave[3], wave[5]]
    lamc=lambda_line[0]
    dlambda_line=lambda_line-lamc
    dldw=-0.0257604

    end
endcase

nw=n_elements(sprofile)
w=findgen(nw)                                         ;;;;;; *(-1.))+(nw-1)
dlambda=(w-wc)*dldw
wl=dlambda_line*0.

; Determine the value of dldw *****************************************************

window, 0
plot, sprofile/max(sprofile)

line=1
mask1 = abs(dlambda -dlambda_line[line]) le 0.15
oplot, mask1, linestyle=line

wtmp1=w[where(mask1)]
ptmp1=sprofile[where(mask1)]

s=where(ptmp1 eq min(ptmp1))

wtmp1=wtmp1[s-4:s+3]
ptmp1=ptmp1[s-4:s+3]

window, 1

c=poly_fit(wtmp1-median(wtmp1), ptmp1, 2)
wl[line]=median(wtmp1)-c[1]/(2*c[2])

plot, wtmp1, ptmp1, linest=1, yst=1, position=[0.05,0.05,0.45,0.95], title='line1'
oplot, wtmp1, poly(wtmp1-median(wtmp1), c), thick=2
oplot, wl[line]+[0,0], [0, 1]*max(ptmp1)

wait, 1
; ----------------------------------------------------------------------------------
line=2
mask2 = abs(dlambda -dlambda_line[line]) le 0.15
wset, 0
plot, mask2, linestyle=line, /noerase

wtmp2=w[where(mask2)]
ptmp2=sprofile[where(mask2)]

s=where(ptmp2 eq min(ptmp2))

wtmp2=wtmp2[s-3:s+4]
ptmp2=ptmp2[s-3:s+4]

c=poly_fit(wtmp2-median(wtmp2), ptmp2, 2)
wl[line]=median(wtmp2)-c[1]/(2*c[2])

wset, 1

plot, wtmp2, ptmp2, linest=1, yst=1, position=[0.55,0.05,0.95,0.95], title='line2', /noerase
oplot, wtmp2, poly(wtmp2-median(wtmp2), c), thick=2
oplot, wl[line]+[0,0], [0, 1]*max(ptmp2)

dldw=(dlambda_line[2]-dlambda_line[1])/(wl[2]-wl[1])

wait, 3
wshow, 1 ,0  ; wshow, 0 or wshow, 0 ,0               ; dldw fitting plot image    If you want to see this window, input "wshow, 1"
wshow, 0, 0  ; wdelete, 0
print, wvname[4], ' -> the value of dldw (A per one pixel) :', dldw
;print, wl
;*********************************************************************************
case strmid(wvband,0,4) of
  '6562': begin
    dlambda=(w-wc)*dldw                    ; dldw -> apply for
    refvalue= lambda_line[1]-lambda_line[0]
    wwc=(-1.)*(refvalue/dldw)+wl[1]
   end

  '8542': begin
    dlambda=(w-wc)*dldw                    ; dldw -> apply for
    refvalue= lambda_line[1]-lambda_line[0]

    ww=(wl*(-1.)+(nw-1.))
    wwc=(refvalue/dldw)+ww[1]
    wwc=(wwc*(-1.))+(nw-1.)
   end
endcase
wvcen=lambda_line[0]
print, 'The line center (Pixel) : ', wwc
end