      program test1
      implicit none
      integer nx,ny,nz,nt
      parameter (nx=33,ny=33,nz=65,nt=65)
      integer i,j,k
      real u(nx,ny,nz)
      real x(nx),y(ny),z(nz)
      real pi,hx,hy,hz
      integer ixdrs
      integer initxdr
      external initxdr
      real r1,r2
      real c1(3),c2(3)

      real xmin,xmax,ymin,ymax,zmin,zmax
      parameter (xmin = -10.0, xmax= 10.0)
      parameter (ymin = -10.0, ymax= 10.0)
      parameter (zmin = -16.0, zmax= 16.0)
      real m1,m2
      parameter (m1 = 1.5, m2=1.0)
      real time
      integer itime
      real dt
      character*16 name(nt),names,names2,namefld(nt)
      
      do i = 1, nt
         names = 'testdata'
         names2 = names
         if (i.le.9) then
            write(names(9:13),'(i1,a4)') i,'.xdr'
            write(names2(9:13),'(i1,a4)') i,'.fld'
         else if (i.le.99) then
            write(names(9:14),'(i2,a4)')i,'.xdr'
            write(names2(9:14),'(i2,a4)')i,'.fld'
         else if (i.le.999) then
            write(names(9:15),'(i3,a4)')i,'.xdr'
            write(names2(9:15),'(i3,a4)')i,'.fld'
         else if (i.le.9999) then
            write(names(9:16),'(i4,a4)')i,'.xdr'
            write(names2(9:16),'(i4,a4)')i,'.fld'
         else
            write(6,*) 'error'
            stop
         end if
         name(i) = names
         namefld(i) = names2
c         write(6,*) namefld(i)
      end do
      
      
      

      pi = 4.0*atan(1.0)
      hx = (xmax-xmin)/(nx-1)
      hy = (ymax-ymin)/(ny-1)
      hz = (zmax-zmin)/(nz-1)
      do i = 1, nx
         x(i) = xmin + (i-1)*hx
      end do
      do j = 1, ny
         y(j) = ymin + (j-1)*hy
      end do
      do k = 1, nz
         z(k) = zmin +  (k-1)*hz
      end do
      dt = 2.0*pi/(nt-1)
C     dt = 0.0 
      do itime = 1 , nt
         time = (itime-1)*dt
         c1(1) = m1*cos(time)
         c1(2) = m1*sin(time)
         c2(1) = -m2*cos(time)
         c2(2) = -m2*sin(time)
         c1(3) = zmin + ( (-zmin/4.0) + 3*time)
         c2(3) = zmax - ((zmax/4.0) +3*time)
         write(6,*) time,c1(1),c1(2),c1(3),m1
         write(6,*) time,c2(1),c2(2),c2(3),m2
         do k = 1, nz
            do j= 1, ny
               do i = 1, nx
                  r1 =sqrt((x(i)-c1(1))**2+(y(j)-c1(2))**2+
     .                 (z(k)-c1(3))**2)
                  r2 =sqrt((x(i)-c2(1))**2+(y(j)-c2(2))**2+
     .                 (z(k)-c2(3))**2)
                  u(i,j,k) = m1/r1 + m2/r2 +m1*m2/(r1*r2)
               end do
            end do
         end do
         write(6,*) 'writing data for ',namefld(itime)
C         open(22,file=namefld(itime))
C         write(22,100) 
 100     format('# AVS field')
C         write(22,*) '# Time = ',time
C         write(22,*) 'ndim = 3'
C         write(22,*) 'dim1 = ',nx
C         write(22,*) 'dim2 = ',ny
C         write(22,*) 'dim3 = ',nz
C         write(22,*) 'nspace= 3'
C         write(22,*) 'veclen = 1'
C         write(22,*) 'data = float'
C         write(22,*) 'field = uniform'
C         write(22,*) 'variable 1 file=',name(itime),' filetype=binary'
C         close(22)
         
         
         ixdrs = initxdr(namefld(itime),'w')
         call xdrfldheaderuniffloat(ixdrs,nx,ny,nz)
         call xdrrmat(ixdrs,nx*ny*nz,u)
         call xdrclose(ixdrs)
      end do
      
      stop
      end
