Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!usc!cs.utexas.edu!uunet!sparky!kent From: allen@viewlogic.com (Dave Allen) Newsgroups: comp.sources.misc Subject: v18i004: planet - planet generation simulator, Part04/04 Message-ID: <1991Apr8.002500.9813@sparky.IMD.Sterling.COM> Date: 8 Apr 91 00:25:00 GMT Sender: kent@sparky.IMD.Sterling.COM (Kent Landfield) Organization: Sterling Software, IMD Lines: 1333 Approved: kent@sparky.imd.sterling.com X-Checksum-Snefru: a04685d2 08008cff ee6a319f 9286e4c8 Submitted-by: Dave Allen Posting-number: Volume 18, Issue 4 Archive-name: planet/part04 Supersedes: tec: Volume 10, Issue 77-78 #! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh src/const.h <<'END_OF_src/const.h' X#define MAXX 100 X#define MAXY 100 X#define PRINTMODE_NONE 0 X#define PRINTMODE_LONG 1 X#define PRINTMODE_SCALE 2 X#define PRINTMODE_SHORT 3 X#define PRINTMODE_GREY 4 X#define PRINTMODE_CLIM 5 X#define DRAW_GREY 1 X#define DRAW_LAND 2 X#define DRAW_CLIM 3 X#define DRAW_TEC 4 X#define LINE_DIAG 1 X#define LINE_CORN 2 X#define LINE_NONE 3 X X#define ABS(xx) ((xx < 0) ? -xx: xx) X#define CMP(x) !strcmp (s, x) X Xextern double cos(), sin(), sqrt(), atof(), asin(), atan2(); END_OF_src/const.h if test 467 -ne `wc -c src/tec.h <<'END_OF_src/tec.h' X#define MAXSPLINTER 10 X#define MAXFRAG 100 X#define MAXPLATE 40 X#define MAXCHANGE 10 X#define REALSCALE 100 X Xstruct plate { int dx, dy, odx, ody, rx, ry, age, area, id, next; }; X Xextern int XSIZE, YSIZE, MAXSTEP; Xextern int MAXLIFE, MAXBUMP, BUMPTOL; Xextern int DRAWEVERY, HYDROPCT, PRINTMODE; Xextern int ZINIT, ZSUBSUME, ZCOAST; Xextern int ZSHELF, ZMOUNTAIN, ZERODE; Xextern int RIFTPCT, DOERODE, ERODERND; Xextern int MAXCTRTRY, RIFTDIST, BENDEVERY; Xextern int BENDBY, SPEEDBASE, SPEEDRNG; Xextern int UNDERSCAN; Xextern double MR[]; END_OF_src/tec.h if test 531 -ne `wc -c src/unix.c <<'END_OF_src/unix.c' X/* This program is Copyright (c) 1991 David Allen. It may be freely X distributed as long as you leave my name and copyright notice on it. X I'd really like your comments and feedback; send e-mail to X allen@viewlogic.com, or send us-mail to David Allen, 10 O'Moore Ave, X Maynard, MA 01754. */ X X#include "const.h" Xextern int XSIZE, YSIZE, MAXSTEP, step; Xint picktype; X Xmain (argc,argv) int argc; char **argv; { X unsigned short seed16v[3]; X X /* Initialize random number generator */ X srand (time ((long *) 0)); /* initialize rand() */ X seed16v[0] = rand(); seed16v[1] = rand(); seed16v[2] = rand(); X seed48 (seed16v); /* initialize lrand48() */ X X init (*++argv); X for (step=0; stepsrc/x.c <<'END_OF_src/x.c' X/* This program is Copyright (c) 1991 David Allen. It may be freely X distributed as long as you leave my name and copyright notice on it. X I'd really like your comments and feedback; send e-mail to X allen@viewlogic.com, or send us-mail to David Allen, 10 O'Moore Ave, X Maynard, MA 01754. X X Based on code written by by George Ferguson (ferguson@cs.rochester.edu), X 28 Apr 1990. That code was intended for an X version of "tec", never X released. */ X X#include X#include X#include X#include X#include X#include X#include X#include X#include "const.h" X#include "clim.h" X X#define SIZE 6 X#define WINDOW_HEIGHT (SIZE*MAXX) X#define WINDOW_WIDTH (SIZE*MAXY) X Xextern int MAXSTEP, XSIZE, YSIZE, step; Xint picktype = M_HEAT; X Xstatic XtAppContext appc; Xstatic Display *Dis; Xstatic Window Win; Xstatic Pixmap Pix; Xstatic Pixmap Pixsave; Xstatic GC gc_def; Xstatic GC Gc[32]; X X Xstatic Arg simple_args[] = { X {XtNfromVert, (XtArgVal) NULL}, X {XtNheight, (XtArgVal) WINDOW_HEIGHT}, X {XtNwidth, (XtArgVal) WINDOW_WIDTH} }; X X Xunsigned char landcols [] = { 11, 28, 3 }, greycols[256], tecols[256]; Xunsigned char climcols [] = { 9, 20, 18, 2, 3, 18, 22, 26, 14, 17, 12 }; Xstatic char *color_names[] = { X "black","white","grey75", "grey50", X "MediumPurple1","purple4","purple3","purple2","purple1", X "blue4","blue3","blue2","blue1", X "DarkGreen","green4","green3","green2","green1", X "DarkGoldenrod4","yellow4","yellow3","yellow2","yellow1", X "orange4","orange3","orange2","orange1", X "brown4","red4","red3","red2","red1" }; X X Xmain (argc,argv) int argc; char **argv; { X Widget toplevel, top_form, surface; X XWindowAttributes wattr; ushort seed16v[3]; X X XtToolkitInitialize (); X appc = XtCreateApplicationContext (); X X Dis = XtOpenDisplay (appc, NULL, NULL, "Demo", NULL, 0, &argc, argv); X toplevel = XtAppCreateShell (NULL, NULL, applicationShellWidgetClass, X Dis, NULL, 0); X top_form = XtCreateManagedWidget ("top form", formWidgetClass, X toplevel, NULL, 0); X surface = XtCreateManagedWidget ("drawing surface", simpleWidgetClass, X top_form, simple_args, XtNumber(simple_args)); X XtRealizeWidget (toplevel); X X Win = XtWindow (surface); X XGetWindowAttributes (Dis, Win, &wattr); X Pix = XCreatePixmap (Dis, Win, WINDOW_WIDTH, WINDOW_HEIGHT, wattr.depth); X Pixsave = XCreatePixmap (Dis, Win, WINDOW_WIDTH, WINDOW_HEIGHT, wattr.depth); X XSelectInput (Dis, Win, X ExposureMask | ButtonPressMask | ButtonReleaseMask | X KeyPressMask | StructureNotifyMask ); X X gc_def = DefaultGC (Dis, DefaultScreen (Dis)); X assign_colors (); X XFlush (Dis); X X /* Initialize random number generator */ X srand (time ((long *) 0)); /* initialize rand() */ X seed16v[0] = rand(); seed16v[1] = rand(); seed16v[2] = rand(); X seed48 (seed16v); /* initialize lrand48() */ X X XFillRectangle (Dis, Pix, Gc[0], 0, 0, MAXX*SIZE, MAXY*SIZE); redraw (); X init (*++argv); X for (step=0; stepsrc/tec1.c <<'END_OF_src/tec1.c' X/* This program is Copyright (c) 1991 David Allen. It may be freely X distributed as long as you leave my name and copyright notice on it. X I'd really like your comments and feedback; send e-mail to X allen@viewlogic.com, or send us-mail to David Allen, 10 O'Moore Ave, X Maynard, MA 01754. */ X X/* This is the file containing all of the important functions except for X trysplit (), which splits a continent into pieces. Also, all of the main X arrays are declared here, even a couple that are only used by functions in X tec2.c. The array declarations are first, followed by the sequencing X function onestep () and some miscellaneous routines including the text X output routines; initialization routines and the routines that do all X the interesting stuff are last. */ X X#include "const.h" X#include "tec.h" X X/* These are the parameters and their default values; each is defined in X params.doc */ Xint XSIZE = 100, YSIZE = 100, MAXSTEP = 100; Xint MAXBUMP = 50, BUMPTOL = 50, UNDERSCAN = 0; Xint HYDROPCT = 70, DRAWEVERY = 0, PRINTMODE = PRINTMODE_NONE; Xint ZINIT = 22, ZSUBSUME = 16, ZCOAST = 16; Xint ZSHELF = 8, ZMOUNTAIN = 48, ZERODE = 16; Xint RIFTPCT = 40, DOERODE = 1, ERODERND = 4; Xint MAXCTRTRY = 50, RIFTDIST = 5, BENDEVERY = 6; Xint BENDBY = 50, SPEEDRNG = 300, SPEEDBASE = 200; X Xdouble MR [] = { 1.0,1.0,1.0,0.7,0.4,0.1,-0.2,-0.5,-0.8,-1.0 }; Xint MAXLIFE = 10; /* Length of MR vector */ Xint change[1]; /* needed for fileio */ X X/* The following arrays are global and most are used by functions in both X source files. The two main ones are m and t. Each is set up to be two X 2-d arrays, where each array is the size of the whole world. M is the X map; elements in m are indices of plates, showing which squares are X covered by which plate. T is the topography; elements in t are altitudes. */ X Xchar m[2][MAXX][MAXY]; unsigned char t[2][MAXX][MAXY]; X X/* Several arrays are used by the binary blob segmenter, segment() in tec2.c. X These include r, which is used to store fragment indices; many fragments X make up one region during a segmentation. Kid is a lookup table; fragment X k belongs to region kid[k] after a segmentation is finished. Karea[k] X is the area of fragment k. */ X Xchar r[MAXX][MAXY], kid[MAXFRAG]; int karea [MAXFRAG]; X X/* The merge routine gets information from the move routine; when the move X routine puts a square of one plate on top of another plate, that information X is recorded in the merge matrix mm. */ X Xchar mm[MAXPLATE][MAXPLATE]; X X/* The erosion routine needs an array to store delta information; during an X erosion, the increases or decreases in elevation are summed in e and then X applied all at once to the topography. */ X Xchar e[MAXX][MAXY]; X X/* Several routines need temporary storage for areas and plate identifiers. */ X Xint tarea[MAXPLATE]; int ids[MAXPLATE]; X X/* The plates in use are stored in this data structure. Dx,dy are the X values to move by THIS STEP ONLY; odx,ody are the permanent move X values; rx,ry are the remainder x and y values used by newdxy() to X determine dx,dy; age is the age of the plate, in steps; area is the X area of the plate, in squares; id is the value in the m array which X corresponds to this plate; next is a pointer to the next occupied X element of the plate array. */ X Xstruct plate p [MAXPLATE]; X X/* The linked list header for available plates and used plates are global, X as is the step counter. */ X Xint pavail, phead, step; X X Xonestep () { X /* This is the sequencing routine called by main once per step. X It just calls the important subfunctions in order: X - trysplit finds a plate to break up, and computes new velocities X - newdxy computes the deltas to move each plate this step X - move moves the plates X - merge determines results when plates rub together X - erode erodes the terrain, adding or subtracting altitude X - draw draw the resulting array once every DRAWEVERY steps X The m and t arrays are double-buffered in the sense that operations go X from m[0] to m[1] or vice-versa; src and dest determine which is which. */ X X int src, dest; X X src = step % 2; dest = 1 - src; X if (rnd (100) < RIFTPCT) trysplit (src); X newdxy (); X move (src, dest); X merge (dest); X if (DOERODE) erode (dest); X draw (DRAW_TEC, LINE_NONE, t[dest], 0); X X if (DRAWEVERY) if (step && !(step % DRAWEVERY)) tecst (dest); X if (!DRAWEVERY && (step == MAXSTEP - 1)) tecst (dest); } X X Xpalloc () { X /* Allocate a plate from the array and return its index. All the fields X of the plate are initialized to 0, except `next'. That field is used to X link together the plate structures in use. */ X X int x; X X if (!pavail) panic ("No more objects"); X x = pavail; pavail = p[x].next; X p[x].next = phead; phead = x; X p[x].area = 0; p[x].age = 0; X p[x].rx = 0; p[x].ry = 0; X p[x].odx = 0; p[x].ody = 0; X p[x].dx = 0; p[x].dy = 0; X return (x); } X X Xpfree (n) int n; { X /* Return a plate array element to the pool of available elements. X To check for infinite loops, the variable guard is incremented X at each operation; if the number of operations exceeds the maximum X possible number, the program panics. */ X X int i, guard = 0; X X if (phead == n) phead = p[n].next; X else { X for (i=phead; p[i].next!=n; i=p[i].next) X if (++guard > MAXPLATE) panic ("Infinite loop in pfree"); X p[i].next = p[n].next; } X p[n].next = pavail; pavail = n; } X X Xmainpar (s) char *s; { X if (CMP ("XSIZE")) getdim (&XSIZE, 0); X else if (CMP ("YSIZE")) getdim (&YSIZE, 0); X else if (CMP ("MOVERATE")) getdvec (&MAXLIFE, MR, 0); X else if (CMP ("MAXSTEP")) getlng (&MAXSTEP, 0); X else if (CMP ("MAXBUMP")) getlng (&MAXBUMP, 0); X else if (CMP ("BUMPTOL")) getlng (&BUMPTOL, 0); X else if (CMP ("DRAWEVERY")) getlng (&DRAWEVERY, 0); X else if (CMP ("PRINTMODE")) getlng (&PRINTMODE, 0); X else if (CMP ("HYDROPCT")) getlng (&HYDROPCT, 0); X else if (CMP ("ZINIT")) getlng (&ZINIT, 0); X else if (CMP ("ZSUBSUME")) getlng (&ZSUBSUME, 0); X else if (CMP ("ZCOAST")) getlng (&ZCOAST, 0); X else if (CMP ("ZSHELF")) getlng (&ZSHELF, 0); X else if (CMP ("ZMOUNTAIN")) getlng (&ZMOUNTAIN, 0); X else if (CMP ("ZERODE")) getlng (&ZERODE, 0); X else if (CMP ("RIFTPCT")) getlng (&RIFTPCT, 0); X else if (CMP ("DOERODE")) getlng (&DOERODE, 0); X else if (CMP ("ERODERND")) getlng (&ERODERND, 0); X else if (CMP ("MAXCTRTRY")) getlng (&MAXCTRTRY, 0); X else if (CMP ("RIFTDIST")) getlng (&RIFTDIST, 0); X else if (CMP ("BENDEVERY")) getlng (&BENDEVERY, 0); X else if (CMP ("BENDBY")) getlng (&BENDBY, 0); X else if (CMP ("SPEEDBASE")) getlng (&SPEEDBASE, 0); X else if (CMP ("SPEEDRNG")) getlng (&SPEEDRNG, 0); X else if (CMP ("UNDERSCAN")) getlng (&UNDERSCAN, 0); X else return (0); X return (1); } X X Xtecst (src) int src; { X /* This function is called whenever map output is called for. It looks X at the parameter `printmode' to decide between long text, simple text, X and PostScript output formats. Note that the default for this X function is no output at all, corresponding to PRINTMODE_NONE. If only X one output map is desired, then move the coastline up or down to meet the X desired hydrographic percentage. */ X X register int i, j, zcoast; int hist[256], goal; unsigned char sk[MAXX][MAXY]; X X if (!PRINTMODE) return (0); X if (!DRAWEVERY) { X /* Create a histogram of the output array */ X for (i=0; i<256; i++) hist[i] = 0; X for (i=0; i0; zcoast--) X if ((i += hist[zcoast]) > goal) break; X X /* If the new coast level is zero, then there wasn't enough land */ X /* to meet the goal, even going right down to the ocean floor. The */ X /* only possible result is to panic since the goal can't be met. */ X if (!zcoast) panic ("Scaled till oceans dried up"); X ZCOAST = zcoast; } X X if (PRINTMODE != PRINTMODE_SHORT) putmat ("LAND", -1, PRINTMODE, t[src], 0); X else { X for (i=0; i ZMOUNTAIN) sk[i][j] = 2; X else sk[i][j] = 1; } X putmat ("LAND", -1, PRINTMODE, sk, 0); } X return (0); } X X Xdouble greyscale (x) int x; { X /* Called by the PostScript print routine, this function simply computes X the intensity from 0-1 corresponding to the altitude 0-255 */ X if (x < ZCOAST) return ((float) 0); X return (1.0 - ((x > 128) ? 128 : x) / 128.0); } X X Xinit (s) char *s; { X /* This is the catchall function that initializes everything. First, X it calls getparams() in fileio.c to allow the user to set parameters. Next, X it links together the plates onto the free list and starts the used list X at empty. The first plate is created by a fractal technique and then X improved. Finally, the fractal is copied to the data array and drawn. X There are two kinds of improvement done here. First, islands are X eliminated by segmenting the blob and erasing all the regions except X for the biggest. Second, oceans inside the blob (holes) are eliminated X by segmenting the _ocean_ and filling in all regions except the biggest. */ X X int besti, x; register int i, j; X X if (s) if (*s) getparams (s); fileinit (); X X for (i=1; i 0) for (i=1; i 0) for (i=1; isrc/tec2.c <<'END_OF_src/tec2.c' X/* This program is Copyright (c) 1991 David Allen. It may be freely X distributed as long as you leave my name and copyright notice on it. X I'd really like your comments and feedback; send e-mail to X allen@viewlogic.com, or send us-mail to David Allen, 10 O'Moore Ave, X Maynard, MA 01754. */ X X/* This file contains the function trysplit(), which is called from X onestep() in tec1.c, and its subfunctions. One of its subfunctions, X segment(), is also called from init() in tec1.c. Trysplit is the X function in charge of splitting one plate into smaller plates. */ X X X#include "const.h" X#include "tec.h" X X#define PI 3.14159 X#define TWOPI 6.28318 X#define TWOMILLIPI 0.00628318 X X/* RIFTMARK is the temporary indicator placed in the arrays to indicate X the squares a rift has just appeared. The function stoprift() puts X them in, and trysplit() takes them out before anybody can see them. */ X#define RIFTMARK -1 X X/* These are all defined in tec1.c */ Xextern char m[2][MAXX][MAXY], r[MAXX][MAXY], kid[MAXFRAG]; Xextern unsigned char t[2][MAXX][MAXY]; Xextern int karea[MAXFRAG], tarea[MAXPLATE], ids[MAXPLATE], step; Xextern struct plate p [MAXPLATE]; X Xtrysplit (src) int src; { X /* Trysplit is called at most once per step in only 40% of the steps. X It first draws a rift on one of the plates, then it segments the result X into some number of new plates and some splinters. If exactly two new X non-splinter plates are found, new plate structures are allocated, new X dx and dy values are computed, and the old plate is freed. If anything X goes wrong, the rift is erased from the array, returning the array to its X previous state. The functions newrift, segment and newplates do most X of the work. */ X X register int i, j, a; int count, old, frag, reg; X X if (newrift (src, &old)) if (segment (src, old, &frag, ®)) if (reg > 1) { X X /* Set tarea[i] to areas of the final segmented regions */ X for (i=0; i MAXSPLINTER) count++; X if (count == 2) { X X /* Compute new dx,dy; update m with the ids of the new plates */ X newplates (src, old); X for (i=0, count=0; i MAXCTRTRY) { *old = 0; return (0); } X x = rnd (XSIZE); y = rnd (YSIZE); X X /* If the location is ocean, try again */ X if (!m[src][x][y]) { tries++; goto getctr; } X X /* Set lx,rx,ty,by to the coordinate values of a box 2*RIFTDIST on a side */ X /* centered on the center. Clip the values to make sure they are on the */ X /* array. Loop through the box; if a point is ocean, try another center. */ X lx = (x < RIFTDIST) ? 0 : x - RIFTDIST; X rx = (x > XSIZE - RIFTDIST - 1) ? XSIZE - 1 : x + RIFTDIST; X ty = (y < RIFTDIST) ? 0 : y - RIFTDIST; X by = (y > YSIZE - RIFTDIST - 1) ? YSIZE - 1 : y + RIFTDIST; X for (i=lx; i TWOPI) tt -= TWOPI; X if (!halfrift (src, x, y, which, tt)) { *old = which; return (0); } X *old = which; return (1); } X X Xhalfrift (src, cx, cy, which, tt) int src, cx, cy, which; double tt; { X /* Draw a rift from cx,cy on plate `which' at angle t. At the beginning, X digitize the angle using Bresenham's algorithm; once in a while thereafter, X modify the angle randomly and digitize it again. For each square travelled, X call stoprift() to see if the rift has left the plate. */ X X int ddx, ddy, rdx, rdy, draw, i, a; double dx, dy, adx, ady; X X checkmouse (); X /* For-loop against SIZE to guard against infinite loops */ X for (i=0; i TWOPI) tt -= TWOPI; if (tt < 0) tt += TWOPI; X X /* Compute dx and dy, scaled so that larger is exactly +1.0 or -1.0 */ X dy = sin (tt); dx = cos (tt); adx = ABS(dx); ady = ABS(dy); X if (adx > ady) { dy = dy / adx; dx = (dx < 0) ? -1.0: 1.0; } X else { dx = dx / ady; dy = (dy < 0) ? -1.0: 1.0; } X X /* Convert to integer value and initialize remainder */ X /* for each coordinate to half value */ X ddx = REALSCALE * dx; ddy = REALSCALE * dy; X rdx = ddx >> 1; rdy = ddy >> 1; } X X /* Main part of loop, draws one square along line. The basic idea */ X /* of Bresenham's algorithm is that if the slope of the line is less */ X /* than 45 degrees, each time you step one square in X and maybe step */ X /* one square in Y. If the slope is greater than 45, step one square */ X /* in Y and maybe one square in X. Here, if the slope is less than 45 */ X /* then ddx == REALSCALE (or -REALSCALE) and the first call to */ X /* stoprift() is guaranteed. If stoprift returns <0, all is ok; */ X /* if zero, the rift ran into the ocean, so stop now; if positive, the */ X /* rift ran into another plate, which is a perverse condition and the */ X /* rift must be abandoned. */ X rdx += ddx; rdy += ddy; X if (rdx >= REALSCALE) { cx++; rdx -= REALSCALE; draw = 1; } X if (rdx <= -REALSCALE) { cx--; rdx += REALSCALE; draw = 1; } X if (draw == 1) { X a = stoprift (src, cx, cy, which); if (a >= 0) return (!a); } X if (rdy >= REALSCALE) { cy++; rdy -= REALSCALE; draw = 2; } X if (rdy <= -REALSCALE) { cy--; rdy += REALSCALE; draw = 2; } X if (draw == 2) { X a = stoprift (src, cx, cy, which); if (a >= 0) return (!a); } } X return (1); } X X Xstoprift (src, x, y, which) int src, x, y, which; { X /* This function is called once for each square the rift enters. It X puts a rift marker into m[src] and decides whether the rift can go on. X It looks at all four adjacent squares. If one of them contains ocean X or another plate, return immediately so that the rift stops (if ocean) X or aborts (if another plate). If none of them do, then check to make X sure the point is within underscan boundaries. If so, return ok. */ X X register int w, a; X X w = which; p[w].area--; m[src][x][y] = RIFTMARK; X a = m[src][x][y+1]; if ((a != w) && (a!= RIFTMARK)) return (a != 0); X a = m[src][x][y-1]; if ((a != w) && (a!= RIFTMARK)) return (a != 0); X a = m[src][x+1][y]; if ((a != w) && (a!= RIFTMARK)) return (a != 0); X a = m[src][x-1][y]; if ((a != w) && (a!= RIFTMARK)) return (a != 0); X if ((x < UNDERSCAN) || (x > XSIZE - UNDERSCAN)) return (1); X if ((y < UNDERSCAN) || (y > YSIZE - UNDERSCAN)) return (1); X return (-1); } X X Xsegment (src, match, frag, reg) int src, match, *frag, *reg; { X /* This routine implements a standard binary-blob segmentation. It looks X at the array m[src]; match is the value of the blob, and everything else X is background. The result is placed into array r and vectors kid and karea. X One 8-connected region can be made up of many fragments; each fragment is X assigned a unique index. Array r contains the frag indices k, while kid[k] X is the region frag k belongs to and karea[k] is the area of frag k. X Variables frag and reg are set on output to the number of fragments and X regions found during the segmentation. The private vector kk provides one X level of indirection for merging fragments; fragment k is merged with X fragment kk[k] where kk[k] is the smallest frag index in the region. */ X X register int i, j, k, k1, k2, k3, l; X char kk [MAXFRAG]; X X /* Initialize all frag areas to zero and every frag to merge with itself */ X for (k=0; k -1) && (m[src][i-1][j] == RIFTMARK)) (dx[a])++; X if ((i+1 < XSIZE) && (m[src][i+1][j] == RIFTMARK)) (dx[a])--; X if ((j-1 > -1) && (m[src][i][j-1] == RIFTMARK)) (dy[a])++; X if ((j+1 < XSIZE) && (m[src][i][j+1] == RIFTMARK)) (dy[a])--; } X X /* For those regions larger than splinters (tarea is set up in trysplit), */ X /* allocate a new plate structure and initialize its area; compute the */ X /* magnitude of the dx dy vector and remember the maximum magnitude; also */ X /* record the total area of new regions */ X for (i=1; i MAXSPLINTER) { X ids[i] = palloc (); p[ids[i]].area = tarea[i]; X totarea += tarea[i]; X a =sqrt ((double) ((dx[i]*dx[i]) + (dy[i]*dy[i]))); X if (a > maxmag) maxmag = a; } X X /* Generate a random speed and predivide so that all speeds computed */ X /* below are less than the random speed. */ X scale = (double) (rnd (SPEEDRNG) + SPEEDBASE) / (maxmag * totarea); X X /* Compute the dx and dy for each new plate; note that the speed the */ X /* plate was moving at before splitting is given by p[old].odx,ody */ X /* but those must be multiplied by MR to get the actual values */ X for (i=1; isrc/tec3.c <<'END_OF_src/tec3.c' X/* This program is Copyright (c) 1991 David Allen. It may be freely X distributed as long as you leave my name and copyright notice on it. X I'd really like your comments and feedback; send e-mail to X allen@viewlogic.com, or send us-mail to David Allen, 10 O'Moore Ave, X Maynard, MA 01754. */ X X#include "const.h" X#include "tec.h" X X/* These are all defined in tec1.c */ Xextern char m[2][MAXX][MAXY], r[MAXX][MAXY], e[MAXX][MAXY], kid[MAXFRAG], X mm[MAXPLATE][MAXPLATE]; Xextern unsigned char t[2][MAXX][MAXY]; Xextern int step, phead, tarea[MAXPLATE], karea[MAXFRAG]; Xextern struct plate p [MAXPLATE]; X X Xmakefrac (src, match) int src, match; { X /* This function uses a very simple fractal technique to draw a blob. Four X one-dimensional fractals are created and stored in array x, then the X fractals are superimposed on a square array, summed and thresholded to X produce a binary blob. Squares in the blob are set to the value in `match'. X A one-dimensional fractal of length n is computed like this. First, X set x [n/2] to some height and set the endpoints (x[0] and x[n]) to 0. X Then do log-n iterations. The first iteration computes 2 more values: X x[n/4] = average of x[0] and x[n/2], plus some random number, and X x[3n/4] = average of x[n/2] and x[n], plus some random number. The second X iteration computes 4 more values (x[n/8], x[3n/8], ...) and so on. The X random number gets smaller by a factor of two each time also. X X Anyway, you wind up with a number sequence that looks like the cross-section X of a mountain. If you sum two fractals, one horizontal and one vertical, X you get a 3-d mountain; but it looks too symmetric. If you sum four, X including the two 45 degree diagonals, you get much better results. */ X X register int xy, dist, n, inc, i, j, k; char x[4][65]; X int xofs, yofs, xmin, ymin, xmax, ymax, bloblevel, hist[256]; X X /* Compute offsets to put center of blob in center of the world, and X compute array limits to clip blob to world size */ X xofs = (XSIZE - 64) >> 1; yofs = (YSIZE - 64) >> 1; X if (xofs < 0) { xmin = -xofs; xmax = 64 + xofs; } X else { xmin = 0; xmax = 64; } X if (yofs < 0) { ymin = -yofs; ymax = 64 + yofs; } X else { ymin = 0; ymax = 64; } X X for (xy=0; xy<4; xy++) { X /* Initialize loop values and fractal endpoints */ X x [xy] [0] = 0; x [xy] [64] = 0; dist = 32; X x [xy] [32] = 24 + rnd (16); n = 2; inc = 16; X X /* Loop log-n times, each time halving distance and doubling iterations */ X for (i=0; i<5; i++, dist>>=1, n<<=1, inc>>=1) X for (j=0, k=0; j>1) + rnd (dist) - inc; } X X /* Superimpose fractals into the output array. x[0] is horizontal, x[1] */ X /* vertical, x[2] diagonal from top left to bottom right, x[3] diagonal */ X /* from TR to BL. While superimposing, create a histogram of the values.*/ X for (i=0; i<256; i++) hist[i] = 0; X for (i=xmin; i> 1] + x[3][(i + j) >> 1]; X if (k < 0) k = 0; if (k > 255) k = 255; X hist[k]++; m[src][i+xofs][j+yofs] = k; } X X /* Pick a threshhold to get as close to the goal number of squares as */ X /* possible, then go back through the array and adjust it */ X bloblevel = XSIZE * YSIZE * (100 - HYDROPCT) / 100; X for (k=255, i=0; k>=0; k--) if ((i += hist[k]) > bloblevel) break; X for (i=xmin; i k) ? match : 0; } X X Xsinglefy (src, match) int src, match; { X /* This is a subfunction of init() which is called twice to improve the X fractal blob. It calls segment() and then interprets the result. If X only one region was found, no improvement is needed; otherwise, the X area of each region must be computed by summing the areas of all its X fragments, and the index of the largest region is returned. */ X X int i, reg, frag, besti, besta; X X segment (src, match, &frag, ®); X if (reg == 1) return (-1); /* No improvement needed */ X X /* Initialize the areas to zero, then sum frag areas */ X for (i=1; i<=reg; i++) tarea[i] = 0; X for (i=1; i<=frag; i++) tarea [kid[i]] += karea [i]; X X /* Pick largest area of all regions and return it */ X for (i=1, besta=0, besti=0; i<=reg; i++) X if (besta < tarea[i]) { besti = i; besta = tarea[i]; } X return (besti); } X X Xnewdxy () { X /* For each plate, compute how many squares it should move this step. X Multiply the plate's basic movement vector odx,ody by the age modifier X MR[], then add the remainders rx,ry from the last move to get some large X integers. Then turn the large integers into small ones by dividing by X REALSCALE and putting the remainders back into rx,ry. This function also X increases the age of each plate, but doesn't let the age of any plate X go above MAXLIFE. This is done to make sure that MR[] does not need to X be a long vector. */ X X register int i, a; X X /* If there is only a single supercontinent, anchor it */ X for (i=phead, a=0; i; i=p[i].next, a++); X if (a == 1) if (p[phead].odx || p[phead].ody) { X p[phead].odx = 0; p[phead].ody = 0; } X X for (i=phead; i; i=p[i].next) { X a = (p[i].odx * MR[p[i].age]) + p[i].rx; X p[i].dx = a / REALSCALE; p[i].rx = a % REALSCALE; X a = (p[i].ody * MR[p[i].age]) + p[i].ry; X p[i].dy = a / REALSCALE; p[i].ry = a % REALSCALE; X if (p[i].age < MAXLIFE-1) (p[i].age)++; } } X X Xmove (src, dest) int src, dest; { X /* This function moves all the plates that are drifting. The amount to X move by is determined in newdxy(). The function simply steps through X every square in the array; if there's a plate in a square, its new location X is found and the topography is moved there. Overlaps between plates are X detected and recorded so that merge() can resolve the collision; mountain X growing is performed. If two land squares wind up on top of each other, X folded mountains are produced. If a land square winds up where ocean was X previously, that square is the leading edge of a continent and grows a X mountain by subsuming the ocean basin. */ X X register int i, j; int a, b, c, x, y; X X /* Clear out the merge matrix and the destination arrays */ X for (i=1; i 0) { X X /* Add the plate's dx,dy to the position to get the square's new */ X /* location; if it is off the map, throw it away */ X x = p[a].dx + i; y = p[a].dy + j; X if ((x >= XSIZE) || (x < 0) || (y >= YSIZE) || (y < 0)) p[a].area--; X else { /* It IS on the map */ X X /* If the destination is occupied, remove the other guy but */ X /* remember that the two plates overlapped; set the new height */ X /* to the larger height plus half the smaller. */ X if (c = m[dest][x][y]) { X (mm[a][c])++; (p[c].area)--; X b = t[src][i][j]; c = t[dest][x][y]; X t[dest][x][y] = (b > c) ? b + (c>>1) : c + (b>>1); } X X /* The destination isn't occupied. Just copy the height. */ X else t[dest][x][y] = t[src][i][j]; X X /* If this square is over ocean, increase its height. */ X if (t[src][x][y] < ZCOAST) t[dest][x][y] += ZSUBSUME; X X /* Plate A now owns this square */ X m[dest][x][y] = a; } } } X X Xmerge (dest) int dest; { X /* Since move has set up the merge matrix, most of the work is done. This X function calls bump once for each pair of plates which are rubbing; note X that a and b below loop through the lower diagonal part of the matrix. X One subtle feature is that a plate can bump with several other plates in X a step; suppose that the plate is merged with the first plate it bumped. X The loop will try to bump the vanished plate with the other plates, which X would be wrong. To avoid this, the lookup table lut is used to provide X a level of indirection. When a plate is merged with another, its lut X entry is changed to indicate that future merges with the vanished plate X should be applied to the plate it has just been merged with. */ X X char lut[MAXPLATE]; int a, aa, b, bb, i; X X for (a=1; a 1.0) rat = 1.0; X X /* Do some math to update the move vectors. This looks complicated */ X /* because a plate's actual movement vector must be multiplied by */ X /* MR[age], and because I have rewritten the equations to maximize */ X /* use of common factors. Trust me, it's just inelastic collision. */ X maa = p[a].area * MR[p[a].age]; mab = p[b].area * MR[p[b].age]; X ta = MR[p[a].age] * area; X p[a].odx = (p[a].odx * maa + p[b].odx * mab * rat) / ta; X p[a].ody = (p[a].ody * maa + p[b].ody * mab * rat) / ta; X tb = MR[p[b].age] * area; X p[b].odx = (p[b].odx * mab + p[a].odx * maa * rat) / tb; X p[b].ody = (p[b].ody * mab + p[a].ody * maa * rat) / tb; X X /* For each axis, compute the remaining relative velocity. If it is */ X /* too large, return without merging the plates */ X if (ABS (p[a].odx*MR[p[a].age] - p[b].odx*MR[p[b].age]) > BUMPTOL) return(0); X if (ABS (p[a].ody*MR[p[a].age] - p[b].ody*MR[p[b].age]) > BUMPTOL) return(0); X X /* The relative velocity is small enough, so merge the plates. Replace */ X /* all references to a with b, free a, and tell merge() a was freed. */ X for (i=0; i 255) z = 255; X X /* If the square just rose above shelf level, look at the four */ X /* adjacent squares. If one is a plate, add the square to that plate */ X if ((z >= ZSHELF) && (t[dest][i][j] < ZSHELF)) { xx = 0; X if (i > 1) if (x = m[dest][i-1][j]) xx = x; X if (i < XSIZE-1) if (x = m[dest][i-1][j]) xx = x; X if (j > 1) if (x = m[dest][i][j-1]) xx = x; X if (j < YSIZE-1) if (x = m[dest][i][j+1]) xx = x; X if (xx) { p[xx].area++; m[dest][i][j] = xx; } } X X /* Add the increment to the old value; if the square is lower than */ X /* shelf level but belongs to a plate, remove it from the plate */ X t[dest][i][j] = z; X if ((z < ZSHELF) && (x = m[dest][i][j])) { X p[x].area--; m[dest][i][j] = 0; } } } X X Xonerode (dest, i, j, di, dj) int dest, i, j, di, dj; { X /* This function is called once per pair of squares in the array. The X amount each square loses to an adjacent but lower square in each step is X one-eighth the difference in altitude. This is coded as a shift right 3 X bits, but since -1 >> 3 is still -1, the code must be repeated to avoid X shifting negative numbers. Also, since erosion is reduced below the X waterline, if an altitude is lower than ZERODE, ZERODE is used instead. */ X X register int z, t1, t2; X X t1 = t[dest][i][j]; t2 = t[dest][i+di][j+dj]; z = 0; X if ((t1 > t2) && (t1 > ZERODE)) { X if (t2 < ZERODE) t2 = ZERODE; X z = ((t1 - t2 + ERODERND) >> 3); } X else if ((t2 > t1) && (t2 > ZERODE)) { X if (t1 < ZERODE) t1 = ZERODE; X z = -((t2 - t1 + ERODERND) >> 3); } X if (z) { e[i][j] -= z; e[i+di][j+dj] += z; } } END_OF_src/tec3.c if test 14070 -ne `wc -c src/globe.c <<'END_OF_src/globe.c' X/* This program is Copyright (c) 1991 David Allen. It may be freely X distributed as long as you leave my name and copyright notice on it. X I'd really like your comments and feedback; send e-mail to X allen@viewlogic.com, or send us-mail to David Allen, 10 O'Moore X Avenue, Maynard, MA 01754. */ X X/* This file contains a spinning globe animation */ X X#include "const.h" X#define N ((MAXX / 2) - 1) X#define AREA (MAXX * MAXX) X#define RADIUS 45 X#define PI 3.14159 X#define TWOPI (2.0 * PI) X#define HALFPI (PI / 2.0) X#define SET(n,x,y,i,j) m[n][(x)+N][(y)+N] = m[0][(i)+N][(j)+N] X Xunsigned char m[21][MAXX][MAXY]; Xint rx[] = {0, 0, 0, 0, 0}, ry[] = {0, N, 0, N, 0}; Xint xofs, prime, xflip, sflip, sxflip, change[1], step = 0; Xint MAXSTEP = 10000, PRINTMODE = PRINTMODE_SCALE, XSIZE = MAXX, YSIZE = MAXX; Xint ZCOAST = 0; X X Xinit (s) char *s; { int i; X fileinit (); X if (!s || !(*s)) panic ("no valid filename"); else getparams (s); X for (i=1, xofs=45; i<6; xofs-=10, i++) { X prime = 16-i; xflip = i+5; sflip = 15+i; sxflip = 6-i; compute (); } } X X Xonestep () { static int i = 0; X draw (DRAW_TEC, LINE_NONE, m[i], 0); i++; step++; X if (i == 5) i = 6; if (i == 16) i = 17; if (i == 21) i = 1; } X X Xmainpar (s) char *s; { int vec [AREA], index = AREA; register int i, j, k; X X if (CMP ("LAND")) { X getlvec (&index, vec, 0); X for (k=0, j=0; j 16) ? vec[k] : 16; X return (1); } X return (0); } X X Xcompute () { X register int i, j, top; int x, y; X X rx[0] = xofs + N; rx[4] = rx[0]; rx[2] = xofs - N; X for (j=0; j HALFPI) { bot++; theta -= HALFPI; } X top = bot + 1; theta = theta / HALFPI; X mapr (bot, r, &botx, &boty); mapr (top, r, &topx, &topy); X *i = botx + theta * (topx - botx) + 0.5; X *j = boty + theta * (topy - boty) + 0.5; X if (*i >= N) *i = N-1; if (*j >= N) *j = N-1; } X X Xmapr (n, r, x, y) int n; double r, *x, *y; { X *x = r * (rx[n] - xofs) + xofs; X *y = r * ry[n]; X if (*x >= N) { *y = *x - N; *x = N; if (n == 4) *y = -(*y); } } X X Xpermute (x, y, i, j) int x, y, *i, *j; { X if ((x >= 0) && (y >= 0)) { *i = N - y - 1; *j = x - N + 0; } X if ((x < 0) && (y >= 0)) { *i = y - N + 0; *j = -N - x - 1; } X if ((x >= 0) && (y < 0)) { *i = y + N + 0; *j = N - x - 1; } X if ((x < 0) && (y < 0)) { *i = -N - y - 1; *j = x + N + 0; } } X X Xgreyscale (x) int x; { } END_OF_src/globe.c if test 3176 -ne `wc -c