/* Next available MSG number is  34 */

/*    

   MOUNTAIN.C

   Copyright (C) 1989, 1990, 1991, 1992, 1993, 1994 by Autodesk, Inc.
 
   Permission to use, copy, modify, and distribute this software in 
   object code form for any purpose and without fee is hereby granted, 
   provided that the above copyright notice appears in all copies and 
   that both that copyright notice and the limited warranty and 
   restricted rights notice below appear in all supporting 
   documentation.
 
   AUTODESK PROVIDES THIS PROGRAM "AS IS" AND WITH ALL FAULTS.  
   AUTODESK SPECIFICALLY DISCLAIMS ANY IMPLIED WARRANTY OF 
   MERCHANTABILITY OR FITNESS FOR A PARTICULAR USE.  AUTODESK, INC.
   DOES NOT WARRANT THAT THE OPERATION OF THE PROGRAM WILL BE 
   UNINTERRUPTED OR ERROR FREE.
 
   Use, duplication, or disclosure by the U.S. Government is subject to 
   restrictions set forth in FAR 52.227-19 (Commercial Computer 
   Software - Restricted Rights) and DFAR 252.227-7013(c)(1)(ii) 
   (Rights in Technical Data and Computer Software), as applicable.
    
   .

      DESCRIPTION:

        Fractal mountain landscape generator

        Designed and implemented in November of 1989 by John Walker

        This  program generates fractal mountain landscapes by Fourier
        filtering of random data.  The landscapes are  represented  as
        pface  meshes  within  an  AutoCAD  drawing.   To  run the
        program, first load the ADS application with the command:

            (xload "mountain")

        and enter the command:

            MOUNTAIN

        You're  asked the to specify the parameters of the mountain as
        follows:

            Mesh size (power of 2) <32>:
            Fractal dimension (typically between 1 and 3) <2.15>:
            Power law scaling exponent <1>:

        The  mountain  range is  modeled as a square  pface  mesh with
        edge size given by the first prompt.  The mesh size must be  a
        power  of  two (this limitation is imposed by the fast Fourier
        transform algorithm; if you specify a mesh size that's  not  a
        power  of two it will be rounded up to the next larger power).

        The  fractal dimension controls the roughness of the generated
        landscape by applying a low-pass spatial filter to the  random
        frequency  domain  data  before performing the inverse Fourier
        transform.  The higher the fractal dimension the  rougher  the
        terrain;  values  between  2  and  3  are typical.  Dimensions
        between 1.7 and  2  generate  Appalachian-type  hills;  values
        between  2 and 2.5 yield Sierra Nevada like peaks; and numbers
        close to 3 result in fantasy landscapes of spires (which  look
        good only at high mesh densities).

        The output of the mountain generation algorithm is a table  of
        heights.   Negative  heights  are  set to zero and coloured as
        water.  You can scale elevations above sea  level  by  raising
        the  raw  values to a power specified by the power law scaling
        exponent.  The default value of 1 simply uses  the  elevations
        as generated.  Powers of less than 1 mimic eroded terrain (try
        0.33--cube root scaling often produces  landscapes  that  look
        glacially sculpted like Yosemite).  Powers greater than 1 make
        increasingly forbidding peaks.  Again, if you use large  power
        law  exponents,  the  results  will  look  best at higher mesh
        densities.

        Every time  you  enter  the  MOUNTAIN  command  you'll  get  a
        different  landscape.   The random number generator is started
        from the time and date, so a given result will normally  never
        recur.   If  you  use  the  command several times, remember to
        ERASE the previous mesh to keep the new mountains  from  being
        drawn  over  top of the old ones.  If you want to generate the
        same terrain repeatedly (for example, to experiment  with  the
        effects of fractal dimension and power law exponent settings),
        you can set the random sequence used to generate  the  terrain
        with the SEED command.  This command prompts:

            Random number seed (0 to reseed) <0xNNNN>:

        where  0xNNNN is the current starting seed in hexadecimal.  If
        you enter a nonzero value it is used as the seed by  the  next
        MOUNTAIN  command.   Subsequent commands advance the generator
        and continue  to  generate  different  terrain,  so  you  must
        re-enter  the SEED before each MOUNTAIN command if you want to
        repeat the last output.  The default for the SEED  command  is
        the  last seed entered, so simply entering a blank response to
        the command resets the seed to the last value, making it  easy
        to  repeatedly generate from the same seed.  If a seed of 0 is
        entered, a new random value is generated  from  the  time  and
        date.

        The  faces  of  the pface mesh  are  coloured  individually to
        represent their elevation.  The WORLD command  allows  you  to
        choose  between  two  different  colour  schemes.  The command
        prompts you:

                World:  0 = Landscape (8 colours)
                        1 = Landscape (256 colours)
                        2 = Clouds
                World type <0>:

        The default world, 0, colours negative elevations as different
        shades of blue (to indicate "water depth"), and terrain  on  a
        scale  from  green through yellows to white for "snowy peaks".
        World type 1 restricts the palette of colours to  the  default
        AutoCAD  8  colour  set;  the  resulting  pictures  are not as
        attractive but will display properly on systems with 8  colour
        capability   or  above.   World  type  2  selects  colours  to
        illustrate affinity between fractal clouds and mountains.  All
        negative  elevations  in  the  clouds world are coloured blue;
        other faces are assigned shades from dark grey to  white  with
        increasing  elevation.   Pictures  generated  with  the clouds
        colouring are  intended  to  be  viewed  in  plan  view,  from
        directly  above,  and  look  best in shaded renderings at high
        resolution.  Cloud colouring requires a  256  colour  display.

        The  mountains are generated on a dark grey (cyan for 8 colour
        landscapes) base, drawn with invisible edges, which serves  to
        obscure the underside of the terrain in hidden line and shaded
        renderings.

        The  MOUNTAIN  command  always  generates  a pface  mesh; this
        mesh  can  be  created  either  by  invocation  of  the  PFACE
        command  through  the  ads_command() facility or directly with
        ads_entmake().  The MAKEMODE command lets you specify how  the
        mesh  is  to  be created, allowing you to try both methods and
        compare their relative performance.  You're prompted:

                Mesh generation:  0 = Directly with ads_entmake()
                                  1 = PFACE command via ads_command()
                Generation type <0>:

        Enter  the  generation  mode  you prefer.  Subsequent MOUNTAIN
        commands will  use  the  mode  you  chose.   If  the  mesh  is
        generated with ads_entmake(), extended entity data is attached
        to it that records the  world  type,  fractal  dimension,  and
        power factor used to create the mesh.

        References:

            Peitgen, H.-O., and Saupe, D. eds., The Science Of Fractal
                Images, New York: Springer Verlag, 1988.

            Press, W. H., Flannery, B. P., Teukolsky, S. A., Vetterling,
                W. T., Numerical Recipes In C, New Rochelle; Cambridge
                University Press, 1988.

*/

