Usenet 1994 October/usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso/misc/volume18/planet/part02

Internet Message Format | 1991-04-10 | 55KB

From: [e-mail] (Dave Allen)
Newsgroups: comp.sources.misc
Subject: REPOST: v18i002:  planet - planet generation simulator, Part02/04
Message-ID: <[e-mail]>
Date: 9 Apr 91 04:17:18 GMT
Approved: [e-mail]
X-Checksum-Snefru: f4e3061e d6b1e950 3b22b131 feaa3b0c

Submitted-by: Dave Allen <[e-mail]>
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 [e-mail] 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   [e-mail], 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   [e-mail], 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   [e-mail], 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: [e-mail]
Sterling Software, IMD           UUCP:     uunet!sparky!kent
Phone:    [phone]         FAX:      [phone]
Please send comp.sources.misc-related mail to [e-mail].