anchor.diff:        s_int = (-b_eq - sqrt(delta))/(two*a_eq)
         s_imp = s_int
   
  -      ! MISCUT first step: P coordinates
  +      ! MISCUT first step: P coordinates (center of curvature of crystalline planes)
         s_P_tmp = (c_rcurv-c_xmax)*sin_mb(-c_miscut)
         x_P_tmp = c_xmax + (c_rcurv-c_xmax)*cos_mb(-c_miscut)
  -!      write(*,*) "P side tmp", s_P_tmp, x_P_tmp
   
         if(s_int < c_rcurv*c_sBend) then
           ! Transform to a new reference system: shift and rotate
  
positive.code:   s_length = c_rcurv*sin_mb(length/c_rcurv)
    L_chan   = length
  
    ! MISCUT second step: fundamental coordinates
  !  s_K = c_rcurv*sin_mb(cry_bend)
  !  x_K = c_rcurv*(1-cos_mb(cry_bend))
  !  s_M = (c_rcurv-c_xmax)*sin_mb(cry_bend)
  !  x_M = c_xmax + (c_rcurv-c_xmax)*(1-cos_mb(cry_bend))
    s_K = c_rcurv*sin_mb(length/c_rcurv)
    x_K = c_rcurv*(1-cos_mb(length/c_rcurv))
    s_M = (c_rcurv-c_xmax)*sin_mb(length/c_rcurv)
    x_M = c_xmax + (c_rcurv-c_xmax)*(1-cos_mb(length/c_rcurv))
    r   = sqrt(s_P**2 + (x-x_P)**2)
  !  write(*,*) "bending var", cry_bend
  !  write(*,*) "bending cal", length/c_rcurv
  !  write(*,*) "Rcurv", c_rcurv
  !  write(*,*) "P", s_P, x_P
  !  write(*,*) "I", s, x
  !  write(*,*) "L", s, c_xmax
  !  write(*,*) "K", s_K, x_K
  !  write(*,*) "M", s_M, x_M
  !  write(*,*) "r", r
  
    ! MISCUT third step: F coordinates on crystal exit face
  !  A_F = tan_mb(cry_bend)**2 + one
  !  B_F = -two*tan_mb(cry_bend)**2*c_rcurv + two*tan_mb(cry_bend)*s_P -two*x_P
  !  C_F = tan_mb(cry_bend)**2*c_rcurv**2 - two*tan_mb(cry_bend)*s_P*c_rcurv + s_P**2 + x_P**2 - r**2
    A_F = (tan_mb(length/c_rcurv))**2 + one
    B_F = ((-two)*(tan_mb(length/c_rcurv))**2)*c_rcurv + (two*tan_mb(length/c_rcurv))*s_P - two*x_P
    C_F = ((tan_mb(length/c_rcurv))**2)*(c_rcurv**2) - ((two*tan_mb(length/c_rcurv))*s_P)*c_rcurv + s_P**2 + x_P**2 - r**2
  
    x_F = (-B_F-sqrt(B_F**2-four*(A_F*C_F)))/(two*A_F)
  !  s_F = -tan_mb(cry_bend)*(x_F-c_rcurv)
    s_F = (-tan_mb(length/c_rcurv))*(x_F-c_rcurv)
  !  write(*,*) "F", s_F, x_F
  !  write(*,*) "F-K", s_F-s_K, x_F-x_K
  !  write(*,*) "F-M", s_F-s_M, x_F-x_M
  
    if(x_F >= x_K .and. x_F <= x_M .and. s_F >= s_M .and. s_F <= s_K) then
  !    write(*,*) "F on exit face"
      ! MISCUT No additional steps required
    else if (c_miscut == 0 .and. abs(x_F-x_K) <= c1m13 .and. abs(s_F-s_K) <= c1m13) then
      ! no miscut, entrance from below: exit point is K
  !    write(*,*) "F on exit face"
      x_F = x_K
      s_F = s_K
    else if (c_miscut == 0 .and. abs(x_F-x_M) <= c1m13 .and. abs(s_F-s_M) <= c1m13) then
      ! no miscut, entrance from above: exit point is M
  !    write(*,*) "F on exit face"
      x_F = x_M
      s_F = s_M
    else
  !    write(*,*) "F NOT on exit face"
  
      ! MISCUT Third step bis: F coordinates on bent side
      if(c_miscut < 0) then
        ! Intersect with bottom side
        alpha_F = (c_rcurv-x_P)/x_P
  
negative.code: 
      Lrefl = xp_rel*r ! Distance of refl. point [m]
      Srefl = sin_mb(xp_rel/two + c_miscut)*Lrefl
  
      if(Lrefl > zero .and. Lrefl < L_chan) then ! VR point inside
  
  