c-------------------------------------------------------------

      subroutine local_move_init(debug)
crc      implicit none

c     Includes
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'  ! Needed by COMMON.LOCAL
      include 'COMMON.GEO'  ! For pi, deg2rad
      include 'COMMON.LOCAL'  ! For vbl
      include 'COMMON.LOCMOVE'

c     INPUT arguments
      logical debug


c     Determine wheter to do some debugging output
      locmove_output=debug

c     Set the init_called flag to 1
      init_called=1

c     The following are never changed
      min_theta=60.D0*deg2rad  ! (0,PI)
      max_theta=175.D0*deg2rad  ! (0,PI)
      dmin2=vbl*vbl*2.*(1.-cos(min_theta))
      dmax2=vbl*vbl*2.*(1.-cos(max_theta))
      flag=1.0D300
      small=1.0D-5
      small2=0.5*small*small

c     Not really necessary...
      a_n=0
      b_n=0
      res_n=0

      return
      end

c-------------------------------------------------------------

      subroutine local_move(n_start, n_end, PHImin, PHImax)
c     Perform a local move between residues m and n (inclusive)
c     PHImin and PHImax [0,PI] determine the size of the move
c     Works on whatever structure is in the variables theta and phi,
c     sidechain variables are left untouched
c     The final structure is NOT minimized, but both the cartesian
c     variables c and the angles are up-to-date at the end (no further
c     chainbuild is required)
crc      implicit none

c     Includes
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.GEO'
      include 'COMMON.CHAIN'
      include 'COMMON.VAR'
      include 'COMMON.MINIM'
      include 'COMMON.SBRIDGE'
      include 'COMMON.LOCMOVE'

c     External functions
      integer move_res
      external move_res
      double precision ran_number
      external ran_number

c     INPUT arguments
      integer n_start, n_end  ! First and last residues to move
      double precision PHImin, PHImax  ! min/max angles [0,PI]

c     Local variables
      integer i,j
      double precision min,max
      integer iretcode


c     Check if local_move_init was called.  This assumes that it
c     would not be 1 if not explicitely initialized
      if (init_called.ne.1) then
        write(6,*)'   ***   local_move_init not called!!!'
        stop
      endif

c     Quick check for crazy range
      if (n_start.gt.n_end .or. n_start.lt.1 .or. n_end.gt.nres) then
        write(6,'(a,i3,a,i3)')
     +       '   ***   Cannot make local move between n_start = ',
     +       n_start,' and n_end = ',n_end
        return
      endif

c     Take care of end residues first...
      if (n_start.eq.1) then
c     Move residue 1 (completely random)
        theta(3)=ran_number(min_theta,max_theta)
        phi(4)=ran_number(-PI,PI)
        i=2
      else
        i=n_start
      endif
      if (n_end.eq.nres) then
c     Move residue nres (completely random)
        theta(nres)=ran_number(min_theta,max_theta)
        phi(nres)=ran_number(-PI,PI)
        j=nres-1
      else
        j=n_end
      endif

c     ...then go through all other residues one by one
c     Start from the two extremes and converge
      call chainbuild
      do while (i.le.j)
        min=PHImin
        max=PHImax
c$$$c     Move the first two residues by less than the others
c$$$        if (i-n_start.lt.3) then
c$$$          if (i-n_start.eq.0) then
c$$$            min=0.4*PHImin
c$$$            max=0.4*PHImax
c$$$          else if (i-n_start.eq.1) then
c$$$            min=0.8*PHImin
c$$$            max=0.8*PHImax
c$$$          else if (i-n_start.eq.2) then
c$$$            min=PHImin
c$$$            max=PHImax
c$$$          endif
c$$$        endif

c     The actual move, on residue i
        iretcode=move_res(min,max,i,c)  ! Discard iretcode
        i=i+1

        if (i.le.j) then
          min=PHImin
          max=PHImax