#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <string.h>
#ifdef UNIX
#include <memory.h>
#endif

#include "adslib.h"

#define AppName  /*MSG1*/"AUTODESK_FRACTAL_MOUNTAINS"

#ifndef M_PI
#define M_PI    3.14159265358979323846
#endif

#if defined(__WATCOMC__) || defined(TURBOC)
#undef max
#undef min
#endif

#define max(a,b) ((a) >  (b) ? (a) : (b))
#define min(a,b) ((a) <= (b) ? (a) : (b))

/* Definitions used to address real and imaginary parts in a two-dimensional
   array of complex numbers as stored by fourn(). */

#define Real(v, x, y)  v[1 + (((x) * n) + (y)) * 2]
#define Imag(v, x, y)  v[2 + (((x) * n) + (y)) * 2]

/*  Data types  */

typedef enum {False = 0, True = 1} Boolean;

/* Definitions  to  wrap  around  submission  of  AutoCAD commands to
   prevent their being echoed.  */

#define Cmdecho  False                /* Make True for debug command output */

#define CommandB()  { struct resbuf rBc, rBb, rBu, rBh; \
        ads_getvar(/*MSG0*/"CMDECHO", &rBc); \
        ads_getvar(/*MSG0*/"BLIPMODE", &rBb); \
        ads_getvar(/*MSG0*/"HIGHLIGHT", &rBh); \
        rBu.restype = RTSHORT; \
        rBu.resval.rint = (int) Cmdecho; \
        ads_setvar(/*MSG0*/"CMDECHO", &rBu); \
        rBu.resval.rint = (int) False; \
        ads_setvar(/*MSG0*/"BLIPMODE", &rBu); \
        ads_setvar(/*MSG0*/"HIGHLIGHT", &rBu)

#define CommandE()  ads_setvar(/*MSG0*/"CMDECHO", &rBc); \
                    ads_setvar(/*MSG0*/"BLIPMODE", &rBb); \
                    ads_setvar(/*MSG0*/"HIGHLIGHT", &rBh); }

/*  Definitions  that permit you to push and pop system variables with
    minimal complexity.  These don't work  (and  will  cause  horrible
    crashes  if  used)  with  string  variables,  but since all string
    variables are read-only, they cannot be saved and restored in  any
    case.  */

#define PushVar(var, cell, newval, newtype) { struct resbuf cell, cNeW; \
        ads_getvar(var, &cell); cNeW.restype = cell.restype;             \
        cNeW.resval.newtype = newval; ads_setvar(var, &cNeW)

#define PopVar(var, cell) ads_setvar(var, &cell); }

/* Set point variable from three co-ordinates */

#define Spoint(pt, x, y, z)  pt[X] = (x);  pt[Y] = (y);  pt[Z] = (z)

/* Copy point from another */

#define Cpoint(d, s)   d[X] = s[X];  d[Y] = s[Y];  d[Z] = s[Z]

/* Definitions for building result buffer chains. */

#define tacky()     struct resbuf *rb, *rbtail, *ri
#define tackrb(t)   ri=ads_newrb(t);  rbtail->rbnext=ri;  rbtail=ri
#define defent(e)   rb=ri=rbtail=ads_newrb(0); ri->resval.rstring = strsave(e)
#define Xed(g)     (1000 + (g))       /* Generate extended entity data code */
#define makent() { int stat=ads_entmake(rb); if (stat!=RTNORM) { \
                   ads_printf(/*MSG2*/"Error creating %s entity\n",       \
                   rb->resval.rstring); } ads_relrb(rb); }
