function make_fake_dark, seq_filnam
;+
;NAME:
;	make_fake_dark
;PURPOSE:
;	Calculate a simple fake dark image from the corners of a FD image
;HISTORY:
;	Written 15-Nov-94 by I. Zayer
;CALLING SEQUENCE:
;	dark=make_fake_dark(seq_filnam)
;INPUT:
;	seq_filnam: file name of the datafile (string)
;OUTPUT:
;	dark: 2-D array of same dimensions as 1st image in seq_filnam
;	      which mimics a dark frame (with gradient calculated from
;	      the corners of the above image assumed to be a FD frame)
;-
;

read_mdi,seq_filnam,0,img & s=size(img) & xmax=s(1)-1 & ymax=s(2)-1

  upperave=total(img(2:21,ymax-21:ymax)/800.)+$
  	   total(img(xmax-21:xmax,ymax-21:ymax)/800.)
  lowerave=total(img(2:21,2:21)/800.)+$
	   total(img(xmax-21:xmax,2:21)/800.)
  dcol=fix(indgen(ymax+1)*(upperave-lowerave)/float(ymax+1) + lowerave)
  for i=0,xmax do img(i,*)=dcol

return,img
end

;*****************************************************************************

pro do_fit_gcross, seq_filnam, images=images, idark=idark, $
	qhr=qhr, qmove_around=qmove_around, nfit=nfit, $
	poutfil=poutfil, spec_code=spec_code, $
	qdisp=qdisp, qplot=qplot, qpause=qpause, qstop=qstop
;
;
; qdisp=1	& qplot=1	& qpause=1
; qdisp=1	& qplot=1	& qpause=0
; qdisp=0	& qplot=0	& qpause=0
;
;+
;NAME:
;	do_fit_gcross
;PURPOSE:
;	Given a series of images, fit a gaussian curve to all of the
;	grid locations.
;HISTORY:
;	Written Feb-94 by M.Morrison
;V2.1	26-Apr-94 (MDM) - Added NFIT keyword option
;V2.2	27-Apr-94 (MDM) - Changed full disk grid separation from 135 to 130 pixels
;			  and  NFIT from 10 to 18.
;			- Change to FIT_GCROSS to make a better first guess.
;			  Resulted in many more good fits
;V2.3	 4-Nov-94 (IZ) - Faked dark gradient for 941012_190252 data...
;V2.4	15-Nov-94 (IZ) - Fake dark calculated from first image (assumed FD)
;-
;
progver = 'DO_FIT_GCROSS  Ver 2.4
;
if (n_elements(spec_code) ne 0) then case spec_code of
    0: begin
	seq_filnam = '/umbradata2/optest/pdist616.'		; 81 images
	images = indgen(40)*2
	idark = 80
	qmove_around = 0
	qhr = 0
       end
    1: begin
	seq_filnam = '/umbradata2/optest/distgrid617.'		;210 images
	images = indgen(104)*2+1
	idark = 209
	qmove_around = 1
	qhr = 0
       end
    2: begin
	seq_filnam = '/umbradata3/optest/940222_190446'		;for Chris
	images = indgen(22)+1
	idark = -1
	qmove_around = 1
	qhr = 1
       end
    3: begin
	seq_filnam = '/umbradata3/optest/940222_203749'		;for Chris
	images = indgen(22)+1
	idark = -1
	qmove_around = 1
	qhr = 1
       end
    4: begin
	seq_filnam = '/umbradata5/optest/log/log/sci160k/940120_134308'
	images = [indgen(40) + 1, indgen(40)+42]
	idark = 82
	qmove_around = 1
	qhr = 1
       end
    5: begin
	seq_filnam = '/umbradata5/optest/log/log/sci160k/940120_141706'
	images = indgen(200) + 1
	idark = 201
	qmove_around = 1
	qhr = 1
       end
endcase
;
read_mdi, seq_filnam, images, code, txtfil, /type_info
outfil = concat_dir(getenv('MDI_CAL_TXT'), txtfil+'.gcross')
if (keyword_set(poutfil)) then outfil = outfil + poutfil	;append special qualifiers (want two files from one series)
;
nimg = n_elements(images)
;
if (n_elements(idark) ne 0) then if (idark ne -1) then $
  read_mdi,seq_filnam,idark,dark,header $
