program dipole

! Program to test numerical methods for integrating particle orbits in
! a magnetic field

  implicit none

  double precision vx0, vy0, vx1, vy1
  double precision x0, y0, x1, y1, Omega
  double precision dt
  integer i

  vx0=1 ;  vy0=0 ; dt=0.0001 ; Omega = 1
  x0=5 ; y0 = 0

  open(unit=20, file="dipole.dat", status="unknown")

  do i=1,1000000

     Omega = 5**1.5/(x0**2+y0**2)**1.5

     vx1 = vx0 + dt*Omega*vy0
     vy1 = vy0 - dt*Omega*vx0
     x1 = x0 + dt*vx0
     y1 = y0 + dt*vy0

     vx0=vx1
     vy0=vy1
     x0=x1
     y0=y1

     if( mod(i,1000) .eq. 0) write(20, *) x1, y1

  enddo

  print *, "Energy = ", vx1**2+vy1**2

end program dipole