Path: csus.edu!wuarchive!julius.cs.uiuc.edu!news.cs.indiana.edu!rutgers!mcnc!uvaarpa!haven!umd5!galileo.cs.jhu.edu!callahan
From: callahan@cs.jhu.edu (Paul Callahan)
Newsgroups: alt.fractals.pictures
Subject: Re: Random Circular Gasket
Message-ID: <callahan.660631086@galileo.cs.jhu.edu>
Date: 8 Dec 90 04:36:25 GMT
References: <callahan.660507999@newton.cs.jhu.edu>
Reply-To: callahan@cs.jhu.edu (Paul Callahan)
Organization: JHU Computer Science Deparment, Baltimore MD
Lines: 794

In article <callahan.660507999@newton.cs.jhu.edu> I write:
>Here's a fractal I generated a while back (I can post source if there's
>sufficient interest). 

I got a reasonable amount of e-mail asking for this, so here it is, in
shar form.  The programs given will generate either Postscript or
an X bitmap (I made the gif using one of the commonly available conversion
filters).  This is also intended to satisfy requests for the Postscript version
of the picture.  However, if you cannot install this program on your system, 
I can still send you a Postscript file.  Again, this file (with 20000 circles) 
is quite large (and doesn't compress particularly well) so I would prefer to 
send it only if necessary.

If there are any problems with installation, feel free to send mail. 

--
Paul Callahan
callahan@cs.jhu.edu

---- cut here ----
# This is a shell archive.  Remove anything before this line,
# then unpack it by saving it in a file and typing "sh file".
#
# Wrapped by newton.cs.jhu.edu!callahan on Fri Dec  7 23:17:54 EST 1990
# Contents:  README f1.inp f2.inp f3.inp hap.c hd.ps makefile plothap.c tl.ps
#	hap2ps
 
echo x - README
sed 's/^@//' > "README" <<'@//E*O*F README//'
To install these programs, type "make" to make the object files.
This will create two binaries, hap and plothap, which perform the following 
functions:

hap:
   Creates a random pattern of non-overlapping circles which fill up 
   a circular disk.

   Input format is as follows:

   line 1:  random seed number (integer)
   line 2:  number of circles to draw (integer)
   line 3:  ratio of actual circle drawn 
            to largest one that will fit (floating point, between 0 and 1)
   line 4:  Postscript/numeric output flag (n or p)

plothap:
   Takes _numeric_ output from hap and produces an X bitmap (.xbm) file.
   Command-line arguments are numeric, and are optional.  
 
   "plothap" produces a bitmap of default size (300x300), with random
   happy-face orientations based on seed determined from the system 
   date function.

   "plothap <size>" produces a bitmap of user-specified dimension with
   seed determined from date.

   "plothap <size> <seed>" produces a bitmap with user-specified 
   dimension and seed number.

The other executable is a _very_ simple shell script that places a Postscript
header and trailer on the Postscript output of hap.

WARNING ABOUT POSTSCRIPT FILES:
-------------------------------
Postscript files containing many thousands of happy faces can take well over an
hour to print on certain Postscript printers, such as Apple laser-writers.
Be sure that you choose a time when others do not need the printer.  (In
spite of the sub-optimal nature of the algorithm used to generate the 
pictures, this turns out to be the main bottleneck in practice).

Examples:
---------
To create a Postscript file type:

    hap <f3.inp | hap2ps >f3.ps

where f3.inp is an input file containing a "p" on the last line.

To create a 500x500 X bitmap type:  
    
    hap <f1.inp | plothap 500 >f1.xbm

where f1.inp is an input file containing a "n" on the last line. 

Both of the above commands should work verbatim using the enclosed
input files.

Overview of Method:
-------------------
This program works by choosing a random sequence of non-overlapping 
circles.  At each step, a random center is chosen as the location of the
circle to be added, and distance from this point to the nearest circle is
determined (if the location is inside of a circle, we pick a new point).
This distance gives the maximum radius of a circle centered at the chosen
location.  In fact, the new circle we actually add has a radius that is
smaller than the maximum by a fixed percentage (the third input parameter).

A search tree is maintained throughout the program to make it possible to
find the nearest circle relatively efficiently.  This tree is not balanced,
however, and in some cases it is necessary to search more than one of the
four subtrees of a given node at some point in the search.  For this 
reason, one cannot guarantee that the time complexity will be better
than O(n^2) in the worst case.  However, the randomness of the data causes the
tree to be relatively bushy, and in practice it seems to take much less
time than the worst case.  For example, 20000 circles is quite feasible on
a low-power workstation of some kind.

It would be nice if the average time were O(n log n) (where n is the number of 
circles).  I suspect that it may be as high as O(n^(3/2)) but I haven't 
collected any empirical data.  I wrote this several years ago, before I had 
studied any computational geometry.  The problem could be treated in terms of 
the on-line generation of Voronoi diagrams of sets of circles (as opposed to 
the more standard point sets).  I have heard that the off-line construction of 
such a Voronoi diagram can be done in O(n log^2 n), but I have no references 
handy.  This line of attack may yield a better method of constructing random 
circular gaskets.  In any case, this program seems reasonably effective over 
a significant range of values of n.
@//E*O*F README//
chmod u=rw,g=r,o=r README
 
echo x - f1.inp
sed 's/^@//' > "f1.inp" <<'@//E*O*F f1.inp//'
63872935
5000
0.9
n
@//E*O*F f1.inp//
chmod u=rw,g=rw,o=rw f1.inp
 
echo x - f2.inp
sed 's/^@//' > "f2.inp" <<'@//E*O*F f2.inp//'
63872935
20000
0.9
p
@//E*O*F f2.inp//
chmod u=rw,g=r,o=r f2.inp
 
echo x - f3.inp
sed 's/^@//' > "f3.inp" <<'@//E*O*F f3.inp//'
65882935
5000
0.9
p
@//E*O*F f3.inp//
chmod u=rw,g=r,o=r f3.inp
 
echo x - hap.c
sed 's/^@//' > "hap.c" <<'@//E*O*F hap.c//'
#include <stdio.h>
#include <math.h>
#define MAXCIRC 30000
#define TRUE 1
#define FALSE 0
struct treenode {
      struct treenode *pt[2][2];  
      int num;
      float mor;
   };

struct cmptreenode {
      struct treenode *tree;
      float dist;
      struct cmptreenode *left;
      struct cmptreenode *right;
   };


int i,n;
int sd1;
float cx[MAXCIRC-1];
float cy[MAXCIRC-1];
float cr[MAXCIRC-1];
float x,y,rmax;
float ratio;
struct treenode *tree;
int inside;
int cnter,try;
struct cmptreenode *cmproot;
struct cmptreenode *freecmp;
char outform;

float rndnum()
{
    return (float)(random()%(10000000))/(float)10000000;
}

freetree(root)
struct cmptreenode *root;
{
   if (root!=NULL) {
      freetree(root->left);
      freetree(root->right);
      root->left = freecmp;
      freecmp = root;
   }
}

getfree(nodeptr)
struct cmptreenode **nodeptr;
{
   if (freecmp==NULL) 
      *nodeptr = (struct cmptreenode *) malloc(sizeof(struct cmptreenode));
   else {
      *nodeptr=freecmp;
      freecmp=freecmp->left;
   }
}

inscomp(root, tree, dist)
struct cmptreenode **root;
struct treenode *tree;
float dist;
{
    if (*root==NULL) {
        getfree(root);
        (*root)->left = NULL;
        (*root)->right = NULL;
        (*root)->tree = tree;
        (*root)->dist = dist;
    }
    else if (dist < (*root)->dist)
        inscomp(&((*root)->left),tree,dist);
    else
        inscomp(&((*root)->right),tree,dist);
}

float max(a,b)
float a,b;
{
   if (a>b) return a;
   else return b;
}

inscirc(i, x, y, r, root)
int i;
float x,y,r;
struct treenode **root;

{
int p1,p2;
float qx,qy;
float orad;

   if (*root==NULL) { 
         *root = (struct treenode *) malloc(sizeof(struct treenode));
         (*root)->num = i;
         cx[(*root)->num] = x;
         cy[(*root)->num] = y;
         cr[(*root)->num] = r;
         (*root)->mor = r;
         (*root)->pt[0][0] = NULL;
         (*root)->pt[0][1] = NULL;
         (*root)->pt[1][0] = NULL;
         (*root)->pt[1][1] = NULL;
   }
   else {
         qx = x-cx[(*root)->num];
         qy = y-cy[(*root)->num];
         orad = sqrt(qx*qx+qy*qy)+r;
	 (*root)->mor = max((*root)->mor,orad);
         if (x<cx[(*root)->num]) p1 = 0; else p1 = 1;
         if (y<cy[(*root)->num]) p2 = 0; else p2 = 1;
	 inscirc(i,x,y,r,&((*root)->pt[p1][p2]));
   }
}


compcirc(x, y, root, rmax, inside, cnter)
float x,y;
struct treenode *root;
float *rmax;
int *inside;
int *cnter;

{
int i;
int p1,p2;
float dist;
float qx,qy;
struct treenode *subtree;

   qx = x-cx[root->num];
   qy = y-cy[root->num];
   dist = sqrt(qx*qx+qy*qy)-cr[root->num];

   if (dist < 0) *inside  =  TRUE; 
   else
      if (dist < *rmax) *rmax = dist;

   (*cnter)++;

   i = 0;
   for (i = 0; (i <= 3) && (!*inside); i++) {
      p1 = (i / 2) % 2;
      p2 = i % 2;
      subtree = root->pt[p1][p2];
      if (subtree!=NULL) { 
         qx = x-cx[subtree->num];
         qy = y-cy[subtree->num];
         dist = sqrt(qx*qx+qy*qy)-subtree->mor;
         if (dist < 0) 
               compcirc(x,y,subtree,rmax,inside,cnter);
         else if (dist < *rmax) 
               inscomp(&cmproot,subtree,dist);
      }
   }
}

recompcirc(x,y,root,rmax,cnter)
float x,y;
struct treenode *root;
float *rmax;
int *cnter;

{
int i;
int p1,p2;
float dist;
float qx,qy;
struct treenode *subtree;

   qx = x-cx[root->num];
   qy = y-cy[root->num];
   dist = sqrt(qx*qx+qy*qy)-cr[root->num];

   if (dist < *rmax) *rmax = dist;

   (*cnter)++;

   for (i = 0; i <= 3; i++) {
      p1 = (i / 2) % 2;
      p2 = i % 2;
      subtree = root->pt[p1][p2];
      if (subtree!=NULL) { 
         qx = x-cx[subtree->num];
         qy = y-cy[subtree->num];
         dist = sqrt(qx*qx+qy*qy)-subtree->mor;
         if (dist < *rmax) 
               recompcirc(x,y,subtree,rmax,cnter);
      }
   }
}

recompall(x,y,rmax,cnter,cmproot)
float x,y;
float *rmax;
int *cnter;
struct cmptreenode *cmproot;
{
    if (cmproot!=NULL) {
       recompall(x,y,rmax,cnter,cmproot->left);
       if (cmproot->dist < *rmax) {
           recompcirc(x,y,cmproot->tree,rmax,cnter);
           recompall(x,y,rmax,cnter,cmproot->right);
       }
    }
}


dumptree(root)
struct cmptreenode *root;
{
    if (root!=NULL) {
       dumptree(root->left);
       printf("%6d %18.10f\n",root->tree->num,root->dist);
       dumptree(root->right);
    }
}

main()
{
    /*get seed number for random number generator*/
    scanf("%d\n",&sd1);
    srandom(sd1);
    /*get number of circles to draw*/
    scanf("%d\n",&n);
    /*get ratio of circle maximum size*/
    scanf("%f\n",&ratio);
    /*get output format (if p then print postscript plot commands,
                        otherwise numerical data is printed)*/
    outform = getchar();
    while (outform != 'p' && outform != EOF) outform = getchar();

    if (outform == 'p') {
        printf("0 0 10000 10000 s\n");
        printf("5000 5000 5000 c\n");
    }
    tree = NULL;
    freecmp = NULL;
    x = rndnum();
    y = rndnum();
    while ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)>0.25) { 
       x = rndnum();
       y = rndnum();
    }

    rmax = 0.5-sqrt((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5));
    i = 1;
    inscirc(i,x,y,rmax*ratio,&tree);
    cnter  =  0;
    try  =  1;
    if (outform == 'p') 
	printf("%d %d %d %d hap\n",
                (int) (cx[i]*10000),
                (int) (cy[i]*10000),
                (int) (cr[i]*10000),
                (int) (360*rndnum()));
    else
        printf("%6d %10.8f %10.8f %10.8f %4d %4d\n",
         	i,cx[i],cy[i],cr[i],cnter,try);
    for (i = 2; i <= n; i++) {
      /*find point not inside a circle and distance from nearest circle*/
      inside = TRUE;
      try = 0;
      while (inside) { 
         x = rndnum();
         y = rndnum();
         while ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)>0.25) {
             x = rndnum();
             y = rndnum();
         }
         rmax = 0.5-sqrt((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5));

         /*find distance from nearest circle*/

         cnter = 0;
         inside = FALSE;
         cmproot = NULL;
         compcirc(x,y,tree,&rmax,&inside,&cnter);

         /*Activate next line to get dump of recomparison tree*/
         /* if (cmproot != NULL) printf("---------");
           dumptree(cmproot);     */
	 if (! inside) 
             recompall(x,y,&rmax,&cnter,cmproot);
         freetree(cmproot);

         try++;
      }
      inscirc(i,x,y,rmax*ratio,&tree);
      if (outform == 'p') 
         printf("%d %d %d %d hap\n", (int) (cx[i]*10000),
                  (int) (cy[i]*10000),
                  (int) (cr[i]*10000),
                  (int) (360*rndnum()));
      else
        printf("%6d %10.8f %10.8f %10.8f %4d %4d\n",
        	i,cx[i],cy[i],cr[i],cnter,try);
   } 
}
@//E*O*F hap.c//
chmod u=rw,g=r,o=r hap.c
 
