anchor.diff:    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
  +      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
   
  -          ! 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
  -
  +          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
   
  
positive.code:       end if
  
        ! Transform back to particle coordinates with opening and 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
  
        ! Last do rotation into collimator frame
        x_in(j)  =  x*cRRot +  z*sRRot
  
negative.code:   real(kind=fPrec), intent(in)    :: cry_length ! Original length (from the db) [m]
  
    integer j,mat,nabs,nhit
    real(kind=fPrec) p0,zlm,x,xp,z,zp,s,p,dpop,x_in0,xp_in0,y_in0,yp_in0,p_in0,s_in0,mirror, &
      tiltangle,cry_bend,cRot,sRot,cRRot,sRRot
  
    ! CRY ---------------------
  