      Subroutine derivs(t, r, drdt, nStts)

      double precision r(*),drdt(*),t
      integer nStts

c     Calculates the derivates for the Newton-Lorentz equation
c     for a charged particle. This version assumes periodic
c     boundary conditions for the em field grids. So, in evaluating
c     the em fields the particle's position is always mapped back
c     onto the 'unit cube'.
c
c     INPUT PARAMETERS
c  
c         t  -- time in sufficiently normalized units
c
c         r  -- a 6-tuple denoting (X,Y,Z,Ux,Uy,Uz) for the ptcl
c
c     OUTPUT PARAMETERS
c
c         drdt -- the 6-tuple denoting the time derivative of r
c
c     The following lines are not needed as it is not needed for this case.

      double precision th0, ph0, rr0, z0
      common /inisph/ th0, ph0, rr0, z0

      integer  field_type, derivative_type
      common /setup/  field_type, derivative_type

      double precision alpha, lambda
      common /sptrm/ nkmin,nkmax,alpha,lambda

c
      logical debug
      common /debug/ debug

c     LOCAL PARAMETER
      double precision x, y, z, xx, yy, zz
c
      double precision bx, by, bz, bb
      double precision bx1, by1, bz1

      double precision alpha3
      common /particle/ alpha3

      double precision ex, ey, ez
      double precision vx, vy, vz
      double precision vt, qom
      double precision check
      double precision vvv,bbb,cosvb

cc    speed of light in AU/min unit.
      double precision vc
c      parameter (vc = 0.1202)
      parameter (vc = 1.0)

      double precision PI
      parameter (PI = 3.14159265358979323846d0)

c     x -> \phi       vx -> v_\phi
c     y -> \theta     vy -> v_\theta
c     z -> r          vz -> v_r
c
c     13 Dec 07 change theta -> lambda 
c               where lambda = PI/2. - theta
c     we gendata as lambda so: y = r(2)
c     and vy = v_lambda = r(5)
c

      x = r(1)
c      y = th0-r(2)
      y = r(2)
      z = r(3)
      
      vx = r(4)
      vy = r(5)
      vz = r(6)

      xx = r(4)
c      yy = th0-r(5)
      yy = r(5)
      zz = r(6)
c
c     Evaluate the EM fields at the particle position
c

      if (field_type .eq. 1) then
c     bz = b0, bx = by = 0 
        call emfldspar1(t,x,y,z,bx,by,bz,nStts)
        if (derivative_type .eq. 1) then
          call emfldspar1(t,xx,yy,zz,bx1,by1,bz1,nStts)
        endif
      endif

      if (field_type .eq. 2) then
c     bz = b0*(r0^2)/r^2, bx = by = 0
        call emfldspar2(t,x,y,z,bx,by,bz,nStts)
        if (derivative_type .eq. 1) then
          call emfldspar2(t,xx,yy,zz,bx1,by1,bz1,nStts)
        endif
      endif

      if (field_type .eq. 3 ) then
c     bz = b0(r0^2)/r^2, bx = (bxslab+bx2d)*r0^2/r^2, 
c     by = (byslab+by2d)*r0^2/r^2
        call emfldspar3(t,x,y,z,bx,by,bz,nStts)
        if (derivative_type .eq. 1) then
          call emfldspar3(t,xx,yy,zz,bx1,by1,bz1,nStts)
        endif
      endif

c      if (field_type .eq. 4 ) then
cc     bz = b0(r0^2)/r^2, bx = bxslab, by = byslab
c         call emfldspar4(t,x,y,z,bx,by,bz,nStts)
c         call emfldspar4(t,xx,yy,zz,bx1,by1,bz1,nStts)
c      endif

c      if (field_type .eq. 5 ) then
cc     bz = b0(r0^2)/r^2, bx = bxslab*r0/r+bx2d, by = byslab*r0/r+by2d
c         call emfldspar5(t,x,y,z,bx,by,bz,nStts)
c         call emfldspar5(t,xx,yy,zz,bx1,by1,bz1,nStts)
c      endif

c      if (field_type .eq. 6 ) then
cc     bz = b0(r0^2)/r^2, bx = bxslab+bx2d, by = byslab+by2d
c         call emfldspar6(t,x,y,z,bx,by,bz,nStts)
c         call emfldspar6(t,xx,yy,zz,bx1,by1,bz1,nStts)
c      endif

      if (field_type .eq. 8 ) then