#define tackint(code, val) tackrb(code); ri->resval.rint = (val)
#define tackreal(code, val) tackrb(code); ri->resval.rreal = (val)
#define tackpoint(code, x, y, z) tackrb(code); \
        Spoint(ri->resval.rpoint, (x), (y), (z))
#define tackvec(code, vec) tackrb(code); Cpoint(ri->resval.rpoint, vec)
#define tackstring(code, s) tackrb(code), ri->resval.rstring = strsave(s)

#define V  (void)

/* Utility definition to get an  array's  element  count  (at  compile
   time).   For  example:

       int  arr[] = {1,2,3,4,5};
       ...
       printf("%d", ELEMENTS(arr));

   would print a five.  ELEMENTS("abc") can also be used to  tell  how
   many  bytes are in a string constant INCLUDING THE TRAILING NULL. */

#define ELEMENTS(array) (sizeof(array)/sizeof((array)[0]))

/*  Display parameters  */

#define COL_BASE    (wtype == 0 ? 4 : 250) /* Baseplate colour */
#define BASE_THICK  0.1               /* Baseplate thickness */

/*  Forward functions  */

long    time _((long *));
void    fourn _((float *, int *, int, int));
void    initgauss _((int));
double  gauss _((void));
void    spectralsynth _((float **, int, double));
void    mtn_srand _((unsigned int));
int     mtn_rand _((void));
int     veryrandom _((void));
void    main _((int, char **));
Boolean funcload _((void));
Boolean initacad _((Boolean));
void    initseed _((void));
int     calccol _((float *, int, int, int, ads_real, ads_real));
void    meshcol _((float *, int, int, int, ads_real, ads_real));
void    genmesh _((float *, int, ads_real, ads_real));
char    *strsave _((char *));
void    makevtx _((ads_point));
void    makeface _((int, int, int, int, int));
void    makemesh _((float *, int, ads_real, ads_real));
void    mountain _((void));
void    seed _((void));
void    world _((void));
void    makemode _((void));

#ifdef DEBUG
void    test _((void));
#endif

/*  Local variables  */

static int nrand;                     /* Gauss() sample count */
static double arand, gaussadd, gaussfac; /* Gaussian random parameters */
static int rseed;                     /* Current random seed */
static Boolean seeded = False;        /* Initial seed computed ? */
static Boolean forceseed = False;     /* Force seed on next computation ? */
static int ccolour = -1;              /* Current face colour */
static ads_real fracdim = 2.15;       /* Fractal dimension */
static ads_real powscale = 1.0;       /* Elevation power factor */
static int wtype = 0;                 /* World type */
static Boolean entmake = True;        /* Use ads_entmake() to create mesh ? */

static Boolean functional;            /* C:command is returning result */

/* The following variables are used by the psuedo-random sequence generator
   veryrandom() and must be unsigned so the right shift works properly */

static unsigned long regA = 1;
static unsigned long regB = 1;
static unsigned long regC = 1;

/*  Command definition and dispatch table.  */

struct {
        char *cmdname;
        void (*cmdfunc)();
} cmdtab[] = {
/*        Name         Function  */
#ifdef DEBUG
{/*MSG3*/"TEST",       test},
#endif
{/*MSG4*/"MAKEMODE",   makemode},
{/*MSG5*/"MOUNTAIN",   mountain},
{/*MSG6*/"SEED",       seed},
{/*MSG7*/"WORLD",      world}
};

/*      FOURN  --  Multi-dimensional fast Fourier transform

        Called with arguments:

           data       A  one-dimensional  array  of  floats  (NOTE!!!   NOT
                      DOUBLES!!), indexed from one (NOTE!!!   NOT  ZERO!!),
                      containing  pairs of numbers representing the complex
                      valued samples.  The Fourier transformed results  are
                      returned in the same array.

           nn         An  array specifying the edge size in each dimension.
                      THIS ARRAY IS INDEXED FROM  ONE,  AND  ALL  THE  EDGE
                      SIZES MUST BE POWERS OF TWO!!!

           ndim       Number of dimensions of FFT to perform.  Set to 2 for
                      two dimensional FFT.

           isign      If 1, a Fourier transform is done, if -1 the  inverse
                      transformation is performed.

        This  function  is essentially as given in Press et al., "Numerical
        Recipes In C", Section 12.11, pp.  467-470.
*/

