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, E0
  integer i, j

  vx0=0.09 ;  vy0=0 ; dt=0.00001 ; Omega = 1
  x0=5 ; y0 = 0
  E0 = vx0**2 + vy0**2  ! initial energy

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

  do j=1,20
     do i=0,10000000

        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,100000) .eq. 0) write(20, *) x1, y1

     enddo
     print *, "finished iteration ", j
  enddo

  print *, "Energy/(Initial Energy) = ", (vx1**2+vy1**2)/E0

end program dipole