pro fitsurf,basepath,outfilname

;  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

;  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 = 4	;  degree of polynomial fit
 
;  start loop over files
for kk = 0,nftot-1 do begin

;  read in the data, make sure overflow etc is compensated
readl1_sbsp,pthnam(kk),dat,hdr
temp = dat(*,*,0)
;nx = sizeof(temp,1) & ny = sizeof(temp,2)
sz = size(temp)
nx = sz(1) & ny = sz(2)
nxm1 = nx - 1

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

if kk eq 0 then begin
;  set up arrays for display
	bstks = bytarr(nx*8,ny)
	window,1,xs=nx*8,ys=ny
	scoef = dblarr(ndeg+1,3,nftot)
endif

;;  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
soveri = fltarr(ny,4)
for jj = 0,ny-2 do begin
	itemp = temp(*,jj)
	whr = where(itemp gt median(itemp))
	avgi = mean_sbsp(itemp(whr))
	for istks = 1,3 do begin
		sdat = 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

;  loop through Stokes QUV parameters, get fit coefficients along slit
yyy = findgen(ny)
for istks = 1,3 do begin
;  bypass 20 points at each end of the slit because of possible vignetting
	scoef(*,istks-1,kk) = poly_fit(yyy(20:ny-20),soveri(20:ny-20,istks),ndeg)
;  recover the fit
	smth = fltarr(nx,ny)
	smth(*,*) = 0.
	for jj = 0,ny-1 do begin
		for mj = 0,ndeg do begin
			smth(*,jj) = smth(*,jj) + $
			scoef(mj,istks-1,kk) *  yyy(jj)^mj
		endfor
	endfor

;  use red end of spectrum where there is more continuum
	xtalk(kk,istks-1) = mean(smth(99:106,10:1000))
	xtalk100(kk,istks-1) = mean(smth(99:106,50:150))
	xtalk500(kk,istks-1) = mean(smth(99:106,450:550))
	xtalk900(kk,istks-1) = mean(smth(99:106,850:950))

;  setup final result display
;  insert the raw data in the first 4 columns
	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)
	xstrt = nx*4 + istks*nx
	bstks(xstrt:xstrt+nxm1,*) = $
	bytscl( (dat(*,*,istks) -smth*dat(*,*,0)),max=mxx,min=mnn)

;  end loop over Stokes parameters
endfor
;window,1,xs=nx*8,ys=ny
tvscl,bstks


;  end loop over files
endfor


;  sort the observations in terms of ascending slit position
isrt = sort(sltpos)
sltpos = sltpos(isrt)
scoef = scoef(*,*,isrt)



;  do boxcar averaging in the slit scan direction
;  optimal smoothing seems to be about 200 scan steps
;  determine smoothing parameter for the present data
nsmoothb = fix( 200.*nftot/(max(sltpos)-min(sltpos)) )
print,' boxcar smoothing parameter in measured slit position steps: ',nsmoothb
;  smooth the fit coefficients

;  bypass smoothing because small changes in high order terms will ruin the fit
goto,jump453
for istks = 0,2 do begin
	for jj = 0,ndeg-1 do begin
			temp = smooth(scoef(jj,istks,*),nsmoothb,/edge_truncate)
			scoef(jj,istks,*) = temp
	endfor
endfor
jump453:



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

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




return
end