Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Path: utzoo!watmath!clyde!burl!ulysses!mhuxr!mhuxn!ihnp4!cbosgd!rlp
From: rlp@cbosgd.UUCP (J.  Knapp)
Newsgroups: net.micro.atari16
Subject: Re: Yet Another Fractal Program
Message-ID: <2078@cbosgd.UUCP>
Date: Sun, 4-May-86 12:05:43 EDT
Article-I.D.: cbosgd.2078
Posted: Sun May  4 12:05:43 1986
Date-Received: Wed, 7-May-86 01:35:49 EDT
References: <2077@cbosgd.UUCP>
Reply-To: rlp@cbosgd.UUCP (J.  Knapp)
Followup-To: net.micro.atari16
Distribution: net
Organization: AT&T Bell Laboratories, Columbus
Lines: 637
Keywords: (but not the Mandelbrot set)


Here is the source code for the uuencoded file 'levy.ttp' posted previously:
----------------------------------------------------------------------------
/*
 * levy.ttp:  Brownian Relief Fractal Generation Program for the Atari 520ST
 *
 * This program is distributed in the public domain with no restrictions on
 * its use.
 *
 * The algorithm implemented herein is described in B. B. Mandelbrot "Fractals:
 * Form, Chance, and Dimension" (Freeman 1977) p. 207. It was first defined
 * by P. Levy "Processus Stochastiques et Mouvement Brownien" (Paris:
 * Gauthier-Villars 1948).
 *
 * program author: J. M. Knapp  AT&T Bell Labs   Columbus OH  (cbosgd!nscs!jmk)
 *
 * Background (from Mandelbrot):
 *
 * [The Brownian model of the landscape] finally provides us with our
 * long-sought example of a curve that (a) is devoid of self intersections,
 * (b) is practically devoid of self contacts, (c) has fractal dimension
 * clearly greater than 1, and (d) is isotropic.
 *
 * More precisely, the dimension of these coastlines is 3/2, which is of
 * course much higher than most of Richardson's values [e.g. 1.26 for the
 * coast of Britain]. The Brownian artificial coastline is therefore
 * limited in its applicability. It recalls northern Canada, Indonesian
 * Islands, perhaps western Scotland, and the Aegean. It is applicable to
 * other examples as well, but certainly not to all.
 *
 * It is a pity, because a relief of dimension D=5/2 and coastlines of
 * dimension 3/2 would have been easy to explain. Indeed, the Brown-Levy
 * function is an excellent approximation of the relief that would have
 * been created by superimposing independent rectilinear faults.  The
 * generative model proceeds simply as follows. Starting with a horizontal
 * plateau, break it along a straight line chosen at random to introduce
 * a kind of cliff, a random difference between the levels of the two
 * sides of the break. Then start all over again, ad infinitum. The
 * process merely generalizes the Poisson process. With no need for
 * mathematical or physical details, we can see that the argument seizes
 * at least one aspect of tectonic evolution...
 *
 * From B.B. Mandelbrot "Fractals: Form, Chance, and Dimension"
 * 
 */

#include <stdio.h>
#include <math.h>
#include <fcntl.h>

/* PUT YOUR ST HEADER FILES HERE */
/* the following are Lattice C header files */
#include <gemlib.h>
#include <osbind.h>

#define TRUE  1
#define FALSE 0

#define NX     320           /* max number of horizontal pixels */
#define NY     200           /* max number of vertical pixels */
#define NCOLOR  15           /* number of colors excluding background */
#define GRMODE   0           /* graphics mode 0 (320x200) */
#define SAVEFILE "levy.pi1"  /* DEGAS default save-file */
#define FSIZE  32034         /* size of DEGAS save-file */

#define SEED   2211955       /* default pseudorandom seed */

/* help */
#define US0 "TTP parameter usage:\n\n"
#define US1 "number of iterations:\n"
#define US1a "   -n<N> must always be specified\n\n"
#define US2 "optionally [default values in ()]:\n"
#define US3 "-s<int> sets random seed (2211955)\n"
#define US4 "-x<1-320> # of horizontal pixels (320)\n"
#define US5 "-y<1-200> # of vertical pixels (200)\n"
#define US6 "-u<N> update screen every N iterations\n"
#define US7 "      (default: only update when done)\n"
#define US8 "-f<filename> save-file name (levy.pi1)\n"
#define US9 "-d<file> display previously-stored file\n"
#define US9a "-i to generate island-dominated relief\n\n"
#define US10 "examples:\n"
#define US11 "   -n10000 -i\n"
#define US12 "   -n100 -s9349490390 -fcoast.pi1\n"
#define US13 "   -dcoast.pi1\n\n"
#define US14 "press any key to abort run in progress\n\n"

