/*
 *   Plasma Cloud and other Fractal Stuff
 *
 *  improvement in drawSpiral, mixer, and cells  by Daniel Cabeza Gras
 */
#ifdef X
#include <X11/X.h>
#include <X11/Xlib.h>
#else
//#include <metawndo.h>
#include <stdlib.h>
//#include "fakex.h"
#endif
#include <stdio.h>
#include <math.h>
#include "trippy.h"
#include "plasma.h"
#include "color.h"
#ifndef M_PI_2
#define M_PI_2 M_PI/2
#endif

#ifndef M_PI_4
#define M_PI_4 M_PI/4
#endif

#ifdef WINDOWS
static HDC hdc;
extern HPALETTE hPal;
#endif


int isqrt(int);
extern foo options;

struct COMPLEX
{
  double real;
  double imag;
};

struct COMPLEX addComplex(struct COMPLEX,struct COMPLEX);
struct COMPLEX multComplex(struct COMPLEX,struct COMPLEX);
struct COMPLEX expComplex(struct COMPLEX);
struct COMPLEX powComplex(struct COMPLEX,struct COMPLEX);
struct COMPLEX divComplex(struct COMPLEX,struct COMPLEX);

/* Plasma code based on OOZE.C, written by Jeff Clough for the IBM/PCs */

static void draw_grid(int,int,int);
static short adjust(int,int,int,int);
static void subdiv (int,int,int,int);

 void project(int *,int *,int*); 
/* short *grid;  grid is a list of the colors */

#ifdef X
extern XEvent event;
#endif

void redrawPlasma(int winno, PAINTSTRUCT* ps)
{
  short cs0,cs1,cs2,cs3;
  hdc = ps->hdc;
  cs0=(short)rndm(numcolors);
  cs1=(short)rndm(numcolors);
  cs2=(short)rndm(numcolors);
  cs3=(short)rndm(numcolors);
  if(GetPixel(hdc,ps->rcPaint.left-1,ps->rcPaint.top-1)==RGB(0,0,0))
    draw_grid(ps->rcPaint.left-1,ps->rcPaint.top-1,cs0);
  if(GetPixel(hdc,ps->rcPaint.right+1,ps->rcPaint.top-1)==RGB(0,0,0))
    draw_grid(ps->rcPaint.right+1,ps->rcPaint.top-1,cs1);
  if(GetPixel(hdc,ps->rcPaint.left-1,ps->rcPaint.bottom+1)==RGB(0,0,0))
    draw_grid(ps->rcPaint.left-1,ps->rcPaint.bottom+1,cs2);
  if(GetPixel(hdc,ps->rcPaint.right+1,ps->rcPaint.bottom+1)==RGB(0,0,0))
    draw_grid(ps->rcPaint.right+1,ps->rcPaint.bottom+1,cs3);

   subdiv(ps->rcPaint.left-1,ps->rcPaint.top-1, ps->rcPaint.right+1,
          ps->rcPaint.bottom+1);
}

/* Generate a plasma cloud */
void draw_plasma(int winno)
{
  int x,y;
  short cs0,cs1,cs2,cs3;

  hdc = GetDC(window[winno]);
  /* make points at (0,0) (0,CY) (CX,0) (CX,CY) of random colors */

/*  if ((grid = (short *)calloc((CX*CY),sizeof(short)))==NULL)
  {
    fprintf(stderr, "Can't Alloc enough memory for the plasma grid\n");
    exit(1);
  }
*/
  x=CX-1;
  y=CY-1;
  cs0=(short)rndm(numcolors);
  cs1=(short)rndm(numcolors);
  cs2=(short)rndm(numcolors);
  cs3=(short)rndm(numcolors);
  draw_grid(0,0,cs0);
  draw_grid(x,0,cs1);
  draw_grid(0,y,cs2);
  draw_grid(x,y,cs3);
  subdiv(0,0,x,y);

  ReleaseDC(window[winno],hdc);
/*  free(grid); */
}

void draw_grid(int x, int y, int colr)
{
/*  XDrawPoint(display,window,color_gcs[colr],x,y);*/
  SelectPalette(hdc,hPal,FALSE);
  SetPixel(hdc,x,y,PALETTEINDEX(colr));
/*  grid[x+y*CX]=colr; */
}

void draw_grid3d(int winno,int x,int y,int colr)
{
  project(&x,&y,&colr);
  XDrawPoint(display,window[winno],color_gcs[colr],x,y);
}

