#define TWOPI 6.2831853
#include <math.h>

simpson(z,x1,x2,y1,y2,n,sr,si)
float z,x1,x2,y1,y2, *si, *sr;  int n;
   { float h,tr,ti,f,x,y,hy,z2,z4,yy,zz,yyzz,t,r,r2; int i,m,j;

     m = (n+1)/2; n = 2*m;  z2=z+z;  z4=z2+z2;  zz=z*z;
     h=(x2-x1)/n;  hy=(y2-y1)/n;
     y=y1; *sr = 0; *si = 0;
     for (j = 0; j <= n; ++j)
        {x=x1; tr = 0; ti = 0; yyzz = y*y + zz;
        for (i = 0; i <= n;  ++i)
          {if ((i == 0) || (i == n)) t = z; else
           if ((i & 1) == 1)  t = z4; else t = z2;
           r2 = yyzz + x*x;
	   r = sqrt((double) r2);
           t = t/(r2*r);   r=TWOPI*r;
           tr = tr + t * cos(r);
           ti = ti + t * sin(r);
           x = x+h; }
        if ((j==0) || (j==n)) t = 1; else
        if ((j & 1) == 1) t = 4; else t = 2;
        *sr = *sr + t*tr;  *si = *si + t*ti;   y = y+hy; }
     t=h*h/(9*TWOPI);   *sr = *sr * t;  *si = *si * t;     
    }

main ()
{
    float   x, y, z, si, sr, Ssr, Ssi;

    int     i,j,k,l,m,n,ro,nn,i1,j1,n1,ii,jj;

    z = 300;  n = 512;  ii = 10;  jj = 20;
    i = ii;  j = jj;
    simpson(z,j/2.0-0.25,j/2.0+0.25,i/2.0-0.25,i/2.0+0.25,8,&Ssr,&Ssi);    
    printf("%d  %e  %e i\n",0,Ssr,Ssi);
    for (l = 1;  l < 1000; ++l)
      {for (i1 = -l; i1 <= l; ++i1)
         {i = i1*n + ii;  j = -l*n + jj;
         simpson(z,j/2.0-0.25,j/2.0+0.25,i/2.0-0.25,i/2.0+0.25,8,&sr,&si);
         Ssr += sr;  Ssi += si;
         j = l*n + jj;
         simpson(z,j/2.0-0.25,j/2.0+0.25,i/2.0-0.25,i/2.0+0.25,8,&sr,&si);
         Ssr += sr;  Ssi += si;}

       for (j1 = -l+1; j1 <= l-1; ++j1)
         {i = -l*n + ii;  j = j1*n + jj;
         simpson(z,j/2.0-0.25,j/2.0+0.25,i/2.0-0.25,i/2.0+0.25,8,&sr,&si);
         Ssr += sr;  Ssi += si;
         i = l*n + ii;
         simpson(z,j/2.0-0.25,j/2.0+0.25,i/2.0-0.25,i/2.0+0.25,8,&sr,&si);
         Ssr += sr;  Ssi += si;}
      printf("%d  %e %e i\n",l,Ssr,Ssi);

      }
}