static void fourn(data, nn, ndim, isign)
  float data[];
  int nn[], ndim, isign;
{
    register int i1, i2, i3;
    int i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
    int ibit, idim, k1, k2, n, nprev, nrem, ntot;
    float tempi, tempr;
    double theta, wi, wpi, wpr, wr, wtemp;

#define SWAP(a,b) tempr=(a); (a) = (b); (b) = tempr

    ntot = 1;
    for (idim = 1; idim <= ndim; idim++)
        ntot *= nn[idim];
    nprev = 1;
    for (idim = ndim; idim >= 1; idim--) {
        n = nn[idim];
        nrem = ntot / (n * nprev);
        ip1 = nprev << 1;
        ip2 = ip1 * n;
        ip3 = ip2 * nrem;
        i2rev = 1;
        for (i2 = 1; i2 <= ip2; i2 += ip1) {
            if (i2 < i2rev) {
                for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2) {
                    for (i3 = i1; i3 <= ip3; i3 += ip2) {
                        i3rev = i2rev + i3 - i2;
                        SWAP(data[i3], data[i3rev]);
                        SWAP(data[i3 + 1], data[i3rev + 1]);
                    }
                }
            }
            ibit = ip2 >> 1;
            while (ibit >= ip1 && i2rev > ibit) {
                i2rev -= ibit;
                ibit >>= 1;
            }
            i2rev += ibit;
        }
        ifp1 = ip1;
        while (ifp1 < ip2) {
            ifp2 = ifp1 << 1;
            theta = isign * 6.28318530717959 / (ifp2 / ip1);
            wtemp = sin(0.5 * theta);
            wpr = -2.0 * wtemp * wtemp;
            wpi = sin(theta);
            wr = 1.0;
            wi = 0.0;
            for (i3 = 1; i3 <= ifp1; i3 += ip1) {
                for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) {
                    for (i2 = i1; i2 <= ip3; i2 += ifp2) {
                        k1 = i2;
                        k2 = k1 + ifp1;
                        tempr = wr * data[k2] - wi * data[k2 + 1];
                        tempi = wr * data[k2 + 1] + wi * data[k2];
                        data[k2] = data[k1] - tempr;
                        data[k2 + 1] = data[k1 + 1] - tempi;
                        data[k1] += tempr;
                        data[k1 + 1] += tempi;
                    }
                }
                wr = (wtemp = wr) * wpr - wi * wpi + wr;
                wi = wi * wpr + wtemp * wpi + wi;
            }
            ifp1 = ifp2;
        }
        nprev *= n;
    }
}
#undef SWAP

/*  INITGAUSS  --  Initialise random number generators.  */

static void initgauss(seed)
  int seed;
{
    nrand = 4;
    arand = pow(2.0, (sizeof(int) * 8) - 1.0) - 1.0;
    gaussadd = sqrt(3.0 * nrand);
    gaussfac = 2 * gaussadd / (nrand * arand);
    mtn_srand(seed);
}

/*  GAUSS  --  Return a Gaussian random number.  */

static double gauss()
{
    int i;
    double sum = 0.0;

    for (i = 1; i <= nrand; i++) {
        sum += mtn_rand();
    }
    return gaussfac * sum - gaussadd;
}

/*  SPECTRALSYNTH  --  Spectrally synthesised fractal motion in two
                       dimensions.  */

static void spectralsynth(x, n, h)
  float **x;
  int n;
  double h;
{
    int i, j, i0, j0, nsize[3];
    double rad, phase;
    float *a;

    a = (float *) malloc((unsigned) (i =
                                     (((n * n) + 1) * 2 * sizeof(float))));
    if (a == NULL) {
        V fprintf(stderr,
                  /*MSG8*/"Cannot allocate %d x %d result array (%d bytes).\n",
                  n, n, i);
        exit(1);
    }
    *x = a;
    V memset((char *) a, 0, i);       /* Clear array to zeroes */

    for (i = 0; i <= n / 2; i++) {
        for (j = 0; j <= n / 2; j++) {
            phase = 2 * M_PI * (mtn_rand() / arand);
            if (i != 0 || j != 0) {
                rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
            } else {
                rad = 0;
            }
            Real(a, i, j) = rad * cos(phase);
            Imag(a, i, j) = rad * sin(phase);
            i0 = (i == 0) ? 0 : n - i;
            j0 = (j == 0) ? 0 : n - j;
            Real(a, i0, j0) = rad * cos(phase);
            Imag(a, i0, j0) = - rad * sin(phase);
        }
    }
    Imag(a, n / 2, 0) = Imag(a, 0, n / 2) = Imag(a, n / 2, n / 2) = 0;
    for (i = 1; i <= n / 2 - 1; i++) {
        for (j = 1; j <= n / 2 - 1; j++) {
            phase = 2 * M_PI * (mtn_rand() / arand);
            rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
            Real(a, i, n - j) = rad * cos(phase);
            Imag(a, i, n - j) = rad * sin(phase);
            Real(a, n - i, j) = rad * cos(phase);
            Imag(a, n - i, j) = - rad * sin(phase);
        }
    }

    nsize[0] = 0;
    nsize[1] = nsize[2] = n;          /* Dimension of frequency domain array */
    fourn(a, nsize, 2, -1);           /* Take inverse 2D Fourier transform */
}

/*  MTN_SRAND  --  Seed the psuedo-random number generator.  */

static void
/*FCN*/mtn_srand(seed)
    unsigned int seed;
{
    regA = regB = regC = seed;
}

/*  MTN_RAND  --  The psuedo-random number generator.  */

static int
/*FCN*/mtn_rand()
{
    int i;
    int rnum = 0;

    for (i=0; i < sizeof(int)*8; i++) {
        rnum |= veryrandom();
        rnum <<= 1;
    }
    return rnum;
}

/*  VERYRANDOM  --  Psuedo-Random Sequence Generator for 32-Bit CPUs,
                    by Bruce Schneier, Dr. Dobb's Journal, February 1992 */

