pro genflt3_sbsp, savnam, fltav, singleframe,shftavg,delshavg, $
	svrr,cc1,cc2,sp1,sp2,sl1,sl2, $
	wc1,wc2,sr1,sr2,nxtnd,nsearch,noplot=noplot, verbose=verbose

;==============================================================================
;+
;
;   procedure: genflt2_sbsp 
;
;   purpose:  Generate gaintable for HINODE SP data.
;
;   history:  modified from the DLSP routine developed by Bruce Lites
;		This routine now accepts the dark-corrected average flat images
;		output from flatgen_sbsp.pro.

;  inputs:  savenam	= string name of file for saving results
;  	    fltav	= average flat images that have been corrected for intensity
;			  variations along the slit and thermal drift
;  	    singleframe	= path to fits file containing one single observation with
;			  the structure along the slit for getting the spectral skew
;		shftavg  = operation-averaged spectral skew variation
;		delshavg = operation-averaged displacement of spectrum between CCDSIDEs
;			This is now used instead of the single frame results
;	    svrr	= array of averaged slit intensity variations for each CCDSIDE
;			  (this array just passed through to the output)
;	    cc1,cc2	= pixel range in data units for averaging continuum intensity
;	    sp1,sp2	= active spectral pixel range for each CCDSIDE, data units
;	    sl1,sl2	= active pixel range along slit for each CCDSIDE, CCD units
;	    wc1,wc2	= spectral pixel range for cross-correlation to get 
;				spectral shifts, data units
;	    sr1,sr2	= pixel range along slit for determination of spectral
;				curvature, CCD units
;	    nxtnd	= number of spectral pixels to extend spectral range
;	    nsearch	= correlation search range
;
;  outputs: mflat	= flat field correction images for each CCDSIDE
;	    shftsl	= skew correction, rel. shift along slit for each CCDSIDE
;	    specshft	= fit of spectral shift along slit between CCDSIDEs
;	    specrfit	= fit to spectral curvature
;  outputs in both IDL save file and SSW .geny file
;
;
;
;==============================================================================
;
;	Check number of parameters.
;
if n_params() ne 18 then begin
	print
	print, "usage:  genflt3_sbsp, savnam,avg_flats,singleframe,svrr "
	print
	print, "Generate gaintable for the Solar-B SP data "
	print
	print, "	Keywords"
	print, "	    noplot  - if set,turn off gaintable plotting	"
	print, "					  (def=do plotting)					"
	print, "	    verbose	- flag to print run-time info			"
	print, "		              0 = minimal printing (def)		"
	print, "		              1 = print everything (/verbose)	"
	print, "		              2 = print all but headers			"
	print
	print
	return
endif
;-
;
	do_plot   = 1 - keyword_set(noplot)
	do_verb   = keyword_set(verbose)

;  get the string descriptor from the file name, use this in the name of the
;  output file
;                spos = strpos(savnam,'_gn')
; 		epos = strpos(savnam,'_20')
;		numm = epos - spos
;		savnam = strmid(avg_darks,spos,numm)
		print,'savnam = ',savnam

;
;	-----------------------------------------
;	Get averaged dark corrected flat spectra.
;	-----------------------------------------

;  routine to generate slit-averaged spectral profiles, then divide
;  into the flat-field spectrum to get a "line-free" image.
;  This process will probably also remove spectral fringes, which
;  will have to be removed in a later process.

;  get sizes of raw array
	nnw = sizeof_sbsp(fltav,1) & nnslt = sizeof_sbsp(fltav,2) 

;  establish multiplicative flat field array for output
	mflat = fltav

;  start loop over two CCDSIDE images
	for kk = 0,1 do begin
		if kk eq 0 then img = fltav(*,*,0) else img = fltav(*,*,1)


;  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)


;  find average profile near center of slit
		avprf = fltarr(nx)
		for jj = ny2-20,ny2+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.
		xxx = findgen(ny)