/* VDI globals */
short contrl[12],
      intin[128],
      intout[128],
      ptsin[128],
      ptsout[128] ;

short oldpalette[16] ;   /* storage for old palette */

/*************/
main(argc,argv) 
int argc ;
char *argv[] ;
{

   short *relief ;   /* pointer to 'relief' data structure */
   long ncyc ;       /* number of iterations of generation process */
   long iseed ;      /* pseudorandom seed */
   int xl, yl ;      /* extent of each axis of screen region */
   short handle ;    /* VDI screen handle */
   short update ;    /* screen update every 'update' iterations */
   char fname[80] ;  /* save-file */
   short loadf ;     /* load flag: if true then load & display pic */
   short island ;    /* if true --> generate island-dominated relief */

   /* VDI palette (increasing altitude) */
   static short palette[NCOLOR+1] = { 0x000, 0x003, 0x004, 0x005,
                                      0x660, 0x550, 0x440, 0x340,
                                      0x040, 0x050, 0x432, 0x543,
                                      0x653, 0x444, 0x555, 0x777 } ;

   short xpalette[NCOLOR+1] ;  /* palette xlated to VDI ordering */

   short cginit(), cgopen() ;

   /* begin */
   conflush() ;   /* clear console input buffer */

   /* process argv[] options */
   options(argc,argv,&xl,&yl,&ncyc,&iseed,&update,fname,&loadf,&island) ;

   if (loadf)  /* load and display DEGAS file */
   {
      handle = cginit(palette,xpalette) ;   /* open screen */
      if (readfrac(fname,palette) < 0)      /* send file data to physbase */
      {
         gemdos(0x1) ;                      /* PAKTC */
         cgclose(handle) ;                  /* file read error */
         exit(1) ;
      }
      Setpallette(palette) ;                /* set new palette */
      gemdos(0x1) ;                         /* PAKTC */
      cgclose(handle) ;                     /* done, restore screen */
      exit(0) ;
   }

   /* generate new relief */

   srand48(iseed) ;  /* randomize */

   /* allocate 'relief' data structure */
   if((relief = (short *)malloc(sizeof(short)*xl*yl)) == 0)
   {
      fprintf(stderr,"Not enough memory\n") ;
      gemdos(0x1) ;  /* PAKTC */
      exit(1) ;
   }

   level(relief,xl,yl) ;  /* initialize to level 'plain' */

   handle = cginit(palette,xpalette) ;  /* initialize graphics device */

   if (cycle(handle,relief,xl,yl,ncyc,update,island))  /* generate relief */
   {
      plotrel(handle,relief,xl,yl,island) ;            /* plot relief */
      gemdos(0x1) ;
   }

   if (writefrac(fname,xpalette) < 0)  /* save screen to DEGAS file */
   {
      gemdos(0x1) ;
      cgclose(handle) ;    /* file write error, clean up & quit */
      exit(1) ;
   }
   cgclose(handle) ;                   /* close graphics device */

   exit(0) ;     /* done */
}

/* fractal generation routine */
int cycle(handle,relief,xl,yl,ncyc,update,island)
short handle ;
short *relief ;
int xl, yl ;
long ncyc ;
short update ;
short island ;
{
   int icyc ;  /* cycle counter */

   printf("working...\n") ;
   for (icyc = 0 ; icyc < ncyc ; icyc++)
   {
      if (drand48() < .5)            /* either inc or dec half-plane */
         inc_hplane(relief,xl,yl) ;
      else
         dec_hplane(relief,xl,yl) ;

      if (update && icyc && !(icyc % update))
         plotrel(handle,relief,xl,yl,island) ;  /* update screen */

      if (Bconstat(2))
      {
         conflush() ;  /* a key was pressed, flush and abort */
         return(FALSE) ;
      }
   }
   return(TRUE) ;
}

