/* * A program to simulate tectonics. * * Copyright 1988 by Mark Isaak. * You may distribute this however you like, as long as you don't sell it * and you keep this notice in it. * * Wish list: * Do it on a sphere, not a square torus * Give plates angular momentum * Make some constants variable * Horsts and grabens * Let 2 plates fuse together and/or 1 plate split apart sometimes * Astroblemes * Take density of rock into account */ /* * Some additional notes, 1998-06-15, Mark J. Stock, mstock@umich.edu * * I added two commands: 'f' will advance five time steps, and 'g' * will write a 16-bit int, ASCII, grayscale PGM image of the land * elevations. Any problems with that code is from my own fault. * I may get to addressing Mark Isaak's wish list, but probably not * for a while. * * Compile with: * cc -O -o plates plates.c -lm */ #include #include /* * ---------- Defined parameters --------------------------------------------- * * Some notes on Earth for comparison: * Mean elevation of land = 840 m. * Kilamanjaro = 5895 m. * Everest = 8848 m. * Dead Sea = -400 m. * Gread Rift Valley = -150 m. * Andes = 6950 m. * Marianas Trench = -11000 m. * Mean depth of sea = -3810 m. */ #if 1 #define NPLATES 20 /* # of plates in world at genesis */ #define MAXPLATES 27 /* maximum # of plates we can have */ #define NCONTINENTS 5 /* # of plates which are continental */ #define NHOTSPOTS 10 /* # of hot spots */ /* This next value scales the resolution up! */ #define CIRCUMF 100 /* # of land cells in one direction */ #else #define NPLATES 6 /* # of plates in world at genesis */ #define MAXPLATES 7 /* maximum # of plates we can have */ #define NCONTINENTS 3 /* # of plates which are continental */ #define NHOTSPOTS 2 /* # of hot spots */ #define CIRCUMF 12 /* # of land cells in one direction */ #endif #define LANDELEV 800 /* starting elevation of land */ #define OCEANELEV -1000 /* starting elevation of oceans */ #define EPOCHMOVE 5 /* max movement during one epoch */ #define FOLDFACTOR 1.2 /* how much to increase ht. of folds */ #define RIFTELEV -400 /* elevation for rifts */ #define SUBSUMEHEIGHT 1200 /* elevation change for coastal ranges */ #define TRENCHHEIGHT -500 /* elev. change for trench by coastal range */ #define HOTSPOTHEIGHT 600 /* average ht. growth of hot spots */ #define RIFTMINERALS 7 /* minerals in rifts */ #define SUBSUMEMINERALS 3 /* minerals in costal ranges */ #define ARCHIPMINERALS 5 /* minerals in archipelagos */ #define HOTSPOTMINERALS 2 /* minerals in underwater hot spots */ #define ERODEFACTOR 0.023 /* how much of a height to erode */ #define SUBMERGEFACTOR 0.2 /* how being underwater slows erosion */ #define SKYELEV 10000 /* no erosion above this */ #define WAVEACTION 400 /* erosion due entirely to waves */ #define SUBSIDEFACTOR 0.05 /* how fast ocean floor sinks */ #define BASINELEV -9000 /* how far ocean floor can sink */ /* * ---------- Structures ----------------------------------------------------- */ typedef char bool; #define TRUE 1 #define FALSE 0 /* * Tectonic plate */ typedef struct { float xvec, yvec; /* direction of plate movement */ short dx, dy; /* how far do move plate per age */ bool continent; /* whether initially continent or ocean */ } Plate_t; /* * Hot spot */ typedef struct { short x, y; /* location of hot spot */ short howhot; /* how much it spews forth */ } Hotspot_t; /* * A piece of land which moves around */ typedef struct rock_s { short elev; /* elevation */ short elevchange; /* change in elevation due to erosion */ Plate_t *plate; /* which plate it belongs to */ bool moved; /* whether we've moved this iteration */ long minerals; /* how many minerals */ struct rock_s *next; /* next in linked list */ } Rock_t; /* * ---------- Global variables ----------------------------------------------- */ Rock_t *surface[CIRCUMF][CIRCUMF]; /* surface of planet */ Plate_t plates[MAXPLATES]; /* descriptions of plates */ Hotspot_t hotspots[NHOTSPOTS]; /* descriptions of hot spots */ int nplates; /* # of tectonic plates */ Rock_t *freeRockList; /* ptr to list of unused rocks */ short maxelev, minelev; /* elevation extremes */ bool verbose; /* how many maps to print */ FILE *outfile; /* output file */ /* * ---------- Utilities ------------------------------------------------------ */ /* * Return a random number from 0 to n-1. */ int Rnd(n) int n; { return (rand() % n); } /* * Wraparound from one side of the world to the other. */ int Wrap(val) register int val; /* x or y coordinate */ { if (val < 0) val += CIRCUMF; else if (val >= CIRCUMF) val -= CIRCUMF; return (val); } /* * Allocate a rock. */ Rock_t * AllocRock(pl) Plate_t *pl; /* which plate it goes on */ { Rock_t *rp; /* ptr to rock */ if (freeRockList) { rp = freeRockList; freeRockList = rp->next; } else rp = (Rock_t *)malloc(sizeof(Rock_t)); rp->plate = pl; rp->elev = (pl->continent ? LANDELEV : OCEANELEV); rp->elevchange = 0; rp->moved = FALSE; rp->minerals = 0; rp->next = NULL; return (rp); } /* * ---------- Initialization ------------------------------------------------- */ /* * All initialization. */ Init() { outfile = stdout; verbose = FALSE; nplates = NPLATES; SeedPlates(); InitSurface(); SeedHotSpots(); } /* * Put plates at random spots on the planet. * Make NCONTINENTS of them land, the rest ocean. */ SeedPlates() { register int i; register Plate_t *pl; /* ptr to plates */ pl = &plates[0]; for (i = 0; i < NPLATES; ++i, ++pl) { pl->dx = Rnd(CIRCUMF); pl->dy = Rnd(CIRCUMF); /* * Check if this spot is already taken. */ if (surface[pl->dx][pl->dy]) { --pl; --i; continue; /* try again */ } surface[pl->dx][pl->dy] = (Rock_t *)pl; pl->xvec = 0.0; pl->yvec = 0.0; pl->continent = (i < NCONTINENTS); } } /* * Initialize the surface. Allocate a rock for each spot on the surface * and assign it to the nearest plate. */ InitSurface() { register short dx, dy; /* x, y distance */ register long dist; /* distance */ register short x, y; /* location on surface */ register Plate_t *pl; /* plate pointer */ register Plate_t *bestpl; /* closest plate */ register long bestdist; /* distance of closest plate */ register int i; /* plate counter */ for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { bestdist = 2 * CIRCUMF * CIRCUMF; for (pl = &plates[0], i = 0; i < NPLATES; ++i, ++pl) { dx = pl->dx - x; if (dx < 0) dx = -dx; if (dx > CIRCUMF/2) dx = CIRCUMF - dx; dy = pl->dy - y; if (dy < 0) dy = -dy; if (dy > CIRCUMF/2) dy = CIRCUMF - dy; dist = dx * dx + dy * dy; if (dist < bestdist) { bestdist = dist; bestpl = pl; } } surface[x][y] = AllocRock(bestpl); } } } /* * Place some hot spots randomly. */ SeedHotSpots() { register int i; /* hot spot counter */ register Hotspot_t *hsp; /* hot spot pointer */ for (hsp = &hotspots[0], i = 0; i < NHOTSPOTS; ++i, ++hsp) { hsp->x = Rnd(CIRCUMF); hsp->y = Rnd(CIRCUMF); hsp->howhot = (Rnd(HOTSPOTHEIGHT) + Rnd(HOTSPOTHEIGHT) + Rnd(HOTSPOTHEIGHT) + Rnd(HOTSPOTHEIGHT) + 2) >> 1; } } /* * ---------- Land movement -------------------------------------------------- */ /* * Give all the plates directions, move them, and see what happens. * They may move a maximum of EPOCHMOVE cells per epoch. */ Epoch() { register int i; AssignDirections(); for (i = 0; i < EPOCHMOVE; ++i) { SetMovementVectors(i); MovePlates(); MountainBuild(); EruptHotSpots(); Erode(); if (verbose) PrintElev(stdout); } } /* * Assign each of the plates a direction and distance (velocity) to move. * Velocity is added to previous value, but isn't allowed to exceed * EPOCHMOVE in any direction. */ AssignDirections() { register int i; /* plate counter */ register Plate_t *pl; /* ptr to plates */ float dist; /* magnitude of vector */ register int dir; /* direction of vector */ #define DEGtoRAD (3.1415926 / 180.0) for (pl = &plates[0], i = 0; i < nplates; ++i, ++pl) { dir = Rnd(360); dist = (float)(Rnd(101) + Rnd(101) + Rnd(101) - 150) * (float)EPOCHMOVE/150.0; pl->xvec += dist * cos(dir * DEGtoRAD); pl->yvec += dist * sin(dir * DEGtoRAD); if ((dist = hypot(pl->xvec, pl->yvec)) > (float)EPOCHMOVE) { dist = (float)EPOCHMOVE / dist; pl->xvec *= dist; pl->yvec *= dist; } } } /* * Set plates' dx,dy to the vector given by (i+1)*v/EPOCHMOVE - i*v/EPOCHMOVE, * where v = plate's (xvec,yvec) vector. This arrangement makes the sum of * (dx,dy) over an epoch equal to (xvec,yvec). */ SetMovementVectors(iter) int iter; /* which iteration */ { register int i; /* plate counter */ register Plate_t *pl; /* ptr to plates */ for (pl = &plates[0], i = 0; i < nplates; ++i, ++pl) { pl->dx = (short)((iter + 1) * pl->xvec / EPOCHMOVE) - (short)(iter * pl->xvec / EPOCHMOVE); pl->dy = (short)((iter + 1) * pl->yvec / EPOCHMOVE) - (short)(iter * pl->yvec / EPOCHMOVE); } } /* * Move each rock according to its plate's movement vector. * Use the rock's 'next' field to keep track of rocks that pile up. */ MovePlates() { register int x, y; /* where in the world */ register int newx, newy; /* rock's new home */ register Rock_t *rp; /* rock pointer */ register Rock_t *prevrp; /* rock before rp */ register Rock_t *nextrp; /* rock after rp */ for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { prevrp = NULL; for (rp = surface[x][y]; rp; prevrp = rp, rp = nextrp) { nextrp = rp->next; if (rp->moved) continue; /* * Figure out where to move rock. */ rp->moved = TRUE; if ((rp->plate->dx == 0) && (rp->plate->dy == 0)) continue; newx = Wrap(x + rp->plate->dx); newy = Wrap(y + rp->plate->dy); /* * Move rock. */ if (prevrp == NULL) surface[x][y] = rp->next; else prevrp->next = rp->next; rp->next = surface[newx][newy]; surface[newx][newy] = rp; } } } } /* * ---------- Mountain building and the like --------------------------------- */ /* * For each spot on the surface, combine rocks which were piled on top of * each other in MovePlates() into mountains, * and create rifts where plates have spread apart leaving nothing. */ MountainBuild() { register int x, y; /* where in the world */ register Rock_t *rp; /* rock pointer */ register Rock_t *rp2; /* another rock pointer */ for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; if (rp == NULL) Rift(x, y); else while ((rp2 = rp->next) != NULL) { if ((rp->elev > 0) && (rp2->elev > 0)) Fold(x, y); /* land + land */ else Subsume(x, y); /* land or ocean + ocean */ } } } } /* * Create land where plates have spread apart. */ Rift(x, y) int x, y; { register Plate_t *pl; /* which plate new land goes with */ register int px, py; /* offset at which to look for plate */ register long tries; /* # of failures to find plate */ register Rock_t *rp; /* new rock */ /* * Find a plate to go with. */ for (tries = 0; tries < 30; ++tries) { px = Wrap(x + Rnd(3) - 1); py = Wrap(y + Rnd(3) - 1); if ((rp = surface[px][py]) != NULL) { pl = rp->plate; break; } } /* * Can't find a plate; start a new one. */ if (pl == NULL) { if (nplates < MAXPLATES) ++nplates; pl = &plates[nplates - 1]; pl->xvec = pl->yvec = 0.0; pl->dx = pl->dy = 0; pl->continent = FALSE; } /* * Make rock. */ rp = AllocRock(pl); rp->elev = RIFTELEV; rp->minerals = RIFTMINERALS; /* should distinguish between ocean and continental rifts!!! */ surface[x][y] = rp; } /* * Continent hits continent; mountains are pushed up. * New height = taller + 0.5 * shorter. The leftover height is * distributed to other rocks in the area. */ Fold(x, y) int x, y; { register Rock_t *rp0; /* first rock on list */ register Rock_t *rp1; /* short rock */ register Rock_t *rp2; /* tall rock */ register short halfshort; /* half the short elevation */ rp0 = surface[x][y]; if (rp0->next->elev > rp0->elev) { rp1 = rp0; rp2 = rp1->next; } else { rp2 = rp0; rp1 = rp2->next; } halfshort = (rp1->elev >> 1) * FOLDFACTOR; ElevNearby(x, y, rp1, (halfshort + 1) >> 1); ElevNearby(x, y, rp2, halfshort >> 1); /* * Combine rocks and free a rock structure. */ rp0->elev = rp2->elev + halfshort; rp0->minerals = (rp1->minerals + rp2->minerals + 1) >> 1; rp0->plate = Rnd(2) ? rp1->plate : rp2->plate; rp1 = rp0->next; rp0->next = rp1->next; rp1->next = freeRockList; freeRockList = rp1; } /* * Raise (or lower) land near x,y by height. * Determine direction from x,y by rp's plate's vector. */ ElevNearby(x, y, rp, height) int x, y; Rock_t *rp; short height; { float h; /* hypotenuse of rp->plate->vec */ int x2, y2; /* rock to raise */ h = hypot(rp->plate->xvec, rp->plate->yvec); if (h < 0.01) { x2 = x; y2 = y; } else { x2 = Wrap(x + (int)(rp->plate->xvec / h + 1.5) - 1); y2 = Wrap(y + (int)(rp->plate->yvec / h + 1.5) - 1); } if (surface[x2][y2]) surface[x2][y2]->elev += height; } /* * Create a costal range where oceanic plate is subsumed by continental * plate. Simply add a constant height and minerals to the continent. */ Subsume(x, y) int x, y; { register Rock_t *rp0; /* first rock on list */ register Rock_t *rp1; /* underwater rock */ register Rock_t *rp2; /* higher rock */ rp0 = surface[x][y]; if (rp0->next->elev > rp0->elev) { rp1 = rp0; rp2 = rp1->next; } else { rp2 = rp0; rp1 = rp2->next; } ElevNearby(x, y, rp2, TRENCHHEIGHT); if (rp2->elev < 0) rp0->minerals = rp2->minerals + ARCHIPMINERALS; else rp0->minerals = rp2->minerals + SUBSUMEMINERALS; rp0->elev = rp2->elev + SUBSUMEHEIGHT; rp0->plate = rp2->plate; rp1 = rp0->next; rp0->next = rp1->next; rp1->next = freeRockList; freeRockList = rp1; } /* * Increase elevation of hot spots. If underwater, increase minerals. */ EruptHotSpots() { register int i; /* hot spot counter */ register Hotspot_t *hsp; /* hot spot pointer */ register Rock_t *rp; /* rock pointer */ for (hsp = &hotspots[0], i = 0; i < NHOTSPOTS; ++i, ++hsp) { rp = surface[hsp->x][hsp->y]; if (rp->elev <= 0) rp->minerals += HOTSPOTMINERALS; rp->elev += hsp->howhot; } } /* * Wear down mountains. Amount of wear is proportional to the average * difference in height between a rock and the surrounding rocks, * except very tall mountains don't erode much (no weather), * there's more erosion around seashores (waves), and underwater * mountains don't erode much. * Also lower sea floors. * Mark everything unmoved. */ Erode() { register int x, y; /* where in the world */ register Rock_t *rp0; /* rock pointer */ register Rock_t *rp1; /* pointer to nearby rock */ register short elev0, elev1; /* elevation of rp0, rp1 */ register int i, j; /* indices to nearby rocks */ register short diff; /* difference in height, roughly */ register long depth; /* average height of area */ int diffplates; /* nearby rocks belonging to other plates */ int deepplates; /* # of nearby rocks under water */ for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp0 = surface[x][y]; elev0 = rp0->elev; depth = elev0; if (elev0 < 0) elev0 = elev0 * SUBMERGEFACTOR; else if (elev0 > SKYELEV) elev0 = SKYELEV; diffplates = 0; deepplates = 0; for (i = -1; i <= 1; ++i) { for (j = -1; j <= 1; ++j) { /* for each surrounding rock */ if ((i == 0) && (j == 0)) continue; rp1 = surface[Wrap(x + i)][Wrap(y + j)]; elev1 = rp1->elev; if (rp1->plate != rp0->plate) ++diffplates; depth += elev1; if (elev1 < 0) { ++deepplates; elev1 = elev1 * SUBMERGEFACTOR; } else if (elev1 > SKYELEV) elev1 = SKYELEV; /* * Erode according to diff between elev0, elev1. */ diff = elev0 - elev1; if ((elev0 <= 0 && elev1 >= 0) || (elev0 >= 0 && elev1 <= 0)) { diff += WAVEACTION; } diff = ERODEFACTOR * diff + 0.5; rp0->elevchange -= diff; rp1->elevchange += diff; } } /* * If 7/8 of the rocks surrounding us belong to different * plates, join the crowd. */ if (diffplates >= 7) { do { i = Rnd(3) - 2; j = Rnd(3) - 2; if (i == 0 && j == 0) continue; rp1 = surface[Wrap(x + i)][Wrap(y + j)]; } while (rp1->plate == rp0->plate); rp0->plate = rp1->plate; } /* * Compute subsidence of ocean basins. */ if (deepplates >= (rp0->elev >= 0 ? 8 : 6)) { rp0->elevchange -= SUBSIDEFACTOR * (depth/9 - BASINELEV); } } } /* * Incorporate elevchange; mark things moved. */ for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp0 = surface[x][y]; rp0->moved = FALSE; rp0->elev += rp0->elevchange; rp0->elevchange = 0; } } } /* * ---------- Output routines ------------------------------------------------ */ /* * Print 2 maps, showing elevations and plates */ PrintElev(fd) FILE *fd; /* output file */ { register int x, y; char line[2 * CIRCUMF + 5]; /* land along a latitude */ register char *cp; /* ptr into line */ register int e; /* elevation */ maxelev = 0; minelev = 0; for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { e = surface[x][y]->elev; if (e > maxelev) maxelev = e; else if (e < minelev) minelev = e; } } for (y = 0; y < CIRCUMF; ++y) { cp = &line[0]; /* * Print elevations. */ for (x = 0; x < CIRCUMF; ++x) { e = surface[x][y]->elev; if (e <= 0) { e = -e * 6 / (-minelev + 1); *cp++ = ":,.. "[e]; } else { /* land */ e = e * 26 / (maxelev + 1); *cp++ = 'A' + e; } } *cp++ = ' '; /* * Print plates. */ for (x = 0; x < CIRCUMF; ++x) { e = surface[x][y]->plate - &plates[0]; *cp++ = 'A' + e; } *cp++ = '\n'; *cp = '\0'; fputs(line, fd); } fprintf(fd, "max elev = %d, min elev = %d\n", maxelev, minelev); putc('\n', fd); } /* * Print Statistics */ PrintStats() { register int x, y; register Rock_t *rp; /* rock pointer */ register long totland; /* total land elevations */ register long totsea; /* total sea elevations */ long totmineral; /* total minerals */ short maxminerals; /* maximum minerals */ long nland; /* # of rocks above water */ maxelev = 0; minelev = 0; maxminerals = 0; totland = totsea = 0; totmineral = 0; nland = 0; for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; if (rp->elev > maxelev) maxelev = rp->elev; else if (rp->elev < minelev) minelev = rp->elev; if (rp->elev > 0) { totland += rp->elev; ++nland; } else totsea += rp->elev; if (rp->minerals > maxminerals) maxminerals = rp->minerals; totmineral += rp->minerals; } } fprintf(outfile, "Maximum elevation = %d, Minimum elevation = %d, %% land = %5.2f\n", maxelev, minelev, (float)(nland * 100)/(float)(CIRCUMF * CIRCUMF)); fprintf(outfile, "Average land elevation = %d, Average sea elevation = %d\n", totland / nland, totsea / ((long)CIRCUMF * CIRCUMF - nland)); fprintf(outfile, "Overall average elevation = %d\n", (totland + totsea) / ((long)CIRCUMF * CIRCUMF)); fprintf(outfile, "Maximum minerals = %d, Average minerals = %4.2f\n", maxminerals, (float)totmineral / ((long)CIRCUMF * CIRCUMF)); } /* * Print a postscript-readable relief map. * Darker colors are higher; water is white. */ PSPrintLand() { float fmaxelev; /* float version of maxelev */ register int x, y; register Rock_t *rp; /* rock pointer */ maxelev = 0; for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; if (rp->elev > maxelev) maxelev = rp->elev; } } fmaxelev = (float)maxelev / 0.94; fprintf(outfile, "%4.1f 45 { dup mul exch dup mul add 1 exch sub } setscreen\n", (float)300/8.0); fprintf(outfile, "/d %4.2f def\t%% diameter of rock\n", 432.0 / (float)CIRCUMF); fprintf(outfile, "/r {\tsetgray moveto\n"); fprintf(outfile, "\td 0 rlineto 0 d rlineto d neg 0 rlineto fill\n} bind def\n"); for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; if (rp->elev > 0) { fprintf(outfile, "%6.2f %6.2f %5.4f r\n", 108.0 + x * 432.0 / (float)CIRCUMF, 612.0 - y * 432.0 / (float)CIRCUMF, 0.94 - (float)(rp->elev)/fmaxelev); } } } fprintf(outfile, "\nshowpage\n\n"); } /* * Print PostScript representation of oceans. * Land is white; deeper is darker. */ PSPrintOcean() { float fminelev; /* float version of minelev */ register int x, y; register Rock_t *rp; /* rock pointer */ minelev = 0; for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; if (rp->elev < minelev) minelev = rp->elev; } } fminelev = (float)minelev / 0.94; fprintf(outfile, "%4.1f 45 { dup mul exch dup mul add 1 exch sub } setscreen\n", (float)300/8.0); fprintf(outfile, "/d %4.2f def\t%% diameter of rock\n", 432.0 / (float)CIRCUMF); fprintf(outfile, "/r {\tsetgray moveto\n"); fprintf(outfile, "\td 0 rlineto 0 d rlineto d neg 0 rlineto fill\n} bind def\n"); for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; if (rp->elev <= 0) { fprintf(outfile, "%6.2f %6.2f %5.4f r\n", 108.0 + x * 432.0 / (float)CIRCUMF, 612.0 - y * 432.0 / (float)CIRCUMF, 0.94 - (float)(rp->elev)/fminelev); } } } fprintf(outfile, "\nshowpage\n\n"); } /* * By now, you should be able to figure out what this routine does. */ PSPrintMinerals() { register short maxmin; /* maximum minerals */ float fmaxmin; /* float version of maxmin */ register int x, y; register Rock_t *rp; /* rock pointer */ maxmin = 0; for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; if (rp->minerals > maxmin) maxmin = rp->minerals; } } fmaxmin = (float)maxmin; fprintf(outfile, "%4.1f 45 { dup mul exch dup mul add 1 exch sub } setscreen\n", (float)300/8.0); fprintf(outfile, "/d %4.2f def\t%% diameter of rock\n", 432.0 / (float)CIRCUMF); fprintf(outfile, "/r {\tsetgray moveto\n"); fprintf(outfile, "\td 0 rlineto 0 d rlineto d neg 0 rlineto fill\n} bind def\n"); for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; if (rp->minerals > 0) { fprintf(outfile, "%6.2f %6.2f %5.4f r\n", 108.0 + x * 432.0 / (float)CIRCUMF, 612.0 - y * 432.0 / (float)CIRCUMF, 1.0 - (float)(rp->minerals)/fmaxmin); } } } fprintf(outfile, "\nshowpage\n\n"); } /* * Print a postscript-readable relief map. * Darker colors are higher; water is white. */ PGMPrintLand() { float fmaxelev; /* float version of maxelev */ register int x, y; register Rock_t *rp; /* rock pointer */ maxelev = 0; for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; if (rp->elev > maxelev) maxelev = rp->elev; } } fprintf(outfile, "P2\n%d %d\n%d\n",CIRCUMF,CIRCUMF,maxelev); for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; if (rp->elev > 0) fprintf(outfile, "%d\n", rp->elev); else fprintf(outfile, "0\n"); } } /* fprintf(outfile, "\nshowpage\n\n"); */ } /* * ---------- Main routine and interpreter ---------------------------------- */ main() { char buf[128]; /* input buffer */ char *cp; /* ptr into buf */ srand(time(0)); Init(); for (;;) { fprintf(stderr, "-> "); fflush(stderr); gets(buf); for (cp = buf; *cp == ' ' || *cp == '\t'; ++cp) ; switch (*cp++) { case 'e': Epoch(); break; case 'f': Epoch(); Epoch(); Epoch(); Epoch(); Epoch(); break; case 'a': GetOutFile(cp); PrintElev(outfile); break; case 's': GetOutFile(cp); PrintStats(); break; case 'l': GetOutFile(cp); PSPrintLand(); break; case 'o': GetOutFile(cp); PSPrintOcean(); break; case 'm': GetOutFile(cp); PSPrintMinerals(); break; case 'g': GetOutFile(cp); PGMPrintLand(); break; case 'i': { register short x, y; register Rock_t *rp; for (x = 0; x < CIRCUMF; ++x) { for (y = 0; y < CIRCUMF; ++y) { rp = surface[x][y]; rp->next = freeRockList; surface[x][y] = 0; freeRockList = rp; } } } Init(); break; case 'v': verbose = !verbose; break; case 'q': exit(0); default: Help(); break; } } } /* * Print a help message. */ Help() { fprintf(stderr, "e\t\tproceed for an epoch\n"); fprintf(stderr, "f\t\tproceed for five epochs\n"); fprintf(stderr, "a [file]\tprint ascii maps of elevations and plates\n"); fprintf(stderr, "s [file]\tprint statistics\n"); fprintf(stderr, "l [file]\tprint PostScript data for land\n"); fprintf(stderr, "o [file]\tprint PostScript data for oceans\n"); fprintf(stderr, "m [file]\tprint PostScript data for minerals\n"); fprintf(stderr, "g [file]\tprint data for land to PGM image file\n"); fprintf(stderr, "i\t\tstart over\n"); fprintf(stderr, "v\t\ttoggle verbosity (prints ascii map after each age)\n"); fprintf(stderr, "q\t\tquit\n"); fprintf(stderr, "If file is not specified, output is to stdout.\n"); } /* * Read name of output file and open it. */ GetOutFile(cp) register char *cp; { register char *cp2; /* * If outfile is already open, close it. */ if (outfile != stdout) fclose(outfile); /* * Find beginning of name. */ while (*cp != ' ' && *cp != '\t' && *cp != '\0') ++cp; /* find blank */ while (*cp == ' ' || *cp == '\t') ++cp; /* skip blanks */ /* * Find end of name, open named file. */ for (cp2 = cp; *cp2 > ' ' && *cp2 <= '~'; ++cp2) ; /* find end of name */ if (*cp <= ' ' || *cp > '~') outfile = stdout; else { *cp2 = '\0'; if ((outfile = fopen(cp, "w")) == NULL) { fprintf(stderr, "Can't open %s\n", cp); outfile = stdout; } } }