c     read from file
         call emfldspar8(t,x,y,z,bx,by,bz,nStts)
         call emfldspar8(t,xx,yy,zz,bx1,by1,bz1,nStts)
      endif

c      if (field_type .eq. 9 ) then
c     spherical harmonic surface
c         call emfldspar9(t,x,y,z,bx,by,bz,nStts)
C         call emfldspar9(t,xx,yy,zz,bx1,by1,bz1,nStts)
c      endif

      if (field_type .eq. 10 ) then
c     bz = b0(r0^2)/r^2, bx = (bxslab+bx2d)*r0^2/r^2, 
c     by = (byslab+by2d)*r0^2/r^2
        call emfldspar3(t,x,y,z,bx,by,bz,nStts)
        if (derivative_type .eq. 1) then
          call emfldspar3(t,xx,yy,zz,bx1,by1,bz1,nStts)
        endif
      endif

      if(nStts .eq. 0) return

c
c     Now form the ODE with spatially interpolated 
c     E(x,y,z,t) and B(x,y,z,t)
c
        
      bb = dsqrt(bx*bx+by*by+bz*bz)


c
c     Spherical coordinate
c
       if (derivative_type .eq. 1) then
c     Created the magnetic fields in spherical
c     bx -> b_\phi
c     by -> b_\theta
c     bz -> b_r
          drdt(1) = bx/bb
          drdt(2) = by/bb
          drdt(3) = bz/bb
          drdt(4) = bx1/bz1
          drdt(5) = by1/bz1
          drdt(6) = 1.0
c          write(*,*) 'bx1,by1,bz1,bb = ',bx1,by1,bz1,bb
       endif

       if (derivative_type .eq. 2) then