static int
/*FCN*/veryrandom()
{
    /* regA is a 32-bit Linear Feedback Shift Register (LFSR) */
    /* regB is a 31-bit LFSR.  regC is a 29-bit LFSR. */
    /* The feedback sequences are chosen to be maximum length. */
    regA = ((((regA>>31)^(regA>>6)^(regA>>4)^(regA>>2)^(regA>>1)^regA)
        & 0x00000001)<<31) | (regA>>1);
    regB = ((((regB>>30)^(regB>>2)) & 0x00000001)<<30) | (regB>>1);
    regC = ((((regC>>28)^(regC>>1)) & 0x00000001)<<28) | (regC>>1);
    return ((regA & regB) | (!regA & regC)) & 0x00000001;
}

/*  MAIN  --   Main ADS transaction processor.  */

void main(argc, argv)
  int argc;
  char *argv[];
{
    int stat, cindex, scode = RSRSLT;
    char errmsg[80];

    ads_init(argc, argv);             /* Initialise the application */

    /* Main dispatch loop. */

    while (True) {

        if ((stat = ads_link(scode)) < 0) {
            V sprintf(errmsg,
                      /*MSG9*/"MOUNTAIN: bad status from ads_link() = %d\n",
                      stat);
#ifdef Macintosh
            macalert(errmsg);
#else
            puts(errmsg);
            fflush(stdout);
#endif /* Macintosh */
            exit(1);
        }

        scode = RSRSLT;               /* Default return code */

        switch (stat) {

        case RQXLOAD:                 /* Load functions.  Called at the start
                                         of the drawing editor.  Re-initialise
                                         the application here. */
            scode = -(funcload() ? RSRSLT : RSERR);
            break;

        case RQSUBR:                  /* Evaluate external lisp function */
            cindex = ads_getfuncode();
            functional = False;
            if (!initacad(False)) {
                ads_printf(/*MSG10*/"\nUnable to initialise application.\n");
            } else {

                /* Execute the command from the command table with
                   the index associated with this function. */

                if (cindex > 0) {
                    cindex--;
                    assert(cindex < ELEMENTS(cmdtab));
                    (*cmdtab[cindex].cmdfunc)();
                }
            }
            if (!functional)
                ads_retvoid();        /* Void result */
            break;

        default:
            break;
        }
    }
}

/* FUNCLOAD  --  Load external functions into AutoLISP */

static Boolean funcload()
{
    char ccbuf[40];
    int i;

    V strcpy(ccbuf, /*MSG0*/"C:");
    for (i = 0; i < ELEMENTS(cmdtab); i++) {
        V strcpy(ccbuf + 2, cmdtab[i].cmdname);
        ads_defun(ccbuf, i + 1);
    }

    return initacad(True);            /* Reset AutoCAD initialisation */
}

/*  INITACAD  --  Initialise the required modes in the AutoCAD
                  drawing.  */

static Boolean initacad(reset)
  Boolean reset;
{
    static Boolean initdone, initok;

    if (reset) {
        initdone = False;
        initok = True;
    } else {
        if (!initdone) {
            struct resbuf *rp;
            static char aname[] = AppName;  /* Static to keep compiler from
                                               making making one copy of
                                               constant string for every
                                               reference below. */

            /* First of all, see if our application is registered
               and register it if not. */

            if ((rp = ads_tblsearch(/*MSG0*/"APPID", aname, False)) == NULL) {
                if (ads_regapp(aname) != RTNORM) {
                    ads_printf(/*MSG11*/"Cannot register application %s.\n",
                               aname);
                    initdone = True;
                    return (initok = False);
                }
            } else {
                ads_relrb(rp);
            }

            /* Reset the program modes to standard values upon
               entry to the drawing editor. */

            initdone = initok = True;

            seeded = False;           /* No seed specified */
            forceseed = False;        /* No seed stuck in its throat */
            fracdim = 2.15;           /* Fractal dimension */
            powscale = 1.0;           /* Elevation power factor */
            wtype = 0;                /* World type */
            entmake = True;           /* Use ads_entmake to create meshes */
        }
    }
    return initok;
}

#ifdef DEBUG
/*  TEST  --  Verify that command invocation works.  */

static void test()
{
    ads_printf(/*MSG12*/"\nTEST command!\n");
#ifdef DEBUGDUMP
    {   ads_name en;
        struct resbuf *rb, *ri;

        ads_entnext(NULL, en);
        rb = ri = ads_entget(en);

        while (ri != NULL) {
            ads_printf(/*MSG13*/"Type %d\n", ri->restype);
            ri = ri->rbnext;
        }
        ads_relrb(rb);
    }
#else
    sleep(30);                        /* Delay to let debugger snag us */
#endif
}
#endif

/*  INITSEED  --  Generate initial random seed, if needed.  */

static void initseed()
{
    if (!seeded) {
        int i;

        mtn_srand(((int) (time((long *) NULL) ^ 0x5B3CF37C) & ((int) ~0)));
        for (i = 0; i < 7; i++)
            V mtn_rand();
        rseed = mtn_rand();
        seeded = forceseed = True;
    }
}

/*  CALCCOL  --  Calculate colour of mesh item from its relative
                 elevation.  */

