#include <cmacs.h>
#include <image.h>

main(argc,argv)
int argc;
char *argv[];
{   float   *R,*I,t,sr,si,dr,di,a,r,rmin,rmax,red,grn,blu,smin,smax;
    int     i,j,k,l,m,n,nn,FIL,OFIL;
IMAGE *Ri,*Gi,*Bi,*Wi;

FIL = -1;
if (argc != 2) printf("Command line must include file name\n"); 
else FIL = open(argv[1],0);
if (FIL == -1) printf("%s - file not found\n",argv[1]);
else
{
printf("creating array\n");
m = 9;
n = 1 << m;   printf("n = %d\n",n);
nn = n * n;
R = (float*) malloc(nn * sizeof(float));
I = (float*) malloc(nn * sizeof(float));
if (R == 0 || I == 0) printf("Oops - couldn't allocate %d words\n",2*nn);
printf("array made\n");
read(FIL,R,nn*sizeof(float));  read(FIL,I,nn*sizeof(float));  close(FIL);
printf("array read\n");

for (i = 0; i < nn; ++i)
   {if ((i % (512*64)) == 0) printf("input %d (line %d)\n",i,i / 512);
    a = ATAN2(R[i],I[i]);    r = HYPOT(R[i],I[i]);
    R[i] = r;  I[i] = a; }

rmin = rmax = HYPOT(R[0],I[0]);
for (i = n/8; i < n/2; ++i) for (j = n/8; j < n/2; ++j)
    {r = R[i*n+j];
    if (r < rmin) rmin = r;     if (r > rmax) rmax = r;}

printf("rmin %e, rmax %e\n\n",rmin,rmax);

rmax=rmax/3.5;

printf("creating pictures\n");
Ri = imgcreat("r.img",0655,8,0,511,0,511);
Gi = imgcreat("g.img",0655,8,0,511,0,511);
Bi = imgcreat("b.img",0655,8,0,511,0,511);
Wi = imgcreat("w.img",0655,8,0,511,0,511);

for (i=0; i<n; ++i)
for (j=0; j<n; ++j) 
   {a = I[i*n+j];  r = (R[i*n+j]-rmin)/(rmax-rmin);
    red = SIN(a);  grn = SIN(a+TWOPI/3.0);  blu = SIN(a+2.0*TWOPI/3.0);
    if (r > 0.5) {smax = 1.0;  smin = (r-0.5)*2.0;}
        else {smax = 2.0*r;  smin=0.0;}
    red=(smax-smin)*(red+1)/2+smin;
    grn=(smax-smin)*(grn+1)/2+smin;
    blu=(smax-smin)*(blu+1)/2+smin;
    l = 255*red; if (l<0) l=0; if (l>255) l=255; setpixel(Ri,i,j,l);
    l = 255*grn; if (l<0) l=0; if (l>255) l=255; setpixel(Gi,i,j,l);
    l = 255*blu; if (l<0) l=0; if (l>255) l=255; setpixel(Bi,i,j,l);
    l = 255*r  ; if (l<0) l=0; if (l>255) l=255; setpixel(Wi,i,j,l);}

imgclose(Ri);
imgclose(Gi);
imgclose(Bi);
imgclose(Wi);
}
}