#include <cmacs.h>
#include <video.h>

integ(z,x1,x2,y1,y2,sr,si)
float z,x1,x2,y1,y2, *sr, *si;
   { 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[200];
    float   x, y, z, sr, si, srs, sis, sa, sm, s, rd, gr, bl;
    int     i,j,k,l,m,n,ord,ro,nn,i1,j1,n1, FIL;
    int     ncpy;
#define SIZE 128

 SELECT(WINDOWSET(-SIZE*1.2,-SIZE,SIZE,SIZE,"RGBW"));

 for (ncpy=1;  ncpy<20;  ++ncpy)
 for (l=1; l<argc; ++l)
  {sscanf(argv[l],"%d",&ord);   sprintf(name,"P%d.%d.img",ord,ncpy);

   z = ord;  printf("(%d) z = %f, %dth copy into %s\n",l,z,ncpy,name);
   s = 0;
     for (i1 = -ncpy*SIZE;  i1 <= ncpy*SIZE;  i1 += SIZE)
     for (j1 = -ncpy*SIZE;  j1 <= ncpy*SIZE;  j1 += SIZE)
       {integ(z,i1-0.25,i1+0.25,j1-0.25,j1+0.25,&sr,&si);
        s += SQRT(sr*sr + si*si); }
   printf("Magnitude %f\n",s);
    for (i = 0; i <= SIZE; ++i)   for (j = 0; j <= i; ++j)
    {srs = sis = 0;
     for (i1 = -ncpy*SIZE;  i1 <= ncpy*SIZE;  i1 += SIZE)
     for (j1 = -ncpy*SIZE;  j1 <= ncpy*SIZE;  j1 += SIZE)
       {integ(z,i1+i/2.0-0.25,i1+i/2.0+0.25,
                j1+j/2.0-0.25,j1+j/2.0+0.25,&sr,&si);
        srs += sr;   sis += si;}
     sm = SQRT(srs*srs + sis*sis);   sa = ATAN2(sis,srs);
     rd = (COS(sa/2)*sm/s);                rd = rd*rd;
     gr = (COS(sa/2+TWOPI*0.33)*sm/s);     gr = gr*gr;
     bl = (COS(sa/2+TWOPI*0.67)*sm/s);     bl = bl*bl;
     PAINT(rd, gr, bl);
     RECTANGLE(i,j,i+1,j+1);      RECTANGLE(j,i,j+1,i+1);
     RECTANGLE(-i,-j,-i-1,-j-1);  RECTANGLE(-j,-i,-j-1,-i-1);
     RECTANGLE(i,-j,i+1,-j-1);    RECTANGLE(-j,i,-j-1,i+1);
     RECTANGLE(-i,j,-i-1,j+1);    RECTANGLE(j,-i,j+1,-i-1);
     }
  
  FILEUP(name);
  }
}