void subdiv(int x1, int y1, int x2, int y2)
{
 /* made these static so we won't need stack space for each iteration */
  COLORREF BLACK = RGB(0,0,0);
  int x,y;
  int horz=x2-x1;
  int vert=y2-y1;
  if (horz<2 && vert<2)
    return;

    x=(x1+x2)>>1;
    y=(y1+y2)>>1;

  if(GetPixel(hdc,x,y1)==RGB(0,0,0))
    draw_grid(x,y1,adjust(x1,y1,x2,y1));
  if(GetPixel(hdc,x2,y)==RGB(0,0,0))
    draw_grid(x2,y,adjust(x2,y1,x2,y2));
  if(GetPixel(hdc,x,y2)==RGB(0,0,0))
    draw_grid(x,y2,adjust(x1,y2,x2,y2));
  if(GetPixel(hdc,x1,y)==RGB(0,0,0))
    draw_grid(x1,y,adjust(x1,y1,x1,y2));

  if(GetPixel(hdc,x,y)==RGB(0,0,0))
  {
    short spooge;
    if(options.opt1) /* Mod 'em */
      spooge=( GetNearestPaletteIndex(hPal,GetPixel(hdc,x1,y1))+
               GetNearestPaletteIndex(hPal,GetPixel(hdc,x2,y2))+
               GetNearestPaletteIndex(hPal,GetPixel(hdc,x1,y2))+
               GetNearestPaletteIndex(hPal,GetPixel(hdc,x2,y1)) )%numcolors;
    else   /* Average 'em */
      spooge=( GetNearestPaletteIndex(hPal,GetPixel(hdc,x1,y1))+
               GetNearestPaletteIndex(hPal,GetPixel(hdc,x1,y2))+
               GetNearestPaletteIndex(hPal,GetPixel(hdc,x2,y1))+
               GetNearestPaletteIndex(hPal,GetPixel(hdc,x2,y2))
               )>>2;

    draw_grid(x,y,spooge);
  }

  subdiv(x1,y1,x,y);
  subdiv(x,y1,x2,y);
  subdiv(x1,y,x,y2);
  subdiv(x,y,x2,y2);
}

short
adjust(int x1, int y1, int x2, int y2)
{
  int c,d;
  int spoo;
  short clr1,clr2,clr3,clr4;
  short horz,vert;
  clr1=GetNearestPaletteIndex(hPal,GetPixel(hdc,x1,y1));
  clr2=GetNearestPaletteIndex(hPal,GetPixel(hdc,x2,y2));
  clr3=GetNearestPaletteIndex(hPal,GetPixel(hdc,x1,y2));
  clr4=GetNearestPaletteIndex(hPal,GetPixel(hdc,x2,y1));

  horz=x2-x1;
  vert=y2-y1;
  d=isqrt(horz*horz+vert*vert);
  /* leave the next 8 lines commented out for mellow colors...*/

  if(((spoo=rndm(100))>options.number)&&(d!=0))
  {
    if(spoo%2==0)
      c=(((clr1+clr2)>>1)+rndm(d))%numcolors;
    else
      c=(((clr1+clr2)>>1)-rndm(d))%numcolors;
  }
  else
     
/*    c=((clr1+clr2+clr3+clr4)>>2); */
      c=(clr1+clr2)>>1;
  
  if(c<0)
    c = abs(c); /* numcolors+c */
  else if(c==0)
    c = 1;

  return c;
}

/* This does a 3D projection... but since we don't use 3D for anything, I'll
   comment it out for now -jj */
  void project(int *x,int *y,int *z) 
  { 
  int zed,xx; /* 'cuz we don't we to change the value of z */
  static int inited=0;
  static double convert0= 0;
  static double convert1= M_PI_4; 
  static double convert2= M_PI_2;

   static double cosc0; 
   static double sinc0;  
   static double cosc1; 
   static double sinc1;
   static double cosc2;
   static double sinc2;
   if(!inited) 
   {
     cosc0=cos(convert0); 
     sinc0=sin(convert0); 
     cosc1=cos(convert1); 
     sinc1=sin(convert1); 
     cosc2=cos(convert2); 
     sinc2=sin(convert2); 
     inited=1; 
   } 

   zed= (*z*CY/numcolors)/2; 
 /*  *x *=(CX/32); 
   *y *=(CY/30); 
*/
 /* these are the 3-D ones  */
   xx=*x; /* save it */
   *x= ((*x*cosc0) + (*y*cosc1) + (zed*cosc2)); 
   *y= ((xx*sinc0) + (*y*sinc1) + (zed*sinc2));
} 