c$$$c     Move the last two residues by less than the others
c$$$          if (n_end-j.lt.3) then
c$$$            if (n_end-j.eq.0) then
c$$$              min=0.4*PHImin
c$$$              max=0.4*PHImax
c$$$            else if (n_end-j.eq.1) then
c$$$              min=0.8*PHImin
c$$$              max=0.8*PHImax
c$$$            else if (n_end-j.eq.2) then
c$$$              min=PHImin
c$$$              max=PHImax
c$$$            endif
c$$$          endif

c     The actual move, on residue j
          iretcode=move_res(min,max,j,c)  ! Discard iretcode
          j=j-1
        endif
      enddo

      call int_from_cart(.false.,.false.)

      return
      end

c-------------------------------------------------------------

      subroutine output_tabs
c     Prints out the contents of a_..., b_..., res_...
      implicit none

c     Includes
      include 'COMMON.GEO'
      include 'COMMON.LOCMOVE'

c     Local variables
      integer i,j


      write(6,*)'a_...'
      write(6,'(8f7.1)')(a_ang(i)*rad2deg,i=0,a_n-1)
      write(6,'(8(2x,3l1,2x))')((a_tab(i,j),i=0,2),j=0,a_n-1)

      write(6,*)'b_...'
      write(6,'(4f7.1)')(b_ang(i)*rad2deg,i=0,b_n-1)
      write(6,'(4(2x,3l1,2x))')((b_tab(i,j),i=0,2),j=0,b_n-1)

      write(6,*)'res_...'
      write(6,'(12f7.1)')(res_ang(i)*rad2deg,i=0,res_n-1)
      write(6,'(12(2x,3l1,2x))')((res_tab(0,i,j),i=0,2),j=0,res_n-1)
      write(6,'(12(2x,3l1,2x))')((res_tab(1,i,j),i=0,2),j=0,res_n-1)
      write(6,'(12(2x,3l1,2x))')((res_tab(2,i,j),i=0,2),j=0,res_n-1)

      return
      end

c-------------------------------------------------------------

      subroutine angles2tab(PHImin,PHImax,n,ang,tab)
c     Only uses angles if [0,PI] (but PHImin cannot be 0.,
c     and PHImax cannot be PI)
      implicit none

c     Includes
      include 'COMMON.GEO'

c     INPUT arguments
      double precision PHImin,PHImax

c     OUTPUT arguments
      integer n
      double precision ang(0:3)
      logical tab(0:2,0:3)


      if (PHImin .eq. PHImax) then
c     Special case with two 010's
        n = 2;
        ang(0) = -PHImin;
        ang(1) = PHImin;
        tab(0,0) = .false.
        tab(2,0) = .false.
        tab(0,1) = .false.
        tab(2,1) = .false.
        tab(1,0) = .true.
        tab(1,1) = .true.
      else if (PHImin .eq. PI) then
c     Special case with one 010
        n = 1
        ang(0) = PI
        tab(0,0) = .false.
        tab(2,0) = .false.
        tab(1,0) = .true.
      else if (PHImax .eq. 0.) then
c     Special case with one 010
        n = 1
        ang(0) = 0.
        tab(0,0) = .false.
        tab(2,0) = .false.
        tab(1,0) = .true.
      else
c     Standard cases
        n = 0
        if (PHImin .gt. 0.) then
c     Start of range (011)
          ang(n) = PHImin
          tab(0,n) = .false.
          tab(1,n) = .true.
          tab(2,n) = .true.
c     End of range (110)
          ang(n+1) = -PHImin
          tab(0,n+1) = .true.
          tab(1,n+1) = .true.
          tab(2,n+1) = .false.
          n = n+2
        endif
        if (PHImax .lt. PI) then
c     Start of range (011)
          ang(n) = -PHImax
          tab(0,n) = .false.
          tab(1,n) = .true.
          tab(2,n) = .true.
c     End of range (110)
          ang(n+1) = PHImax
          tab(0,n+1) = .true.
          tab(1,n+1) = .true.
          tab(2,n+1) = .false.
          n = n+2
        endif
      endif

      return
      end

c-------------------------------------------------------------

      subroutine minmax_angles(x,y,z,r,n,ang,tab)
