subroutine cart2int(cart,qint)
      implicit none
!     This version merges both coordinate transformation routines into
!     one.  JTmod's sscales(2:3) are ignored.
!     This is the first version to be compatible with one of my proper 6D fits
!     Time-stamp: <2024-10-22 13:52:59 dwilliams>

!     Input (cartesian, in Angström)
!     cart(:,1):   N
!     cart(:,1+i): Hi
!     Output
!     qint(i):     order defined in JTmod.
!     Internal Variables
!     no(1:3):     NO distances 1-3
!     pat_in:      temporary coordinates
!     axis:        main axis of NO3
      include 'nnparams.incl'
      include 'JTmod.incl'


      real*8 cart(3,4),qint(maxnin)

      real*8 no(3), r1, r2, r3
      real*8 v1(3), v2(3), v3(3)
      real*8 n1(3), n2(3), n3(3), tr(3)
      real*8 ortho(3)
      real*8 pat_in(maxnin)
      logical ignore_umbrella,dbg_umbrella
      logical dbg_distances

!..   Debugging parameters
!..   set umbrella to 0
      parameter (ignore_umbrella=.false.)
!     parameter (ignore_umbrella=.true.)

!..   break if umbrella is not 0
      parameter (dbg_umbrella=.false.)
!      parameter (dbg_umbrella=.true.)

!..   break for tiny distances
      parameter (dbg_distances=.false.)
!      parameter (dbg_distances=.true.)

      integer k

!..   get N-O vectors and distances:
      do k=1,3
         v1(k)=cart(k,2)-cart(k,1)
         v2(k)=cart(k,3)-cart(k,1)
         v3(k)=cart(k,4)-cart(k,1)
      enddo
      no(1)=norm(v1,3)
      no(2)=norm(v2,3)
      no(3)=norm(v3,3)

!..   temporarily store displacements
      do k=1,3
         pat_in(k)=no(k)-offsets(1)
      enddo

      do k=1,3
         v1(k)=v1(k)/no(1)
         v2(k)=v2(k)/no(2)
         v3(k)=v3(k)/no(3)
      enddo

!..   compute three normal vectors for the ONO planes:
      call xprod(n1,v1,v2)
      call xprod(n2,v2,v3)
      call xprod(n3,v3,v1)

      do k=1,3
         tr(k)=(n1(k)+n2(k)+n3(k))/3.d0
      enddo
      r1=norm(tr,3)
      do k=1,3
         tr(k)=tr(k)/r1
      enddo

!     rotate trisector
      call rot_trisec(tr,v1,v2,v3)

!..   determine trisector angle:
      if (ignore_umbrella) then
         pat_in(7)=0.0d0
      else
         pat_in(7)=pi/2.0d0 - acos(scalar(v1,tr,3))
         pat_in(7)=sign(pat_in(7),cart(1,2))
      endif

!..   molecule now lies in yz plane, compute projected ONO angles:
      v1(1)=0.d0
      v2(1)=0.d0
      v3(1)=0.d0
      r1=norm(v1,3)
      r2=norm(v2,3)
      r3=norm(v3,3)
      do k=2,3
         v1(k)=v1(k)/r1
         v2(k)=v2(k)/r2
         v3(k)=v3(k)/r3
      enddo

!     make orthogonal vector to v3
      ortho(1)=0.0d0
      ortho(2)=v3(3)
      ortho(3)=-v3(2)

!..   projected ONO angles in radians
      pat_in(4)=get_ang(v2,v3,ortho)
      pat_in(5)=get_ang(v1,v3,ortho)

      pat_in(6)=dabs(pat_in(5)-pat_in(4))

!..   account for rotational order of atoms
      if (pat_in(4).le.pat_in(5)) then
         pat_in(5)=2*pi-pat_in(4)-pat_in(6)
      else
         pat_in(4)=2*pi-pat_in(5)-pat_in(6)
      endif

      pat_in(4)=rad2deg*pat_in(4)-offsets(2)
      pat_in(5)=rad2deg*pat_in(5)-offsets(2)
      pat_in(6)=rad2deg*pat_in(6)-offsets(2)
      pat_in(7)=rad2deg*pat_in(7)

      call genANN_ctrans(pat_in)

      qint(:)=pat_in(:)

      contains

