!> functions to subtract from dipole
!! Long description
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> bxqq0(x,y,z)
!!
        function bxqq0(x,y,z)
!
! ... Local variables ..................................................
!
#include "dipole.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::bxqq0(x,y,z)"
#endif
        rsq = x**2 + y**2 + z**2
        r = sqrt(rsq)
        if ( r .le. rinner ) then
          bxqq0 = -(geoqmu/rsq)*(3.*x*z/rsq)/r
        else if ( r .ge. router ) then
          bxqq0 = 0.
        else
          rprim = amax1(rinner,amin1(router,r))
          bxqq0 = -(geoqmu/rsq)*(3.*x*z/rsq)/r*(a1 + b1*rprim)
        endif
!
! ... End function bxqq0(x,y,z).........................................
!
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> byqq0(x,y,z)
!!
        function byqq0(x,y,z)
!
! ... Local Variables ..................................................
!
#include "dipole.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::byqq0(x,y,z)"
#endif
        rsq = x**2 + y**2 + z**2
        r = sqrt(rsq)
        if ( r .le. rinner ) then
          byqq0 = -(geoqmu/rsq)*(3.*y*z/rsq)/r
        else if ( r .ge. router ) then
          byqq0 = 0.
        else
          rprim = amax1(rinner,amin1(router,r))
          byqq0 = -(geoqmu/rsq)*(3.*y*z/rsq)/r*(a1 + b1*rprim)
        endif
!
! ... End function byqq(x,y,z) .........................................
!
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> bzqq0
!!
      function bzqq0(x,y,z)
!
! ... Local variables...................................................
!
#include "dipole.inc"
! 
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::bzqq0(x,y,z)"
#endif

      rsq = x**2 + y**2 + z**2
      r = sqrt(rsq)
      if ( r .le. rinner ) then
         bzqq0 = -(geoqmu/rsq)*((2. - 3.*(x**2+y**2)/rsq)/r)
      else if ( r .ge. router ) then
         bzqq0 = 0.
      else
         rprim = amax1(rinner,amin1(router,r))
         bzqq0 = -(geoqmu/rsq)*(2.-3.*(x**2+y**2)/rsq)
     $        /(r)*(a1 + b1*rprim)
      endif
!
! ... End function bzqq0(x,y,z).........................................
!   
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> ax(x,y,z)
!! vector potentials for initial load
!!
        function ax(x,y,z)
!
! ... Local variables ..................................................
!
#include "dipole.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::ax(x,y,z)"
#endif
        r = sqrt(x**2+y**2+z**2)
        s = sqrt( ((x-xq_off)*aq)**2 + (y**2+z**2)*bq**2)
        if ( s .lt. sinner) then
          AX = (geoqmu/r**2)*(y/r)
        elseif ( s .lt. 1.0) then
          AX = (geoqmu/r**2)*(y/r)*(1.-s)/(1.-sinner)
        else
           AX = 0.
        endif
!
! .. End function ax(x,y,z) ...........................................
!
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> ay(x,y,z)
!! vector potential for initial load
!!
        function ay(x,y,z)
! ... Local varaibles...................................................
#include "dipole.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::ay(x,y,z)"
#endif
        r = sqrt(x**2+y**2+z**2)
        s = sqrt( ((x-xq_off)*aq)**2 + (y**2+z**2)*bq**2)
        if ( s .lt. sinner) then
          AY = -(geoqmu/r**2)*(x/r)
        elseif ( s .lt. 1.0) then
          AY = -(geoqmu/r**2)*(x/r)*(1.-s)/(1.-sinner)
        elseif ( s .lt. souter ) then
          AY = -Bxq*(z*cos_tilt-x*sin_tilt)*(s-1.0)/(souter-1.0)
        else
          AY = -Bxq*(z*cos_tilt-x*sin_tilt)
        endif
!
! ... End function ay(x,y,z) ...........................................
!
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> az(x,y,z)
!!
      function az(x,y,z)
! ... Variables ........................................................
#include "dipole.inc"
!
! ... Begin ............................................................
!
c      for N-S dipoles this is set to zero, otherwise all 3 Vec pots
c      neeed to include dipole direction vector
c
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::az(x,y,z)"
#endif
      az = 0.
!
! ... End function az(x,y,z) ...........................................
!
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> bxDIP(x,y,z)
!!
        function bxDIP(x,y,z)
! ... Variables ........................................................
#include "dipole.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::bxDIP(x,y,z)"
#endif
        rsq = x**2 + y**2 + z**2
        r = sqrt(rsq)
          bxDIP = -(geoqmu/rsq)*(3.*x*z/rsq)/r
!
! ... End function bxDIP(x,y,z).........................................
!
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> byDIP(x,y,z)
!!
        function byDIP(x,y,z)
! ... Variables ........................................................
#include "dipole.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::byDIP(x,y,z)"
#endif
        rsq = x**2 + y**2 + z**2
        r = sqrt(rsq)
          byDIP = -(geoqmu/rsq)*(3.*y*z/rsq)/r
!
! ... End function byDIP(x,y,z) ........................................
!
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> bzDIP(x,y,z)
!!
        real*4 function bzDIP(x,y,z)
! ... Variables ........................................................
#include "dipole.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::bzDIP(x,y,z)"
#endif
        rsq = x**2 + y**2 + z**2
        r = sqrt(rsq)
          bzDIP = -(geoqmu/rsq)*((2. - 3.*(x**2+y**2)/rsq)/r)
!
! ... End function bzDIP(x,y,z) ........................................
!
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#if defined(USE_MIX) && !defined(DISABLE_RCM)
!> ax_dip(x,y,z)
!!  vector potentials for dipole field
!!  VGM added for RCM coupling
        function ax_dip(x,y,z)
! ... Variables ........................................................
#include "dipole.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::ax_dip(x,y,z)"
#endif
        r = sqrt(x**2+y**2+z**2)
        s = sqrt( ((x+xq_off)*aq)**2 + (y**2+z**2)*bq**2)
          AX_dip = (geoqmu/r**2)*(y/r)
!
! ... End function ax_dip(x,y,z)........................................
!
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> ay_dip 
!! Vector potential for dipole field
!! VGM added for RCM coupling
        function ay_dip(x,y,z)
! ... Variables ........................................................
#include "dipole.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in dipole.F::ay_dip(x,y,z)"
#endif
        r = sqrt(x**2+y**2+z**2)
        s = sqrt( ((x+xq_off)*aq)**2 + (y**2+z**2)*bq**2)
          AY_dip = -(geoqmu/r**2)*(x/r)
! ... End function ay_dip(x,y,z) .......................................
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#endif