echo x - hd.ps
sed 's/^@//' > "hd.ps" <<'@//E*O*F hd.ps//'
%!
%%Creator: gondor.cs.psu.edu:callahan(Paul B. Callahan,22W or 23W,31943,2340258)
%%Title: Unix plot file
%%CreationDate: Mon Mar  7 19:05:52 1988
%%DocumentFonts: Courier
%%EndComments
% lib/psplot.pro -- included prolog for plot files
% Copyright (c) 1984 Adobe Systems, Inc.  All Rights Reserved.
save 50 dict begin /psplot exch def
/StartPSPlot
   {newpath 0 0 moveto 20 setlinewidth 0 setgray 1 setlinecap
    /imtx matrix currentmatrix def /dmtx matrix defaultmatrix def
    /fnt /Courier findfont def /smtx matrix def fnt 8 scalefont setfont}def
/solid {{}0}def
/dotted	{[2 nail 10 nail ] 0}def
/longdashed {[10 nail] 0}def
/shortdashed {[6 nail] 0}def
/dotdashed {[2 nail 6 nail 10 nail 6 nail] 0}def
/min {2 copy lt{pop}{exch pop}ifelse}def
/max {2 copy lt{exch pop}{pop}ifelse}def
/len {dup mul exch dup mul add sqrt}def
/nail {0 imtx dtransform len 0 idtransform len}def