;	dont fit the last 10 points, avoiding the blank top row
		coef = poly_fit(xxx(0:ny-10), shft(0:ny-10), 5)
		shfit = poly(xxx, coef)

;  now find slit-averaged profile over more of slit and do
;	cross-correlation again
		avprf(*) = 0.
		yst = 50 & ynd = ny-50
		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.
		shft = shft - mean_sbsp(shft)
;
;       Fit the shifts with a fourth-order polynomial.
		xxx = findgen(ny)
;	dont fit the last 10 points, avoiding the blank top row
		coef = poly_fit(xxx(0:ny-10), shft(0:ny-10), 5)
		shfit = poly(xxx, coef)







;  the following insert required because spectral response
;  seems to vary along the slit, or there are residual shifts after
;  averaging
;  find image of average profile over 100 slit pixels, running mean 
;  along slit
apimg = img	; define image of running mean slit average profiles
avprg = avprf
for jj = 0,ny1 do begin
	avprg(*) = 0.
	yst = max([0,jj-50])
	ynd = yst + 99
	if ynd gt ny1 then begin
		ynd = ny1
		yst = ny1 - 99
	endif
	for jjk = yst,ynd do begin
		avprg = avprg + fshft_sbsp(img(*,jjk),-shfit(jjk),0,nx1, nxtnd)
	endfor
	apimg(*,jj) = avprg/float(ynd-yst+1)
endfor
;  remove the average spectral slope of these running-mean profiles
;  remove the linear trend of the average profile, since we want
;  the flat field correction applied to data to compensate for this
;  linear trend.  Fit all points above the median value with a 
;  linear function.  Median value isolates continuum pretty well.
;  get the slit averaged profile without shifting
	avprg(*) = 0.
	for jj = 0,ny1 do begin
		avprg(*) = avprg(*) + apimg(*,jj)
	endfor
	sp1g = sp1(kk) & sp2g = sp2(kk)
	avprgg = avprg(sp1g:sp2g)
	whr = where(avprgg gt median(avprgg))
	xxx = findgen(n_elements(avprg))
	xx2 = xxx(whr) + sp1g
	yyy = avprg(whr)
	clin = linfit(xx2,yyy)
	linft = clin(0) + clin(1)*xxx
	linft = linft(0)/linft
;  correct the average profile by the linear fit, then normalize
	avprf = avprf*linft
	avprf = avprf/median(avprf)

;  remove the spectral slope from the image of average profile
	for jj=0,ny1 do apimg(*,jj) = linft(*)*apimg(*,jj)
;  find average intensity variation along slit, then remove that
        sltv = fltarr(ny)
        for jj = 0,ny1 do sltv(jj) = mean_sbsp(apimg(*,jj))
	sltv = sltv/mean_sbsp(sltv)
	for jj = 0,ny1 do apimg(*,jj) = apimg(*,jj)/sltv(jj)
;  restore shifts to average profiles
	for jj = 0,ny1 do begin
		tmp = fshft_sbsp(apimg(*,jj),shfit(jj),0,nx1,nxtnd)
		apimg(*,jj) = tmp(*)
	endfor
;  normalize by median of profile image, then divide flat image
;  by this average shifted profile
	apimg = apimg/median(apimg)
	clrsp = img/apimg
;  end of modification inspired by lores






;       Generate flat-field image
		if kk eq 0 then begin
			clrsp1 = clrsp
			avclr1 = mean_sbsp(clrsp)
			fltimg1 = avclr1/clrsp
			avprf1 = avprf
			shft1 = shft
			shfit1 = shfit
			mflat(*,*,kk) = fltimg1
		endif else begin
			clrsp2 = clrsp
			avclr2 = mean_sbsp(clrsp)
