pro specross_sbsp,dat,sr1,sr2,nxtnd,specurv,specrfit,avprf

;  routine to get the spectral curvature of CCDSIDE0 by
;  cross-correlation of the slit-averaged spectral profile.  Use
;  the gain-corrected flat field image.

;  INPUTS:	dat		= gain-corrected flat field Stokes I
;		sr1,sr2		= pixel range along slit for fitting curvature,
;					data units
;		nxtnd		= number of spectral pixels to extend spectral range

;  OUTPUTS:	specurv		= measured spectral shifts along the slit, pixels
;		specrfit 	= polynomial fits to spectral shifts
;		avprf		= new average spectral profile


;;  include hardwired parameters sr1, sr2, wl1, wl2
;@params_sbsp.com

;  get gain-corrected CCDSIDE0 image
img = dat(*,*,0)

;  dimensions of image
nx = sizeof_sbsp(img,1) & ny = sizeof_sbsp(img,2)
nx1 = nx-1 & ny1 = ny-1 & ny2 = ny/2
shft = fltarr(ny)

;  start image 1, find average profile near center of slit
avprf = fltarr(nx)
mid_y = ny/2
for jj = mid_y-20,mid_y+20 do avprf = avprf + img(*,jj)
avprf = avprf/41.

; loop through slit length finding cross-correlation of the lines
for jj = 0, ny1 do begin
        tempr = img(*, jj)
        shft(jj) = corshft_sbsp(avprf,tempr)
endfor

;       Subtract out the average to get relative shifts.
shft = shft - mean_sbsp(shft)

;       Fit the shifts with a fourth-order polynomial.
x1 = findgen(ny)
coef = poly_fit(x1(sr1:sr2), shft(sr1:sr2), 5)
shfit = poly(x1, coef)


;  now find slit-averaged profile over more of slit and do cross-correlation again
avprf(*) = 0.
yst = mid_y-ny2+25 & ynd = mid_y+ny2-25
for jj = yst,ynd do begin
	avprf = avprf + fshft_sbsp(img(*,jj),-shfit(jj),0,nx1,nxtnd)
endfor
avprf = avprf/float(ynd-yst+1)

; loop through slit length finding cross-correlation of the lines
for jj = 0, ny1 do begin
        tempr = img(*, jj)
        shft(jj) = corshft_sbsp(avprf,tempr)
endfor

;       Subtract out the average to get relative shifts.
specurv = fltarr(ny)
specurv(sr1:sr2) = shft(sr1:sr2) - mean_sbsp(shft(sr1:sr2))

;       Fit the shifts with a fourth-order polynomial.
x1 = findgen(ny)
coef = poly_fit(x1(sr1:sr2), specurv(sr1:sr2), 5)
specrfit = poly(x1, coef)
window,xs=700,ys=500,/free
plot,x1,specrfit,title='spectral curvature vs slit length'
oplot,x1(sr1:sr2), specurv(sr1:sr2),psym=3


return

end