pro ltf, filnam1, img1, filnam2, img2, filnamdc, imgdc, label, $
		imin, imax, qplot, method, $
		xst, xen0, yst, yen0, nx, ny, xstep, ystep, $
		areadnoise2, aavg_off, max_gainc, $
		gainc_lin, gainc_log, gainc_hist, $
		avg_off_lin, readnoise2_lin, b_lin, fwhm, $
		bls_noise, avg_off_img, $
		err_linfit, err_logfit, err_histfit, maxoff, avgoff
;
;							V1.0  MDM 21-Feb-89
;							V1.01 MDM 22-Feb-89
;							V2.00 MDM 15-Mar-89
;							V2.01 MDM 16-Mar-89
;							V3.00 MDM 17-Mar-89
;							V3.01 MDM 20-Mar-89
;							V3.02 MDM  5-May-89
;							V3.03 MDM  2-Jan-90
;							V3.04 MDM  2-Jan-90
;							V4.00 MDM  8-Jan-90
strver = 					'MDM LTF Ver 4.0 (8-Jan-90)
;
;LTF Ver 4.0	- Used a dark current image to adjust for no DC offset in
;		  signal (aavg_off should be zero if DC image is used)
;		  Changed plotting information and format slightly
;		  Log fit does not subtract "areadnoise2"
;LTF Ver 3.04	- Modified "gainc_hist_arr" to discard the last bin
;		  because "<max_gainc" causes all values greater to go into that
;		  bin
;LTF Ver 3.03	- Changed "areadnoise2" from being overwritten (to use the
;		  "BLS_noise" calculated value) to the value passed by the
;		  calling routine
;		- "img" parameter is "strver" if filnam1 = '****'
;		- Put a 4 element box car smooth on gain histogram vector
;		  for fitting (found spikes elsewise)
;LTF Ver 3.02	- Modified "gainc_hist_arr" to only use the "avg" and
;		  "var" values that are within the specified region (imax)
;
;
;input:	filnam1 - file name 1 (used by header)
;	img1	- image array 1
;	filnam2 - file name 2 (used by header)
;	img2	- image array 2
;	imax	- image max intensity to fit line for gain to
;	qplot	- what to plot 	bit 0 (value=1) - Fit to data
;				bit 1 (value=2) - Histogram
;				bit 2 (value=4) - Log-Log fit
;				bit 7 (value=128) - Pause between plots
;	method	- 
;output-gainc	- gain constant of data given
;
if (filnam1 eq '****') then begin
    img1 = strver
    return
end
;
;----- initialize certain variables
;
siz=size(img1)
nximg = siz(1)
nyimg = siz(2)
;
if (xen0 lt 0) then xen = nximg - abs(xen0) - 1 else xen=xen0
if (yen0 lt 0) then yen = nyimg - abs(yen0) - 1 else yen=yen0
;
nx_step = fix ((xen - xst -nx) / xstep) + 1
ny_step = fix ((yen - yst -ny) / ystep) + 1
ngainc = nx_step*ny_step 
;
; ----- get noise and variance data
;
ltf_get_avgvar, img1, img2, imgdc, $
	xst, xen, yst, yen, nx, ny, xstep, ystep, ngainc, $
	avg, var, maxoff, avgoff
;
x1 = 5.
x2 = 14.
y1 = 25.
y2 = nyimg-25.
n = (x2-x1+1)*(y2-y1+1)
bls_noise = total( (long(img1(x1:x2,y1:y2)) - long(img2(x1:x2,y1:y2)))^2/2 )/n
;
x1 = 35.
x2 = nximg-35
y1 = 5
y2 = 14
n = (x2-x1+1)*(y2-y1+1)
avg_off_img = total( long(img1(x1:x2,y1:y2)) + long(img2(x1:x2,y1:y2)) )/2.0/n
;
; ----- Do the gain calculations for a line fit
;
;	(photon_noise^2 + read_noise^2) = gain*(photon + DC_offset)
;			y		=    m*x
;
;Assume signal=0 ==> photon_noise^2 = 0
;		read_noise^2 = gain*DC_offset
;
if (imax eq 0.) then imax = 3500.
valid = where ((avg lt imax) and (avg gt imin))

favg = avg(valid)
fvar = var(valid)
ss=sort(favg)
favg=favg(ss)
fvar=fvar(ss)

lin_coeff = poly_fit(favg, fvar, 1, lin_fit, yband, sigma, a)	; linear fit
gainc_lin = 1./lin_coeff(1)
b_lin = lin_coeff(0)	;should be zero (need to change poly_fit to chi^2??)
;
;;p=favg			;minimize chi squared??		predicted
;;o=fvar			;				observed
;;gain = total( p ) / total( p^2/o )
;;gainc_lin = 1/gain
;;lin_fit = gain*p
;
;		read_noise^2 = gain*DC_offset	;"B" "fix" is not understood?
avg_off_lin = (areadnoise2 - lin_coeff(0))/lin_coeff(1)
readnoise2_lin = lin_coeff(0) - lin_coeff(1)*aavg_off 
;need to shift y-intercept over to correct for DC offset
err_linfit = total(abs(lin_fit-fvar))/n_elements(fvar)	;normalize for # points
;
; ----- Get the gain for each pair of data points
;
;areadnoise2 = readnoise2_lin	;!!! TODO REmove??
;areadnoise2 = bls_noise	;!!! TODO REmove??
;gainc_hist_arr = (avg-aavg_off) / (var - areadnoise2)	;taken out 5-May-89
gainc_hist_arr = (favg-aavg_off) / (fvar - areadnoise2)
;
xmin = 0.
xmax = max_gainc
;;xmax = max(gainc_hist_arr)
;;xmax = xmax < max_gainc	;because 4x4 has some wierdness
fact = 100./(xmax-xmin) 	;stretch out or compress gainc range (integerizing)
h = histogram((gainc_hist_arr<max_gainc)*fact)
n = n_elements(h)
h = h(0:n-2)			;discard last bin
xx = findgen(n_elements(h))/fact
;;dummy = max(h(25:*))
;;peak = (!c+25)/fact
peak = 0.