/m {newpath moveto}def
/n {lineto currentpoint stroke moveto}def
/p {newpath moveto gsave 1 setlinecap solid setdash
    dmtx setmatrix .4 nail setlinewidth
    .05 0 idtransform rlineto stroke grestore}def
/l {moveto lineto currentpoint stroke moveto}def
/t {smtx currentmatrix pop imtx setmatrix show smtx setmatrix}def
/a {gsave newpath /y2 exch def /x2 exch def 
    /y1 exch def /x1 exch def /yc exch def /xc exch def
    /r x1 xc sub dup mul y1 yc sub dup mul add sqrt
       x2 xc sub dup mul y2 yc sub dup mul add sqrt add 2 div def
    /ang1 y1 yc sub x1 xc sub atan def
    /ang2 y2 yc sub x2 xc sub atan def
    xc yc r ang1 ang2 arc stroke grestore}def
/hap {gsave newpath /ori exch def /r exch def /y exch def /x exch def
      x y r 0 360 arc stroke
      /re r 0.5 mul def
      /xe1 45 ori add cos re mul x add def
      /ye1 45 ori add sin re mul y add def
      /xe2 135 ori add cos re mul x add def
      /ye2 135 ori add sin re mul y add def
      xe1 ye1 r 0.04 mul 0 360 arc stroke 
      xe2 ye2 r 0.04 mul 0 360 arc stroke
      x y 0.67 r mul ori 202 add ori 338 add arc stroke grestore}def
