#include <cmacs.h>
#include <math.h>

integ(z,x1,x2,y1,y2,sr,si)
float z,x1,x2,y1,y2, *si, *sr;
   { float x,y,zz,t,r,ra,x1x1,y1y1,x2x2,y2y2,t11,t22,t21,t12,tp,tn;

     x = (x1+x2)/2;   y = (y1+y2)/2;   zz = z*z;
     x1x1 = x1*x1;  y1y1 = y1*y1;     x2x2 = x2*x2;  y2y2 = y2*y2;
     t22 = x2*y2/(z*sqrt(x2x2+y2y2+zz));  t11 = x1*y1/(z*sqrt(x1x1+y1y1+zz));
     tp = (t22+t11)/(1+t22*t11);
     t12 = x2*y1/(z*sqrt(x2x2+y1y1+zz));  t21 = x1*y2/(z*sqrt(x1x1+y2y2+zz));
     tn = (t12+t21)/(1+t12*t21);

     t = (tp-tn)/(1+tp*tn);  t = ATAN(t)/TWOPI;  t = ABS(t);  t = SQRT(t);
     
     r = SQRT(x*x + y*y + zz);   ra = TWOPI*r;
     *sr = t * COS(ra);  *si = t * SIN(ra);
    }

main(argc,argv)
int argc;
char *argv[];
{   char name[50];
    float   x, y, z, si, sr;
    int     i,j,k,l,m,n,ord,ro,nn,i1,j1,n1, FIL;
 for (l=1; l<argc; ++l)
  {
   sscanf(argv[l],"%d",&ord);
   sprintf(name,"/usr/hpm/x/wave/PC%d.9",ord);

   z = ord;  printf("(%d) z = %f into %s\n",l,z,name);

    FIL = creat(name,0666);
    for (i = 0; i <= 512; ++i)   for (j = 0; j <= i; ++j)
    {integ(z,i/2.0-0.25,i/2.0+0.25,j/2.0-0.25,j/2.0+0.25,&sr,&si);
     write(FIL,&sr,sizeof(float));   write(FIL,&si,sizeof(float)); }
    close(FIL); }
}