function smei_theta2radial(thetay, thetax)

	    real	thetay
	    real	thetax

	bad = BadR4()

	ty = thetay
	if (abs(thetay) .gt. 30.0) ty = 0.0

	tx = thetax*thetax
	smei_theta2radial = 9.36e-4*tx-1.044e-7*tx*tx+thetay*
     &		(-18.0226+3.6d-5*tx+thetay*(0.1664-1d-5*tx+thetay*(-7.44e-3+ty*(-1.5e-4))))

	return
	end

	function smei_theta2delta_r(thetay)

	    real	thetay

	save		dr, tx

	smei_theta2delta_r = smei_theta2radial(thetay, tx)-dr

	return

	entry smei_theta2delta_r_set(delta_r, thetax)

	smei_theta2delta_r_set = 0.0	! Have to return something

	dr = delta_r
	tx = thetax

	return
	end

	function smei_radial2theta(delta_r, thetax)

	    real	delta_r
	    real	thetax

	logical		nrZbrac
	real		nrZbrent

	tol = 1.0d-3

	! The first estimate of the zero point is the lowest order solution
	! (it may be possible to improve on this by including one more higher order term,
	! but it probably won't make any difference).

	tmp = -delta_r/18.0226e0

	! Set up the initial range of values close to the zero.

	ty1 = tmp-(1.0/18.0226)/2.0
	ty2 = tmp+(1.0/18.0226)/2.0

	! Function smei_theta2delta_r needs values of thetax and delta_r
	! to be able to find the correct zero

	tmp = smei_theta2delta_r_set(delta_r, thetax)

	! The above range does not necessarily bracket the zero.
	! Expand the range outward until it does.

	if (nrZbrac('smei_theta2delta_r', ty1, ty2, 0)) then

	    ! The range ty1, ty2 should now bracket the zero. Warnings are displayed by
	    ! nrZBrent if it doesn't.

	    thetay = nrZBrent('smei_theta2delta_r', ty1, ty2, tol)

	else
	    thetay = BadR4()

	end if

	smei_radial2theta = thetay

	return
	end