Subroutine MFXTR(RTO,XF,PF,ANGF,IERR)
C   =======================================
C
C Author: G.Bruni
C Steering of a RK extrapolation to surface R=RTO
C
C
C    ierr=1 error in Runge Kutta
C         2 prediction not reached
C         3 Input momentum below cut (PinMin)
C
C
C       ps(1:3)  starting momentum components
C
C       xf       (x,y,z)-final
C       angf(1)  theta  -final
C            2   phi    -final
C       pf(1)    px     -final
C          2     py     -final
C          3     pz     -final
C
C
      Implicit None
C
#include "partap.inc"
#include "mfrtz.inc"
C
#include "mfflg.inc"
C
      Integer ierr
      Real xf(3),angf(2),pf(3),rto
C
      Logical first /.true./
      Save    first
      Integer i,c1,c2,iht,j,nsteps
      Real    pcell,Pi,x(3),ddx(3),xx(3),ddxs(3),Pin,Step,H,DH,ParStart,
     .        ParNow,ParTarg,dels,ddxe(3),xe(3),Dphi,TwoPi,Eloss,Phf,Thf
      Logical Ok
C
      Real ParEnd
      Integer Ient/4/
      Real ST     /0.10/
      Real PinMin /0.75/
      Real TarDif /0.20/
      Real dEdx  /1./
      Save ST,PinMin,Ient,iht,dEdx,TarDif
C
      if(first) then
        first=.false.
        Pi=ACOS(-1.)
        TwoPi=2.*Pi
        DPhi=TwoPi/60.
        if(Mdebug.ge.2) then
          write(6,*) ' ============================ '
          write(6,*) ' MFXTR  parameters             '
          write(6,*) ' ============================ '
          write(6,*) ' STEP   : ',st
          write(6,*) ' PinMin : ',PinMin
          write(6,*) ' TarDif : ',TarDif
          write(6,*) ' dE/dX  : ',dedx
          write(6,*) ' ============================ '
        endif
      endif
C
      ParEnd=Rto*0.01
C
      ddx(3)=1./Sqrt(1.+mfrtz_x(3)**2+mfrtz_x(4)**2)
      ddx(1)=mfrtz_x(3)*ddx(3)
      ddx(2)=mfrtz_x(4)*ddx(3)
      ddxs(1)=ddx(3)
      ddxs(2)=ddx(1)
      ddxs(3)=ddx(2)
      xx(1)=mfrtz_x(6)/100.
      xx(2)=mfrtz_x(1)/100.
      xx(3)=mfrtz_x(2)/100.
      Pin=1./mfrtz_x(5)
      if(Abs(Pin).lt.PinMin) then
        ierr=3
        goto 9900
      endif
      ParStart=Sqrt(xx(2)**2+xx(3)**2)
c      ParStart=xx(1)
      Step=ST
      H=ParEnd-ParStart
      if(Abs(H).lt.Step) then
       nsteps=1
      else
       nsteps=Abs(H)/Step+1
      endif
      DH=Sign(1.,H)*Step
C
      if(Mdebug.ge.2) then
        write(6,*)
        write(6,*) ' =====>  MFXTR <===== '
        write(6,*) ' ParStart    ',ParStart
        write(6,*) ' ParEnd      ',ParEnd
        write(6,*) ' Nsteps      ',Nsteps
        write(6,*) ' x_start     ',xx
        write(6,*) ' dx_start    ',ddxs
        write(6,*) ' Pin         ',Pin
      endif
C
C
      Do 6 j=1,nsteps
       if(j.eq.nsteps) then
         DH=abs(H)-step*(nsteps-1)
       endif
       Call MAMRKUT(xx,ddxs,pin,DH,Ient,xe,ddxe,dels,ierr)
       if(Mdebug.ge.3) then
        write(6,*) 'Step/ierr      ',j,ierr
        write(6,*) 'xe             ',xe
        write(6,*) 'dxe            ',ddxe
       endif
       if(ierr.gt.0) then
         ierr=1
         goto 9900
       endif
       Call UCOPY(xe,xx,3)
       Call UCOPY(ddxe,ddxs,3)
C... Take approx. into account energy loss.
       if(xx(1).lt.6.173.and.xx(1).gt.5.or.xx(1).gt.0.225.and.
     .    xx(1).lt.3.525) then
         Pin=Pin+Sign(1.,Pin)*Abs(DH)*dEdx
         if(Mdebug.ge.3) then
           write(6,*) 'Pin adjusted by dedx... new Pin ',Pin
         endif
       endif
    6 Continue
      ParNow=Sqrt(xe(2)**2+xe(3)**2)
      thf=acos(ddxe(1))
      phf=atan2(ddxe(3),ddxe(2))
c      phf=atan2(xe(3),xe(2))
      if(thf.lt.0.) thf=thf+pi
      if(phf.lt.0.) phf=phf+twopi
      if(Mdebug.ge.2) then
        write(6,*) 'PARNOW,PAREND,TARDIF ',parnow,parend,tardif
      endif
      if(abs(ParNow-ParEnd).gt.TarDif) then
       Write(6,*) 'MFXT. Prediction not reached, skip.'
       ierr=2
       goto 9900
      endif
      angf(1)=thf
      angf(2)=phf
      xf(1)=xe(2)*100.
      xf(2)=xe(3)*100.
      xf(3)=xe(1)*100.
      pf(1)=Pin*ddxe(2)
      pf(2)=Pin*ddxe(3)
      pf(3)=Pin*ddxe(1)
C
 9900 Continue
C
      Return
      End
C