c     When solutions do not exist, assume all angles
c     are acceptable - i.e., initial geometry must be correct
      implicit none

c     Includes
      include 'COMMON.GEO'
      include 'COMMON.LOCMOVE'

c     Input arguments
      double precision x,y,z,r

c     Output arguments
      integer n
      double precision ang(0:3)
      logical tab(0:2,0:3)

c     Local variables
      double precision num, denom, phi
      double precision Kmin, Kmax
      integer i


      num = x*x + y*y + z*z
      denom = x*x + y*y
      n = 0
      if (denom .gt. 0.) then
        phi = atan2(y,x)
        denom = 2.*r*sqrt(denom)
        num = num+r*r
        Kmin = (num - dmin2)/denom
        Kmax = (num - dmax2)/denom

c     Allowed values of K (else all angles are acceptable)
c     -1 <= Kmin <  1
c     -1 <  Kmax <= 1
        if (Kmin .gt. 1. .or. abs(Kmin-1.) .lt. small2) then
          Kmin = -flag
        else if (Kmin .lt. -1. .or. abs(Kmin+1.) .lt. small2) then
          Kmin = PI
        else
          Kmin = acos(Kmin)
        endif

        if (Kmax .lt. -1. .or. abs(Kmax+1.) .lt. small2) then
          Kmax = flag
        else if (Kmax .gt. 1. .or. abs(Kmax-1.) .lt. small2) then
          Kmax = 0.
        else
          Kmax = acos(Kmax)
        endif

        if (Kmax .lt. Kmin) Kmax = Kmin

        call angles2tab(Kmin, Kmax, n, ang, tab)

c     Add phi and check that angles are within range (-PI,PI]
        do i=0,n-1
          ang(i) = ang(i)+phi
          if (ang(i) .le. -PI) then
            ang(i) = ang(i)+2.*PI
          else if (ang(i) .gt. PI) then
            ang(i) = ang(i)-2.*PI
          endif
        enddo
      endif

      return
      end

c-------------------------------------------------------------

      subroutine construct_tab
c     Take a_... and b_... values and produces the results res_...
c     x_ang are assumed to be all different (diff > small)
c     x_tab(1,i) must be 1 for all i (i.e., all x_ang are acceptable)
      implicit none

c     Includes
      include 'COMMON.LOCMOVE'

c     Local variables
      integer n_max,i,j,index
      logical done
      double precision phi


      n_max = a_n + b_n
      if (n_max .eq. 0) then
        res_n = 0
        return
      endif

      do i=0,n_max-1
        do j=0,1
          res_tab(j,0,i) = .true.
          res_tab(j,2,i) = .true.
          res_tab(j,1,i) = .false.
        enddo
      enddo

      index = 0
      phi = -flag
      done = .false.
      do while (.not.done)
        res_ang(index) = flag

c     Check a first...
        do i=0,a_n-1
          if ((a_ang(i)-phi).gt.small .and.
     +         a_ang(i) .lt. res_ang(index)) then
c     Found a lower angle
            res_ang(index) = a_ang(i)
c     Copy the values from a_tab into res_tab(0,,)
            res_tab(0,0,index) = a_tab(0,i)
            res_tab(0,1,index) = a_tab(1,i)
            res_tab(0,2,index) = a_tab(2,i)
c     Set default values for res_tab(1,,)
            res_tab(1,0,index) = .true.
            res_tab(1,1,index) = .false.
            res_tab(1,2,index) = .true.
          else if (abs(a_ang(i)-res_ang(index)).lt.small) then
c     Found an equal angle (can only be equal to a b_ang)
            res_tab(0,0,index) = a_tab(0,i)
            res_tab(0,1,index) = a_tab(1,i)
            res_tab(0,2,index) = a_tab(2,i)
          endif
        enddo
c     ...then check b
        do i=0,b_n-1
          if ((b_ang(i)-phi).gt.small .and.
     +         b_ang(i) .lt. res_ang(index)) then
c     Found a lower angle
            res_ang(index) = b_ang(i)
