anchor.diff:    else
   
       if(x < zero) then ! Crystal hit from below
  +!      write(*,*) 'Crystal hit from below'
   !      write(*,*) 'CHECK', x
         xp_tangent = sqrt((-(two*x)*c_rcurv + x**2)/c_rcurv**2)
   !      write(*,*) 'tangent', xp_tangent, xp
  
positive.code: 
      else
  
        s = c_rcurv*sin_mb(cry_length/c_rcurv)
        x = x + s*xp
        z = z + s*zp
  
negative.code:         s      = zero
  
          tilt_int = s_int/c_rcurv
  !        write(*,*) "tilt_int", tilt_int
          xp       = xp-tilt_int
  
          ! MISCUT first step (bis): transform P in new reference system
          ! Translation
          s_P_tmp = s_P_tmp - s_int
          x_P_tmp = x_P_tmp - x_int
  !        write(*,*) "P side after translation", s_P_tmp, x_P_tmp
          ! Rotation
          s_P = s_P_tmp*cos_mb(tilt_int) + x_P_tmp*sin_mb(tilt_int)
          x_P = -s_P_tmp*sin_mb(tilt_int) + x_P_tmp*cos_mb(tilt_int)
  !        write(*,*) "P side after rotation", s_P, x_P
  
          call cry_interact(mat,x,xp,z,zp,p,cry_length-(tilt_int*c_rcurv),s_P,x_P)
          s   = c_rcurv*sin_mb(cry_bend - tilt_int)
  