; preserve the overall relative scaling of the two sides here
; so that gets corrected in the final flat-field
			;fltimg2 = avclr2/clrsp
			fltimg2 = avclr1/clrsp
			avprf2 = avprf
			shft2 = shft
			shfit2 = shfit
			mflat(*,*,kk) = fltimg2
		endelse

;  top row of image has low intensity, set flat field to unity
	mflat(*,ny1,*) = 1.

;  end loop over CCDSIDEs
	endfor


;  use raw gain image to get the shifts in the direction along the slit 
;  as a function of spectral pixel: shftsl
;  Use a single frame with lots of structure along the slit for this
;  determination.  Expectation is that these shifts will not change in time
	slitshftnp_sbsp,singleframe,cc1,cc2,sp1,sp2,sl1,sl2, $
		ssum,ssl1,ssl2,nsearch,shftsl,delsh

;  recover the operation-averaged shftsl, delsh instead of using the
;  single-frame values
shftsl = shftavg
delsh = delshavg

dumm = ' '
print,' spatial shifts along the slit for two CCDSIDES vs spectrum'
;read,' enter any character to continue: ',dumm

;  set range along slit for both curvature and spectral shift
;  according to summing along slit
	ssr1 = sr1/ssum
	ssr2 = sr2/ssum

;  use gain-corrected flat image to get the spectral shifts between
;  the two CCDSIDEs
	dat = fltav*mflat
	;specshft_sbsp,dat,ssl1,ssl2,wc1,wc2,specshft
	specshft_sbsp,dat,ssr1,ssr2,wc1,wc2,specshft

print,' spectral shifts between sides vs. distance along the slit'
;read,' enter any character to continue: ',dumm

;  use gain-corrected flat image to get the spectral curvature
;;  set range for curvature according to summing along slit
;	ssr1 = sr1/ssum
;	ssr2 = sr2/ssum
	specross_sbsp,dat,ssr1,ssr2,nxtnd,specurv,specrfit,avprf

print,' spectral curvature vs length along the slit, CCDSIDE0'
;read,' enter any character to continue: ',dumm

;  set gain to unity for blank row (at top in reformatted images)
;  at top of images, set the gain in the vignetted area to that of
;  last useable row.  Because of offset between the two CCDSIDEs,
;  CCDSIDE1 starts 3 pixels before CCDSIDE0
	for jj = ssl2-1,ny1 do mflat(*,jj,0) =  $
		mean_sbsp(mflat(*,ssl2-2,0))
	for jj = ssl2-4,ny1 do mflat(*,jj,1) =  $
		mean_sbsp(mflat(*,ssl2-5,1))
;  Also correct the last lines of svrr
	for iside = 0,1 do svrr(ssl2+1:*,iside) = svrr(ssl2,iside)

;  Display the final flat fielded image
	mflat2 = bytarr(nnw*2,nnslt)
	medn = median(mflat(*,*,0))
	mflat2(0:nnw-1,*) = bytscl((mflat(*,*,0)>0.8*medn)<1.2*medn)
	medn = median(mflat(*,*,1))
	mflat2(nnw:*,*) = bytscl((mflat(*,*,1)>0.8*medn)<1.2*medn)
	window,xs=224,ys=1024,/free
	tvscl, mflat2

;  save data relevant to flat fielding, merging, etc.

;	Set output savefile name including the path.
	;strsum = strcompress(string(ssum),/remove_all)
	;savefile = strcompress('gaintable' + savnam + $
	;	'_sum' + strsum + '.save',/remove_all)
	savefile = strcompress('gaintable' + savnam + '.save',/remove_all)
	print,'The flat results will be saved in    ',savefile

	save,filename=savefile,mflat,shftsl, $
		delsh,specshft,specrfit,avprf1,shft1,shfit1,clrsp1, $
		avprf2,shft2,shfit2,clrsp2,fltav,svrr

;  save the needed data in ssw geny format
	savegenx,file='gaintable'+savnam+'.geny',mflat,shftsl,specrfit,specshft,svrr

	return
end
