anchor.diff:  
       Lrefl = xp_rel*r ! Distance of refl. point [m]
       Srefl = sin_mb(xp_rel/two + c_miscut)*Lrefl
  -!    write(*,*) 'Lrefl', Lrefl
  -!    write(*,*) 'L_chan', L_chan
   
       if(Lrefl > zero .and. Lrefl < L_chan) then ! VR point inside
   
  
positive.code:       A_F = alpha_F**2 + one
        B_F = two*(alpha_F*beta_F) - two*c_rcurv
        C_F = beta_F**2 - c_xmax*(c_xmax-two*c_rcurv)
  !      write(*,*) alpha_F, beta_F, A_F, B_F, C_F
      endif
      
      x_F = (-B_F-sqrt(B_F**2-four*(A_F*C_F)))/(two*A_F)
      s_F = alpha_F*x_F + beta_F
  !    write(*,*) B_F**2, A_F*C_F, four*A_F*C_F, B_F**2-four*A_F*C_F
  !    write(*,*) 'new F', x_F, s_F
      
    endif
  
    ! MISCUT fourth step: deflection and length
    a = sqrt(s_F**2+(x-x_F)**2)
    tP = acos_mb((2*(r**2)-a**2)/(2*(r**2)))
    tdefl = asin_mb((s_f-s_P)/r)
    L_chan = r*tP
  !  write(*,*) "Results", tdefl, L_chan
  
    xp_rel = xp - c_miscut
  
  
negative.code: 
    else
  
      if(x < zero) then ! Crystal can be hit from below
        xp_tangent = sqrt((-(two*x)*c_rcurv + x**2)/(c_rcurv**2))
      else              ! Crystal can be hit from above
        xp_tangent = asin_mb((c_rcurv*(one - c_cBend) - x)/sqrt(((two*c_rcurv)*(c_rcurv - x))*(one - c_cBend) + 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
  