pro even_light_corrections, hhvv, spot_centroid, sym_spot_nr, good_spot, bins_axis,	$
	bins_step, img_nr, center0, dcenter, optaxis0, doptaxis, binwidth, dbinwidth,	$
	y_run, bmp=bmp, paper=paper, outdir=outdir, stopq=stopq

;+
; NAME:
;	even_light_corrections
; PURPOSE:
; CATEGORY:
; INPUTS:
;	hhvv				identifies frame sequence
;	spot_centroid		x-y coordinates of spot centroids (relative to center0)
;	good_spot			list of spot indices to be used for obtaining corrections
;	bins_axis
;	bins_step
;	img_nr
;	center0				center of fov arc (pix)
; OUTPUTS:
;	center0
;	dcenter
;	optaxis0
;	doptaxis
;	binwidth
;	dbinwidth
;	y_run
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; MODIFICATION HISTORY:
;	???-????, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, bmp  , /key
InitVar, paper, /key
InitVar, stopq, /key

; Note that the first image is not used so image 0 is file with number 01
; Correction for center for 'neutral position' (images 0....4)

rr = spot_centroid[*,good_spot,*]

xycenter = lsqCircleFit(rr)				; [3,*]: x, y, R

prcenter = xycenter[1:2,*]				; [2,*]: y, R (y is overwritten by azimuth later)
xycenter = xycenter[0:1,*]				; [2,*]: x, y

nspot = n_elements(good_spot)
nimg  = n_elements(img_nr)

rr = rr-SuperArray(xycenter, nspot, after=1)

sz = size(rr)

rr = reform(rr, sz[1], sz[sz[0]+2]/sz[1])
rr = cv_coord(from_rect=rr, /to_polar)
rr = reform(rr, sz[1:sz[0]], /overwrite)

rr = reform(rr[0,*,*])
rr = lsqQuadraticFit(sym_spot_nr[good_spot], rr)

prcenter[0,*] = rr[0,*]					; Azimuth for optical axis
binw  = reform(rr[1,*])
dbinw = reform(rr[2,*])

; Decide whether this is a horizontal or vertical run

vbeg = total(prcenter[*, 0: 4],2)/5		; Start values for azimuth, radius of optical axis
vend = total(prcenter[*,44:46],2)/3 	; End values for azimuth, radius of optical axis

doptaxis = (vend-vbeg)/40				; Rate of change per frame

; y_run should evaluate to [1,0] for a horizontal run and [0,1] for a vertical run.
; The rate of change in azimuth for a horizontal run is ~10^-4 radians/frame.
; The rate of change in radius for a vertical run is 0.125 pixels/frame.

y_run = [abs(doptaxis[0]) gt 0.00005, abs(doptaxis[1]) gt 0.04]

if total(y_run) ne 1 then message, 'shoot'

message, /info, hhvv+' : '+(['Horizontal','Vertical'])[y_run[1]]+' run'

doptaxis = y_run*doptaxis						; Only keep one rate of change, set the other to zero

old_center = center0							; The xycenter[0:1,*] values are relative to old_center


; To first order we set all parameters to averages over the whole sequence

center0	  = total(xycenter,2)/nimg				; Center of fov arc
optaxis0  = total(prcenter,2)/nimg				; Azimuth, radius of optical axis
binwidth  = mean(binw)							; Binwidth
dbinwidth = mean(dbinw)							; Change in binwidth along box

; We take two systematic drifts into account:
; the motion of the light box, and the drift in the center of the fov arc

doptaxis  = [0d0,0d0]							; Drift in optical axis
dcenter   = [0d0,0d0]							; Drift in center of fov arc

; Modify some of the parameter to track these drifts

vrun = y_run[1]									; 0: horizontal, 1: vertical run

; Update the X or Y coordinate, and rate of change
; (only X is updated for horizontal run; only Y is updated for vertical run)

vbeg = mean(xycenter[vrun, 0: 4])
vend = mean(xycenter[vrun,44:46])

center0[vrun] = vbeg							; Start value
dcenter[vrun] = (vend-vbeg)/40					; Rate of change

; Update azimuth or R of optical axis, and rate of change
; (only azimuth is updated for horizontal run; only R is updated for vertical run)

vbeg = mean(prcenter[vrun, 0: 4])
vend = mean(prcenter[vrun,44:46])

optaxis0[vrun] = vbeg
doptaxis[vrun] = (vend-vbeg)/40

; For a horizontal run we set dbinwidth to its start value.
; The drift in dbinwidth is taken into account using an analytic
; expression, so we don't need a rate of change.

vbeg = mean(dbinw[0:4])
if 1-vrun then dbinwidth = vbeg


center0 = old_center+center0



dpr = 180d0/!dpi

if paper then			$
	set_page, /portrait	$
else					$
	set_page, /portrait, bmp=bmp, filepath(root=outdir,hhvv+'_geometry.bmp')

if IsDisplay() then begin
	twin, /show, xsize=750, ysize=750
	chars = 2.2
endif else		$
	chars = 1.3


!p.multi = [0,2,3]

xtitle = 'frame number'
Rmean  = mean(prcenter[1,*])

plot, img_nr, xycenter[0,*], /ynozero, xtitle=xtitle, ytitle='X!D0!N-'+strcompress(old_center[0],/rem), chars=chars
plot, img_nr, xycenter[1,*], /ynozero, xtitle=xtitle, ytitle='Y!D0!N-'+strcompress(old_center[1],/rem), chars=chars
plot, img_nr, prcenter[0,*]*dpr+90 , /ynozero, xtitle=xtitle, ytitle='!9P!X!D0!N+90 (deg)', chars=chars
plot, img_nr, prcenter[1,*]-Rmean, /ynozero, xtitle=xtitle, ytitle='R!D0!N-'+strcompress(Rmean,/rem), chars=chars

plot, img_nr, (binw*dpr-1.5)*1e3, /ynozero, xtitle=xtitle, ytitle='binw-1.5 (mdeg)', chars=chars
plot, img_nr, dbinw*dpr*1e3, /ynozero, xtitle=xtitle, ytitle='d(binw) (mdeg/spot)', chars=chars
;plot, img_nr, 2*(binw*dpr-1.5)*1e3, /ynozero, xtitle=xtitle, ytitle='!7Du!X-3 (mdeg)', chars=chars
;plot, img_nr, 4*dbinw*dpr*1e3, /ynozero, xtitle=xtitle, ytitle='!7du!X (mdeg)', chars=chars

!p.multi = 0

if bmp or paper then	$
	get_page			$
else					$
	if stopq then qset_page, /spit

return  &  end