else dark=make_fake_dark(seq_filnam)
;
if (n_elements(lun) ne 0) then free_lun, lun
openw, lun, outfil, /get_lun
printf, lun, progver
printf, lun, 'Program Run: ' + !stime
printf, lun, ' '
printf, lun, 'File Base: ', seq_filnam
printf, lun, ' '
printf, lun, 'Image   Init Guess  Gauss Fit Loc           X-Fit              Fit-Dev                 Y-Fit
printf, lun, ' #       X    Y     X        Y       Amp       Width   Offset   (y-fit)  Amp       Width   Offset   (y-fit)
             ;0000    158  236  160.373  236.251 -756.282    3.958  855.141   28.696 -737.560    4.179  899.344   15.817
printf, lun, ' '
for iimg=0,nimg-1 do begin
    images0 = images(iimg)
    read_mdi, seq_filnam, images0, img, header
;    if (n_elements(idark) ne 0) then if (idark ne -1) then img = img - dark
    if (n_elements(idark) ne 0) then img = img - dark
    ;
    nx = n_elements(img(*,0))
    ny = n_elements(img(0,*))
    if (not keyword_set(qhr)) then begin	;full disk case
	n = 200
	if (keyword_set(qmove_around)) then begin
	    dummy = max(img(20:nx-20,*), imax)
	    x0 = ((imax mod (nx-39))+20) < (nx-n)
	    y0 = (imax/(nx-39)) < (ny-n)
	end else begin
	    x0 = 400
	    y0 = 400
	end
	subpart = img(x0:(x0+n-1)<(nx-15), y0:y0+n-1)
	xarr = total(subpart,2)/n
	yarr = total(subpart,1)/n
	;
	xmin = min(xarr, ixmin)
	ymin = min(yarr, iymin)
	avg = total(subpart)/n_elements(subpart)
	;
	approx_x = x0+ixmin
	approx_y = y0+iymin

	;steps = 135
	steps = 130
	if (n_elements(nfit) eq 0) then nfit = 18
    end else begin
	xarr = total(img,2)/n_elements(img(0,*))
	yarr = total(img,1)/n_elements(img(*,0))

	n = fix(0.90*nx)
	x0 = fix(0.05*nx)	;drop first few lines/columns
	y0 = fix(0.05*ny)
	xmin = min(xarr(x0:x0+n-1), ixmin)
	ymin = min(yarr(x0:x0+n-1), iymin)
	avg = total(xarr(x0:x0+n-1) + yarr(x0:x0+n-1))/ (n*2)
	;
	approx_x = x0+ixmin
	approx_y = y0+iymin

	steps = 440
	if (n_elements(nfit) eq 0) then nfit = 20
    end
    ;
    if (keyword_set(qdisp)) then begin
	mag = 1
	if (nx gt 1000) then mag = 0.5
	if (nx le 256) then mag = 2.0
	erase
	img2 = rebin(img, nx*mag, ny*mag)
	tvscl, img2
	for i=approx_x mod steps,nx, steps do for j=approx_y mod steps,nx,steps do plots, i*mag, j*mag, psym=2, /dev
	if (keyword_set(qpause)) then pause
    end
    ;
    for ypos=approx_y mod steps,nx,steps do for xpos=approx_x mod steps,ny, steps do begin
	fit_gcross, img, nfit, avg, xpos, ypos, 0, coeff0, dev0, qplot=qplot, qpause=qpause
	fit_gcross, img, nfit, avg, xpos, ypos, 1, coeff1, dev1, qplot=qplot, qpause=qpause
	ss = [0,2,3]
	fmt = '(i4.4, 2x, 2i5, 2f9.3, 10f9.3)'
	str = string(images0, xpos, ypos, coeff0(1), coeff1(1), $
		coeff0(ss), dev0, coeff1(ss), dev1, format=fmt)
	print, str
	printf, lun, str
	flush, lun
    end
end
;
if (keyword_set(qstop)) then stop
free_lun, lun
end