C+
C NAME:
C	smei_axis_cam
C PURPOSE:
C	Convert polar CCD coordinates to sky coordinates
C	in camera frame, and v.v
C CALLING SEQUENCE:
	subroutine smei_axis_cam(id,icam,mode,dr,phi,rx,ry,rz,tx,ty)
C INPUTS:
C	id	integer			id=0: CCD pixel coord to camera sky coord
C					id=1: camera sky coord to CCD pixel coord
C	icam	integer			camera id (1,2,3)
C	mode	integer			mode (0,1,2)
C
C	id=0:
C	dr	double precision	radial distance to optical axis in
C					units of engineering mode pixels
C	phi	double precision	angle from optical axis in degrees
C
C	id=1:
C	rx	double precision	x-component of unit vector
C	ry	double precision	y-component of unit vector
C	rz	double precision	z-component of unit vector (always positive)
C OUTPUTS:
C	id=0:
C	rx	double precision	x-component of unit vector
C	ry	double precision	y-component of unit vector
C	rz	double precision	z-component of unit vector (always positive)
C
C	id=1:
C	dr	double precision	radial distance to optical axis in
C					units of engineering mode pixels
C	phi	double precision	angle from optical axis in degrees
C
C	tx	double precision	theta-x = arcsin(rx)
C	ty	double precision	theta-y = arcsin(ry)
C CALLS:
C	BadR8
C RESTRICTIONS:
C	The input values for x,y need to be sensible, i.e. not too far
C	outside the field of view.
C RESTRICTIONS:
C	If a unit vector is specified (id=1) then the conversion
C	is only done if the location is within a reasonable
C	distance of the fov (currently abs(tx) < 45 and abs(ty) < 10
C	degrees. For locations outside this range dr and phi are
C	set to BadR8() on return.
C PROCEDURE:
C	(rx,ry,rz) is a unit vector representing the location
C	on the sky of (x,y) on the CCD, i.e. rx,ry,rz are the
C	direction cosines.
C	Quadratic coefficients from Andys July 2001 memo. In addition some
C	ad-hoc stretching of tx and ty is done
C MODIFICATION HISTORY:
C	FEB-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer		id
	    integer		icam
	    integer		mode
	    double precision	dr
	    double precision	phi
	    double precision	rx
	    double precision	ry
	    double precision	rz
	    double precision	tx
	    double precision	ty

	! Ad-hoc correction to tx and ty scale

	double precision	tx_scale1(3)	/ 0.9995d0, 0.99700d0, 1.00480d0/
	double precision	tx_scale2(3)	/ 0.0000d0, 0.00000d0,-0.00007d0/
	double precision	ty_scale (3)	/-0.0001d0,-0.00005d0,-0.00016d0/

	double precision	a
	double precision	b
	double precision	c
	double precision	tx2

	double precision	dsind		! Needed for GNU compiler
	double precision	datan2d

	double precision	BadR8

	nbin = 2**mode				! Binning factor

	if (id .eq. 0) then

	    tx  = (tx_scale1(icam)+tx_scale2(icam)*phi)*phi
	    tx2 = tx*tx

	    a  = 0.1664d0 - 0.00001d0*tx2
	    b  = 0.000036d0*tx2 - 18.0226d0
	    c  = 0.000936d0*tx2 - dr

	    ty = (-b-dsqrt(b*b-4.0d0*a*c))/(2.0d0*a)
	    ty = ty+ty_scale(icam)*tx2

	    rx = dsind(tx)
	    ry = dsind(ty)
	    rz = dsqrt(1.0d0-(rx*rx+ry*ry))

	else

	    tx = datan2d(rx,dsign(1.0d0,rz)*dsqrt(1.0d0-rx*rx))
	    ty = datan2d(ry,dsign(1.0d0,rz)*dsqrt(1.0d0-ry*ry))

	    ! Apply the conversion to CCD coordinates over a range of tx,ty
	    ! that generously covers the FOV.

	    if (dabs(tx) .lt. 45.0d0 .and. dabs(ty) .lt. 10.0d0) then

		tx2 = tx*tx

		dr = ty-ty_scale(icam)*tx2
		dr = (0.1664d0-0.00001d0*tx2)*dr*dr+(0.000036d0*tx2-18.0226d0)*dr+0.000936d0*tx2

		phi = (tx-tx_scale2(icam)*tx*tx)/tx_scale1(icam)

	    else

		dr  = BadR8()
		phi = dr

	    end if

	end if

	return
	end