c     Copy the values from b_tab into res_tab(1,,)
            res_tab(1,0,index) = b_tab(0,i)
            res_tab(1,1,index) = b_tab(1,i)
            res_tab(1,2,index) = b_tab(2,i)
c     Set default values for res_tab(0,,)
            res_tab(0,0,index) = .true.
            res_tab(0,1,index) = .false.
            res_tab(0,2,index) = .true.
          else if (abs(b_ang(i)-res_ang(index)).lt.small) then
c     Found an equal angle (can only be equal to an a_ang)
            res_tab(1,0,index) = b_tab(0,i)
            res_tab(1,1,index) = b_tab(1,i)
            res_tab(1,2,index) = b_tab(2,i)
          endif
        enddo

        if (res_ang(index) .eq. flag) then
          res_n = index
          done = .true.
        else if (index .eq. n_max-1) then
          res_n = n_max
          done = .true.
        else
          phi = res_ang(index)  ! Store previous angle
          index = index+1
        endif
      enddo

c     Fill the gaps
c     First a...
      index = 0
      if (a_n .gt. 0) then
        do while (.not.res_tab(0,1,index))
          index=index+1
        enddo
        done = res_tab(0,2,index)
        do i=index+1,res_n-1
          if (res_tab(0,1,i)) then
            done = res_tab(0,2,i)
          else
            res_tab(0,0,i) = done
            res_tab(0,1,i) = done
            res_tab(0,2,i) = done
          endif
        enddo
        done = res_tab(0,0,index)
        do i=index-1,0,-1
          if (res_tab(0,1,i)) then
            done = res_tab(0,0,i)
          else
            res_tab(0,0,i) = done
            res_tab(0,1,i) = done
            res_tab(0,2,i) = done
          endif
        enddo
      else
        do i=0,res_n-1
          res_tab(0,0,i) = .true.
          res_tab(0,1,i) = .true.
          res_tab(0,2,i) = .true.
        enddo
      endif
c     ...then b
      index = 0
      if (b_n .gt. 0) then
        do while (.not.res_tab(1,1,index))
          index=index+1
        enddo
        done = res_tab(1,2,index)
        do i=index+1,res_n-1
          if (res_tab(1,1,i)) then
            done = res_tab(1,2,i)
          else
            res_tab(1,0,i) = done
            res_tab(1,1,i) = done
            res_tab(1,2,i) = done
          endif
        enddo
        done = res_tab(1,0,index)
        do i=index-1,0,-1
          if (res_tab(1,1,i)) then
            done = res_tab(1,0,i)
          else
            res_tab(1,0,i) = done
            res_tab(1,1,i) = done
            res_tab(1,2,i) = done
          endif
        enddo
      else
        do i=0,res_n-1
          res_tab(1,0,i) = .true.
          res_tab(1,1,i) = .true.
          res_tab(1,2,i) = .true.
        enddo
      endif

c     Finally fill the last row with AND operation
      do i=0,res_n-1
        do j=0,2
          res_tab(2,j,i) = (res_tab(0,j,i) .and. res_tab(1,j,i))
        enddo
      enddo

      return 
      end

c-------------------------------------------------------------

      subroutine construct_ranges(phi_n,phi_start,phi_end)
c     Given the data in res_..., construct a table of 
c     min/max allowed angles
      implicit none

c     Includes
      include 'COMMON.GEO'
      include 'COMMON.LOCMOVE'

c     Output arguments
      integer phi_n
      double precision phi_start(0:11),phi_end(0:11)

c     Local variables
      logical done
      integer index


      if (res_n .eq. 0) then
c     Any move is allowed
        phi_n = 1
        phi_start(0) = -PI
        phi_end(0) = PI
      else
        phi_n = 0
        index = 0
        done = .false.
        do while (.not.done)
c     Find start of range (01x)
          done = .false.
          do while (.not.done)
            if (res_tab(2,0,index).or.(.not.res_tab(2,1,index))) then
              index=index+1
            else
              done = .true.
              phi_start(phi_n) = res_ang(index)
            endif
            if (index .eq. res_n) done = .true.
          enddo
