integer function r2s(r1,r2,r3,s1,s2,s3,unsym)
      !
      ! For a given r1,r2,r3 this routine 
      ! gives s1,s2,s3 such that
      !
      ! 0 < s3 < s/(2*sqrt(6))
      ! 0 < s2 < s/(2*sqrt(2))
      ! 0 < s3 < sqrt(1/3) s2
      !
      ! So that (s2,s3) point lies in the
      ! region (1) defined in the notes.
      ! On output, the integer between 1 and 6
      ! indicates the region where the r1,r2,r3
      ! point was before the transformations.
      ! 
      implicit none
      integer i
      logical unsym
      real*8 s1,s2,s3,r1,r2,r3,s,q
      real*8 s1o,s2o,s3o
      logical debug
      data debug /.false./

      r2s = 0      
      q = dsqrt(1d0/3d0)
      s = r1+r2+r3
      s1o=(r1+r2+r3)/dsqrt(3.0d0)
      s2o=(r2-r3)/dsqrt(2.0d0)
      s3o=(2.0d0*r1-r2-r3)/dsqrt(6.0d0)

      s1 = s1o 
      s2 = s2o 
      s3 = s3o 

      if(unsym) return   ! return if we do not need to symmetrize

      if(s3o.ge.q*s2o.and.s3o.ge.-q*s2o ) then
!          if(s3.ge.-dsqrt(1d0/3d0)*s2) then
              if(s2o.ge.0d0) then
                  ! region 1
                  ! (123)->(123)                                    
                  i=1
              else
                  ! region 6
                  ! (23)->(32)                                    
                  i=6
                  s2=(r3-r2)/dsqrt(2.0d0)
              end if
      elseif(s3o.ge.q*s2o.and.s3o.lt.-q*s2o ) then
              ! region 5
              ! (123)->(312)                                
              i = 5              
              s2=(r1-r2)/dsqrt(2.0d0)
              s3=(2.0d0*r3-r1-r2)/dsqrt(6.0d0)
      elseif(s3o.lt.q*s2o.and.s3o.lt.-q*s2o ) then
              if(s2o.gt.0d0) then
                  ! region 3
                  ! (123)->(231)                  
                  i=3
                  s2=(r3-r1)/dsqrt(2.0d0)
                  s3=(2.0d0*r2-r3-r1)/dsqrt(6.0d0)                  
              else
                  ! region 4
                  i=4
                  ! (13)->(31)
                  s2=(r2-r1)/dsqrt(2.0d0)
                  s3=(2.0d0*r3-r2-r1)/dsqrt(6.0d0)                  
              end if
      elseif(s3o.lt.q*s2o.and.s3o.ge.-q*s2o) then
              ! region 2
              ! (12)->(21)
              i = 2              
              s2=(r1-r3)/dsqrt(2.0d0)
              s3=(2.0d0*r2-r1-r3)/dsqrt(6.0d0)
      end if
      if(debug) then
      write(*,'(a5,6f12.6,a16,i3)') 'Point',s1o,s2o,s3o,r1,r2,r3,
     $ ' lies in region ',i
      end if
      r2s = i
      
      return
      end