C+
C NAME:
C	rotate_dcm2quat
C PURPOSE:
C	Convert direction cosine matrix to quaternion
C CATEGORY:
C	gen/for/lib
C CALLING SEQUENCE:
	subroutine rotate_dcm2quat(dcm,qq)
C INPUTS:
C	dcm(3,3)	double precision	direction cosine matrix
C OUTPUTS:
C	q(4)		double precision	quaternion
C PROCEDURE:
C MODIFICATION HISTORY:
C	JUL-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    double precision	dcm(3,3)
	    double precision	qq(4)

	double precision	dcm_trace
	double precision	tmp
	double precision	xx
	double precision	yy
	double precision	zz

	dcm_trace = dcm(1,1)+dcm(2,2)+dcm(3,3)

	! The trace of the DCM matrix is 3*q0*q0-(qx*qx+qy*qy+qz*qz) = 4*q0*q0-1
	! If this is positive then abs(q0) > 1/2.

	if (dcm_trace .gt. 0.0d0) then

            ! Take the positive solution for q0

	    tmp = dsqrt(1.0d0+dcm_trace)
	    qq(4) = 0.5d0*tmp
	    tmp = 0.5d0/tmp

	    ! Difference of terms on opposite sides of diagonal

	    qq(1) = (dcm(3,2)-dcm(2,3))*tmp
    	    qq(2) = (dcm(1,3)-dcm(3,1))*tmp
	    qq(3) = (dcm(2,1)-dcm(1,2))*tmp

	else

	    xx = dcm(1,1)
	    yy = dcm(2,2)
	    zz = dcm(3,3)

	    if (xx .gt. yy .and. xx .gt. zz) then
		i = 1
	    else if (yy .gt. zz) then
		i = 2
	    else
		i = 3
	    end if

	    j = mod(i,3)+1	
	    k = mod(j,3)+1

	    tmp = dsqrt(1.0d0+(dcm(i,i)-(dcm(j,j)+dcm(k,k))))
	    qq(i) = 0.5d0*tmp
	    tmp = 0.5d0/tmp

	    qq(4) = (dcm(k,j)-dcm(j,k))*tmp
	    qq(j) = (dcm(j,i)+dcm(i,j))*tmp
	    qq(k) = (dcm(k,i)+dcm(i,k))*tmp

	end if

	return
	end