static int calccol(a, n, x, y, rmax, rmin)
  float *a;
  int n, x, y;
  ads_real rmax, rmin;
{
    ads_real h1 = Real(a, x, y),
             h2 = Real(a, x + 1, y),
             h3 = Real(a, x, y + 1),
             h4 = Real(a, x + 1, y + 1),
             h;
    int rcol;

    static struct {
       ads_real cthresh;
       int tcolour;
    } ctab[] = {
       {0.25, 80},
       {0.6, 90},
       {0.75, 50},
       {0.80, 40},
       {0.90, 131},
       {1.00, 255}
      },
      etab[] = {
       {0.50, 3},
       {0.60, 2},
       {0.80, 1},
       {0.90, 6},
       {1.00, 7}
      },
      clouds[] = {
       {0.1666, 250},
       {0.3333, 251},
       {0.5000, 252},
       {0.6666, 253},
       {0.8333, 254},
       {1.0000, 255}
      };

    h = max(h1, max(h2, max(h3, h4)));
    if (h <= 0) {
        int ih;

        if (wtype == 0 || wtype == 2) {
            rcol = 5;                 /* Flat blue */
        } else {
            h = (h1 + h2 + h3 + h4) / 4;
            if (rmin == 0)
                rmin = 1;
            h = h / rmin;
            ih = (h * 3) + 0.5;
            rcol = 140 + min(3, ih) * 10;
        }
    } else {
        int i;

        if (rmax == 0)
            rmax = 1;
        h = h / rmax;
        if (wtype == 0) {
            for (i = 0; i < ELEMENTS(etab); i++) {
                if (h <= etab[i].cthresh)
                    break;
            }
            assert(i < ELEMENTS(etab));
            rcol = etab[i].tcolour;
        } else if (wtype == 2) {
            for (i = 0; i < ELEMENTS(clouds); i++) {
                if (h <= clouds[i].cthresh)
                    break;
            }
            assert(i < ELEMENTS(clouds));
            rcol = clouds[i].tcolour;
        } else {
            for (i = 0; i < ELEMENTS(ctab); i++) {
                if (h <= ctab[i].cthresh)
                    break;
            }
            assert(i < ELEMENTS(ctab));
            rcol = ctab[i].tcolour;
        }
    }
    return rcol;
}


/*  MESHCOL  --  Set colour of mesh item from its relative elevation.  */

static void meshcol(a, n, x, y, rmax, rmin)
  float *a;
  int n, x, y;
  ads_real rmax, rmin;
{
    int rcol = calccol(a, n, x, y, rmax, rmin);

    if (ccolour != rcol) {
        ccolour = rcol;
        ads_command(RTSTR, /*MSG31*/"_.Colour", RTSHORT, ccolour, RTNONE);
    }
}

/*  GENMESH  --  Generate mesh from elevation array.  This version
                 creates the mesh by submitting a PFACE command
                 with the ads_command function.  */

static void genmesh(a, n, rmax, rmin)
  float *a;
  int n;
  ads_real rmax, rmin;
{
    int i, j;

    CommandB();
    ads_command(RTSTR, /*MSG32*/"_.Pface", RTNONE);

    /* Send the vertex table. */

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            ads_point p;

            Spoint(p, i, j,
                   BASE_THICK + (n / (rmax * 2)) * max(0, Real(a, i, j)));
            ads_command(RT3DPOINT, p, RTNONE);
        }
    }

    /* Send the baseplate vertices.  We compute and store the
       four edges of the baseplate, addressed by edge number from
       0 to 3 and vertex number along the edge. */

    for (i = 0; i < n; i++) {
        ads_point p1, p2, p3, p4;

        Spoint(p1, i, 0, 0);
        Spoint(p2, i, n - 1, 0);
        Spoint(p3, 0, i, 0);
        Spoint(p4, n - 1, i, 0);
        ads_command(RT3DPOINT, p1, RT3DPOINT, p2,
                    RT3DPOINT, p3, RT3DPOINT, p4, RTNONE);
    }

    ads_command(RTSTR, "", RTNONE);

    /* Send the mesh tiles. */

    ccolour = -1;
#define Vtx(x, y) ((((x) * n) + (y)) + 1)     /* Address mesh vertex */
#define Bpx(e, x) (((n * n) + 1) + (((x) * 4) + (e)))
    for (i = 0; i < (n - 1); i++) {
        for (j = 0; j < (n - 1); j++) {
            meshcol(a, n, i, j, rmax, rmin);
            ads_command(RTSHORT, Vtx(i, j),
                        RTSHORT, Vtx(i, j + 1) * (j == n - 2 ? 1 : -1),
                        RTSHORT, Vtx(i + 1, j + 1) * (i == n - 2 ? 1 : -1),
                        RTSHORT, Vtx(i + 1, j),
                        RTSTR, "", RTNONE);
        }
    }


    /* Draw the baseplate */

    ads_command(RTSTR, /*MSG33*/"_.Colour", RTSHORT, COL_BASE, RTNONE);
    ads_command(RTSHORT, -Bpx(0, 0), RTSHORT, -Bpx(1, 0),
                RTSHORT, -Bpx(1, n - 1), RTSHORT, -Bpx(0, n - 1),
                RTSTR, "", RTNONE);

    /* Build the "walls" that connect the baseplate to the edges
       of the model. */

    for (i = 0; i < n - 1; i++) {
        ads_command(RTSHORT, -Bpx(0, i),
                    RTSHORT, -Vtx(i, 0),
                    RTSHORT, -Vtx(i + 1, 0),
                    RTSHORT, -Bpx(0, i + 1),
                    RTSTR, "",

                    RTSHORT, -Bpx(1, i),
                    RTSHORT, -Vtx(i, n - 1),
                    RTSHORT, -Vtx(i + 1, n - 1),
                    RTSHORT, -Bpx(1, i + 1),
                    RTSTR, "",

                    RTSHORT, -Bpx(2, i),
                    RTSHORT, -Vtx(0, i),
                    RTSHORT, -Vtx(0, i + 1),
                    RTSHORT, -Bpx(2, i + 1),
                    RTSTR, "",

                    RTSHORT, -Bpx(3, i),
                    RTSHORT, -Vtx(n - 1, i),
                    RTSHORT, -Vtx(n - 1, i + 1),
                    RTSHORT, -Bpx(3, i + 1),
                    RTSTR, "", RTNONE);
    }

