pro icrossqu_sbsp, stks, sr1, sr2
;+
;
;  procedure:  icrossqu_sbsp
;
;  purpose:  Determine residual i -> Q,U crosstalk from the continuum
;	     spectral region, then apply the inverse of that crosstalk to
;	     the data.  Crosstalk determined and applied on a row-by-row basis
;	     by finding a linear trend of the crosstalk along the slit.
;  THIS ROUTINE CORRECTS ONLY Q,U SPECTRA, SINCE V SPECTRUM HAS NO APPARENT
;  RESIDUAL FEED-THROUGH OF GRANULATION PATTERN IN CALIBRATED DATA
;
;==============================================================================
;
;	Check number of parameters.
;
;;if n_params() ne 10 then begin
if n_params() ne 3 then begin
	print
	print, "usage:  icross_sbsp, stks"
	print
	print, "	Determine residual i -> Q,U,V crosstalk from the"
	print, "	continuum spectral region, then apply the inverse of"
	print, "	that crosstalk to the data.  Crosstalk determined and"
	print, "	applied on a row-by-row basis."
	print
	print, "	Arguments"
	print, "		stks	- spectral images (all 4 are input"
	print, "			  but only qa,ua,va are modified)"
	print, "			  first index is along slit, second"
	print, "			  is wavelength"
	print, "	sr1,sr2		- range of pixels along the slit for "
	print, "			  fitting the crosstalk (avoid ends)"
	print
	print, "	Keywords"
	print, "		(none)"
	print
	return
endif
;-


;
;	Get dimensions of arrays.
;

nw = sizeof_sbsp(stks, 1)
ny = sizeof_sbsp(stks, 2)
nw1 = nw-1
ny1 = ny-1

;	Get median of Stokes vector in continuum window.
;
pmed = fltarr(4,ny)		; zero array

;  get average of spectrum along the slit
spect = fltarr(nw)
for jj = sr1,sr2 do spect = spect + stks(*,jj,0)

;  get median of spectrum, find continuum wavelengths
spect = spect/float(sr2-sr1+1)
mdn = median(spect)
whr = where(spect gt mdn)
for kk = 0,2 do begin
	for jj = 0, ny1 do pmed(kk,jj) =  $
	    mean_sbsp(stks(whr,jj,kk))
endfor

;	Calculate crosstalk levels.
for kk = 1, 2 do  pmed(kk,*) = pmed(kk,*)/pmed(0,*)

;	Subtract crosstalk correction from the polarization spectra
for kk = 1,2 do begin
	for jj = 0, ny1 do begin
;  check for very low intensities (off limb), where no correction will
;  be applied because of the possibility of intensity values very close to
;  zero
		if pmed(0,jj) gt 100. then begin
			stks(*,jj,kk) = stks(*,jj,kk) $
				- pmed(kk,jj)*stks(*,jj,0)
		endif
	endfor
endfor


return
end