anchor.diff:      if(x < zero) then ! Crystal hit from below
   !      write(*,*) 'Crystal potentially hit from below'
   !      write(*,*) 'CHECK', x
  -      xp_tangent = sqrt((-(two*x)*c_rcurv + x**2)/c_rcurv**2)
  +      xp_tangent = sqrt((-(two*x)*c_rcurv + x**2)/(c_rcurv**2))
   !      write(*,*) 'tangent', xp_tangent, xp
       else ! Crystal hit from above
   !      write(*,*) 'Crystal potentially hit from above'
  
positive.code:       alpha_F = (c_rcurv-x_P)/x_P
        beta_F = -(r**2-s_P**2-x_P**2)/(two*s_P)
        A_F = alpha_F**2 + one
        B_F = two*alpha_F*beta_F - two*c_rcurv
        C_F = beta_F**2
      else
        ! Intersect with top side
        alpha_F = (c_rcurv-x_P)/s_P
        beta_F = -(r**2+c_xmax*(c_xmax-two*c_rcurv)-s_P**2-x_P**2)/(two*s_P)
        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
  
negative.code: !  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
  