#undef Vtx
#undef Bpx
    ads_command(RTSTR, "", RTNONE);
    CommandE();
}

/*  STRSAVE  --  Allocate a duplicate of a string.  */

static char *strsave(s)
  char *s;
{
    char *c = malloc((unsigned) (strlen(s) + 1));

    if (c == NULL)
        ads_abort(/*MSG15*/"Out of memory");
    V strcpy(c, s);
    return c;
}

/*  MAKEVTX  --  Append vertex entity to the database. */

static void makevtx(p)
  ads_point p;
{
    tacky();

    defent(/*MSG0*/"VERTEX");
    tackvec(10, p);                   /* Vertex point */
    tackint(70, 128 + 64);            /* Flags = Pface vertex point */
    makent();
}

/*  MAKEFACE  --  Append vertex entity representing a face in the mesh. */

static void makeface(colour, p1, p2, p3, p4)
  int colour, p1, p2, p3, p4;
{
    tacky();

    defent(/*MSG0*/"VERTEX");
    tackpoint(10, 0, 0, 0);           /* Vertex point */
    tackint(62, colour);              /* Face colour */
    tackint(70, 128);                 /* Vertex flags = Pface mesh face */
    tackint(71, p1);                  /* Vertex 1 */
    tackint(72, p2);                  /* Vertex 2 */
    tackint(73, p3);                  /* Vertex 3 */
    tackint(74, p4);                  /* Vertex 4 */
    makent();
}

/*  MAKEMESH  --  Create the mesh using the ads_entmake mechanism.  */

static void makemesh(a, n, rmax, rmin)
  float *a;
  int n;
  ads_real rmax, rmin;
{
    int i, j;
    tacky();

    defent(/*MSG0*/"POLYLINE");
    tackint(66, 1);                   /* Complex entity flag */
    tackint(70, 64);                  /* Polyline type flags: pface mesh */
    tackint(71, n * n + 4 * n);       /* Vertex count = Mesh plus
                                                        baseplate edges */
    /*      Total faces = n^2 +                   mesh tiles
                                1 +               baseplate
                                    4 * (n - 1)   wall tiles. */
    tackint(72, n * n + 1 + 4 * (n - 1));

    /* Tack on extended entity data that records the generation
       parameters of the mesh. */

    tackrb(-3);                       /* Emplace start of application data */

    tackstring(Xed(1), AppName);      /* Application name */
    tackint(Xed(70), wtype);          /* World type */
    tackreal(Xed(40), fracdim);       /* Fractal dimension */
    tackreal(Xed(40), powscale);      /* Elevation power factor */
    makent();

    /* Send the vertex table. */

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            ads_point p;

            Spoint(p, i, j,
                   BASE_THICK + (n / (rmax * 2)) * max(0, Real(a, i, j)));

            makevtx(p);
        }
    }

    /* Send the baseplate vertices.  We compute and store the
       four edges of the baseplate, addressed by edge number from
       0 to 3 and vertex number along the edge. */

    for (i = 0; i < n; i++) {
        ads_point p;

        Spoint(p, i, 0, 0);
        makevtx(p);
        Spoint(p, i, n - 1, 0);
        makevtx(p);
        Spoint(p, 0, i, 0);
        makevtx(p);
        Spoint(p, n - 1, i, 0);
        makevtx(p);
    }

    /* Send the mesh tiles. */

#define Vtx(x, y) ((((x) * n) + (y)) + 1)     /* Address mesh vertex */
#define Bpx(e, x) (((n * n) + 1) + (((x) * 4) + (e)))
    for (i = 0; i < (n - 1); i++) {
        for (j = 0; j < (n - 1); j++) {
            makeface(calccol(a, n, i, j, rmax, rmin),
                     Vtx(i, j),
                     Vtx(i, j + 1) * (j == n - 2 ? 1 : -1),
                     Vtx(i + 1, j + 1) * (i == n - 2 ? 1 : -1),
                     Vtx(i + 1, j));
        }
    }

    /* Draw the baseplate */

    makeface(COL_BASE,  -Bpx(0, 0), -Bpx(1, 0),
             -Bpx(1, n - 1), -Bpx(0, n - 1));

    /* Build the "walls" that connect the baseplate to the edges
       of the model. */

    for (i = 0; i < n - 1; i++) {
        makeface(COL_BASE,
                 -Bpx(0, i),
                 -Vtx(i, 0),
                 -Vtx(i + 1, 0),
                 -Bpx(0, i + 1));

        makeface(COL_BASE,
                 -Bpx(1, i),
                 -Vtx(i, n - 1),
                 -Vtx(i + 1, n - 1),
                 -Bpx(1, i + 1));

        makeface(COL_BASE,
                 -Bpx(2, i),
                 -Vtx(0, i),
                 -Vtx(0, i + 1),
                 -Bpx(2, i + 1));

        makeface(COL_BASE,
                 -Bpx(3, i),
                 -Vtx(n - 1, i),
                 -Vtx(n - 1, i + 1),
                 -Bpx(3, i + 1));
    }
