FUNCTION smei_theta2radial, thetay, thetax

;+
; NAME:
;	smei_theta2radial
; PURPOSE:
;	Converts from angle thetay in the sky to radial offset
;	Rccd-R0 on the CCD
; CATEGORY:
;	camera
; CALLING SEQUENCE:
;	R = smei_theta2radial(thetay [, thetax])
; INPUTS:
;	thetay		array[*]; type: any
;					thetay argument for Eq. 3
; OPTIONAL INPUT PARAMETERS:
;	thetax		array[*] or array[2,*]; type: same as thetay
;					if thetax is absent then Equation 2 from the memo is
;					used (= Eq. 3 with thetax=0)
;
;					CASE 1: array[*]: if thetax is an array with the same structure
;						as thetay then the radial offset from Eq. 3 is returned
;					CASE 2: array[2,*]: this is used by smei_radial2theta to determine the
;						zero for Eq 3 for given thetax and radial offset.
;						thetax[0,*]: stores the thetax values
;						thetax[1,*]: stores the radial offset values
; OUTPUTS:
;	R			array[*]; type: same type as thetax
;					CASE 1: radial offset from Eq 3 in pixels
;					CASE 2: radial offset from Eq 3, minus radial offset specified in
;						thetax[1,*] (the value of thetay that makes this difference zero is
;						the value that smei_radial2theta is after).
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	IsType, SubArray
; RESTRICTION:
;	Thetay varies across the field of view from about -1.5 to +1.5 degrees. Since Eq. 3 is
;	an expansion for small thetay it doesn't make much sense to use Eq. 3 outside this range.
;	An arbitrary cutoff at thetay=+/-30 degrees is used to apply Eq. 3. Outside this range
;	only the linear term is applied (have to return something you know).
;	The full expression for the radial offset has a second zero near thetay = -80 degrees
;	(in addition to thetay=0, due to the thetay^4 term). So the cutoff also avoids that
;	near-zero offsets are returned for large negative thetay (suggesting that it's inside
;	the field of view !!)
; PROCEDURE:
;	See Andy's memo "Defining the field of view for the SMEI cameras" (11-July-2001)
;	This function implements Equation 3. All angles are in degrees.
; MODIFICATION HISTORY:
;	FEB-2003, Paul Hick (UCSD/CASS)
;	JAN-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Modification to intercept large thetay values
;-

ty = thetay
rr = where( abs(thetay) GT 30.0 )
IF rr[0] NE -1 THEN ty[rr] = 0.0

CASE IsType(thetax, /defined) OF

0: rr = thetay*(-18.0226 + thetay*(0.1664 + thetay*(-0.00744 + ty*(-0.00015))))

1: BEGIN

	dr_present = n_elements(thetax) EQ 2*n_elements(thetay) 

	CASE dr_present OF
	0: tx = thetax
	1: BEGIN
		tx = SubArray(thetax, element=0)
		dr = SubArray(thetax, element=1)
	END
	ENDCASE

	tx = tx*tx
	rr = 0.000936*tx - 1.044E-7*tx*tx + thetay*(-18.0226 + 3.6E-5*tx + thetay*(0.1664 - 1.0E-5*tx+thetay*(-0.00744 +ty*(-0.00015))))

	IF dr_present THEN rr = rr-dr

END

ENDCASE

RETURN, rr  &  END