int isqrt(int n)
{
  /*  x*x-n=0 */
  int x,xsqu;
  int xh=256;
  int xl=4;
/* do the quick ones first */
  if (n<=0)
    return 0;
  else if (n<=2)
    return 1;
  else if(n<=6)
    return 2;
  else if (n<=12)
    return 3;
  else
    do
    {
      if(xh==xl)
        return xh;
      x=(xh+xl)>>1;
      xsqu=x*x;
      if(xsqu>n)
      {
        if(xh==x)
          return x;
        else
          xh=x;
      }
      else if (xsqu<n)
      {
        if(xl==x)
          return x;
        else
          xl=x;
      }
      else
      {
        return x;
      }
    }  while(xsqu!=n);
  return x;
}

int ihypot(int x, int y)
{
  int temp;
  temp = (x * x) + (y * y);
  return isqrt(temp);
}

XPoint PolartoRect(float radius, float angle)
{
  XPoint temp;
  
  temp.x=(int) (radius*cos(angle));
  temp.y=(int) (radius*sin(angle));
  
  return temp;
}


void draw_spiral(int winno)
{
  /* A spiral has the formula rad=ang */ 
  
  XPoint plotme;
  int midX,midY;
  float ang=0.0;
  short clr=numcolors-1;
  float end=hypot(CX,CY);
  midX= CX>>1;
  midY= CY>>1;
  do
  {
    plotme=PolartoRect(ang/3,ang);
    ang+=(float)M_PI/737;

/* a 3d spiral looks kinda neat, but not enough to keep */
/*    draw_grid3d(plotme.x+midX,plotme.y+midY,clr); */

    XDrawPoint(display,window[winno],color_gcs[clr],
               plotme.x+midX,plotme.y+midY);
    clr--;
    if(clr<=0)
      clr=numcolors-1;
  }  while((2*ang/3)<end);
} 

void
draw_mandel(int winno,int l,int t,int r, int b, RECT* rt )
{
  /* do funky shit with complex numbers */
  int x,y,i;
  double left,right,top,bottom;
  struct COMPLEX C;
  struct COMPLEX Z;

  top=(double)t/(CY>>1)-1.25;
  bottom=(double)b/(CY>>1)-.75;
  left=(double)l/(CX>>1)-1.875;
  right=(double)r/(CX>>1)-1.375;
/*   printf("Drawing Mandel from (%lf,%lf) to (%lf,%lf)\n", */
/*       ((right-left)*0/CX)+left, */
/*       ((bottom-top)*0/CY)+top, */
/*       ((right-left)*CX/CX)+left, */
/*       ((bottom-top)*CY/CY)+top); */

  for(y=rt->top;y<=rt->bottom;y++)
  {
#ifdef X
    if (XCheckMaskEvent(display,~0L,&event)==True) 
      handle_event(&event);
#else
  #ifdef WINDOWS
    if(CheckforMessage(window[winno]))
      return;
  #endif
#endif
    for(x=rt->top;x<rt->bottom;x++)
    {
      Z.real=0;  Z.imag=0;
      C.real=(double)((right-left)*x/CX)+left;
      C.imag=(double)((bottom-top)*y/CY)+top;
      for  (i=0;i<(options.number*10);i++)
      {
        Z=multComplex(Z,Z);
        Z=addComplex(Z,C);
        if((Z.real*Z.real)+(Z.imag*Z.imag)>4.0)
          break;
      }
      XDrawPoint(display,window[winno],
                 color_gcs[(i%numcolors)?(i%numcolors):1],x,y);
    }
  }
}

/* julia.. the fractal set, not the person */
void
draw_julia(int winno,RECT* r)
{
  /* do funky shit with complex numbers */
  int x,y,i; struct COMPLEX C;
  struct COMPLEX Z;
  for(y=r->top;y<r->bottom;y++)
  {
#ifdef X
    if(XCheckMaskEvent(display,~0L,&event)==True)       handle_event(&event);
#else
  #ifdef WINDOWS
    if(CheckforMessage(window[winno]))
      return;
  #endif
#endif
    for(x=r->left;x<r->right;x++)
    {
      Z.real=(double)(-(CX/2)+x)/(double)(CX/2.5);
      Z.imag=(double)(-(CY/2)+y)/(double)(CY/2.5);
      C.real =0.3;
      C.imag =0.6;
      /* F(z)=z^2+c */
      for(i=0;i<(options.number*10);i++)
      {
        Z=multComplex(Z,Z);
        Z=addComplex(Z,C);
        if((Z.real*Z.real)+(Z.imag*Z.imag)>4.0)
          break;
      }
      XDrawPoint(display,window[winno],color_gcs[(i%numcolors)+1],x,y);
    }
  }
}

