program integral
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      write(6,*)' N, ax, ay, az, nx, ny, nz:'
      read(5,*)N,ax,ay,az,nx,ny,nz
      pi=4.0d0*atan(1.0d0)
      akx=pi*nx/ax
      aky=pi*ny/ay
      akz=pi*nz/az
      ak2=akx**2+aky**2+akz**2
      dx=ax/dble(N)
      dy=ay/dble(N)
      dz=az/dble(N)
      sr=0.0d0
      si=0.0d0
      x=-dx/2.0d0
      do 1 ix=1,N
      x=x+dx
      xx=x**2
      skx=akx*x
      y=-dy/2.0d0
      do 1 iy=1,N
      y=y+dy
      yy=xx+y**2
      sky=skx+aky*y
      z=-dz/2.0d0
      do 1 iz=1,N
      z=z+dz
      zz=yy+z**2
      skz=sky+akz*z
      dd=1.0d0/dsqrt(zz)
      si=si-sin(skz)*dd
1     sr=sr+cos(skz)*dd
      sr=sr*dx*dy*dz
      si=si*dx*dy*dz
      v=ax*ay*az
      write(6,*)'Real:',sr,sr/(4.0d0*pi/ak2)
      write(6,*)'Imag:',si,si/(4.0d0*pi/ak2)
c
c     composite
c     first cube  0..dx etc
      N1=10
      d1x=dx/dble(N1)
      d1y=dy/dble(N1)
      d1z=dz/dble(N1)
      sr1=0.0d0
      si1=0.0d0
      x=-d1x/2.0d0
      do 2 ix=1,N1
      x=x+d1x
      xx=x**2
      skx=akx*x
      y=-d1y/2.0d0
      do 2 iy=1,N1
      y=y+d1y
      yy=xx+y**2
      sky=skx+aky*y
      z=-d1z/2.0d0
      do 2 iz=1,N1
      z=z+d1z
      zz=yy+z**2
      skz=sky+akz*z
      dd=1.0d0/dsqrt(zz)
      si1=si1-sin(skz)*dd
2     sr1=sr1+cos(skz)*dd
      sr1=sr1*d1x*d1y*d1z
      si1=si1*d1x*d1y*d1z
c
c     rest dx..ax etc
      sr2=0.0d0
      si2=0.0d0
      x=-dx/2.0d0
      do 3 ix=1,N
      x=x+dx
      xx=x**2
      skx=akx*x
      y=-dy/2.0d0
      do 3 iy=1,N
      y=y+dy
      yy=xx+y**2
      sky=skx+aky*y
      z=-dz/2.0d0
      do 3 iz=1,N
      z=z+dz
      if(ix.ne.1.or.iy.ne.1.or.iz.ne.1)then
       zz=yy+z**2
       skz=sky+akz*z
       dd=1.0d0/dsqrt(zz)
       si2=si2-sin(skz)*dd
       sr2=sr2+cos(skz)*dd
      endif
3     continue
      sr2=sr2*dx*dy*dz
      si2=si2*dx*dy*dz

      sr=sr1+sr2
      si=si1+si2
      write(6,*)'Composite:'
      write(6,*)'Real:',sr
      write(6,*)'Imag:',si
      stop
      end