function SkyStereoStates, TT, Loc, Dir, Fov, Tilt, stereo=Stereo, degrees=Degrees
;+
; NAME:
;	SkyStereoStates
; PURPOSE:
;	Takes a single viewing state and splits it up into two states
;	representing a stereo pair
; CATEGORY:
;	Number strangling
; CALLING SEQUENCE:
;	R = SkyStereoStates(TT, Loc, Dir, Fov, Tilt, stereo=Stereo, /degrees)
; INPUTS:
;	/stereo
;	stereo=[Angle,Distance]
;			/stereo		if set, stereo pairs are produced using default
;						settings for stereo angle (6 degrees) and
;						distance (Sun-viewer distance)
;			stereo=[Angle, Distance]
;						if specified as a two-element array, the first
;						element is used as stereo angle, the second as
;						stereo distance. If the distance is set to zero
;						then the viewer-origin distance Loc[2,*] is used.
;	TT[n], Loc[3,n], Dir[2,n], Fov[n], Tilt[n]
;			arrays describing n viewing state (see SkyImage)
; OPTIONAL INPUT PARAMETERS:
;	/degrees	if set, all angles are in degrees (default: radians)
; OUTPUTS:
;	R			scalar; type: integer
;						0 if no stereo states were calculated
;						1 if stereo states were calculate
;
; OPTIONAL OUTPUT PARAMETERS:
;	(Only calculated if return value R = 1)
;	TT[2*n], Loc[3,2*n], Dir[2,2*n], Fov[2*n], Tilt[2*n]
;				The stereo pairs for the n input states are stored as
;				TT[2*i  ]		right eye view for state i=0,n-1
;				TT[2*i+1]		left  eye view for state i=0,n-1
;				(Note that TT, Fov and Tilt values will be the
;				same for left and right view; Loc and Dir will be
;				slightly different from the input Loc and Dir)
; CALLS:
;	ToRadians, EulerRotate
; PROCEDURE:
; >	Called by href=SkyImage= to set up stereo pairs
;	In stereo mode each state is split into a stereo pair.
; >	Coordinate systems are defined in document SkyGrid.doc (probably in D:\work\doc).
; >	The stereo pair states are defined in the X,Y,Z coordinate system (+Y is the
;	horizontal direction toward the right). In this coordinate system the viewer
;	location is at the origin and the viewing direction is the positive Z-axis.
; >	The stereo pair is defined by two viewing locations close to the origin and
;	two viewing direction close to the +Z axis.
;	Take a point P on the +Z axis at distance StereoDist. Take viewing locations
;	R on the +Y axis and L on -Y axis at equal distance to the origin such that
;	the lines RP and YP connecting these points to P make an angle StereoAngle
;	L and R are the positions of left and right eye, respectively.
;	LP and RP are the viewing directions from left and right eye, respectively.
;	The separation LR is 2*Dist*tan(StereoAngle/2.)
; MODIFICATION HISTORY:
;	FEB-1998, Paul Hick (UCSD/CASS)
;	MAY-2000, Paul Hick (UCSD/CASS, pphick@ucsd.edu); converted from procedure to function
;-
on_error, 2

COMPILE_OPT hidden

if n_elements(Stereo) ne 2 and not keyword_set(Stereo) then return, 0B

n = Stereo
if n_elements(n) ne 2 then if Degrees then n = [6.,0.] else n = [6./!radeg,0.]
StereoAngle = n[0]
StereoDist  = n[1]

rpd = ToRadians(degrees=Degrees)

n  = n_elements(TT)					; # input states
n2 = n*2							; # output states (left/right pairs)

; Convert Loc from spherical to Cartesian coordinates

Loc = cv_coord(from_sphere=Loc, /to_rect, degrees=Degrees)

; Double the state arrays by interleaving states for left and right views

TT	 = reform([TT]  ,1,n)
Fov  = reform([Fov] ,1,n)
Tilt = reform([Tilt],1,n)

TT   = reform( [TT  , TT  ]   , n2)
Loc  = reform( [Loc , Loc ], 3, n2)
Dir  = reform( [Dir , Dir ], 2, n2)
Fov  = reform( [Fov , Fov ]   , n2)
Tilt = reform( [Tilt, Tilt]   , n2)

tanAngle = tan(StereoAngle/2.*rpd)

print, 'skystereostates
for i=0,n-1 do begin
	i2 = i*2				    		; Use viewer-origin distance for stereo distance
    if StereoDist gt 0 then Dist = StereoDist else Dist = sqrt(total(Loc[*,i2]*Loc[*,i2]))

	dOff = Dist*tanAngle
;			  	Right view		Left view
	Off = [ [0., dOff,  0.], [0.,-dOff,  0.] ,	$	; Positions R,L on +Y/-Y axis
			[0.,-dOff,Dist], [0., dOff,Dist]	]	; Viewing directions RP, LP in YZ plane

; The Euler angles ABG transform the vectors Off defined in the X,Y,Z coordinate
; system to the x',y',z' coordinate system. A rotation over -Tilt transforms from
; X,Y,Z system to x'',y'',z'' system
print, off
print, dir[1,i2],dir[0,i2]
    Off = EulerRotate([-Tilt[i2]*rpd, !pi/2.-Dir[1,i2]*rpd, !pi-Dir[0,i2]*rpd], Off, /rectangular)
print, off
	Loc[*,i2:i2+1] = Loc[*,i2:i2+1]+Off[*,0:1]		; Positions R,L in x,y,z coordinate system
													; (at this point Loc contains Cartesian coordinates)
													; Convert viewing directions RP, LP to spherical coordinates
    Off = ( cv_coord(from_rect=Off[*,2:3], /to_sphere, degrees=Degrees) )

    Dir[*,i2:i2+1] = Off[0:1,*]						; Store the direction angles back into Dir

endfor

Loc = cv_coord(from_rect=Loc, /to_sphere, degrees=Degrees)	; Convert Loc back to spherical coordinates

return, 1B  &  end