/c {gsave newpath 0 360 arc stroke grestore}def
/cf {gsave newpath 0 360 arc fill grestore}def
/e {gsave copypage copypage showpage grestore newpath 0 0 moveto}def
/f {load exec setdash}def
/s {/ury exch def /urx exch def /lly exch def /llx exch def
    imtx setmatrix newpath clippath pathbbox newpath
    /dury exch def /durx exch def /dlly exch def /dllx exch def
    /md durx dllx sub dury dlly sub min def
    /Mu urx llx sub ury lly sub max def
    dllx dlly translate md Mu div dup scale llx neg lly neg translate}def
/EndPSPlot {clear psplot end restore}def
% end fixed prolog

StartPSPlot
%%EndProlog
%%Page: 1 1
@//E*O*F hd.ps//
chmod u=rw,g=r,o=r hd.ps
 
echo x - makefile
sed 's/^@//' > "makefile" <<'@//E*O*F makefile//'

all: hap plothap

hap: hap.c
	cc hap.c -o hap -lm

plothap: plothap.c
	cc plothap.c -o plothap -lm

@//E*O*F makefile//
chmod u=rw,g=r,o=r makefile
 
echo x - plothap.c
sed 's/^@//' > "plothap.c" <<'@//E*O*F plothap.c//'
#include <stdio.h>
#include <math.h>
#define BMAX 1200
#define DEFBSIZE 300
#define A1 (double)(3.14159/4)
#define A2 (double)(3.14159*3/4)