#undef Vtx
#undef Bpx

    /* Finally, tack on the SEQEND entity that triggers generation
       of the whole shebang. */

    defent(/*MSG0*/"SEQEND");
    makent();
}

#undef defent
#undef tackrb

/*  MOUNTAIN  --  Make a mountain.  */

static void mountain()
{
    float *a;
    int i, j, uds;
    static int n = 32;
    char tbuf[80];
    ads_real rmin = 1e50, rmax = -1e50;

    ads_initget(2 + 4, NULL);
    V sprintf(tbuf, /*MSG16*/"\nMesh size (power of 2) <%d>: ", n);
    uds = ads_getint(tbuf, &n);
    if (uds == RTCAN)
        return;
    else if (uds != RTNONE) {
        if (n > 1)
            for (uds = n; (uds & 1) == 0; uds >>= 1) ;
        else
            n = 2;
        if (uds != 1) {
            for (uds = 2; uds < n; uds <<= 1) ;
            n = uds;
        }
    }

    while (True) {
        int uds;

        ads_initget(2 + 4, NULL);
        V sprintf(tbuf, /*MSG17*/"\n\
Fractal dimension (typically between 1 and 3) <%g>: ",
                  fracdim);
        uds = ads_getreal(tbuf, &fracdim);
        if (uds == RTCAN)
            return;
        else if (uds == RTNONE)
            break;
        if (fracdim >= 0.0 && fracdim <= 4.0) {
            break;
        }
        ads_printf(/*MSG18*/"Fractal dimension must be between 0 and 4\n");
    }

    ads_initget(2, NULL);
    V sprintf(tbuf, /*MSG19*/"\nPower law scaling exponent <%g>: ", powscale);
    uds = ads_getreal(tbuf, &powscale);
    if (uds == RTCAN)
        return;

    initseed();
    if (forceseed) {
        initgauss(rseed);
        forceseed = False;
    }

    spectralsynth(&a, n, 3.0 - fracdim);

    /* Apply power law scaling if non-unity scale is requested. */

    if (powscale != 1.0) {
        for (i = 0; i < n; i++) {
            for (j = 0; j < n; j++) {
                if (Real(a, i, j) > 0) {
                    Real(a, i, j) = pow((double) Real(a, i, j), powscale);
                }
            }
        }
    }

    /* Compute extrema for autoscaling. */

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            rmin = min(rmin, Real(a, i, j));
            rmax = max(rmax, Real(a, i, j));
        }
    }

    if (entmake)
        makemesh(a, n, rmax, rmin);
    else
        genmesh(a, n, rmax, rmin);

    free((char *) a);
}

/*  SEED  --  Set seed for random numbers.  */

static void seed()
{
    int uds;
    char tbuf[80], ibuf[134];

    initseed();

    V sprintf(tbuf, /*MSG20*/"Random number seed (0 to reseed) <0x%X>: ",
              rseed);
    uds = ads_getstring(False, tbuf, ibuf);
    if (uds == RTNORM) {
        if (strlen(ibuf))
            rseed = atoi(ibuf);
        if (rseed == 0) {
            seeded = False;
            initseed();
            ads_printf(/*MSG21*/"New seed is 0x%X.\n", rseed);
        }
        forceseed = True;
    }
}

/*  WORLD  --  Set world we're operating in.  */

static void world()
{
    int i;
    char tbuf[80];

    ads_textpage();
    ads_printf(/*MSG22*/"\nWorld:  0 = Landscape (8 colours)");
    ads_printf(/*MSG23*/"\n        1 = Landscape (256 colours)");
    ads_printf(/*MSG24*/"\n        2 = Clouds");
    V sprintf(tbuf, /*MSG25*/"\nWorld type <%d>: ", wtype);
    if (ads_getint(tbuf, &i) == RTNORM) {
        if (i < 0 || i > 2) {
            ads_printf(/*MSG26*/"Invalid world type.\n");
        } else {
            wtype = i;
        }
    }
}

/*  MAKEMODE  --  Set mesh generation method.  */

static void makemode()
{
    int i;
    char tbuf[80];

    ads_printf(
       /*MSG27*/"\nMesh generation:  0 = Directly with ads_entmake()");
    ads_printf(
       /*MSG28*/"\n                  1 = PFACE command via ads_command()");
    V sprintf(tbuf, /*MSG29*/"\nGeneration type <%d>: ", entmake ? 0 : 1);
    if (ads_getint(tbuf, &i) == RTNORM) {
        if (i < 0 || i > 1) {
            ads_printf(/*MSG30*/"Invalid mesh generation mode.\n");
        } else {
            entmake = i == 0 ? True : False;
        }
    }
}



#ifdef  HIGHC

/*  Early versions of High C put abort() in the same module with exit();
    ADS defines its own exit(), so we have to define our own abort().  */

static void abort()
{
    ads_abort("");
}

#endif  /* HIGHC */
