pro FSMT_scan_check,path,lambda0,lambda1,step
;Scan check of Huairou Full_disk_center_data
;path : the path where the data is stored
;lambda0 ; star position
;lambda1 ; end position
;step   the scan step (from - to + step>0; from + to - step <0)
wclear

if N_elements(path) eq 0 then path='./'
if N_elements(show) eq 0 then show=0
if N_elements(lambda0) eq 0 then lambda0=-0.1
if N_elements(lambda1) eq 0 then lambda1=0.5
if N_elements(step) eq 0 then step=0.02


cd,path


;some setting for the spline fit line center
;you can changed it depends the situation you met. 
if N_elements(wd) eq 0 then wd=4


file=file_search(path,'*L5.fit')

nn=N_elements(file)

Ic=fltarr(992,992,nn)
mag=fltarr(992,992,nn)
judge=intarr(nn)
judge[*]=1

goto,Recheck

Aligndata:
;---------Align sun disk center-------------------


for i=0,nn-1 do begin
print,'==========================',i,nn,file[i]
l=readfits(file[i],hl)
sww=size(l)
nxo=sww[1]
nyo=sww[2]
nzo=sww[3]
if nyo eq 992 then begin
l=readfits(file[i],hl)
endif


;不同阶段的CCD处理不同
if nyo eq 1004 and nxo ne 1 then begin
l=swap_endian(l)
l=l[0:991,6:997,*]
endif

;没有存储左右旋的数据
if nyo eq 1004 and sww[0] eq 2  then begin
temp=readfits(file[i],hl)
l=fltarr(992,1004,2)
l[*,*,0]=temp
l[*,*,1]=temp
l=l[0:991,6:997,*]
endif
help,l
;print,hl

l1=fltarr(992,992) & l2=fltarr(992,992)


time=sxpar(hl,'startobs')
time=ht2sswt(time)

l1=l(*,*,0) &  ls=congrid(l1,992/2.,992/2.)  

l2=l(*,*,1)

vmax=max(0.5*(l1+l2))

sv=10000.*(l1-l2)/(l1+l2)

im=0.5*(l1+l2)

;rotate
im=rotate(im,5)
sv=rotate(sv,5)

p={nx:992,ny:992,xd:980,yd:980}
center_sobel4,im,p,para

;para=[496,496,975,975]

;first
if para[2] LT 0 or para[2] LT 900 then begin
para[2]=980&para[3]=980
ss=where(l1 lt 0.08*vmax)
fxaddpar,hl,'STATE','unkown'
endif else begin
;imc=im
dist_circle,imc,992,para[0],para[1], /DOUBLE
ss=where(imc GE 0.25*(para(2)+para(3)))
sv(ss)=0.0
im(ss)=0.0
sv=shift(sv,992/2.0-para[0],992/2.0-para[1])
im=shift(im,992/2.0-para[0],992/2.0-para[1])
fxaddpar,hl,'STATE','success'
endelse

fxaddpar,hl,'para0',para[0]
fxaddpar,hl,'para1',para[1]

print,'para0=',para

DATE=strjoin(strsplit(time,'-',/extract))
DATE=strmid(DATE,0,8)
Time0=strmid(time,11,8)

sunparam,'GMT',DATE,time0,rsun,p_angle,b0,cl0,DOY
print,'p_angle',p_angle

sv=rot(sv,p_angle[0],/cubic,/PIVOT)
im=rot(float(im),p_angle[0],/cubic,/PIVOT)

center_sobel4,im,p,para
dist_circle,imc,992,para[0],para[1], /DOUBLE
ss=where(imc GE 0.25*(para(2)+para(3)))
sv(ss)=0.0
im(ss)=0.0
;atv,im
sv=shift(sv,992/2.0-para[0],992/2.0-para[1])
im=shift(im,992/2.0-para[0],992/2.0-para[1])


if sww[0] eq 2 and nyo eq 1004 then begin
sv=l[*,*,0]
para=[sxpar(hl,'CRPIX1'),sxpar(hl,'CRPIX2'),sxpar(hl,'RADIUS')*2,sxpar(hl,'RADIUS')*2]
sv=shift(sv,992/2.0-para[0],992/2.0-para[1])
im=shift(im,992/2.0-para[0],992/2.0-para[1])

