       subroutine rkf45(dsi,yold,fold,ynew,fnew,ds4,ds5)
       
c***********************************************************************
c      Runge-Kutta-Fehlberg integrator
c      Written by Orlando Camargo Rodriguez
c      orodrig@ualg.pt
c      Universidade do Algarve
c      Physics Department
c      Signal Processing Laboratory
c      Faro, 17/12/2009 at 14:30
c****&******************************************************************

c      Include global (common) variables:

       include 'global.for'
         
c      Define local variables:

       integer*4 j
 
       real*8 yold(4),ynew(4),yrk4(4),yrk5(4),yk(4)
       real*8 fold(4),fnew(4) 
       
       real*8 k1(4),k2(4),k3(4),k4(4),k5(4),k6(4)

       real*8 es(2)
       
       real*8 a1,a3,a4,a5
       real*8 b1,b3,b4,b5,b6
       real*8 ri,zi
       real*8 ci,cc
       real*8 cri,czi
       real*8 crri,czzi,crzi
       real*8 sigmai,sigmar,sigmaz
       real*8 dsi,ds4,ds5
       real*8 sri,szi
       real*8 dummy
       
c***********************************************************************
       
c***********************************************************************

       a1 =   25.0/216.0
       a3 = 1408.0/2565.0
       a4 = 2197.0/4101.0
       a5 =   -1.0/5.0
       
       b1 =    16.0/135.0
       b3 =  6656.0/12825.0
       b4 = 28561.0/56430.0
       b5 =    -9.0/50.0
       b6 =     2.0/55.0

c***********************************************************************
       
       ri = yold(1)
       zi = yold(2) 
       
       call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)

       do j = 1,4

          k1(j) = fold(j)
          yk(j) = yold(j) + 0.25*dsi*k1(j)
 
       end do

           ri = yk(1)
           zi = yk(2)
       sigmar = yk(3)
       sigmaz = yk(4)
       
       sigmai = sqrt( sigmar**2 + sigmaz**2 )

       es(1) = sigmar/sigmai
       es(2) = sigmaz/sigmai
              
       call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)
       
       k2(1) = es(1)
       k2(2) = es(2)
       k2(3) = sri
       k2(4) = szi
        
       do j = 1,4
	
	  yk(j) = yold(j) + dsi*(3.0/32.0*k1(j)+9.0/32.0*k2(j))
	
       end do
        
           ri = yk(1)
           zi = yk(2)
       sigmar = yk(3)
       sigmaz = yk(4)
           
       sigmai = sqrt( sigmar**2 + sigmaz**2 )
       
       es(1) = sigmar/sigmai
       es(2) = sigmaz/sigmai
              
       call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)
       
       k3(1) = es(1)
       k3(2) = es(2)
       k3(3) = sri
       k3(4) = szi
       
       do j = 1,4
	
	  yk(j) = yold(j) + dsi*(1932.0/2197.0*k1(j)
     &      -7200.0/2197.0*k2(j)+7296.0/2197.0*k3(j))
	
       end do
        
           ri = yk(1)
           zi = yk(2)
       sigmar = yk(3)
       sigmaz = yk(4)
           
       sigmai = sqrt( sigmar**2 + sigmaz**2 )
       
       es(1) = sigmar/sigmai
       es(2) = sigmaz/sigmai
              
       call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)

       k4(1) = es(1)
       k4(2) = es(2)
       k4(3) = sri
       k4(4) = szi
       
       do j = 1,4
	
	  yk(j) = yold(j) + dsi*(439.0/216.0*k1(j)-8.0*k2(j)
     &        +3680.0/513.0*k3(j)-845.0/4104*k4(j))
	
       end do
        
           ri = yk(1)
           zi = yk(2)
       sigmar = yk(3)
       sigmaz = yk(4)
           
       sigmai = sqrt( sigmar**2 + sigmaz**2 )
       
       es(1) = sigmar/sigmai
       es(2) = sigmaz/sigmai
              
       call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)

       k5(1) = es(1)
       k5(2) = es(2)
       k5(3) = sri
       k5(4) = szi
             
       do j = 1,4
	
	  yk(j) = yold(j) + dsi*(2.0*k2(j)-8.0/27.0*k1(j)
     &   +3544.0/2565.0*k3(j)+1859.0/4104*k4(j)-11.0/40.0*k5(j))
	
       end do
        
           ri = yk(1)
           zi = yk(2)
       sigmar = yk(3)
       sigmaz = yk(4)
           
       sigmai = sqrt( sigmar**2 + sigmaz**2 )
       
       es(1) = sigmar/sigmai
       es(2) = sigmaz/sigmai
              
       call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)
       
       k6(1) = es(1)
       k6(2) = es(2)
       k6(3) = sri
       k6(4) = szi
             
       do j = 1,4
	
	  yrk4(j)=yold(j)+dsi*(a1*k1(j)+a3*k3(j)+a4*k4(j)+a5*k5(j))
	  yrk5(j)=yold(j)+dsi*(b1*k1(j)+b3*k3(j)+b4*k4(j)+b5*k5(j)
     &           +b6*k6(j))
	  
	  ynew(j) = yrk5(j)
	  
       end do
       
           ri = ynew(1)
           zi = ynew(2)
       sigmar = ynew(3)
       sigmaz = ynew(4)
       
       dr = yrk4(1) - yold(1)
       dz = yrk4(2) - yold(2)
           
       ds4 = sqrt(dr*dr+dz*dz)       

       dr = yrk5(1) - yold(1)
       dz = yrk5(2) - yold(2)
       
       ds5 = sqrt(dr*dr+dz*dz)
           
       sigmai = sqrt( sigmar**2 + sigmaz**2 )
       
       es(1) = sigmar/sigmai
       es(2) = sigmaz/sigmai
              
       call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)

       fnew(1) = es(1)
       fnew(2) = es(2)
       fnew(3) = sri
       fnew(4) = szi
  
c*********************************************************************** 
c      Back to tracer...
c*********************************************************************** 
       
       return 
       
       end 

