allen@viewlogic.com (Dave Allen) (04/09/91)
Submitted-by: Dave Allen <allen@viewlogic.com>
Posting-number: Volume 18, Issue 2
Archive-name: planet/part02
Supersedes: tec: Volume 10, Issue 77-78
#! /bin/sh
# into a shell via "sh file" or similar. To overwrite existing files,
# type "sh file -c".
# The tool that generated this appeared in the comp.sources.unix newsgroup;
# send mail to comp-sources-unix@uunet.uu.net if you want that tool.
# Contents: ./doc/params.clim ./src/fileio.c ./src/tec2.c ./src/tec3.c
# Wrapped by kent@sparky on Mon Apr 8 22:39:14 1991
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
echo If this archive is complete, you will see the following message:
echo ' "shar: End of archive 2 (of 4)."'
if test -f './doc/params.clim' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'./doc/params.clim'\"
else
echo shar: Extracting \"'./doc/params.clim'\" \(9407 characters\)
sed "s/^X//" >'./doc/params.clim' <<'END_OF_FILE'
XIf you give clim a filename as a command line argument, it will read
Xparameters from that file. If you hand it '-' as an argument, it will
Xread them from stdin. For example, you could type "clim - < foo", or
X"clim foo". The '-' option is handy for pipes.
X
XThe parameter file is optional; all of the parameters you can change
Xhave defaults. Thus, the parameter file need contain only the parameters
Xyou want to change. A parameter file looks like LISP; for example,
Xto change the XSIZE parameter, the file would have the one line
X
X(XSIZE 40)
X
XParameters can also be vectors; for example
X
X(RAINCUT (40 60 110 160 180))
X
XPrinting and sizing parameters:
X
XXSIZE - Horizontal size of arrays. Default 60.
XYSIZE - Vertical size of arrays. Default 30.
XBSIZE - Number of seasons in animation and computation. Default 4. Note
X that the storage arrays are statically sized with the constant
X MAXB (in const.h); to increase BSIZE, you must increase MAXB and
X recompile.
XPRINTMODE - Default 0. The value 4 produces PostScript output; any other
X value produces short summary output.
XPRINTEMP - Default 0. Any nonzero value produces one text map for each
X season, where the entries represent scaled temperature by a single
X hex digit (0..F). A key is also printed, containing the Farenheit
X temperature equivalents for each digit.
XPRINTPR - Default 0. Any nonzero value produces one text map for each
X season, where the entries are 1 for a region of low pressure,
X 2 for a region of high pressure, and 3 for the heat equator
X (a zone of low pressure).
XPRINTWIND - Default 0. Any nonzero value produces one text map for each
X season, where the entries are a single hex digit interpreted as
X a bitmap. Northbound winds are indicated by bit 0, southbound
X by bit 1, eastbound by bit 2 and westbound by bit 3. For example,
X 5 indicates NE winds.
XPRINTRAIN - Default 0. Any nonzero value produces one text map for each
X season, where the entries range from 0..F. Zero indicates no
X rainfall, and F indicates a lot.
XPRINTCLIM - Default 0. Any nonzero value produces one text map for each
X season, where each entry is in the range 0..A. Here is a key:
X 0 Ocean 7 Savannah
X 3 Icebergs 8 Deciduous forest
X 4 Tundra 9 Jungle
X 5 Steppe A Swamp
X 6 Desert
X
XSo to produce a text file showing the climate and rainfall, you could use:
X
X(PRINTCLIM 1) (PRINTRAIN 1)
X
X
XThe rest of the parameters deal with the computational models.
X
XParameters affecting temperature
X
XTILT - Default 23.0. This is the tilt of the planet with respect to
X its plane of orbit, in degrees. Smaller numbers produce less
X seasonality; numbers above 45 violate some of the assumptions
X of the models used.
XECCENT - Default 0.0. The eccentricity of the planet's orbit; this
X parameter affects seasonality as well. Numbers above 0.5 are
X probably unrealistic.
XECCPHASE - Default 0.0. This parameter describes the phase offset of the
X eccentricity with respect to the tilt, in radians. You can
X produce climates with complicated seasonality by varying this.
XLCONST - Default 275.0. The basic temperature for land squares, assuming
X no tilt, eccentricity, or nearby ocean.
XLCOS - Default 45.0. The amount by which land temperatures should vary
X from north pole to equator. Land temperature, ignoring ocean
X effects, varies from LCONST - LCOS/2 at the poles to LCONST +
X LCOS/2 at the equator.
XLTILT - Default 1.0. The fraction of the tilt parameter that should be
X applied to temperature adjustment for land. Typically, land
X temperatures vary more from season to season than the ocean
X temperatures do, so LTILT should be higher than OTILT.
XLSMOOTH - Default 0.6. One equation governs the effect of land on ocean
X temperatures and vice versa. The equation involves LSMOOTH,
X LDIV, OSMOOTH and ODIV. Given the land and sea temperatures, and
X the number of land squares in a 11 x 5 box around the square,
X the final temperature is a weighted sum of the two temperatures.
X The weights are related to LSMOOTH and OSMOOTH, and the importance
X of nearby land is diminished by increasing LDIV or ODIV.
XLDIV - Default 180.0. See above.
XOCONST - Default 275.0. Same as LCONST, only for the ocean.
XOCOS - Default 30.0. Same as LCOS, only for the ocean.
XOTILT - Default 0.2. See LTILT.
XOSMOOTH - Default 0.2. See LSMOOTH.
XODIV - Default 250.0. See LSMOOTH.
X
XParameters affecting pressure
X
XOLTHRESH - Default 1. Ocean pressure zones essentially ignore land masses
X whose radius is equal to or less than this number, like islands.
XOOTHRESH - Default 5. Ocean pressure zones must be at least this many
X squares away from the nearest (non-ignored) land.
XOLMIN - Default 40. If the unscaled temperature of an ocean square is
X greater than OLMIN and less than OLMAX, then that square is a
X low pressure zone.
XOLMAX - Default 65. See above.
XOHMIN - Default 130. If the unscaled temperature of an ocean square is
X greater than OHMIN and less than OHMAX, then that square is a
X high pressure zone.
XOHMAX - Default 180. See above.
XLOTHRESH - Default 3. Land pressure zones essentially ignore ocean bodies
X whose radius is less than or equal to this number, like lakes.
XLLTHRESH - Default 7. Land pressure zones must be at least this many
X squares away from the nearest (non-ignored) ocean.
XLLMIN - Default 220. If the unscaled temperature of a land square is
X greater than LLMIN and less than LLMAX, then that square is a
X low pressure zone.
XLLMAX - Default 255. See above.
XLHMIN - Default 0. If the unscaled temperature of a land square is
X greater than LHMIN and less than LHMAX, then that square is a
X high pressure zone.
XLHMAX - Default 20. See above.
X
XParameters affecting wind
X
XBARSEP - Default 16. Winds are determined from pressure; a smooth
X pressure map ranging from 0..255 is built by interpolating between
X highs and lows. Wind lines are contour lines on this map, and
X BARSEP indicates the pressure difference between lines on the map.
X
XParameters affecting rainfall
X
XMAXFETCH - Default 5. Fetch is the term that describes how many squares a
X given wind line travels over water. A high fetch indicates a
X moist wind. This number is the maximum depth for the tree walking
X algorithm which finds fetch; the effect of wind in one square
X can travel at most this number of squares before stopping.
XRAINCONST - Default 32. This is the base amount of rainfall in each square.
XLANDEL - Default -10. This is the amount by which rainfall is increased
X in every land or mountain square; that is, rainfall goes down.
XMOUNTDEL - Default 32. For each unit of fetch which is stopped by a mountain,
X rainfall in the mountain square increases by this amount.
XFETCHDEL - Default 4. The amount of rainfall in a square is increased by
X this number for each unit of fetch in the square.
XHEQDEL - Default 32. The amount of rainfall in a square is increased by
X this amount if the square is on the heat equator.
XNRHEQDEL - Default 24. The amount of rainfall in a square is increased by
X this amount if the square is next to a square on the heat equator.
XFLANKDEL - Default -24. The amount of rainfall in a square is increased by
X this amount if the square is on the "flank" of a circular wind
X pattern. This happens when the wind blows south.
XNRFDEL - Default -3. The amount of rainfall in a square is increased by
X this amount for each adjacent square which is on a "flank".
X
XParameters affecting climate determination
X
XICEBERGK - Default 263. If an ocean square is below this temperature
X (measured in deg Kelvin) all year round, then the ocean square
X is icebergs.
XTEMPCUT - Default is the vector (0 40 90 120). The climate array found
X in climate.c/climkey is 4 x 5; the first index is based on
X average annual temperature. The temperature is relative, based
X on the range 0..255; this vector determines the cutoff points.
X For example, with the default vector, a scaled temperature of 20
X falls into the first "bin" and 121 falls into the fourth.
XRAINCUT - Default is the vector (40 60 110 160 180). The second index of
X the climate array is based on average annual rainfall, scaled into
X the range 0..255. This vector determines the cutoff points. For
X example, rainfall of 35 falls into the first "bin".
XMTDELTA - Default 20. This is the amount, in degrees Farenheit, by which
X temperature in the mountains is decreased before the climate
X lookup is performed.
END_OF_FILE
if test 9407 -ne `wc -c <'./doc/params.clim'`; then
echo shar: \"'./doc/params.clim'\" unpacked with wrong size!
fi
# end of './doc/params.clim'
fi
if test -f './src/fileio.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'./src/fileio.c'\"
else
echo shar: Extracting \"'./src/fileio.c'\" \(12090 characters\)
sed "s/^X//" >'./src/fileio.c' <<'END_OF_FILE'
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 all of the file input and output routines; putmat()
X is the output routine. The function getparam() is called to start the
X parameter file input; the routines to read in the parameters follow. */
X
X#include "const.h"
X#include "clim.h"
X#include <stdio.h>
X
X/* The input is handled by one function, gettoken(), which returns one
X of the following codes each time it is called. */
X#define TOK_EOF 0
X#define TOK_OPEN 1
X#define TOK_CLOSE 2
X#define TOK_WORD 3
X
X
X/* These parameters are needed by every source file; they are defined
X in main.c and described in params.doc */
X
Xextern int XSIZE, YSIZE, BSIZE, ZCOAST, PRINTMODE;
X
X/* The file pointer is used to read the parameters from a file. Linecount
X is incremented each time a newline is read from the file, so error
X messages can refer to a line number. Getbuf is the buffer each word
X is read into. Getfname is a private copy of the filename used, so it
X can be reread if needed. The two lookup tables are used by putmat();
X scalelut is set up by init(). Change is set by any of the parameter
X input routines which change a parameter. */
Xstatic FILE *fp;
Xstatic int linecount = 1;
Xstatic char getbuf [256];
Xchar scalelut[256], shortlut[] = "0123456789ABCDEF";
Xextern int change[];
X
X
Xusermessage (s) char *s; { fprintf (stderr, "%s\n", s); }
X
X
Xfileinit () { int i; for (i=0; i<256; i++) scalelut[i] = shortlut[i>>4]; }
X
X
Xputmat (s, buf, mode, cra, lra)
X char *s; int buf, mode; unsigned char cra[MAXX][MAXY], lra[MAXX][MAXY]; {
X register int i, j, k;
X
X if (mode == PRINTMODE_CLIM) { psprint (cra, lra, 1); return (0); }
X else if (PRINTMODE == PRINTMODE_GREY) {
X psprint (cra, lra, 0); return (0); }
X
X if (buf > -1) printf ("(%s %d (\n", s, buf);
X else printf ("(%s (\n", s);
X
X if (PRINTMODE == PRINTMODE_LONG)
X for (j=0, k=0; j<YSIZE; j++) for (i=0; i<XSIZE; i++) {
X printf ("%4d", cra[i][j]); if (++k == 15) { printf ("\n"); k = 0; } }
X else if (mode == PRINTMODE_SHORT) for (j=0; j<YSIZE; j++) {
X for (i=0; i<XSIZE; i++) printf ("%c", shortlut[cra[i][j]]);
X if (j < YSIZE-1) printf ("\n"); }
X else if (mode == PRINTMODE_SCALE) for (j=0; j<YSIZE; j++) {
X for (i=0; i<XSIZE; i++) printf ("%c", scalelut[cra[i][j]]);
X if (j < YSIZE-1) printf ("\n"); }
X printf ("))\n"); return (0); }
X
X
Xgetparams (s) char *s; {
X /* This function is called either by init() or picktype(), in main.c.
X In either case, the argument is a string which is supposed to be either
X "-", in which case stdin is opened, or a filename. If the file doesn't
X exist, the program exits via panic(). Once the file is open, the loop
X expects to find paren-delimited pairs of (name value). Value could be
X a vector, like (name (val1 val2)). After a name is found, params(),
X above, is called to recognize the parameter. If there is a syntax
X error or an unrecognized parameter, the program exits via panic().
X Each computation source file has its own parameter function; if all of
X them fail the parameter name is unknown. */
X
X int x;
X
X if (!strcmp (s, "-")) fp = stdin;
X else if (!(fp = fopen (s, "r"))) panic ("Can't find input file");
X while ((x = gettoken (getbuf)) != TOK_EOF) {
X if (x != TOK_OPEN) geterr ("expected open paren");
X if (gettoken (getbuf) != TOK_WORD) geterr ("expected param name");
X if (!mainpar (getbuf)) geterr ("unknown parameter name");
X if (gettoken (getbuf) != TOK_CLOSE) geterr ("expected close paren"); }
X fclose (fp); }
X
X
Xgettoken (s) char *s; {
X /* This is the only function which actually reads from the file. It
X maintains a one-character private unget buffer. It ignores initial
X whitespace, but counts newlines in order to maintain an accurate linecount.
X Open or close parentheses count as tokens, and so does any amount of text
X with no white space or parens. The return code indicates whether an
X open or close paren, end of file, or a word was found. If a word was
X found, it is copied into the string pointed to by s. When a word is found,
X the character which terminated it (whitespace, EOF, or paren) is put into
X the unget buffer to be read the next time gettoken() is called. */
X
X static char buf = 0; char c;
X
X white: if (buf) { c = buf; buf = 0; } else c = getc (fp);
X switch (c) {
X case '\n': linecount++;
X case '\t':
X case ' ' : goto white;
X case EOF : return (TOK_EOF);
X case '(' : return (TOK_OPEN);
X case ')' : return (TOK_CLOSE); }
X text: if ((c==' ')||(c=='\t')||(c=='\n')||(c==')')||(c=='(')||(c==EOF)) {
X buf = c; *s = 0; return (TOK_WORD); }
X else { *s++ = c; c = getc (fp); goto text; } }
X
X
Xgeterr (s) char *s; {
X /* Use panic() to print a nicely formatted message, containing the input
X file line number, describing why the program just spat up. */
X sprintf (getbuf, "Par error: %s on line %d", s, linecount);
X fclose (fp);
X panic (getbuf); }
X
X
Xgetlng (x, chg) int *x, chg; { int n;
X /* Called after reading a name associated with a single int, this
X function just reads a int from the input file and uses atoi() to
X convert it to a int, stored at the address passed in. Note that any
X text, including a string or a real, will be read in; no type checking
X is performed. If the value is different from the previous value of the
X parameter, change[chg] is set to one. */
X
X if (gettoken (getbuf) != TOK_WORD) geterr ("expected int");
X n = atoi (getbuf); if (n != *x) change[chg] = 1; *x = n; }
X
X
Xgetdbl (x, chg) double *x; int chg; { double n;
X /* This function reads in a double format number, like getlng above. If
X the value is different from the previous value, change[chg] is set to 1. */
X
X if (gettoken (getbuf) != TOK_WORD) geterr ("expected double");
X n = atof (getbuf); if (n != *x) change[chg] = 1; *x = n; }
X
X
Xgetmat (dim1, dim2, chg, m)
X int *dim1, *dim2, chg; unsigned char m[MAXX][MAXY]; {
X /* This function expects to find a vector of strings, where each string
X contains only 0-9, A-Z. The characters are converted into the numbers
X 0..36 and stored in the character array whose address is passed into the
X function. If one of the strings is too long, or if there are too many
X strings, the function spits up. Before exiting, the function sets the
X array limits passed in; dim1 is set to the length of the longest string,
X while dim2 is set to the number of strings. If any of the characters,
X or either dimension, is changed, then change[chg] is set to 1. */
X
X int x, i = 0, j = 0, imax = -1; char c;
X
X if (gettoken (getbuf) != TOK_OPEN) geterr ("expected open paren");
X while ((x = gettoken (getbuf)) == TOK_WORD) {
X if (j == *dim2) geterr ("matrix too high");
X for (i=0; getbuf[i]; i++) {
X if ((c = getbuf[i] - '0') > 9) c = 10 + getbuf[i] - 'A';
X if (c != m[i][j]) change[chg] = 1; m[i][j] = c; }
X if (i > *dim1) geterr ("string too long");
X if (i > imax) imax = i; j++; }
X if (x != TOK_CLOSE) geterr ("expected close paren");
X if ((*dim1 != imax) || (*dim2 != j)) change[chg] = 1;
X *dim1 = imax; *dim2 = j; }
X
X
Xgetlvec (dim, v, chg) int *dim, *v, chg; {
X /* This function expects to find a list of ints, starting with an open
X paren and ending with a close paren. The parameter dim should contain the
X compiled-in array limit; if the input file contains more than this number
X of entries, the function spits up. No type checking is done; atoi() is
X used to convert any text to an int. On exit, the dimension is set to
X the right value. If either the dimension, or any of the elements, have
X changed from the previous values, change[chg] is set to one. */
X
X int x, i = 0, n;
X
X if (gettoken (getbuf) != TOK_OPEN) geterr ("expected open paren");
X while ((x = gettoken (getbuf)) == TOK_WORD) {
X if (i == *dim) geterr ("vector too long");
X n = atoi (getbuf);
X if (n != v[i]) change[chg] = 1; v[i++] = n; }
X if (x != TOK_CLOSE) geterr ("expected close paren");
X if (i != *dim) change[chg] = 1; *dim = i; }
X
X
Xgetdim (x, chg) int *x, chg; { int y;
X /* Called right after reading a name associated with an array bound,
X this function reads one int from the input file; if the value is
X greater than the default value assigned above, the user is trying
X to exceed a compiled-in limit. That's an error. */
X
X if (gettoken (getbuf) != TOK_WORD) geterr ("expected int");
X y = atoi (getbuf);
X if (y > *x) geterr ("dimension exceeds reserved space");
X *x = y; change[chg] = 1; }
X
X
Xgetdvec (dim, v, chg) int *dim; double *v; int chg; {
X /* Called right after reading a name associated with a one-dimensional
X array of doubles, this function expects to find a list of doubles
X starting with an open paren and ended by a closing paren. The parameter
X dim contains the compiled-in array limit; if the input file contains
X more than this number of entries, the function complains. No type checking
X is done; atof() is used to convert any text to a double. On exit, the
X dimension is set to the correct value. */
X
X int x, i = 0;
X
X if (gettoken (getbuf) != TOK_OPEN)
X geterr ("expected open paren");
X while ((x = gettoken (getbuf)) == TOK_WORD) {
X if (i == *dim) geterr ("vector is too long");
X v[i++] = atof (getbuf); }
X if (x != TOK_CLOSE) geterr ("expected close paren");
X *dim = i; change[chg] = 1; }
X
X
X#define D ((double) 648 / XSIZE)
X#define XX(x) (72 + ((x) * D))
X#define YY(y) (72 + ((y) * D))
X
Xstatic char *psdef[] = {
X "/b { setgray moveto d 0 rlineto 0 d rlineto d neg 0 rlineto fill",
X "/v { moveto d 0 rlineto stroke",
X "/h { moveto 0 d rlineto stroke",
X "/x { moveto d d rlineto d neg 0 rmoveto d d neg rlineto stroke",
X "/p { moveto e 0 rmoveto 0 d rlineto e neg dup rmoveto d 0 rlineto stroke",
X 0 };
Xstatic double climlut[] = { 0.00, 0.00, 0.00, 0.95, 0.95, 0.85, 0.70, 0.55,
X 0.40, 0.70, 0.70 };
Xextern double greyscale ();
X
X
Xpsprint (cra, lra, doclim)
X unsigned char cra[MAXX][MAXY], lra[MAXX][MAXY]; int doclim; {
X register int i, j, k; double z;
X static int firstpage = 1;
X
X if (firstpage) {
X printf ("%%!PS-Adobe-1.0\n/d %4.2f def /e %4.2f def\n", D, D/2);
X firstpage = 0;
X for (i=0; psdef[i]; i++) printf ("%s } bind def\n", psdef[i]); }
X printf ("%6.2f %6.2f moveto\n", XX(-0.25), YY(-0.25));
X printf ("%6.2f %6.2f lineto\n", XX(YSIZE+0.25), YY(-0.25));
X printf ("%6.2f %6.2f lineto\n", XX(YSIZE+0.25), YY(XSIZE+0.25));
X printf ("%6.2f %6.2f lineto\n", XX(-0.25), YY(XSIZE+0.25));
X printf ("%6.2f %6.2f lineto\nstroke\n", XX(-0.25), YY(-0.25));
X
X if (doclim) {
X for (j=0; j<YSIZE;j++) for (i=0; i<XSIZE; i++) {
X z = climlut[cra[i][j]];
X if (z > 0.0) printf ("%6.2f %6.2f %5.4f b\n", XX(j), YY(i), z); }
X printf ("0 setgray\n");
X for (j=0; j<YSIZE;j++) for (i=0; i<XSIZE; i++) {
X if (cra[i][j] == C_JUNGLE)
X printf ("%6.2f %6.2f p\n", XX(j), YY(i));
X if (cra[i][j] == C_SWAMP)
X printf ("%6.2f %6.2f x\n", XX(j), YY(i)); } }
X
X else for (j=0; j<YSIZE;j++) for (i=0; i<XSIZE; i++) {
X z = greyscale (cra[i][j]);
X if (z > 0.0) printf ("%6.2f %6.2f %5.4f b\n", XX(j), YY(i), z); }
X
X printf ("0 setgray\n");
X if (lra) for (j=0; j<YSIZE;j++) for (i=0; i<XSIZE; i++) {
X if (lra[i][j] & LINE_0V) printf ("%6.2f %6.2f v\n", XX(j), YY(i));
X if (lra[i][j] & LINE_0H) printf ("%6.2f %6.2f h\n", XX(j), YY(i)); }
X
X printf ("1 setgray\n");
X if (lra) for (j=0; j<YSIZE;j++) for (i=0; i<XSIZE; i++) {
X if (lra[i][j] & LINE_1V) printf ("%6.2f %6.2f v\n", XX(j), YY(i));
X if (lra[i][j] & LINE_1H) printf ("%6.2f %6.2f h\n", XX(j), YY(i)); }
X
X printf ("\nshowpage\n"); }
END_OF_FILE
if test 12090 -ne `wc -c <'./src/fileio.c'`; then
echo shar: \"'./src/fileio.c'\" unpacked with wrong size!
fi
# end of './src/fileio.c'
fi
if test -f './src/tec2.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'./src/tec2.c'\"
else
echo shar: Extracting \"'./src/tec2.c'\" \(14644 characters\)
sed "s/^X//" >'./src/tec2.c' <<'END_OF_FILE'
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<MAXPLATE; i++) tarea[i] = 0;
X for (i=1; i<=frag; i++) tarea[kid[i]] += karea[i];
X
X /* Give up unless exactly two regions are large enough */
X for (i=1, count=0; i<=reg; i++) if (tarea[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<XSIZE; i++) for (j=0; j<YSIZE; j++) {
X if (a = r[i][j]) m[src][i][j] = ids[kid[a]];
X if (m[src][i][j] == RIFTMARK) {
X m[src][i][j] = 0; t[src][i][j] = 0; } }
X pfree (old); return (0); } }
X
X /* If execution reaches here, the split operation failed; remove rift */
X if (old) for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++)
X if (m[src][i][j] == RIFTMARK) {
X m[src][i][j] = old; p[old].area++; }
X return (0); }
X
X
Xnewrift (src, old) int src, *old; {
X /* This function randomly picks a center for a new rift, and draws in
X a curving line until the line hits either the coast or another plate.
X If another plate is hit, the rift is invalid and the function returns 0.
X To find a center, the function generates random x,y values until it
X finds one that is at least RIFTDIST squares from any ocean square. If a
X center is found, a random angle is generated; the rift will pass through
X the center at that angle. Next, halfrift() is called twice. Each call
X generates the rift leaving the center in one direction. If everything
X works out, the function returns the id of the plate the rift is on. */
X
X int x, y, lx, rx, ty, by, i, j, tries = 0, which; double tt;
X
X /* Generate a random x, y value */
X getctr: if (tries > 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<rx; i++) for (j=ty; j<by; j++)
X if (!m[src][i][j]) { tries++; goto getctr; }
X
X /* Found a good center, on plate `which'. Put a rift indicator in the */
X /* center. Generate a random angle t, which is really an integer in the */
X /* range 0-499 multiplied by 2 PI / 1000. Call halfrift once for each */
X /* half of the rift; t is the initial angle for the first call, and */
X /* t + PI is the initial angle for the second call. If halfrift() */
X /* returns zero, abort and return 0; otherwise, return the plate id. */
X which = m[src][x][y]; m[src][x][y] = RIFTMARK;
X tt = rnd (500) * TWOMILLIPI;
X if (!halfrift (src, x, y, which, tt)) { *old = which; return (0); }
X tt += PI; if (tt > 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<XSIZE; i++) {
X
X /* If first square or 1/6 chance at each step, digitize */
X if (!i || !rnd (BENDEVERY)) {
X
X /* If not first step, modify angle a little */
X if (i) tt = tt + (rnd (BENDBY<<1) * TWOMILLIPI) - (BENDBY*TWOMILLIPI);
X if (tt > 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<MAXFRAG; k++) { kk[k] = k; karea[k] = 0; }
X checkmouse ();
X
X /* Look at every point in the array */
X for (k=0, i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) {
X
X /* If too many fragments, give up */
X if (k == MAXFRAG) return (0);
X
X /* If this square isn't part of the blob, try the next square */
X if (m[src][i][j] != match) { r[i][j] = 0; goto bottom; }
X
X /* It is part of the blob. Set k1 to the frag id of the square to */
X /* its left, and set k2 to the frag id of the square above it. Note */
X /* that because of the for-loop direction, both of these squares have */
X /* already been processed. */
X k1 = i ? kk [r [i-1] [j]] : 0; k2 = j ? kk [r [i] [j-1]] : 0;
X
X /* If k1 and k2 are both background, start a new fragment */
X if (!k1 && !k2) { r[i][j] = ++k; karea[k]++; goto bottom; }
X
X /* If k1 and k2 are part of the same frag, add this square to it */
X if (k1 && (k1 == k2)) { r[i][j] = k1; karea[k1]++; goto bottom; }
X
X /* If k1 and k2 belong to different frags, merge them by finding */
X /* all the frags merged with max(k1,k2) and merging them instead */
X /* with min(k1,k2). Add k to that fragment as well. */
X if (k1 && k2) {
X if (k2 < k1) { k3 = k1; k1 = k2; k2 = k3; }
X for (l=1; l<=k; l++) if (kk[l] == k2) kk[l] = k1;
X r[i][j] = k1; karea[k1]++; goto bottom; }
X
X /* Default case is that one of k1,k2 is a fragment and the other is */
X /* background. Add k to the fragment. */
X k3 = (k1) ? k1 : k2; r[i][j] = k3; karea[k3]++;
X bottom: continue; }
X
X /* Set up vector kid to map from fragments to regions by using i to count */
X /* unique groups of fragments. A unique group of fragments is */
X /* characterized by kk[k] == k; otherwise, frag k is merged with some */
X /* other fragment. */
X for (i=0, j=1; j<=k; j++) {
X if (j == kk[j]) kid[j] = ++i;
X else kid[j] = kid [kk [j]]; }
X
X /* Make sure the id of the background is zero; set up return values */
X kid[0] = 0; *frag = k; *reg = i; return (1); }
X
X
X
Xnewplates (src, old) int src, old; {
X /* Compute new dx and dy values for plates right after fragmentation. This
X function looks at the rift markers in m[src]; variable old is the index of
X the plate from which the new plates were created. For each plate adjacent
X to the rift, this function subtracts the number of plate squares to the left
X of the rift from the number to the right; this gives some indication of
X whether the plate should move left or right, and how fast. The same is done
X for squares above and below the rift. The results are put into dx[] and
X dy[]. At this point some unscaled movement vector is available for both of
X the new plates. The vectors are then scaled by the relative sizes of the
X plates. The idea is that if one plate is much larger than the other, the
X small one should move faster. New plate structures are allocated for the
X new plates, and the computed dx and dy values are put in them. */
X
X int dx[MAXPLATE], dy[MAXPLATE];
X register int i, j, a; int totarea=0, maxmag=0; double scale, b;
X
X for (i=1; i<MAXPLATE; i++) { dx[i] = 0; dy[i] = 0; ids[i] = 0; }
X checkmouse ();
X
X /* For every point in the array, set a to the region id (kid is the */
X /* lookup table and r contains frag indices); if a is nonzero and */
X /* the rift is adjacent, adjust counters appropriately */
X for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) if (a = kid[r[i][j]]) {
X if ((i-1 > -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<MAXPLATE; i++) if (tarea[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; i<MAXPLATE; i++) if (ids[i]) {
X b = scale * (totarea - tarea[i]);
X p[ids[i]].odx = p[old].odx * MR [p[old].age] + dx[i] * b;
X p[ids[i]].ody = p[old].ody * MR [p[old].age] + dy[i] * b; } }
END_OF_FILE
if test 14644 -ne `wc -c <'./src/tec2.c'`; then
echo shar: \"'./src/tec2.c'\" unpacked with wrong size!
fi
# end of './src/tec2.c'
fi
if test -f './src/tec3.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'./src/tec3.c'\"
else
echo shar: Extracting \"'./src/tec3.c'\" \(14070 characters\)
sed "s/^X//" >'./src/tec3.c' <<'END_OF_FILE'
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<n; j++, k+=dist)
X x[xy][k+inc] = ((x[xy][k]+x[xy][k+dist])>>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<xmax; i++) for (j=ymin; j<ymax; j++) {
X k = x[0][i] + x[1][j] + x[2][(i - j + 64) >> 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<xmax; i++) for (j=ymin; j<ymax; j++)
X m[src][i+xofs][j+yofs] = (m[src][i+xofs][j+yofs] > 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<MAXPLATE; i++) for (j=1; j<MAXPLATE; j++) mm[i][j] = 0;
X for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) {
X m[dest][i][j] = 0; t[dest][i][j] = 0; }
X
X checkmouse ();
X /* Look at every square which belongs to a plate */
X for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) if ((a = m[src][i][j]) > 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<MAXPLATE; a++) lut[a] = a;
X for (a=2; a<MAXPLATE; a++) for (b=1; b<a; b++) if (mm[a][b] || mm[b][a]) {
X aa = lut [a]; bb = lut[b];
X if (aa != bb) if (bump (dest, aa, bb)) {
X lut[aa] = bb;
X for (i=1; i<MAXPLATE; i++) if (lut[i] == aa) lut[i] = bb; } } }
X
X
Xbump (dest, a, b) int dest, a, b; {
X /* Plates a and b have been moved on top of each other by some amount;
X alter their movement rates for a slow collision, possibly merging them.
X The collision "strength" is a ratio of the area overlap (provided by
X move ()) to the total area of the plates involved. A fraction of each
X plate's current movement vector is subtracted from the movement vector
X of the other plate. If the two vectors are now within some tolerance
X of each other, they are almost at rest so merge them with each other. */
X
X double maa, mab, ta, tb, rat, area; register int i, j, x;
X
X checkmouse();
X /* Find a ratio describing how strong the collision is */
X x = mm[a][b] + mm[b][a]; area = p[a].area + p[b].area;
X rat = x / (MAXBUMP + (area / 20)); if (rat > 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<XSIZE; i++) for (j=0; j<YSIZE; j++)
X if (m[dest][i][j] == a) m[dest][i][j] = b;
X p[b].area += p[a].area; pfree (a);
X return (a); }
X
X
Xerode (dest) int dest; {
X /* This function takes the topography in t[dest] and smooths it, lowering
X mountains and raising lowlands and continental margins. It does this by
X stepping across the entire array and doing a computation once for each
X pair of 8-connected pixels. The computation is done by onerode(), below.
X The computation result for a pair is a small delta for each square, which
X is summed in the array e. When the computation is finished, the delta
X is applied; if this pushes an ocean square high enough, it is added to
X an adjacent plate if one can be found. Also, if a land square is eroded
X low enough, it is turned into ocean and removed from its plate. */
X
X register int i, j, x, z, xx;
X
X /* Zero out the array for the deltas first */
X for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) e[i][j] = 0;
X checkmouse();
X
X /* Step across the entire array; each pixel is adjacent to 8 others, and */
X /* it turns out that if four pairs are considered for each pixel, each */
X /* pair is considered exactly once. This is important for even erosion */
X for (i=1; i<XSIZE; i++) for (j=1; j<YSIZE; j++) {
X onerode (dest, i, j, 0, -1);
X onerode (dest, i, j, -1, -1);
X onerode (dest, i, j, -1, 0);
X if (i < XSIZE-1) onerode (dest, i, j, -1, 1); }
X
X /* Now go back across the array, applying the delta values from e[][] */
X for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) {
X z = t[dest][i][j] + e[i][j]; if (z < 0) z = 0; if (z > 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_FILE
if test 14070 -ne `wc -c <'./src/tec3.c'`; then
echo shar: \"'./src/tec3.c'\" unpacked with wrong size!
fi
# end of './src/tec3.c'
fi
echo shar: End of archive 2 \(of 4\).
cp /dev/null ark2isdone
MISSING=""
for I in 1 2 3 4 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 4 archives.
rm -f ark[1-9]isdone
else
echo You still must unpack the following archives:
echo " " ${MISSING}
fi
exit 0
exit 0 # Just in case...
--
Kent Landfield INTERNET: kent@sparky.IMD.Sterling.COM
Sterling Software, IMD UUCP: uunet!sparky!kent
Phone: (402) 291-8300 FAX: (402) 291-4362
Please send comp.sources.misc-related mail to kent@uunet.uu.net.