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
   
  
positive.code:           x     = x + xp*Srefl
            y     = y + yp*Srefl
  
            Dxp = (((-three*Ang_rms)*xp_rel)/(two*xpcrit) + Ang_avr) + ((three*Ang_rms)*(L_chan/c_rcurv))/(two*xpcrit)
            xp  = xp + Dxp
            x   = x + (half*xp)*(s_length-Srefl)
            y   = y + (half*yp)*(s_length-Srefl)
  
negative.code:   else
  
      if(x < zero) then ! Crystal hit from below
        write(*,*) 'CHECK', x
        xp_tangent = sqrt((-(two*x)*c_rcurv + x**2)/c_rcurv**2)
        write(*,*) 'tangent', xp_tangent, xp
      else ! Crystal hit from above
        write(*,*) 'CHECK', x
        xp_tangent = asin_mb((c_rcurv*(one - c_cBend) - x)/sqrt(((two*c_rcurv)*(c_rcurv - x))*(one - c_cBend) + x**2))
        write(*,*) 'tangent', xp_tangent, xp
      end if
  
      ! If the hit is below, the angle must be greater or equal than the tangent,
  