pro terp_xmat_sbsp,slpos,calpos,skx,istrt,iend,xinv

;  routine to read interpolate the Solar_B SP calibration X-matrix for a given
;  slit position, CCDSIDE, Spectral ROI.  This routine delivers a calibration
;  matrix for the entire slit length.  If there is any binning or ROI selection
;  along the slit, the output will need to be adjusted accordingly externally

;  input:	slpos	= the Solar_B SP slit position, scan step number (-1024
;                         to 1024)
;		calpos	= array of Solar_B SP slit positions for which fitted
;			  parameters exist
;		skx	= array (4,4,4,4,2,9) of fit coefficients for
;			  inverse X-matrix
;		istrt	= the spectral ROI starting index (measured from 0)
;		iend	= the spectral ROI ending index 

;  output:	xinv	= the Solar_B SP inverse X-matrix array expanded to the
;			  original data size, for each CCDSIDE

calparams_file=file_search(get_logenv('$SOT_SP_CALIBRATION'),'calparams_sbsp*.com')
if not file_exist(calparams_file) then begin 
   box_message,'Cannot find $SOT_SP_CALIBRATION/calparams_sbsp*'
   return
endif 

pedata=strnocomment(strarrcompress(rd_tfile(last_nelem(calparams_file))),com=';')
pstat=1
for i=0,n_elements(pedata)-1 do begin 
   estat=execute(pedata(i))
   pstat=pstat and estat
endfor

;  define output array
nwl = iend-istrt+1
xinv = fltarr(nwl,nyraw,4,4,2)


nx = xdim & ny = ydim
smth = fltarr(nx,ny)		; smoothed output data temporary array

;  array of smoothed output data expanded to original coordinates
smthbig = fltarr(nx*nxstr,ny*nystr)

;  number of calibration slit positions
npos = n_elements(calpos)


;  do a linear interpolation to get fit coefficients between 
;  discrete cal observation positions
if slpos le calpos(0) then begin
	k1 = 0 
endif else if slpos gt calpos(npos-1) then begin
	k1 = npos-2 
endif else begin
	for kk = 0,npos-2 do begin
	if slpos gt calpos(kk) and slpos le calpos(kk+1) then k1 = kk
	endfor
endelse
k2 = k1 + 1
skk = ((calpos(k2)-slpos)*skx(*,*,*,*,*,k1) +  $
	(slpos-calpos(k1))*skx(*,*,*,*,*,k2))/(calpos(k2)-calpos(k1))


;  now recover smoothed variations from these interpolated coefficients
xpx = dblarr(nx,ny) & ypx = dblarr(nx,ny)

;  begin restoration loop over CCDSIDEs
for iside = 0,1 do begin
;  begin loop over matrix elements
	for kk = 0,3 do for ll = 0,3 do begin
		smth(*,*) = 0.
		xpx(*,*) = findgen(nx) # replicate(1., ny) ;X values at each point
		ypx(*,*) = replicate(1.,nx) # findgen(ny)
;  begin loop over coefficients
		for mj = 0,ndeg do for mi = 0,ndeg do begin
			smth = smth +  skk(mj,mi,ll,kk,iside) * xpx^mi * ypx^mj
		endfor

;  rebin the smoothed fits back up by a factor of nxstr (4) in x, nystr (16) in y
		smthbig(*,*) = rebin(smth,nx*nxstr,ny*nystr)
		xinv(*,*,ll,kk,iside) = smthbig(istrt:iend,*)

;  end loop over matrix elements
	endfor

;  end loop over CCDSIDEs
endfor


return
end