function TMO_skymotion, rr, p_axis, frame, degrees=degrees, fov=fov, vv=vv,	$
	plotx=plotx, south=south, north=north

;+
; NAME:
;	TMO_skymotion
; PURPOSE:
;	Corrects positions on CCD for skymotion during transit
; CATEGORY:
;	SMEI prototype camera: TMO
; CALLING SEQUENCE:
;	rr_out = TMO_skymotion(rr [, p_axis, frame, /degrees, to_zenith=to_zenith,	$
;		fov=fov, vv=vv, plotx=plotx
; INPUTS:
;	rr				array[2,n]; type: float
;						positions on ccd as pairs off azimuthal angle (along long
;						dimension of fov), and radius
;	p_axis			array[2]; type: float; default: [0,0]
;						position of optical axis on ccd: azimuthal angle and radius
;						Currently only p_axis[0] is used.
;	frame			array[m]; type: any; default: 1
;						list of frame numbers (see PROCEDURE)
; OPTIONAL INPUT PARAMETERS:
;	/south			if set, the south position (19.5 south of zenith) is assumed
;	/north			if set, the north position (10.0 north of zenith) is assumed
;	/degrees		if set, all angles are in degrees (default: radians)
;	/plotx			plot results
; OUTPUTS:
;	rr_out			array[2,n,m]; type: float
;						positions rr, corrected for motion across the sky
; OPTIONAL OUTPUT PARAMETERS:
;	fov=fov			array[*]; CCD azimuth angle (-30,..,+30) of points along
;						long dimension of fov
;	vv=vv			array[*]; type: float
;						skymotion at positions 'fov' in pixels/minute
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; RESTRICTIONS:
;	No check is made whether the locations rr are actually inside the fov
; CALLS:
;	gridgen, ToRadians, BadValue, SubArray, SuperArray, SyncDims
; PROCEDURE:
; >	The SMEI fov is assumed to be along the local meridian
; > The motion of stars during transit are calculated in pixels/minute
;	under the assumption that the 3 degree fov is covered by 27.16 pixels, and
;	that the geographic latitude of TMO is 34.382 degree.
; > The 'frame' array is used to define the time axis in units of minutes.
;	(TMO data were taken at a 1 minute cadence). Usually 'frame' will be an
;	integer array of frame numbers, but the effective cadence time can be adjusted
;	by multiplication of the 'frame' array with the cadence time in minutes.
; >	For each location 'rr' on the CCD inside the fov, the corrected position is
;	the position of 'rr' after a time 'frame' minutes has elapsed.
; >	Only the radial component r[1,*] is corrected, i.e. stars are assumed to
;	transit on the CCD in the radial direction.
; MODIFICATION HISTORY:
;	JUL-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

rpm = ToRadians(degrees=degrees)
InitVar, p_axis, [0,0]
InitVar, frame , 1
InitVar, plotx , /key

case 1 of
south: begin
	to_zenith = -19.5
	message, /info, 'pointing south'+strcompress(to_zenith)+' degrees'
end
north: begin
	to_zenith = +10.0
	message, /info, 'pointing north'+strcompress(to_zenith)+' degrees'
end
else: begin
	to_zenith =   0.0*rpm*!radeg
	message, /info, 'pointing to zenith'
end
endcase

to_zenith = to_zenith/!radeg

colat_tmo = (90-34.382)/!radeg			; Geographic latitude of TMO is 34.382 degrees
halfw_deg = 1.5/!radeg					; Full width FOV is 3 degrees
halfw_pix = 27.16 ;26.6					; Halfwidth FOV in pixels ???
ns_offset = 0/!radeg					; North-south offset (~ 2 degrees??)
										; 'fov' covers field of view south to north relative to zenith
fov	  = gridgen(65,range=32/!radeg*[-1,1])
colat = colat_tmo-fov-to_zenith*rpm 	; Co-latitude along FOV from south to north

gamma  = (90+ns_offset)/!radeg

sindha = sin(halfw_deg)*sin(gamma)
cosdha = cos(halfw_deg)*sin(colat)-sin(halfw_deg)*cos(colat)*cos(gamma)

dha    = atan(sindha, cosdha)			; Hour angle from center to edge of FOV (increases from S to N)

dt     = dha/(2*!pi)*24*60				; Time it takes for a star to move from center to edge (minutes)
vv     = -halfw_pix/dt					; # pixels/minute of motion for stars across fov

nr     = n_elements(rr)/2				; # positions

rr_out = BadValue(0.0)					; Return value if something goes wrong

if nr ne 0 then begin					; Coordinates present
	vr = SubArray(rr)					; Azimuth subarray (angle along long dim of fov)
	nf = where(finite(vr))				; Pick up finite coordinates

	if nf[0] ne -1 then begin			; Finite coordinates present
		sz = size(rr)
		rr = reform(rr, 2, nr, /overwrite)

		nimg = n_elements(frame)

		; Currently only radial motion is incorporated

		vr = interpol(vv, fov, (vr[nf]-p_axis[0])*rpm)
		message, /info, 'average radial motion'+strcompress(mean(vr))+' pixels/frame'

		vr = SuperArray(vr, nimg, /trail)*SuperArray(frame, n_elements(vr), /lead)

		rr_out = SuperArray(rr, nimg, /trail)
		rr_out[1,nf,*] = rr_out[1,nf,*]+vr	; Add radial displacements

		SyncDims, rr, size=sz
		SyncDims, rr_out, size=[sz[0]+1,sz[1:sz[0]],nimg,sz[sz[0]+1],sz[sz[0]+2]*nimg]
	endif
endif

fov = fov/rpm

if plotx then begin 					; Plot test case
	plot, fov*rpm*!radeg, vv, xstyle=1, xtitle='Zenith distance', ytitle='Pixels/minute'
	oplot, [30.66,15.23,-4.1,-16.15]-to_zenith*rpm*!radeg , -[1.877,2.93,3.86,4.2], psym=2
endif

return, rr_out  &  end