if (n_elements(xx) lt 10) then begin
    f_xx = xx
    f_h  = h
end else begin
    f_xx = xx(5:*)		;avoid "spike" at zero
    temp = smooth(h, 4)		;added 2-Jan-90
    f_h  = temp(5:*)
end
hist_coeff = fltarr(4)
hist_fit = gaussfit(f_xx,f_h,hist_coeff)
gainc_hist = hist_coeff(1)
fwhm = hist_coeff(2)
;
err_histfit = total(abs(hist_fit-f_h))/n_elements(f_h)	;normalize for # points
;
; ----- Do the gain calculations for a log-log line fit
;
if (imax eq 0.) then imax = 3500.
valid = where ((avg lt imax) and (avg gt imin))

lavg = alog10(avg - aavg_off)		;need to subtract assumed DC offset
;lvar = alog10(var - areadnoise2)	;need to subtract readnoise
lvar = alog10(var)

flavg = lavg(valid)
flvar = lvar(valid)
ss=sort(flavg)
flavg=flavg(ss)
flvar=flvar(ss)

;!!!!! slope might not be one (which it should be)!!!!!
;;log_coeff = poly_fit(flavg, flvar, 1, log_fit, yband, sigma, a)	; linear fit
;;log_gainc = 10.0 ^ ( -log_coeff(0)/log_coeff(1) )

;	var = gain * DN
;	log(var) = log(gain) + log(DN)
;	log(var) = -log(gain_constant) + log(DN)	gain = 1/gain_constant
;	  y      =     b  +                x 
;							log(gain) = -log(gain_cont)
;
p=flavg				;minimize chi squared??		predicted
o=flvar				;				observed
lgain = - total( (p-o)/o ) / total( 1/o )
gainc_log = 10.0 ^ ( -lgain )
log_fit = p + lgain
;
err_logfit = total(abs(log_fit-flvar))/n_elements(flvar) ;normalize for # points
;
flavg   = [-lgain, flavg]
log_fit = [0, log_fit]
;
; ----- now plot up results
;
bits, qplot, b
in=' '
;
if (b(0) eq 1) then begin			;Stern LTF plot with fit
    pmin = 0
    pmax = 0		;display everything, overlay plot of fit
    qplot0 = 1
    ltf_plot, avg, var, favg, lin_fit, pmin, pmax, qplot0, gainc_lin
    ltf_plotinfo, filnam1, filnam2, filnamdc, $
	imin, imax, label, xst, xen, yst, yen, $
	nx, ny, xstep, ystep, maxoff, avgoff, strver
    if (b(7) eq 1) then read,'Pause',in
end

;if ((b(0) eq 1) and (!cymax gt 3*imax)) then begin	;Stern LTF plot with fit
;    pmin = 0
;    pmax = imax
;    qplot0 = 1
;    ltf_plot, avg, var, favg, lin_fit, pmin, pmax, qplot0, gainc_lin
;    ltf_plotinfo, filnam1, filnam2, filnamdc, $
;	imin, imax, label, xst, xen, yst, yen, $
;	nx, ny, xstep, ystep, maxoff, avgoff, strver
;    if (b(7) eq 1) then read,'Pause',in
;end

clearplot

if (b(1) eq 1) then begin			;gain histogram plot
    ltf_plothist, avg, var, xx, h, f_xx, hist_fit, peak, gainc_hist, fwhm
    ltf_plotinfo, filnam1, filnam2, filnamdc, $
	imin, imax, label, xst, xen, yst, yen, $
	nx, ny, xstep, ystep, maxoff, avgoff, strver
    ltf_plotinfo2, gainc_hist, areadnoise2, aavg_off
    if (b(7) eq 1) then read,'Pause',in
end
;
clearplot
;
if (b(2) eq 1) then begin			;JPL (log) LTF plot
    pmin = 1
    pmax = 10000.		;display everything, overlay plot of fit
    qplot0 = 1
    ltf_plotlog, 10.^lavg, 10.^lvar, 10.^flavg, 10.^log_fit, pmin, pmax, $
		qplot0, gainc_log
    ltf_plotinfo, filnam1, filnam2, filnamdc, $
	imin, imax, label, xst, xen, yst, yen, $
	nx, ny, xstep, ystep, maxoff, avgoff, strver
    ;ltf_plotinfo2, gainc_log, areadnoise2, aavg_off
    if (b(7) eq 1) then read,'Pause',in
end

clearplot
;
return
END