/* newton.. the fractal set, not the person */
void
draw_newton(int val,int winno,RECT *rt)
{ 
  /*do funky shit with complex numbers */ 
  int x,y,i;
  struct COMPLEX Z;
  struct COMPLEX minus1={-1,0};
  struct COMPLEX expon;
  struct COMPLEX expDown;

  expon.imag=expDown.imag=0;
  expon.real=val; expDown.real= val-1;

  for(y=rt->top;y<rt->bottom;y++)
  {
#ifdef X
    if (XCheckMaskEvent(display,~0L,&event)==True) 
      handle_event(&event);
#else
  #ifdef WINDOWS
    if(CheckforMessage(window[winno]))
      return;
  #endif
#endif

    for(x=rt->left;x<rt->right;x++)
    {
      Z.real=(double)(-(CX/2)+x)/(double)(CX/2);
      Z.imag=(double)(-(CY/2)+y)/(double)(CY/2);
      /* z^val-1 = 0 */

      for (i=0;i<(options.number*10);i++)
      {
        double lastR=Z.real; double lastI=Z.imag;
/*      Z(j+1)=Z(i)-( (Z(i)^3-1)/(3*Z(i)^2)) */
/* Whoa! this is REAL MESSY! */
        Z= addComplex(Z,multComplex(minus1,(divComplex(
                                (addComplex(powComplex(Z,expon),minus1)),
                                (multComplex(expon,powComplex(Z,expDown)))))));

        if((fabs(Z.real-1.0)<0.01)&&(fabs(Z.imag)<0.01))
        {
          /*  printf("Breaking out "); */
          break;
        }
/*
  fprintf(stderr,"Iter[%d] Z.real = %lf \t Z.imag= %lf\n",
  i,Z.real,Z.imag);
  */
        if(Z.real==lastR && Z.imag==lastI) /* we're stuck.. quit now */
        {
          i=10*options.number;
          break;
        }
      }
      XDrawPoint(display,window[winno],color_gcs[(i%numcolors)+1],x,y); 
    }
  }
}

/* from the sci.fractals FAQ 
   Q9a: How does complex arithmetic work?
   A9a: It works mostly like regular algebra with a couple additional formulas:
   (note: a,b are reals, x,y are complex, i is the square root of -1)
   Powers of i: i^2 = -1
   Addition: (a+i*b)+(c+i*d) = (a+c)+i*(b+d)
   Multiplication: (a+i*b)*(c+i*d) = a*c-b*d + i*(a*d+b*c)
   Division: (a+i*b)/(c+i*d) = (a+i*b)*(c-i*d)/(c^2+d^2)
   Exponentiation: exp(a+i*b) = exp(a)(cos(b)+i*sin(b))
   Sine: sin(x) = (exp(i*x)-exp(-i*x))/(2*i)
   Cosine: cos(x) = (exp(i*x)+exp(-i*x))/2
   Magnitude: |a+i*b| = sqrt(a^2+b^2)
   Log: log(a+i*b) = log(|a+i*b|)+i*arctan(b/a)  (Note: log is multivalued.)
   Complex powers: x^y = exp(y*log(x))
   DeMoivre's theorem: x^a = r^a * [cos(a*theta) + i * sin(a*theta)]
*/


struct COMPLEX divComplex(struct COMPLEX foo, struct COMPLEX bar)
{
  struct COMPLEX temp;
  temp.real = ((foo.real*bar.real)+(foo.imag*bar.imag))/(bar.real*bar.real+
                                                         bar.imag*bar.imag);
  temp.imag = ((foo.imag*bar.real)-(foo.real*bar.imag))/(bar.real*bar.real+
                                                         bar.imag*bar.imag);
  return temp;
}

struct COMPLEX logComplex(struct COMPLEX foo)
{
  struct COMPLEX temp;
  double gaaah=fabs(foo.real+foo.imag);
  if(gaaah<=0.0)
    temp.real=-10;
  else
    temp.real = log10(gaaah);
  temp.imag = atan2(foo.imag,foo.real);
  return temp;
}