c     If a start was found (index < res_n), find the end of range (x10)
c     It may not be found without wrapping around
          if (index .lt. res_n) then
            done = .false.
            do while (.not.done)
              if ((.not.res_tab(2,1,index)).or.res_tab(2,2,index)) then
                index=index+1
              else
                done = .true.
              endif
              if (index .eq. res_n) done = .true.
            enddo
            if (index .lt. res_n) then
c     Found the end of the range
              phi_end(phi_n) = res_ang(index)
              phi_n=phi_n+1
              index=index+1
              if (index .eq. res_n) then
                done = .true.
              else
                done = .false.
              endif
            else
c     Need to wrap around
              done = .true.
              phi_end(phi_n) = flag
            endif
          endif
        enddo
c     Take care of the last one if need to wrap around
        if (phi_end(phi_n) .eq. flag) then
          index = 0
          do while ((.not.res_tab(2,1,index)).or.res_tab(2,2,index))
            index=index+1
          enddo
          phi_end(phi_n) = res_ang(index) + 2.*PI
          phi_n=phi_n+1
        endif
      endif

      return
      end

c-------------------------------------------------------------

      subroutine fix_no_moves(phi)
      implicit none

c     Includes
      include 'COMMON.GEO'
      include 'COMMON.LOCMOVE'

c     Output arguments
      double precision phi

c     Local variables
      integer index
      double precision diff,temp


c     Look for first 01x in gammas (there MUST be at least one)
      diff = flag
      index = 0
      do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
        index=index+1
      enddo
      if (res_ang(index) .le. 0.D0) then ! Make sure it's from PHImax
c     Try to increase PHImax
        if (index .gt. 0) then
          phi = res_ang(index-1)
          diff = abs(res_ang(index) - res_ang(index-1))
        endif
c     Look for last (corresponding) x10
        index = res_n - 1
        do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
          index=index-1
        enddo
        if (index .lt. res_n-1) then
          temp = abs(res_ang(index) - res_ang(index+1))
          if (temp .lt. diff) then
            phi = res_ang(index+1)
            diff = temp
          endif
        endif
      endif

c     If increasing PHImax didn't work, decreasing PHImin
c     will (with one exception)
c     Look for first x10 (there MUST be at least one)
      index = 0
      do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
        index=index+1
      enddo
      if (res_ang(index) .lt. 0.D0) then ! Make sure it's from PHImin
c     Try to decrease PHImin
        if (index .lt. res_n-1) then
          temp = abs(res_ang(index) - res_ang(index+1))
          if (res_ang(index+1) .le. 0.D0 .and. temp .lt. diff) then
            phi = res_ang(index+1)
            diff = temp
          endif
        endif
c     Look for last (corresponding) 01x
        index = res_n - 1
        do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
          index=index-1
        enddo
        if (index .gt. 0) then
          temp = abs(res_ang(index) - res_ang(index-1))
          if (res_ang(index-1) .ge. 0.D0 .and. temp .lt. diff) then
            phi = res_ang(index-1)
            diff = temp
          endif
        endif
      endif

c     If it still didn't work, it must be PHImax == 0. or PHImin == PI
      if (diff .eq. flag) then
        index = 0
        if (res_tab(index,1,0) .or. (.not.res_tab(index,1,1)) .or.
     +       res_tab(index,1,2)) index = res_n - 1
c     This MUST work at this point
        if (index .eq. 0) then
          phi = res_ang(1)
        else
          phi = res_ang(index - 1)
        endif
      endif

      return
      end

c-------------------------------------------------------------

      integer function move_res(PHImin,PHImax,i_move)
c     Moves residue i_move (in array c), leaving everything else fixed
c     Starting geometry is not checked, it should be correct!
c     R(,i_move) is the only residue that will move, but must have
c     1 < i_move < nres (i.e., cannot move ends)
c     Whether any output is done is controlled by locmove_output
crc      implicit none

c     Includes
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CHAIN'
      include 'COMMON.GEO'
      include 'COMMON.LOCMOVE'

c     External functions
      double precision ran_number
      external ran_number

