pro genresidcross_sbsp,basepath,outfilname

;  Routine adapted from fitsurf.pro, the previous version.  Name
;  of program was no longer appropriate

;  routine modified from previous versions of fitsurf, 15 Aug 2011
;  This routine fits the variation of the continuum crosstalk in
;  the direction along the slit only.  Present version fits with 
;  a 4th-order polynomial (ndeg = 4)

;  Smooths the variations of the fit parameters with a running
;  boxcar using the IDL smooth proceedure

;  INPUT:	
;	basepath	= pathname to level1 data calibrated without
;					the residual crosstalk correction
;   outfilname  = the name for the curent output file, usually named
;       after the date of the flat field operation used to generate
;       the crosstalk data, for example residxtalk_20120702_114140

; revised version to use the 2012 eclipse season flat results

;  MUST BE RUN UNDER SSW!

;  OUTPUT:	
;	scoef	 = array of coefficients(3 Stokes,ndeg,nfiles)
;	outputs .geny file to be used by SP_PREP

;  This routine avoids the locations around the lines where
;  the polarization might bias the fit
;  Also modified Nov. 2012 to avoid fitting vignetted regions at upper
;  and lower ends of slit
;  Modified for smoothing the raw residual crosstalk image for 2013
;  data because amplitude of residual crosstalk is growing rapidly

;  determines surface fit coefficients of Q/I, U/I, V/I
;  for the full array

;  get the list of input files
derf = file_search(basepath + '*.fits')
nftot = n_elements(derf)
print,'directory: ',basepath,'; number of files = ',nftot
pthnam = derf

;  output array containing slit positions of data
;  The variables below are for checking the results only
sltpos = fltarr(nftot)
xtalk = fltarr(nftot,3)
xtalk100 = xtalk & xtalk500 = xtalk & xtalk900 = xtalk

ndeg = 7	;  degree of polynomial fit


;  first cycle through the images to determine the slit scan position
;  order
for kk = 0,nftot-1 do begin
;  read in the data, make sure overflow etc is compensated
	readl1_sbsp,pthnam(kk),dat,hdr

;  define crosstalk image
	if kk eq 0 then begin
		temp = reform(dat(*,*,0))
		sz = size(temp)
		nx = sz(1) & ny = sz(2)
		nxm1 = nx - 1
;  set up arrays for display
		bstks = bytarr(nx*8,ny)
		scoef = dblarr(ndeg+1,3,nftot)
;  define working arrays
		xtalkimg = fltarr(nftot,ny,3)
		xtalkraw = fltarr(nftot,ny,3)
		xtalkfit = fltarr(nftot,ny,3)
		soveri = fltarr(ny,4)
	endif

;  get the slit position of observation from the header
	sltpos(kk) = sxpar(hdr,'SLITPOS')
endfor

;  sort the observations in terms of ascending slit position
isrt = sort(sltpos)
;  test to see if sequence of files is in correct order of ascending slitpos
iseq = lindgen(nftot)
iorder = 0
idiff = isrt-iseq
whrdif = where(isrt-iseq ne 0L,countdd)
if countdd gt 0 then begin
	iorder = 1
	sltpos = sltpos(isrt)
	print,'file order not in ascending slit position.  Sequence re-ordered'
endif


;  start processing loop over files
for ksl = 0,nftot-1 do begin
	if iorder ne 0 then kk = isrt(ksl) else kk = iseq(ksl)

;  read in the data, make sure overflow etc is compensated
	readl1_sbsp,pthnam(kk),dat,hdr
	temp = reform(dat(*,*,0))

;  display raw data first
	bstks(0:nxm1,*) = bytscl(dat(*,*,0))                 

;  now make variations of Q/I, U/I, V/I along slit for continuum
;  pixels
	for jj = 0,ny-2 do begin
		itemp = reform(temp(*,jj))
		whr = where(itemp gt median(itemp))
		avgi = mean_sbsp(itemp(whr))
		for istks = 1,3 do begin
			sdat = reform(dat(*,jj,istks))
			sdat = sdat(whr)
			soveri(jj,istks) = mean_sbsp(sdat)/avgi
		endfor
;  put mean intensity into the first slot
		soveri(jj,0) = avgi
	endfor

	for istks = 1,3 do begin
;  load the xtalk images
		for jj = 0,ny-1 do begin
			xtalkraw(kk,jj,istks-1) = soveri(jj,istks)
		endfor
	endfor

;  end loop over files
endfor



;  smooth the raw results, the fit with 1-D polynomial at each slit position
;  use smoothing width of 10 columns in scan direction, 50 rows along slit
yyy = findgen(ny)
smth = fltarr(nftot,ny,3)
xtcorr = fltarr(nx,ny)
for istm1 = 0,2 do begin
	istks = istm1+1
	temp = reform(xtalkraw(*,*,istm1))
	stemp = smooth(temp,[10,50])
	xtalkimg(*,*,istm1) = stemp
	for kk = 0,nftot-1 do begin
;  use red end of spectrum where there is more continuum
		xtalk100(kk,istm1) = mean(xtalkimg(99:106,50:150,istm1))
		xtalk500(kk,istm1) = mean(xtalkimg(99:106,450:550,istm1))
		xtalk900(kk,istm1) = mean(xtalkimg(99:106,850:950,istm1))
	endfor
endfor


;  cycle through the data again to review the results of the correction
xtalkr = fltarr(nx,ny)
window,1,xs=nx*8,ys=ny
for ksl = 0,nftot-1 do begin
	if iorder ne 0 then kk = isrt(ksl) else kk = iseq(ksl)
	readl1_sbsp,pthnam(kk),dat,hdr

;  setup final result display
;  insert the raw data in the first 4 panels
	for istks = 0,3 do begin
		xstrt = istks*nx
		mxx = max(dat(*,*,istks)) & mnn = min(dat(*,*,istks))
		if istks gt 0 then mxx = 100. & mnn = -100.
		bstks(xstrt:xstrt+nxm1,*) = bytscl(dat(*,*,istks),max=mxx,min=mnn)
	endfor
;  insert the corrected data in the last 4 panels
	for istks = 0,3 do begin
		xstrt = nx*4 + istks*nx
		if istks eq 0 then begin
			bstks(xstrt:xstrt+nxm1,*) = bstks(0:nxm1,*)
			goto,jump37
		endif else begin
			for jj = 0,ny-1 do xtalkr(*,jj) = xtalkimg(kk,jj,istks-1)
			bstks(xstrt:xstrt+nxm1,*) = $
				bytscl( (dat(*,*,istks) -xtalkr*dat(*,*,0)),max=mxx,min=mnn)
		endelse
		jump37:
	endfor
	tvscl,bstks
	wait,0.4
endfor

;  saving original crosstalk image (xtalkraw), the spatially smoothed image
;  (xtalksm), and the image reconstructed from polynomial fits along the slit (smth)

save,filename='tempjunk.save',xtalkraw,xtalkimg,smth,braw,bsm,bsmth, $
	ndeg,sltpos,nftot

;  output the results
;save,filename='coef_xtalk_yy.save',nftot,scoef,sltpos,xtalk,ndeg
save,filename='residcrossamples.save',xtalk100, xtalk500, xtalk900,xtalkimg, $
	xtalkraw,sltpos,nftot

;  create the .geny file
savegenx,file=outfilname,xtalkimg,sltpos,nftot



return
end