/* increment half-plane */
inc_hplane(relief,xl,yl)
short *relief ;
int xl, yl ;
{
   register short *endcolp, *ap ;  /* end-of-column pointer, array pointer */
   register short iy, ryl ;        /* row index, reg yl holder */
   int m, b ;                      /* slope & intercept of random line */

   ryl = yl ;             /* put yl in register */
   rndmb(xl,yl,&m,&b) ;   /* get random line parameters m & b */

   ap = relief + yl * (xl - 1) ;  /* point to beginning of last column */
   endcolp = ap + yl ;            /* point to end of last column */

   /* increment all elements below  random line */
   while (xl--)  /* for all columns */
   {
      iy = ((m * xl) >> 8) + b ;  /* starting y-value to increment */
      if(iy < 0) iy = 0 ;         /* can't start before 0th element! */
      if (iy < yl)                /* no action if past last element in col */
      {
         ap += iy ;                  /* point to 1st value to inc */
         while(ap < endcolp)         /* increment rest of column */
            (*(ap++))++ ;
         endcolp -= ryl ;            /* point to end of previous column */
         ap -= (ryl << 1) ;          /* point to beginning of previous column */
      }
      else
      {
         ap -= ryl ;          /* point to end of previous column */
         endcolp -= ryl ;     /* point to beginning of previous column */
      }
   }   
}

/* decrement half-plane */
dec_hplane(relief,xl,yl)
short *relief ;
int xl, yl ;
{
   register short *endcolp, *ap ;  /* end-of-column pointer, array pointer */
   register short iy, ryl ;        /* row index, reg yl holder */
   int m, b ;                      /* slope & intercept of random line */

   ryl = yl ;             /* put yl in register */
   rndmb(xl,yl,&m,&b) ;   /* get random line parameters m & b */

   ap = relief + yl * (xl - 1) ;  /* point to beginning of last column */
   endcolp = ap + yl ;            /* point to end of last column */

   /* decrement all elements to right of random line */
   while (xl--)  /* for all columns */
   {
      iy = ((m * xl) >> 8) + b ;  /* starting y-value to decrement */
      if(iy < 0) iy = 0 ;         /* can't start before 0th element! */
      if (iy < yl)                /* no action if past last element in col */
      {
         ap += iy ;                  /* point to 1st value to dec */
         while(ap < endcolp)         /* decrement rest of column */
            (*(ap++))-- ;
         endcolp -= ryl ;            /* point to end of previous column */
         ap -= (ryl << 1) ;          /* point to beginning of previous column */
      }
      else
      {
         ap -= ryl ;          /* point to end of previous column */
         endcolp -= ryl ;     /* point to beginning of previous column */
      }
   }   
}

/* set all elements to 0 (level plane) */
level(relief,xl,yl)
short *relief ;
int xl, yl ;
{
   register short *ap, *endp ;   /* array pointer, end pointer */

   ap = relief ;
   endp = relief + xl*yl ;
   while (ap < endp)
      *(ap++) = 0 ;
}

/* return slope (m) and intercept (b) of random line through plane */
/* (slope is scaled by 256) */
rndmb(xl,yl,m,b)
int xl, yl ;
int *m, *b ;
{
   float dx, dy ;  /* delta-x, delta-y */
   int x0, y0 ;  /* random point on plane */
   double drand48() ;

   dx = 0 ; dy = 0 ;
   while((dx = 2.0*drand48() - 1.0) == 0) ; /* prevent division by 0 */
   dy = drand48() ;
   x0 = xl*drand48() ;
   y0 = yl*drand48() ;       /* (x0,y0) random point */

   *m = (256.0*dy)/dx ;
   *b =  y0 - (*m*x0)/256 ;     /* calculate y-intercept */
}

/* initialize screen device */
short cginit(palette,xpalette)
short *palette, *xpalette ;
{
   /* translation from color index to palette index  PI = ctrans[CI] */
   static short ctrans[NCOLOR+1] =
      { 0, 2, 3, 6, 4, 7, 5, 8, 9, 10, 11, 14, 12, 15, 13, 1 } ;
   short colind ;     /* color index */
   short handle ;     /* screen handle */
   short cgopen() ;

   handle = cgopen() ;

   /* save old palette */
   for (colind = 0 ; colind < NCOLOR+1 ; colind++)
      oldpalette[colind] = Setcolor(colind,-1) ;

   /* initialize color look-up table */
   for (colind = 0 ; colind < NCOLOR+1 ; colind++)
      xpalette[colind] = palette[ctrans[colind]] ;
   Setpallette(xpalette) ;

   return(handle) ;
}

