anchor.diff:    z  = z - zp*s_rot
     s  = zero
   
  -!  write(*,*) 'x after reference change', x, xp
  -
     ! Check that particle hit the crystal
     if(x >= zero .and. x < c_xmax) then
   
  -!    write(*,*) "front face hit"
  -    ! MISCUT first step: P coordinates
  +    ! MISCUT first step: P coordinates (center of curvature of crystalline planes)
       s_P = (c_rcurv-c_xmax)*sin_mb(-c_miscut)
       x_P = c_xmax + (c_rcurv-c_xmax)*cos_mb(-c_miscut)
  -!    write(*,*) "P front", s_P, x_P
   
       call cry_interact(mat,x,xp,z,zp,p,cry_length,s_P,x_P)
       s   = c_rcurv*c_sBend
  
positive.code:           x = x + (half*xp)*(s_length - Srefl)
            y = y + (half*yp)*(s_length - Srefl)
          else
  !          write(*,*) 'TRAM'
            iProc = proc_TRAM
            x = x + xp*Srefl
            y = y + yp*Srefl
  
negative.code:       s_int = (-b_eq - sqrt(delta))/(two*a_eq)
        s_imp = s_int
  
        ! 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)
  
        if(s_int < c_rcurv*c_sBend) then
          ! Transform to a new reference system: shift and rotate
  