C+
C NAME:
C	smei_ccd2cam
C PURPOSE:
C	Convert CCD coordinates to sky coordinates
C CALLING SEQUENCE:
	subroutine smei_ccd2cam(icam,r0,x,y,rx,ry,rz,tx,ty)
C INPUTS:
C	icam	integer			camera id (1,2,3)
C	r0	double precision	distance center of fov curvature to optical axis (in pixels)
C	x	double precision	hor distance to center of fov curvature (in pixels)
C	y	double precision	vert distance to center of fov curvature (in pixels)
C OUTPUTS:
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	tx	double precision	theta-x = arcsin(rx)
C	ty	double precision	theta-y = arcsin(ry)
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 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	NOV-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer		icam
	    double precision	r0
	    double precision	x
	    double precision	y
	    double precision	rx
	    double precision	ry
	    double precision	rz
	    double precision	tx
	    double precision	ty

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

	double precision	dsind		! Needed for GNU compiler
	double precision	dcosd
	double precision	datan2d

	double precision	BadR8

	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/

	a = x		! In case input x,tx or y,ty are same variable
	b = y

	tx = datan2d(a,-b)	! b < 0 inside fov
	if (icam .ne. 0) tx = tx_scale1(icam)*tx+tx_scale2(icam)*tx*tx

	tx2 = tx*tx
	r = dsqrt(a*a+b*b)-r0

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

	ty = (-b-dsqrt(b*b-4d0*a*c))/(2d0*a)
	if (icam .ne. 0) ty = ty+ty_scale(icam)*tx2

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

	return
C+
C NAME:
C	smei_cam2ccd
C PURPOSE:
C	Convert sky coordinates in camera fram to CCD coordinates
C CALLING SEQUENCE:
	entry smei_cam2ccd(icam,r0,x,y,rx,ry,rz,tx,ty)
C INPUTS:
C	icam	integer			camera id (1,2,3)
C	r0	double precision	distance center of fov curvature to optical axis (in pixels)
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
C OUTPUTS:
C	x	double precision	hor distance to center of fov curvature (in pixels)
C	y	double precision	vert distance to center of fov curvature (in pixels)
C	tx	double precision	theta-x = arcsin(rx)
C	ty	double precision	theta-y = arcsin(ry)
C RESTRICTIONS:
C	The conversion is only performed for sensible (rx,ry,rz), i.e. for unit
C	vectors that project close to the field of view.
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	JAN-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	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. 31.0d0 .and. dabs(ty) .lt. 3.0d0) then

	    tx2 = tx*tx

	    r = ty
	    if (icam .ne. 0) r = r-ty_scale(icam)*tx2
	    r = r0+(0.1664d0-0.00001d0*tx2)*r*r+(0.000036d0*tx2-18.0226d0)*r+0.000936d0*tx2

	    x = tx
	    if (icam .ne. 0) x = (x-tx_scale2(icam)*x*x)/tx_scale1(icam)

	    y = -r*dcosd(x)
	    x =  r*dsind(x)

	else

	    x = BadR8()
	    y = x

	end if

	return
	end