c     Input arguments
      double precision PHImin,PHImax
      integer i_move

c     RETURN VALUES:
c     0: move successfull
c     1: Dmin or Dmax had to be modified
c     2: move failed - check your input geometry


c     Local variables
      double precision X(0:2),Y(0:2),Z(0:2),Orig(0:2)
      double precision P(0:2)
      logical no_moves,done
      integer index,i,j
      double precision phi,temp,radius
      double precision phi_start(0:11), phi_end(0:11)
      integer phi_n

c     Set up the coordinate system
      do i=0,2
        Orig(i)=0.5*(c(i+1,i_move-1)+c(i+1,i_move+1)) ! Position of origin
      enddo

      do i=0,2
        Z(i)=c(i+1,i_move+1)-c(i+1,i_move-1)
      enddo
      temp=sqrt(Z(0)*Z(0)+Z(1)*Z(1)+Z(2)*Z(2))
      do i=0,2
        Z(i)=Z(i)/temp
      enddo

      do i=0,2
        X(i)=c(i+1,i_move)-Orig(i)
      enddo
c     radius is the radius of the circle on which c(,i_move) can move
      radius=sqrt(X(0)*X(0)+X(1)*X(1)+X(2)*X(2))
      do i=0,2
        X(i)=X(i)/radius
      enddo

      Y(0)=Z(1)*X(2)-X(1)*Z(2)
      Y(1)=X(0)*Z(2)-Z(0)*X(2)
      Y(2)=Z(0)*X(1)-X(0)*Z(1)

c     Calculate min, max angles coming from dmin, dmax to c(,i_move-2)
      if (i_move.gt.2) then
        do i=0,2
          P(i)=c(i+1,i_move-2)-Orig(i)
        enddo
        call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),
     +       P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),
     +       P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),
     +       radius,a_n,a_ang,a_tab)
      else
        a_n=0
      endif

c     Calculate min, max angles coming from dmin, dmax to c(,i_move+2)
      if (i_move.lt.nres-2) then
        do i=0,2
          P(i)=c(i+1,i_move+2)-Orig(i)
        enddo
        call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),
     +       P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),
     +       P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),
     +       radius,b_n,b_ang,b_tab)
      else
        b_n=0
      endif

c     Construct the resulting table for alpha and beta
      call construct_tab()

      if (locmove_output) then
        print *,'ALPHAS & BETAS TABLE'
        call output_tabs()
      endif

c     Check that there is at least one possible move
      no_moves = .true.
      if (res_n .eq. 0) then
        no_moves = .false.
      else
        index = 0
        do while ((index .lt. res_n) .and. no_moves)
          if (res_tab(2,1,index)) no_moves = .false.
          index=index+1
        enddo
      endif
      if (no_moves) then
        if (locmove_output) print *,'   ***   Cannot move anywhere'
        move_res=2
        return
      endif

c     Transfer res_... into a_...
      a_n = 0
      do i=0,res_n-1
        if ( (res_tab(2,0,i).neqv.res_tab(2,1,i)) .or.
     +       (res_tab(2,0,i).neqv.res_tab(2,2,i)) ) then
          a_ang(a_n) = res_ang(i)
          do j=0,2
            a_tab(j,a_n) = res_tab(2,j,i)
          enddo
          a_n=a_n+1
        endif
      enddo

c     Check that the PHI's are within [0,PI]
      if (PHImin .lt. 0. .or. abs(PHImin) .lt. small) PHImin = -flag
      if (PHImin .gt. PI .or. abs(PHImin-PI) .lt. small) PHImin = PI
      if (PHImax .gt. PI .or. abs(PHImax-PI) .lt. small) PHImax = flag
      if (PHImax .lt. 0. .or. abs(PHImax) .lt. small) PHImax = 0.
      if (PHImax .lt. PHImin) PHImax = PHImin
c     Calculate min and max angles coming from PHImin and PHImax,
c     and put them in b_...
      call angles2tab(PHImin, PHImax, b_n, b_ang, b_tab)
c     Construct the final table
      call construct_tab()

      if (locmove_output) then
        print *,'FINAL TABLE'
        call output_tabs()
      endif