int bsize=DEFBSIZE;

main(argc, argv)
int argc; 
char *argv[];
{
static char bitmap[BMAX][BMAX];
char inpline[100];
int j1,j2,j3,smrad,sd;
double x,y,r,angle;

if (argc>2) {
  sscanf(argv[2],"%d",&sd);
  srandom (sd);
}  
else srandom ((int) time (NULL));

if (argc>1) {
  sscanf(argv[1],"%d",&bsize);
  if (bsize>BMAX) { 
    printf("Bitmap size too large (max is %d)\n",BMAX);
    exit(-1);
  }
}  

clearmap(bitmap);

while(gets(inpline)) {
   angle=((double)(random()%628318))/100000;
   sscanf(inpline,"%d %lf %lf %lf %d %d",&j1,&x,&y,&r,&j2,&j3);
   plotc(bitmap,(int)(bsize*x),(int)(bsize*y),(int)(bsize*r));
   smrad=(int)(bsize*r*0.67);
   plotsm(bitmap,(int)(bsize*x),(int)(bsize*y),smrad,angle);
   if (smrad>30) plotsm(bitmap,(int)(bsize*x),(int)(bsize*y),smrad-1,angle);
   plotcf(bitmap,(int)(bsize*x+bsize*r*cos(A1+angle)*0.5),
                (int)(bsize*y+bsize*r*sin(A1+angle)*0.5),(int)(bsize*r*0.08));
   plotcf(bitmap,(int)(bsize*x+bsize*r*cos(A2+angle)*0.5),
                (int)(bsize*y+bsize*r*sin(A2+angle)*0.5),(int)(bsize*r*0.08));
}

printmap(bitmap);
exit(0);
}