/* plot relief */
plotrel(handle,relief,xl,yl,island)
short handle ;
short *relief ;
int xl, yl ;
short island ;
{
   int iy, ix, icol ;   /* row index, column index, color index */
   short minval, maxval ;  /* min, max altitude */
   register short *ap ;    /* array pointer */

   /* extents of color regions (increasing altitude): */

   /* continental relief */
   static float cext_con[NCOLOR+1] = { 0, .15, .25, .35,
                                          .40, .45, .50,
                                          .55, .65, .70,
                                          .75, .80, .85,
                                          .90, .95, 1.00 } ;

   /* island relief */
   static float cext_ile[NCOLOR+1] = { 0, .30, .55,  .65,
                                          .68, .71,  .74,
                                          .77, .85,  .88,
                                          .91, .93,  .95,
                                          .97, .99, 1.00 } ;

   float *cext ;

   /* integer (absolute altitude) extents of color regions */
   short icext[NCOLOR+1] ;

   /* set extents pointer */
   if (island)
      cext = cext_ile ;
   else
      cext = cext_con ;

   minmax(relief,xl,yl,&maxval,&minval) ;  /* find max and min altitude */

   /* calculate integer extents */
   for (icol = 0 ; icol < NCOLOR+1 ; icol++)
      icext[icol] = minval + cext[icol] * (maxval - minval) ;

   icext[NCOLOR]++ ;  /* tricky: bump last value by one */

   ap = relief ;  /* point to beginning */
   for (ix= 0 ; ix < xl ; ix++)  /* for all columns */
   {
      register short colm[NY] ;  /* color indices for screen column */
      for (iy = 0 ; iy < yl ; iy++)   /* for all rows */
      {
         register short *cp, colind ;  /* color pointer, color index */
         cp = icext ;                  /* point to color region 0 */
         colind = 0 ;
         while(*ap >= *(cp++)) colind++ ;
         colm[iy] = colind ;
         ap++ ;
      }
      column(handle,colm,ix,yl) ;  /* plot column */
   }
}

/* find maximum & minimum values of altitude */
minmax(relief,xl,yl,maxval,minval)
short *relief ;
int xl, yl ;
short *maxval, *minval ;
{
   register short *ap, *endp ;   /* array pointer, end pointer */

   *maxval = *relief ; *minval = *relief ;
   ap = relief ;
   endp = ap + xl*yl ;
   while (ap < endp)
   {
      if (*ap > *maxval) *maxval = *ap ;
      if (*ap < *minval) *minval = *ap ;
      ap++ ;
   }
}

/* flush any characters in console input buffer */
conflush()
{
   char ch ;
   while(Bconstat(2))
      ch = Bconin(2) ;
}

/* process argv options */
options(argc,argv,xl,yl,ncyc,iseed,update,fname,loadf,island)
int argc ;
char *argv[] ;
int *xl, *yl ;
long *ncyc ;
long *iseed ;
short *update ;
char *fname ;
short *loadf ;
short *island ;
{
   int help ;
   int optind ;

   /* set defaults */
   help = FALSE ;
   *iseed = SEED ;
   *ncyc = 0 ;
   *xl = NX ;
   *yl = NY ;
   *update = 0 ;
   sprintf(fname,"%s",SAVEFILE) ;
   *loadf = FALSE ;
   *island = FALSE ;

   optind = 1 ;
   while ((optind < argc) && (*argv[optind] == '-'))
   {
      switch(*(argv[optind]+1))
      {
         case 's':
         case 'S':
           *iseed = atol(argv[optind]+2);
           break;

         case 'n':
         case 'N':
           *ncyc = atol(argv[optind]+2);
           break;

         case 'u':
         case 'U':
            *update = atoi(argv[optind]+2) ;
            break ;

         case 'f':
         case 'F':
            sprintf(fname,"%s",argv[optind]+2) ;
            break ;

         case 'd':
         case 'D':
            sprintf(fname,"%s",argv[optind]+2) ;
            *loadf = TRUE ;
            break ;

         case 'y':
         case 'Y':
           *yl = atoi(argv[optind]+2) ;
           break ;

         case 'x':
         case 'X':
           *xl = atoi(argv[optind]+2) ;
           break ;

         case 'i':
         case 'I':
            *island = TRUE ;
            break ;

         default:
           help = TRUE ;
           break ;
      }
      optind += 1 ;
   }

   if((optind != argc) || (argc == 1 )) help = TRUE ;

   if( help )
   {
      fprintf(stderr,"%s", US0);
      fprintf(stderr,"%s", US1);
      fprintf(stderr,"%s", US1a) ;
      fprintf(stderr,"%s", US2);
      fprintf(stderr,"%s", US3);
      fprintf(stderr,"%s", US4);
      fprintf(stderr,"%s", US5);
      fprintf(stderr,"%s", US6);
      fprintf(stderr,"%s", US7);
      fprintf(stderr,"%s", US8);
      fprintf(stderr,"%s", US9);
      fprintf(stderr,"%s", US9a) ;
      fprintf(stderr,"%s", US10);
      fprintf(stderr,"%s", US11);
      fprintf(stderr,"%s", US12);
      fprintf(stderr,"%s", US13);
      fprintf(stderr,"%s", US14) ;
      fprintf(stderr,"press any key to continue") ;
      gemdos(0x1) ;
      exit(1);
   }
}