c     Check that there is at least one possible move
      no_moves = .true.
      if (res_n .eq. 0) then
        no_moves = .false.
      else
        index = 0
        do while ((index .lt. res_n) .and. no_moves)
          if (res_tab(2,1,index)) no_moves = .false.
          index=index+1
        enddo
      endif

      if (no_moves) then
c     Take care of the case where no solution exists...
        call fix_no_moves(phi)
        if (locmove_output) then
          print *,'   ***   Had to modify PHImin or PHImax'
          print *,'phi: ',phi*rad2deg
        endif
        move_res=1
      else
c     ...or calculate the solution
c     Construct phi_start/phi_end arrays
        call construct_ranges(phi_n, phi_start, phi_end)
c     Choose random angle phi in allowed range(s)
        temp = 0.
        do i=0,phi_n-1
          temp = temp + phi_end(i) - phi_start(i)
        enddo
        phi = ran_number(phi_start(0),phi_start(0)+temp)
        index = 0
        done = .false.
        do while (.not.done)
          if (phi .lt. phi_end(index)) then
            done = .true.
          else
            index=index+1
          endif
          if (index .eq. phi_n) then
            done = .true.
          else if (.not.done) then
            phi = phi + phi_start(index) - phi_end(index-1)
          endif
        enddo
        if (index.eq.phi_n) phi=phi_end(phi_n-1) ! Fix numerical errors
        if (phi .gt. PI) phi = phi-2.*PI

        if (locmove_output) then
          print *,'ALLOWED RANGE(S)'
          do i=0,phi_n-1
            print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
          enddo
          print *,'phi: ',phi*rad2deg
        endif
        move_res=0
      endif

c     Re-use radius as temp variable
      temp=radius*cos(phi)
      radius=radius*sin(phi)
      do i=0,2
        c(i+1,i_move)=Orig(i)+temp*X(i)+radius*Y(i)
      enddo

      return
      end

c-------------------------------------------------------------

      subroutine loc_test
crc      implicit none

c     Includes
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.LOCMOVE'

c     External functions
      integer move_res
      external move_res

c     Local variables
      integer i,j
      integer phi_n
      double precision phi_start(0:11),phi_end(0:11)
      double precision phi
      double precision R(0:2,0:5)

      locmove_output=.true.

c      call angles2tab(30.*deg2rad,70.*deg2rad,a_n,a_ang,a_tab)
c      call angles2tab(80.*deg2rad,130.*deg2rad,b_n,b_ang,b_tab)
c      call minmax_angles(0.D0,3.8D0,0.D0,3.8D0,b_n,b_ang,b_tab)
c      call construct_tab
c      call output_tabs

c      call construct_ranges(phi_n,phi_start,phi_end)
c      do i=0,phi_n-1
c        print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
c      enddo

c      call fix_no_moves(phi)
c      print *,'NO MOVES FOUND, BEST PHI IS',phi*rad2deg

      R(0,0)=0.D0
      R(1,0)=0.D0
      R(2,0)=0.D0
      R(0,1)=0.D0
      R(1,1)=-cos(28.D0*deg2rad)
      R(2,1)=-0.5D0-sin(28.D0*deg2rad)
      R(0,2)=0.D0
      R(1,2)=0.D0
      R(2,2)=-0.5D0
      R(0,3)=cos(30.D0*deg2rad)
      R(1,3)=0.D0
      R(2,3)=0.D0
      R(0,4)=0.D0
      R(1,4)=0.D0
      R(2,4)=0.5D0
      R(0,5)=0.D0
      R(1,5)=cos(26.D0*deg2rad)
      R(2,5)=0.5D0+sin(26.D0*deg2rad)
      do i=1,5
        do j=0,2
          R(j,i)=vbl*R(j,i)
        enddo
      enddo
      i=move_res(R(0,1),0.D0*deg2rad,180.D0*deg2rad)
      print *,'RETURNED ',i
      print *,(R(i,3)/vbl,i=0,2)

      return
      end

c-------------------------------------------------------------