struct COMPLEX powComplex(struct COMPLEX foo, struct COMPLEX bar)
{
  struct COMPLEX temp;
  return (expComplex(multComplex(bar,logComplex(foo))));
}

struct COMPLEX expComplex(struct COMPLEX foo)
{
  struct COMPLEX temp;
  temp.real = exp(foo.real)*(cos(foo.imag));
  temp.imag = exp(foo.real)*sin(foo.imag);
  return temp;
}

struct COMPLEX addComplex(struct COMPLEX foo, struct COMPLEX bar)
{
  struct COMPLEX temp;
  temp.real = foo.real + bar.real;
  temp.imag = foo.imag + bar.imag;
  return temp;
}

struct COMPLEX multComplex(struct COMPLEX foo ,struct COMPLEX bar)
{
  struct COMPLEX temp;
  temp.real = (foo.real * bar.real) - (foo.imag * bar.imag);
  temp.imag = (foo.imag * bar.real) + (foo.real * bar.imag);
  return temp;
}

void wander(int winno,int *x, int *y)
{
  int spoo;

  XDrawPoint(display,window[winno],color_gcs[HC[0]],*x,*y);
  if((spoo=rndm(3))==1)
    *x+=1;
  else if (spoo==2)
    *x-=1;
  if((spoo=rndm(3))==1)
    *y+=1;
  else if (spoo==2)
    *y-=1;
  if (*x<=0)
    *x=CX;
  else if (*x>=CX)
    *x=0;
  if (*y<=0)
    *y=CY;
  else if (*y>=CY)
    *y=0;
  if(options.multi)
    HC[0]=(HC[0]+1)%numcolors;
}

#ifndef X
void draw_munch(int spoo,int winno)
{
  int x=0;
  int foo;
  for(x=0;x<CX;x++)
  {
    HDC hdc = GetDC(window[winno]);
    foo=(GetNearestPaletteIndex(hPal,(GetPixel(hdc,x,(x^spoo)%CY)<<1))+1)%numcolors;
    ReleaseDC(window[winno],hdc);
    if(foo==0||foo>numcolors) foo=1;
    XDrawPoint(0,window[winno],x,(x^spoo++)%CY,foo);
  }
}
#endif

void draw_funky(int winno)
{
  /* pick a point and a line length
     rotate the line around the midpoint
     decrease the line lenth by 1
     increase the color
     until the length = 0
   */
  int x1,y1;
  int howfunky=options.number;
  XPoint plotme,plotme2;
  float length=hypot(CX,CY); /*rndm(100);*/
  int x=CX/2; /* rndm(CX); */
  int y=CY/2;/* rndm(CY); */
  float ang=rndm(16)*M_PI/8;
  int clr=rndm(numcolors);
  while(length>0)
  {
    plotme=PolartoRect(length/2,ang);

    ang+=(float)M_PI/1024;
    if(howfunky==1)
    {
      plotme2=PolartoRect(length/2,ang+M_PI_2);
      XDrawLine(display,window[winno],clr,x-plotme.x,y-plotme.y,
                x-plotme2.x,y-plotme2.y);
      XDrawLine(display,window[winno],clr,x-plotme2.x,y-plotme2.y,
                x+plotme.x,y+plotme.y);
      XDrawLine(display,window[winno],clr,x+plotme.x,y+plotme.y,
                x+plotme2.x,y+plotme2.y);
      XDrawLine(display,window[winno],clr,x+plotme2.x,y+plotme2.y,
                x-plotme.x,y-plotme.y);
      length-=.5;
    }
    else if (howfunky==2)
    {
      XPoint plotme3;
      plotme2=PolartoRect(length/2,ang+(M_PI/3));
      plotme3=PolartoRect(length/2,ang+(2*(M_PI/3)));
      XDrawLine(display,window[winno],clr,CX-(x-plotme.x),CY-(y-plotme.y),
                CX-(x-plotme2.x),CY-(y-plotme2.y));
      XDrawLine(display,window[winno],clr,CX-(x-plotme2.x),CY-(y-plotme2.y),
                CX-(x-plotme3.x),CY-(y-plotme3.y));
      XDrawLine(display,window[winno],clr,CX-(x-plotme3.x),CY-(y-plotme3.y),
                CX-(x-plotme.x),CY-(y-plotme.y));

      XDrawLine(display,window[winno],clr,x-plotme.x,y-plotme.y,
                x-plotme2.x,y-plotme2.y);
      XDrawLine(display,window[winno],clr,x-plotme2.x,y-plotme2.y
                ,x-plotme3.x,y-plotme3.y);
      XDrawLine(display,window[winno],clr,x-plotme3.x,y-plotme3.y,
                x-plotme.x,y-plotme.y);

      length-=.2;
    }
    else
    {
      XDrawLine(display,window[winno],clr,x-plotme.x,y-plotme.y,
                x+plotme.x,y+plotme.y);
      length-=.05;
    }
    clr=(clr+1)%numcolors;
    if(clr<1)
      clr=1;
  }
  return;
}