cccccc
c       Spherical Coordinate :
c       bx -> b_phi
C       by -> b_theta
c       bz -> b_r
c       x -> \phi   -> r(1)
c       y -> \theta -> r(2)
c       z -> r      -> r(3)
c       vx -> r(4) -> v_{\phi}
c       vy -> r(5) -> v_{\theta}
c       vz -> r(6) -> v_{r}
c       drdt(1) ==> \dot{\phi} = \frac{v_{\phi}}{r\sin{\theta}}
c       drdt(2) ==> \dot{\theta} = \frac{v_{\theta}}{r}
c       drdt(3) ==> \dot{r} = v_{r}
c
c     ==> drdt(4) --> \dot{v}_{\phi}
c        \dot{v}_{\phi} = \frac{q}{m\gamma}(v_{r}B_{\theta} - v_{\theta}B_{r})
c                      - \frac{v_{r}v_{\phi}}{r}
c                      - \frac{v_{\theta}v_{\phi}}{r\tan{\theta}}
c     ==> drdt(5) --> \dot{v}_{\theta}
c        \dot{v}_{\theta} = \frac{q}{m\gamma}(v_{\phi}B_{r} - v_{r}*B_{\phi})
c                        - \frac{v_{r}v_{\theta}}{r}
c                        + \frac{v_{\phi}^{2}}{r\tan{\theta}}
c     ==> drdt(6) --> \dot{v}_{r}
c        \dot{v}_{r} = \frac{q}{m\gamma}(v_{\theta}B_{\phi} - v_{\phi}B_{\theta})
c                   + \frac{v^2_{\theta}}{r}
c                   + \frac{v^2_{\phi}{r}
cccccc
         qom = alpha3

         drdt(1) = r(4)
         drdt(2) = r(5)
         drdt(3) = r(6)
         drdt(4) = qom*(r(5)*bz-r(6)*by)
         drdt(5) = qom*(r(6)*bx-r(4)*bz)
         drdt(6) = qom*(r(4)*by-r(5)*bx)
                  
         
c         drdt(1) = r(4)/(r(3)*dsin(r(2)))
c         drdt(2) = r(5)/r(3)
c         drdt(3) = r(6) 

c          vt=dsqrt(r(4)**2.+r(5)**2.+r(6)**2.)
c         gamma_v = dsqrt(1.0-((vt**2.)/(vc**2.)))

c         drdt(4) = qom*(vz*by - vy*bz) - (vz*vx/z)
c     &             - (vy*vx/(z*dtan(y)))
c         drdt(5) = qom*(vx*bz - vz*bx) - (vz*vy/z)
c     &             + (vx*vx/(z*dtan(y)))
c         drdt(6) = qom*(vy*bx - vx*by) + (vy*vy/z)
c     &             + (vx*vx/z)

     
cc
cc    check v||b  
cc
c         vvv = dsqrt(vx*vx+vy*vy+vz*vz)
c         bbb = dsqrt(bx*bx+by*by+bz*bz)  
c         cosvb = (vx*bx + vy*by + vz*bz)
c     &           /(vvv*bbb)
c         write(*,*) 'cosvb= ',cosvb,x,y,z,bx,by,bz    
    
c         check=vz/vt
c         write(*,*)'bx,by,bz = ',bx,by,bz
c         write(*,*)'vt,mu = ',vt,check
c         check=(1.0-(r(6)**2./vt**2.))/dsqrt(bx*bx+by*by+bz*bz)
c         write(*,*)'vt,(1-mu^2)/B = ',vt,check
c         write(*,*)'bphi,btheta,br = ',bx,by,bz

      endif

       if (derivative_type .eq. 3) then
cccccc
c       Spherical Coordinate :
c       bx -> b_phi
c       by -> b_theta
c       bz -> b_r
c       x -> \phi   -> r(1)
c       y -> \theta -> r(2)
c       z -> r      -> r(3)
c       vx -> r(4) -> v_{\phi}
c       vy -> r(5) -> v_{\theta}
c       vz -> r(6) -> v_{r}
c       drdt(1) ==> \dot{\phi} = \frac{v_{\phi}}{r\sin{\theta}}
c       drdt(2) ==> \dot{\theta} = \frac{v_{\theta}}{r}
c       drdt(3) ==> \dot{r} = v_{r}
c
c     ==> drdt(4) --> \dot{v}_{\phi}
c        \dot{v}_{\phi} = \frac{q}{m\gamma}(v_{r}B_{\theta} - v_{\theta}B_{r})
c                      - \frac{v_{r}v_{\phi}}{r}
c                      - \frac{v_{\theta}v_{\phi}}{r\tan{\theta}}
c     ==> drdt(5) --> \dot{v}_{\theta}
c        \dot{v}_{\theta} = \frac{q}{m\gamma}(v_{\phi}B_{r} - v_{r}*B_{\phi})
c                        - \frac{v_{r}v_{\theta}}{r}
c                        + \frac{v_{\phi}^{2}}{r\tan{\theta}}
c     ==> drdt(6) --> \dot{v}_{r}
c        \dot{v}_{r} = \frac{q}{m\gamma}(v_{\theta}B_{\phi} - v_{\phi}B_{\theta})
c                   + \frac{v^2_{\theta}}{r}
c                   + \frac{v^2_{\phi}{r}
cccccc
         qom = alpha3

         drdt(1) = r(4)
         drdt(2) = r(5)
         drdt(3) = r(6) 

          vt=dsqrt(r(4)**2.+r(5)**2.+r(6)**2.)
c         gamma_v = dsqrt(1.0-((vt**2.)/(vc**2.)))

         drdt(4) = qom*((by*vz/(z*dsin(y)))-(bz*vy/dsin(y)))
     &         -(2.0*vz*vx/z)-(2.0*vy*vx/dtan(y))
         drdt(5) = qom*((dsin(y)*vx*bz)-(vz*bx/z))
     &         -(2.0*vz*vy/z)+(vx*vx*dsin(y)*dcos(y))
         drdt(6) = qom*((z*vy*bx)-(z*dsin(y)*by*vx))
     &         +(vy*vy*z)+(vx*vx*z*dsin(y)*dsin(y))


c         check=vz/vt
c         write(*,*)'bx,by,bz = ',bx,by,bz
c         write(*,*)'vt,mu = ',vt,check
c         check=(1.0-(r(6)**2./vt**2.))/dsqrt(bx*bx+by*by+bz*bz)
c         write(*,*)'vt,(1-mu^2)/B = ',vt,check
c         write(*,*)'bphi,btheta,br = ',bx,by,bz
      endif

cccccc
c     Spherical coordinate
c
c     13 Dec 07 change theta -> lambda   
c
c     trace magnetic field in terms of (phi, lambda, r)          
c     using derivative_type = 4
c
c     trace particle trajectories in terms of (phi, lambda, r)
c     using derivative_type = 5
c
c

       If (derivative_type .eq. 4) then
c     Created the magnetic fields in spherical
c     bx -> b_phi
c     by -> b_lambda
c     bz -> b_r
c
c     equations for tracing magnetic field lines are:
c     dr/b_r = r d_lambda / b_lambda = r cos(lambda) d_phi / b_phi
c
c     d_phi / dr    = b_phi / (b_r r cos(lambda) )
c     d_lambda / dr = b_lambda / (b_r r)  
c

c
c    use these lines if you gendata in degrees
c
c          drdt(1) = bx/(bb*z*dcos((y*PI/180.0d0)))
c          drdt(2) = by/(z*bb)
c          drdt(3) = bz/bb
c          drdt(4) = bx1/(bz1*zz*dcos((yy*PI/180.0d0)))
c          drdt(5) = by1/(zz*bz1)
c          drdt(6) = 1.0

c          write(*,*) 'bx1,by1,bz1,bb = ',bx1,by1,bz1,bb

c
c    use these lines if you gendata in radian
c
          drdt(1) = bx/(bb*z*dcos(y))
          drdt(2) = by/(z*bb)
          drdt(3) = bz/bb
          Drdt(4) = bx1/(bz1*zz*dcos(yy))
          drdt(5) = by1/(zz*bz1)
          drdt(6) = 1.0


       endif

       if (derivative_type .eq. 5) then
cccccc
c       Spherical Coordinate :
c       bx -> b_phi
c       by -> b_lambda
c       bz -> b_r
c       x -> \phi   -> r(1)
c       y -> \lambda -> r(2)
c       z -> r      -> r(3)
c       vx -> r(4) -> v_{\phi}
c       vy -> r(5) -> v_{\lambda} = -v{\theta}
c       vz -> r(6) -> v_{r}
c       drdt(1) ==> \dot{\phi} = \frac{v_{\phi}}{r\sin{\theta}}
c       drdt(2) ==> \dot{\lambda} = \frac{v_{\lambda}}{r}
c       drdt(3) ==> \dot{r} = v_{r}
c
c     ==> drdt(4) --> \dot{v}_{\phi}
c        \dot{v}_{\phi} = \frac{q}{m\gamma}(v_{\lambda}B_{r} - v_{r}B_{lambda})
c                      - \frac{v_{r}v_{\phi}}{r}
c                      + \frac{v_{\lambda}v_{\phi}}{r\cot{\lambda}}
c     ==> drdt(5) --> \dot{v}_{\lambda}
c        \dot{v}_{\lambda} = \frac{q}{m\gamma}(v_{r}B{{\phi} - v_{\phi}*B_{r})
c                        - \frac{v_{r}v_{\lambda}}{r}
c                        - \frac{v_{\phi}^{2}}{r\cot{\lambda}}
c     ==> drdt(6) --> \dot{v}_{r}
c        \dot{v}_{r} = \frac{q}{m\gamma}(v_{\phi}B_{\lambda} - v_{\lambda}B_{\phi})
c                   + \frac{v^2_{\lambda}}{r}
c                   + \frac{v^2_{\phi}{r}
c
c
cccccc
         qom = alpha3

         drdt(1) = r(4)/(r(3)*dcos(r(2)))
         drdt(2) = r(5)/r(3)
         drdt(3) = r(6) 

c          vt=dsqrt(r(4)**2.+r(5)**2.+r(6)**2.)
c         gamma_v = dsqrt(1.0-((vt**2.)/(vc**2.)))

         drdt(4) = qom*(vy*bz-vz*by) - (vz*vx/z)
     &             + (vy*vx*dtan(y)/z)
         drdt(5) = qom*(vz*bx-vx*bz) - (vz*vy/z)
     &             - (vx*vx*dtan(y)/z)
         drdt(6) = qom*(vx*by-vy*bx) + (vy*vy/z)
     &             + (vx*vx/z)

cc
cc    check v||b  
cc
c         vvv = dsqrt(vx*vx+vy*vy+vz*vz)
c         bbb = dsqrt(bx*bx+by*by+bz*bz)  
c         cosvb = (vx*bx + vy*by + vz*bz)
c     &           /(vvv*bbb)
c         write(*,*) 'cosvb= ',cosvb,x,y,z,bx,by,bz    
    
c         check=vz/vt
c         write(*,*)'bx,by,bz = ',bx,by,bz
c         write(*,*)'vt,mu = ',vt,check
c         check=(1.0-(r(6)**2./vt**2.))/dsqrt(bx*bx+by*by+bz*bz)
c         write(*,*)'vt,(1-mu^2)/B = ',vt,check
c         write(*,*)'bphi,btheta,br = ',bx,by,bz

      endif
      return
      end