sxaddpar,hl,'CRPIX1',992/2.0
sxaddpar,hl,'CRPIX2',992/2.0

endif

fxaddpar,hl,'P_ANGLE',0.0
fxaddpar,hl,'B0',b0[0]
fxaddpar,hl,'FNDLMBXC',992./2.0
fxaddpar,hl,'FNDLMBYC',992./2.0
fxaddpar,hl,'FNDLMBMI',0.25*(para[2]+para[3])
fxaddpar,hl,'CDELT1',rsun[0]/(0.25*(para[2]+para[3]))
fxaddpar,hl,'CDELT2',rsun[0]/(0.25*(para[2]+para[3]))
fxaddpar,hl,'STARTOBS',time[0]
fxaddpar,hl,'t_start',time[0]
fxaddpar,hl,'HRFILE',file_basename(file[i])
print,'resolution=',rsun[0]/(0.25*(para[2]+para[3]))

resolution=rsun[0]/(0.25*(para[2]+para[3]))

if abs(resolution) GE 1.5*2.0 then judge[i]=0


Ic[*,*,i]=im
mag[*,*,i]=sv

if i eq 0 and judge[i] eq 1 then begin
FNDLMBMI=0.25*(para[2]+para[3])
endif 

if i eq 0 and judge[i] eq 0 then begin
FNDLMBMI=480.
;index=fitshead2struct(hl)
endif

endfor

save,judge,FNDLMBMI,Ic,filename='Ic.sav'

Recheck:
if file_exist('Ic.sav') ne 1 then goto,Aligndata

restore,'Ic.sav'
sw=size(Ic)
nx=sw(1)
ny=sw(2)
nz=sw(3)

lincen=fltarr(nx,ny)

y0=linspace(0,nz-1,1.0)
wav=linspace(lambda0,lambda1,step)
num=linspace(0,nz-1,1.0)

LL=where(judge eq 1)

y0=y0[LL]
wav=wav[LL]
num=num[LL]
Ic=Ic[*,*,LL]

;for i=289,289 do begin
;for j=406,406 do begin
for i=496,496 do begin

for j=496,496 do begin
Ic0=Ic[i,j,*]

LMin=where(Ic0 eq min(Ic0))

x0=deriv(num,Ic0)


wd=wd
if N_elements(Lmin) ne 1 or $
         Lmin[0]-wd LT 0 or $
         Lmin[0]+wd GT nz-1 or $
         sqrt((i-496.)^2+(j-496.)^2) GT 0.05*FNDLMBMI then continue
                ;0.05*index[0].FNDLMBMI 
x=x0[Lmin-wd:Lmin+wd]
y=y0[Lmin-wd:Lmin+wd]

pos=interpol(y,x,0.0,/spline)
lincen[i,j]=interpol(wav,num,pos)
;if show eq 1 then begin
;window,/free
;plot,y0,Ic0,psym=-2
;
;window,/free
;plot,y0,x0,psym=-2
;window,/free
;plot,y,x,psym=2
;
;y2=interpol(y,x,x,/spline)
;
;;print,lincen[i,j]
;oplot,y2,x,psym=-6,symsize=2.0
;endif


endfor
endfor

loadct,0
device,decomposed=0
window,/free
erase,255
plot,wav,Ic[496,496,*]/1e+5,psym=-2,xtitle='Wavelength(A)',$
        ytitle='Disk Center Intensity 10!E5',ycharsize=1.5,$
        Title=Strmid(String(lambda0,'(F5.2)'),0,5)+'_'$
             +Strmid(String(lambda1,'(F5.2)'),0,5)+'_'$
             +Strmid(String(step,'(F5.2)'),0,5)+'  Lincen='+$
              Strmid(String(lincen[496,496],'(F5.3)'),0,5),$
             color=0,back=255
oplot,[lincen[496,496],lincen[496,496]],[-100,100],linestyle=1,thick=1.2,color=0
write_jpeg,'scanprof.jpg',tvrd(true=1),true=1,quality=100

print,'Line center position=',lincen[496,496]
end