!-------------------------------------------------------------------
!     compute vector product n1 of vectors v1 x v2
      subroutine xprod(n1,v1,v2)
      implicit none

      real*8 n1(3), v1(3), v2(3)

      n1(1) = v1(2)*v2(3) - v1(3)*v2(2)
      n1(2) = v1(3)*v2(1) - v1(1)*v2(3)
      n1(3) = v1(1)*v2(2) - v1(2)*v2(1)

      end subroutine


!-------------------------------------------------------------------
!     compute scalar product of vectors v1 and v2:
      real*8 function scalar(v1,v2,n)
      implicit none
      integer i, n
      real*8 v1(*), v2(*)

      scalar=0.d0
      do i=1,n
         scalar=scalar+v1(i)*v2(i)
      enddo

      end function

!-------------------------------------------------------------------
!     compute norm of vector:
      real*8 function norm(x,n)
      implicit none
      integer i, n
      real*8 x(*)

      norm=0.d0
      do i=1,n
         norm=norm+x(i)**2
      enddo
      norm=sqrt(norm)

      end function

!-------------------------------------------------------------------

      subroutine rot_trisec(tr,v1,v2,v3)
      implicit none

      real*8 tr(3),v1(3),v2(3),v3(3)

      real*8 vrot(3)
      real*8 rot_ax(3)
      real*8 cos_phi,sin_phi

!     evaluate cos(-phi) and sin(-phi), where phi is the angle between
!     tr and (1,0,0)
      cos_phi=tr(1)
      sin_phi=dsqrt(tr(2)**2+tr(3)**2)

      if (sin_phi.lt.1.0d-12) then
         return
      endif

!     determine rotational axis
      rot_ax(1) = 0.0d0
      rot_ax(2) = tr(3)
      rot_ax(3) = -tr(2)
!     normalize
      rot_ax=rot_ax/sin_phi

!     now the rotation can be done using Rodrigues' rotation formula
!     v'=v*cos(p) + (k x v)sin(p) + k (k*v) (1-cos(p))
!     for v=tr k*v vanishes by construction:

!     check that the rotation does what it should
      call rodrigues(vrot,tr,rot_ax,cos_phi,sin_phi)
      if (dsqrt(vrot(2)**2+vrot(3)**2).gt.1.0d-12) then
         write(6,*) "ERROR: BROKEN TRISECTOR"
         stop
      endif
      tr=vrot

      call rodrigues(vrot,v1,rot_ax,cos_phi,sin_phi)
      v1=vrot
      call rodrigues(vrot,v2,rot_ax,cos_phi,sin_phi)
      v2=vrot
      call rodrigues(vrot,v3,rot_ax,cos_phi,sin_phi)
      v3=vrot


      end subroutine

!-------------------------------------------------------------------

      subroutine rodrigues(vrot,v,axis,cos_phi,sin_phi)
      implicit none
      real*8 vrot(3),v(3),axis(3)
      real*8 cos_phi,sin_phi

      real*8 ortho(3)

      call xprod(ortho,axis,v)
      vrot = v*cos_phi + ortho*sin_phi
     >        + axis*scalar(axis,v,3)*(1-cos_phi)

      end subroutine

!-------------------------------------------------------------------

      real*8 function get_ang(v,xaxis,yaxis)
      implicit none
!     get normalized [0:2pi) angle from vectors in the yz plane
      real*8 v(3),xaxis(3),yaxis(3)

      real*8 phi

      real*8 pi
      parameter (pi=3.141592653589793d0)

      phi=atan2(scalar(yaxis,v,3),scalar(xaxis,v,3))
      if (phi.lt.0.0d0) then
         phi=2*pi+phi
      endif
      get_ang=phi

      end function

      end subroutine cart2int