anchor.diff:      lint(j)   = -one
       indiv(j)  = -one
   
  -    ! CRY ---------------------
  -    nabs = 0
  -    s    = zero
  -    ! CRY ---------------------
  -
  -    x   = x_in(j)
  -    xp  = xp_in(j)
  -    z   = y_in(j)
  -    zp  = yp_in(j)
  -    p   = p_in(j)
  +    x  = x_in(j)
  +    xp = xp_in(j)
  +    z  = y_in(j)
  +    zp = yp_in(j)
  +    p  = p_in(j)
   
       ! CRY ---------------------
       dpop   = (p-p0)/p0
  
positive.code:   else
  
      if(x < zero) then ! Crystal hit from below
        XP_tangent = sqrt((-(two*x)*c_rcurv + x**2)/(c_rcurv**2))
        if(xp >= XP_tangent) then
          ! if it hits the crystal, calculate in which point and apply the
          ! transformation and drift to that point
          a_eq  = (one + xp**2)
          b_eq  = (two*xp)*(x-c_rcurv)
          c_eq  = -(two*x)*c_rcurv + x**2
          Delta = b_eq**2 - four*(a_eq*c_eq)
          S_int = (-b_eq - sqrt(Delta))/(two*a_eq)
  
          if(S_int < c_rcurv*sin_mb(cry_bend)) then
            !  transform to a new ref system:shift and rotate
            x_int  = xp*s_int + x
            xp_int = xp
            z      = z + zp*s_int
            x      = zero
            s      = zero
  
            tilt_int = s_int/c_rcurv
            xp       = xp-tilt_int
  
            call cryst(mat,x,xp,z,zp,p,(cry_length-(tilt_int*c_rcurv)),j)
            s   = c_rcurv*sin_mb(cry_bend-tilt_int)
            zlm = c_rcurv*sin_mb(cry_bend-tilt_int)
            if(iProc /= proc_out) then
              x_rot    = x_int
              s_rot    = s_int
              xp_rot   = xp_int
              s_shift  =  s_rot*c_cnTilt + x_rot*c_snTilt
              x_shift  = -s_rot*c_snTilt + x_rot*c_cnTilt
              xp_shift = xp_rot + cry_tilt
  
              if(cry_tilt < zero) then
                s_impact  = s_shift
                x0  = x_shift+shift
                xp0 = xp_shift
              else
                x0  = x_shift
                s_impact  = s_shift
                xp0 = xp_shift
              end if
  
              nhit = nhit + 1
              lhit(j) = ie
              lhit_turn(j) = iturn
              impact(j) = x0
              indiv(j) = xp0
            end if
  
            ! un-rotate
            x_temp  = x
            s_temp  = s
            xp_temp = xp
            s       =  s_temp*cos_mb(-tilt_int) + x_temp*sin_mb(-tilt_int)
            x       = -s_temp*sin_mb(-tilt_int) + x_temp*cos_mb(-tilt_int)
            xp      = xp_temp + tilt_int
  
            ! 2nd: shift back the 2 axis
            x = x+x_int
            s = s+s_int
  
          else
  
            s = c_rcurv*sin_mb(cry_length/c_rcurv)
            x = x + s*xp
            z = z + s*zp
  
          end if
  
        else
  
  
negative.code: 
        end if
  
      else
  
        s = c_rcurv*sin_mb(cry_length/c_rcurv)
        x = x + s*xp
        z = z + s*zp
  
      end if
  
    end if
  
    ! transform back from the crystal to the collimator reference system
  