anchor.diff:        end if
   
         ! Transform back to particle coordinates with opening and offset
  -      x   = x + c_aperture/2 + mirror*c_offset
  +      x  = x + c_aperture/2 + mirror*c_offset
   
         ! Now mirror at the horizontal axis for negative X offset
  -      x   = mirror * x
  -      xp  = mirror * xp
  +      x  = mirror*x
  +      xp = mirror*xp
   
         ! Last do rotation into collimator frame
         x_in(j)  =  x*cRRot +  z*sRRot
  
positive.code:     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)
  
      ! CRY ---------------------
      dpop   = (p-p0)/p0
  
negative.code:   else
  
      if(x < zero) then ! Crystal hit from below
        xp_tangent = sqrt((-(two*x)*c_rcurv + x**2)/(c_rcurv**2))
      else ! Crystal hit from above
        xp_tangent = asin_mb((c_rcurv*(one-cos_mb(cry_bend))-x)/sqrt(two*c_rcurv*(c_rcurv-x)*(one-cos_mb(cry_bend))+x**2))
      end if
  
      ! If the hit is below, the angle must be greater or equal than the tangent,
      ! or if the hit is above, the angle must be smaller or equal than the tangent
      if((x < zero .and. xp >= xp_tangent) .or. (x >= zero .and. 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 reference 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
              x0  = x_shift + shift
              xp0 = xp_shift
            else
              x0  = x_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
  
  