C+
C NAME:
C	smei_axis_ccd
C PURPOSE:
C	Convert rectangular CCD coordinates to polar coordinates
C	relative to optical axis and v.v.
C CALLING SEQUENCE:
	subroutine smei_axis_ccd(id,icam,mode,x,y,dr,phi,r)
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	x,y	double precision	pixel coordinates for specified input mode
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					(usually dr,phi are output from smei_axis_cam)
C OUTPUTS:
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					(usually dr,phi are output from smei_axis_cam)
C	id=0:
C	x,y	double precision	pixel coordinates for specified input mode
C
C	r	double precision	radial distance to origin of polar
C					transformation, i.e. r,phi are polar coordinates
C					on the CCD, and r = r0+dr, where r0 is the
C					distance of the optical axis from the origin.
C SEE ALSO:
C	smei_axis_cam
C PROCEDURE:
C	Binned (mode 1,2) and unbinned (mode 0, engineering mode) pixel
C	indices are related by (nbin = 2^mode):
C	   (unbinned index) = nbin*(binned index)-0.5*(nbin-1)
C MODIFICATION HISTORY:
C	FEB-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer		id
	    integer		icam
	    integer		mode
	    double precision	x
	    double precision	y
	    double precision	dr
	    double precision	phi
	    double precision	r

	!r0	distance center of fov curvature to optical axis (in pixels)
	!x0,y0	center of fov curvature (origin of polar coordinates)

	double precision	r0(3)		/1193.5d0, 1188.5d0, 1198.5d0/
	double precision 	x0(3)		/ 634.0d0, 626.7d0, 637.7d0/
	double precision 	y0(3)		/1242.1d0,1239.1d0,1244.9d0/

	!double precision 	affine(6)	/6*0.0d0/

	double precision	a
	double precision	b
	double precision	dbin

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

	nbin = 2**mode				! Binning factor
	dbin = 0.5d0*(nbin-1)

	if (mod(id,2) .eq. 0) then

	    a = nbin*x-dbin			! Convert to mode 0 (engineering mode) pixel index
	    b = nbin*y-dbin			! unbinned index = nbin*(binned index)-0.5*(nbin-1)

	    a = a-x0(icam)			! Shift to origin of polar coordinates
	    b = b-y0(icam)

	    ! Affine transformation (plate scale changes?)

	    !r = a
	    !a = affine(1)+(1.0d0+affine(2))*r+       affine(3) *b
	    !b = affine(4)+       affine(5) *r+(1.0d0+affine(6))*b

	    r   = dsqrt(a*a+b*b)		! Convert to polar coordinates
	    phi = datan2d(a,-b)			! b < 0 inside arc of fov

	    dr  = r-r0(icam)			! Radial distance to optical axis on CCD

	else

	    r =  r0(icam)+dr
	    a =  r*dsind(phi)			! Convert to rectangular coordinates
	    b = -r*dcosd(phi)			! relative to origin of polar coordinates

	    a = x0(icam)+a			! Shift origin
	    b = y0(icam)+b			! Engineering mode pixels

	    x = (a+dbin)/nbin			! Converto binned pixel coordinates
	    y = (b+dbin)/nbin			! binned index = (unbinned index + 0.5*(nbin-1))/nbin

	end if

	if (id .ge. 2) r = r/r0(icam)

	return
	end