void
draw_static(int winno,RECT* r)
{
  int x,y;

  for(y=r->top;y<=r->bottom;y++)
    for(x=r->left;x<r->right;x++)
    {
      XDrawPoint(display,window[winno],color_gcs[rndm(numcolors-1)+1],x,y);
    }
}


void
draw_max(int winno)
{
/* Point 0 is the center of the screen, initially */
/* Point 1-3 are the 3 points of the cube */
  static XPoint **pts;
  static int initd=0;
  static unsigned int *length;
  
  if(!initd)
  {
    int i;
      
    pts=(XPoint **)malloc(1*sizeof(XPoint*));
    length=(unsigned int *)malloc(1*sizeof(unsigned int));
    for(i=0;i<1;i++)
    {
      double ang;
      pts[i]=(XPoint *)calloc(4,sizeof(XPoint));
      pts[i][0].x=CX>>1;
      pts[i][0].y=CY>>1;
      pts[i][1].x=CX;
      pts[i][1].y=CY;
      pts[i][2].x=0;
      pts[i][2].y=CY;
      pts[i][3].x=CX>>1;
      pts[i][3].y=0;
      length[i]=hypot(pts[i][1].x-pts[i][0].x,
                      pts[i][1].y-pts[i][0].y);
      ang=atan2(pts[i][1].y-pts[i][0].y, pts[i][1].x-pts[i][0].x);
/*      pts[i][2]=PolartoRect(length[i],ang+M_PI*.6667);  PI*2/3 */
/*      pts[i][3]=PolartoRect(length[i],ang+M_PI*1.3333);  PI*4/3 */
    }
    initd=1;
  }
  
  XDrawLine(display,window[winno],color_gcs[1],
            pts[winno][0].x, pts[winno][0].y,
            pts[winno][1].x, pts[winno][1].y);
  
  XDrawLine(display,window[winno],color_gcs[numcolors/3],
            pts[winno][0].x, pts[winno][0].y,
            pts[winno][2].x, pts[winno][2].y);
  
  XDrawLine(display,window[winno],color_gcs[numcolors/2],
            pts[winno][0].x, pts[winno][0].y,
            pts[winno][3].x, pts[winno][3].y);

  if(pts[winno][0].x++>CX)
    pts[winno][0].x=0;
  else if(pts[winno][0].x<0)
    pts[winno][0].x=CX-1;
  if(pts[winno][0].y++>CY)
    pts[winno][0].y=0;
  else if (pts[winno][0].y<0)
    pts[winno][0].y=CY-1;

}

/* The paint mixer */
void
draw_mixer(int winno,RECT* r)
{
  int x,y,X,Y,Color;
  long HYP;

  HYP = (CX*CX+CY*CY)>>2;
  for(y=r->top;y<r->bottom;y++)
  {
    Y=y-(CY/2);
//      if (XCheckMaskEvent(display,~0L,&event)==True)
//      handle_event(&event);
    for(x=r->left;x<r->right;x++)
    {
      X=x-(CX/2);
      Color=((X==0)&&(Y==0)) ? 0 :
        (((int)floor((((atan2((double)Y,(double)X)/M_PI)+1)/2+
                      (1+cos(M_PI*(double)(X*X+Y*Y)/HYP*2))/1)*
                     (numcolors-1)))%(numcolors-1))+1;

      XDrawPoint(display,window[winno],color_gcs[Color],x,y); 
    }
  }
}

/* The cells */
void
draw_cells(int winno,RECT* r)
{
  int x,y,Color;

  for(y=r->top;y<r->bottom;y++)
    {
//      if (XCheckMaskEvent(display,~0L,&event)==True)
//      handle_event(&event);
      for(x=r->left;x<r->bottom;x++)
      {
        Color=(((int)((sin(M_PI*(double) x/80) *sin(M_PI*(double) y/80) + 1)/2*
                      (numcolors-1)))%(numcolors-1))+1;
        
        XDrawPoint(display,window[winno],color_gcs[Color],x,y); 
      }
    }
  return;
}