clearmap(bitmap)
char bitmap[][BMAX];
{
int i,j;

for (i=0; i<bsize; i++)
   for (j=0; j<bsize; j++) bitmap[i][j]='-';
}

printmap(bitmap)
char bitmap[][BMAX];
{
int i,j,k;
unsigned int byte;

printf("#define _width %d\n",bsize);
printf("#define _height %d\n\n",bsize);

printf("static char _bits[] = {");

for (j=0; j<bsize; j++) {
   for (i=0; i<=(int)(bsize/8)*8-8; i+=8) {
       byte=0;
       for (k=7; k>=0; k--) byte=(byte<<1)+(bitmap[i+k][j]=='#');
       printbyte(byte);
   }
   if (i<bsize) {
       byte=0;
       for (i=bsize-1; i>=(int)(bsize/8)*8; i--) 
              byte=(byte<<1)+(bitmap[i][j]=='#');
       printbyte(byte);
   }
}

printf(" };\n");
}

plot8p(bitmap,cx,cy,x,y)
char bitmap[][BMAX];
int cx,cy,x,y;
{
    bitmap[(cx+x+bsize)%bsize][(cy+y+bsize)%bsize]='#';
    bitmap[(cx-x+bsize)%bsize][(cy+y+bsize)%bsize]='#';
    bitmap[(cx+x+bsize)%bsize][(cy-y+bsize)%bsize]='#';
    bitmap[(cx-x+bsize)%bsize][(cy-y+bsize)%bsize]='#';
    bitmap[(cx+y+bsize)%bsize][(cy+x+bsize)%bsize]='#';
    bitmap[(cx-y+bsize)%bsize][(cy+x+bsize)%bsize]='#';
    bitmap[(cx+y+bsize)%bsize][(cy-x+bsize)%bsize]='#';
    bitmap[(cx-y+bsize)%bsize][(cy-x+bsize)%bsize]='#';
}


plotc(bitmap,cx,cy,r)
char bitmap[][BMAX];
int cx,cy,r;
{
int x,y,p;

x=0;
y=r;
p=3-2*r;

while (x<y) {
    plot8p(bitmap,cx,cy,x,y);
    if (p<0) p=p+4*x+6;
    else {
      p=p+4*(x-y)+10;
      y=y-1;
    }
  x=x+1;
}
if (x==y) plot8p(bitmap,cx,cy,x,y);
}

plot8pf(bitmap,cx,cy,x,y)
char bitmap[][BMAX];
int cx,cy,x,y;
{
int yi;
for(yi=y; yi>=0; yi--) {
    bitmap[(cx+x+bsize)%bsize][(cy+yi+bsize)%bsize]='#';
    bitmap[(cx-x+bsize)%bsize][(cy+yi+bsize)%bsize]='#';
    bitmap[(cx+x+bsize)%bsize][(cy-yi+bsize)%bsize]='#';
    bitmap[(cx-x+bsize)%bsize][(cy-yi+bsize)%bsize]='#';
}
for(yi=x; yi>=0; yi--) {
    bitmap[(cx+y+bsize)%bsize][(cy+yi+bsize)%bsize]='#';
    bitmap[(cx-y+bsize)%bsize][(cy+yi+bsize)%bsize]='#';
    bitmap[(cx+y+bsize)%bsize][(cy-yi+bsize)%bsize]='#';
    bitmap[(cx-y+bsize)%bsize][(cy-yi+bsize)%bsize]='#';
}
}