/* close graphics screen and enable cursors */
cgclose(handle)
short handle ;
{
   v_show_c(handle) ;
   Cursconf(1,0) ;
   Setpallette(oldpalette) ;
   v_clsvwk(handle) ;
   appl_exit() ;
}

/* open graphics screen and disable cursors */
short cgopen()
{
   short handle ;
   short work_in[12],
         work_out[57] ;
   int iwk, set_type ;

   appl_init() ;

   for (iwk = 0 ; iwk < 10 ; work_in[iwk++] = 1) ;
   work_in[10] = 2 ; /* raster coordinates */

   v_opnvwk(work_in,&handle,work_out) ;
   v_hide_c(handle) ;  /* bye mouse */
   Cursconf(0,0) ;     /* bye cursor */
   v_clrwk(handle) ;   /* clear screen */
   set_type = vsm_type(handle,1) ;  /* marker type = point */
   return(handle) ;
}

/* plot column on graphics screen */
column(handle,colm,ix,npoints)
short handle ;
register short *colm ;
int ix ;
int npoints ;
{
   register short pxy[2] ;

   *pxy = ix ;
   *(pxy+1) = npoints ;
   while((*(pxy+1))--)
   {
      vsm_color(handle,colm[*(pxy+1)]) ;
      v_pmarker(handle,1,pxy) ;
   }
}

/* save screen data to disk in DEGAS format */
int writefrac(fname,palette) 
char *fname ;
short *palette ;
{
   char *scrnp ;
   int fd, n ;
   short gmode[1] ;

   scrnp = (char *)Physbase() ;  /* get pointer to display frame */
   *gmode = GRMODE ;

   if ((fd = open(fname, O_RAW | O_CREAT | O_WRONLY)) < 0) {
      fprintf(stderr,"open of '%s' failed\n",fname) ;
      return(-1) ;
      }

   n = write(fd, gmode, 2) ;      /* write graphics mode */
   n += write(fd, palette, 32) ;  /* write palette */
   if ((n += write(fd, scrnp, 32000)) != FSIZE)  /* write screen */
   {
      fprintf(stderr,"file write error: only %d bytes written\n",n) ;
      close(fd) ;
      return(-1) ;
   }
   close(fd) ;
   return(0) ;
}

/* read file (DEGAS format) */
int readfrac(fname,palette)
char *fname ;
short *palette ;
{
   int fd, n ;
   short gmode[1] ;
   char *scrnp ;

   if ((fd = open(fname, O_RAW | O_RDONLY)) < 0) {
      fprintf(stderr,"open of '%s' failed\n",fname) ;
      return(-1) ;
      }

   scrnp = (char *)Physbase() ;
   n = read(fd, gmode, 2) ;      /* read graphics mode */
   n += read(fd, palette, 32) ;  /* read palette */
   if ((n += read(fd, scrnp, 32000)) != FSIZE)  /* read screen */
   {
      fprintf(stderr,"file write error: only %d bytes read\n",n) ;
      close(fd) ;
      return(-1) ;
   }
   close(fd) ;
   return(0) ;
}