plotcf(bitmap,cx,cy,r)
char bitmap[][BMAX];
int cx,cy,r;
{
int x,y,p;

x=0;
y=r;
p=3-2*r;

while (x<y) {
    plot8pf(bitmap,cx,cy,x,y);
    if (p<0) p=p+4*x+6;
    else {
      p=p+4*(x-y)+10;
      y=y-1;
    }
  x=x+1;
}
if (x==y) plot8pf(bitmap,cx,cy,x,y);
}

plot8ps(bitmap,cx,cy,x,y,s1,s2,c)
char bitmap[][BMAX];
int cx,cy,x,y,s1,s2,c;
{
 if ( s1*x+s2*y<=c) bitmap[(cx+x+bsize)%bsize][(cy+y+bsize)%bsize]='#';
 if (-s1*x+s2*y<=c) bitmap[(cx-x+bsize)%bsize][(cy+y+bsize)%bsize]='#';
 if ( s1*x-s2*y<=c) bitmap[(cx+x+bsize)%bsize][(cy-y+bsize)%bsize]='#';
 if (-s1*x-s2*y<=c) bitmap[(cx-x+bsize)%bsize][(cy-y+bsize)%bsize]='#';
 if ( s1*y+s2*x<=c) bitmap[(cx+y+bsize)%bsize][(cy+x+bsize)%bsize]='#';
 if (-s1*y+s2*x<=c) bitmap[(cx-y+bsize)%bsize][(cy+x+bsize)%bsize]='#';
 if ( s1*y-s2*x<=c) bitmap[(cx+y+bsize)%bsize][(cy-x+bsize)%bsize]='#';
 if (-s1*y-s2*x<=c) bitmap[(cx-y+bsize)%bsize][(cy-x+bsize)%bsize]='#';
}

plota(bitmap,cx,cy,r,a1,a2)
char bitmap[][BMAX];
int cx,cy,r;
double a1,a2;
{
int x,y,p,x1,y1,x2,y2,s1,s2,c;

x=0;
y=r;
p=3-2*r;
x1=(int)((double)r*cos(a1));
y1=(int)((double)r*sin(a1));
x2=(int)((double)r*cos(a2));
y2=(int)((double)r*sin(a2));
s1=y1-y2;
s2=x2-x1;
c= -x1*y2+y1*x2;

while (x<y) {
    plot8ps(bitmap,cx,cy,x,y,s1,s2,c);
    if (p<0) p=p+4*x+6;
    else {
      p=p+4*(x-y)+10;
      y=y-1;
    }
  x=x+1;
}
if (x==y) plot8ps(bitmap,cx,cy,x,y,s1,s2,c);
}

plotsm(bitmap,cx,cy,r,angle)
char bitmap[][BMAX];
int cx,cy,r;
double angle;
{
   plota(bitmap,cx,cy,r,angle+(202*6.28318/360),angle+(338*6.28318/360));
}


printbyte(x)
unsigned x;
{
static int col=12;
static char delim=' ';

   if (++col<12) printf(", 0x%02x",x);
   else {
      printf("%c\n0x%02x",delim,x);
      col=0;
      delim=',';
   }
}
@//E*O*F plothap.c//
chmod u=rw,g=r,o=r plothap.c
 
echo x - tl.ps
sed 's/^@//' > "tl.ps" <<'@//E*O*F tl.ps//'
e
%%Trailer
EndPSPlot
@//E*O*F tl.ps//
chmod u=rw,g=r,o=r tl.ps
 
echo x - hap2ps
sed 's/^@//' > "hap2ps" <<'@//E*O*F hap2ps//'
#!/bin/sh

cat <hd.ps
cat
cat <tl.ps

@//E*O*F hap2ps//
chmod u=rwx,g=rx,o=rx hap2ps
 
exit 0