terrain evolution simulator
1998 (updated : 2026.02.13)
The software is designed to model the evolution of a square piece of terrain. It uses a fractal pattern to generate the terrain, then uses some basic geographic concepts to model erosion due to rainfall.
If you're a Geography or Geology student, perhaps you'll find this useful.
I wrote this in the early 90's for the DOS platform. It required the fastgraf code library. The gods only know where you'll find a compatible version after all these years. In 2026, I had Claude.AI look at it, and what do you know it turned out to mostly work because the bugs cancelled each other out. You can see the new code in action at etherlabs.
/*
THE TERRAIN EVOLUTION SIMULATOR v 1.33
by Michael Werneburg
Originally v 1.21B, 1994-98
v 1.30-1.33 corrected 2026 with the following bug fixes:
1. river_a changed from unsigned char to unsigned int. [v1.30]
The original's unsigned char (0-255) silently overflowed for watersheds
larger than 255 cells, wrapping back to 0-255 and accidentally capping
erosion rates. Now uses unsigned int (0-65535) for accurate watershed
tracking, with an explicit cap of 255 on the effective value used in
erosion formulas to preserve the original's intended erosion rates.
2. Terrain clamping added after erosion loop. [v1.30]
The original's 16-bit int (-32768 to 32767) silently overflowed during
erosion, accidentally preventing runaway feedback. With corrected
arithmetic, an explicit clamp to +/-32000 is added after each erosion
step to serve the same stabilizing purpose.
3. Watershed tracing while-loop reordered. [v1.30]
The original checked f_model[] and delta_z[] BEFORE checking array bounds,
causing out-of-bounds memory access when flow traced off the grid edge.
Also fixed stale vg_row (based on pre-update X1) in the while condition.
Bounds checks now come first, exploiting C's short-circuit evaluation.
4. Deposition neighbour check corrected. [v1.30]
The original's condition (local_DX != 0 && local_DX != 0) checked DX
twice instead of DX and DY, incorrectly excluding cells where DX==0.
Corrected to (local_DX != 0 || local_DY != 0) to exclude only center.
Also reordered so bounds check (SANITY_TEST) precedes array access.
5. Lake ponding bounds check added. [v1.30]
The inner tdx/tdy loop in the ponding code accessed
lakes[(X+DX+tdx)*127+(Y+DY+tdy)] without checking whether the indices
were within the array. Added bounds checking.
6. Wash/rill deposition direction corrected. [v1.31]
The original deposited eroded material at entrain[vg_row+Y1+DY], which
only offset in Y, missing the X (row) component. For east/west flow
(DY==0), material was deposited back on the same cell, cancelling the
erosion entirely. For diagonal flow, it went to the wrong cell. This
caused a strong directional bias in erosion patterns. Now correctly
deposits at entrain[(X1+DX)*127+(Y1+DY)].
7. Lake detection restricted to true depressions. [v1.31]
The original flagged lakes for any cell with f_model==0 (no downhill
direction), which includes both genuine depressions AND flat high-
elevation plateaus. Since river_a starts at 1, any flat area where
adjacent cells differed by <=1mm triggered lake detection, causing
most of the terrain to be incorrectly flagged as lake. Now requires
at least one neighbour to be strictly higher than the cell before
lake detection proceeds.
8. Smoothing divisor changed from unsigned to signed int. [v1.32]
In smooth_4stats(), the divisor variable was declared as unsigned int.
When dividing long sum by unsigned int divisor, C's promotion rules
convert the signed long to unsigned if the types are the same width.
On the original 16-bit DOS target (long=4 bytes, unsigned int=2 bytes)
the unsigned int is promoted to long and the division is correct. But
on any 32-bit platform where long and unsigned int are both 4 bytes
(including WASM32), a negative sum is reinterpreted as a huge positive
number before division, producing wildly wrong smoothed values.
Changing divisor to signed int fixes this portability bug.
9. Tilt applied uniformly after terrain generation. [v1.33]
The original applied the tilt offset inside the fractal generation
loop using the variable 'aplace', which at that point held the index
of whichever neighbour was last written in the dx/dy loop. If the
last direction failed the bounds check, 'aplace' was stale from an
earlier iteration or even a previous source cell. The result was
that tilt was applied to one arbitrary cell per source rather than
systematically across the terrain. Now applied as a uniform linear
ramp to all cells after the fractal loop completes.
This program develops a model terrain surface through the use of a fractal
(self-similar at all scales) framework. It also then models the pattern of
drainage which would develop on that surface, and simulates the effects of
the flowing water on the terrain surface. This results in an evolution of
the terrain surface.
This evolution is enacted through the use of two erosional procedures,
one of which accounts for the combined effects of overland flow and rilling,
the other which accounts for the effects of channel erosion. The first of
these is a function of surface slope and slope length, while the second is a
function of slope and runoff. The equations on which the coded implementations
are based may be found in the accompanying text.
The modelling process involves an artificial terrain 127 x 127 elements
in size. For reasons of accuracy with regards to stream modelling (see text),
the ground dimensions of these positions has been set to 10 metres by 10 metres,
resulting in an array position ground area of 100m�. If the ground scale is to
represent real values, then the elveation values used in the model must also
represent measurable real values. The scale of the vertical range was
determined by calculating the average slope found by establishing a sound
range of values which could be consistently modelled using an integer data
structure to hold elevation values.
This was accomplished as follows. The integer data type meant that
elevation values could range from -32535 to 32535. It was deemed that a
greater range in elevation was desirable as a greater vertical range would
result in a finer measurement of real elevation values: this finer measurement
was imperative as the rates of denudation were in the order of a few
millimetres per year. Since the model used to create the terrain elevation
values was stochastic, the range of values thus created could not be known
until a number of runs offered a good estimate (although this could be roughly
controlled by setting the vertical range of possible values for each random
generation.
The ability to control the results of the vertical scale generation was
used with the goal of achieving an average local slope of 4o throughout the
terrain model. The value of 4o was judged as a reasonable approximation of
global average slope.
Through the fiddling of the vertical range generation to result in an
average slope of 4o, a real elevation value of 1 mm could be assigned to each
vertical unit. This vertical value, described hereafter as the terrain
evolution simulator (TES) scale, means that a range of real elevations of
50 metres is possible.
The process by which the evolution model was implemented is discussed
below, at the head of the declartions of those functions which implement each
step. */
/////////////////////////////////////////////////////////////////////////////
// The following inclusions compile standard C library headers.
#include
//#include
#include
#include
#include
#include
#include
#include "tes.h"
#define toupper(c) ((c)+'A'-'a')
#define SANITY_TEST ((Y+local_DY)<127)&&((X+local_DX)<127)&&((X+local_DX)>=0)&&((Y+local_DY)>=0)
#define SANITY ((Y+(D_xy*dy))<127)&&((X+(D_xy*dx))<127)&&((X+(D_xy*dx))>=0)&&((Y+(D_xy*dy))>=0)
#define SANITY_CHECK ((Y+DY)<127)&&((X+DX)<127)&&((X+DX)>=0)&&((Y+DY)>=0)
#define XY (X*127)+Y
// The following functions are described in the code below where they are defined.
void intro(),
setup_graphics(),
terrain_gen(),
terrain_rec(),
erode(),
legend(int green_or_grey),
ldraw(int x, int y, int col),
draw(int x, int y, int col);
int vari8(unsigned char scale,unsigned char midX,unsigned char midY),
initialize(),
values_gen(unsigned char isit_first),
get_terrain(),
smooth_4stats();
unsigned char huge *f_model; // The use of these data structures is outlined below
/* [v1.30 FIX 1] Changed from unsigned char to unsigned int.
The original silently overflowed at 255, accidentally capping erosion. */
unsigned int huge *river_a;
int huge *terrain;
int huge *delta_z;
int huge *entrain;
int huge *lakes;
unsigned int X, // The X position within the arrays.
Y; // The Y position within the arrays.
unsigned char ask=1, // holder for command-line driven request for no DEM output.
rivers=1;
int highcount = -32000, // This is used to determine the highest terrain elevation.
lowcount = 32000, // This is used to determine the lowest terrain elevation.
old_vmode,
lines,
A, // These four variables are used as limits for windowing procedures which are performed during
B, // smoothing, slope calculation, aspect generation, and watershed calculation.
C,
D,
tilt=0, // the terrain will be titled....
tiltz=0,
maxx, // width of screen
maxy; // height of screen
long avecount; // This is a holder for average slope and elevation values.
char
D_type = 'O', // holder for data type recorded to and read from Idrisi DEMs on disk
noise, // holder for user key hits, so they're out of the way
output_h[30], // string for output of highest elevation, slope values
output_a[30], // string for output of average elevation, slope values
output_l[30]; // string for output of lowest elevation, slope values
float hold_v; // holder for output values.
////////////////////////////////////////MAIN()//FUNCTION///////////////////
/* This function calls the others in order of operation. It first
* initializes the arrays f_model[] and terrain[] through the initialize()
* function. It then initiates the graphics mode in run_bgi(), and generates
* the fractal surface and terrain model with the terrain_gen() function.
* These models are simultaneously draWn to screen. Following this, the
* terrain model is smoothed within the smooth_4stats() function (the name of
* this function comes from its output of the elevation high, low, and
* average values). The slope/aspect generation function values_gen() is
* then performed to create arrays holding the slope and aspect values at
* each point in the model. The initial watershed generation function
* river_gen() is then performed; the model is now ready for the erosion
* simulation. This simulation is carried out through the use of the erode()
* function. Then the graphics mode is exited and the user is alerted to the
* completion of the program. */
void main(int argc, char *argv[]) {
char yesno;
int initmode;
erase();
textcolor(12);
initmode = fg_getmode();
fg_setmode(initmode);
fg_locate(0,26);
fg_setcolor(12);
fg_text("TERRAIN EVOLUTION SIMULATOR", 27);
fg_locate(3,30);
fg_text("M. Werneburg 1994-98", 20);
fg_locate(6,0);
/* These next few lines determine whether the user wishes to view the optional
* introduction, export the generated terrain models at various stages
* throughout the program) to Idrisi-format DEMs, or display rivers during
* the draWing of the terrain model to screen. */
for (yesno = 1; yesno < 5; yesno++) {
if (argv[argc-yesno][0] == ':') intro();
if (argv[argc-yesno][0] == '-') {
if (argv[argc-yesno][1] == 'o') ask = 0;
else if (argv[argc-yesno][1] == 'r') rivers = 0;
}
}
textcolor(15);
fg_locate(15,0);
fg_setcolor(15);
fg_text ("This model is capable of generating tilted terrains. This has the effect of", 76);
fg_locate(16,0);
fg_text ("adding a trend to the drainage pattern.", 39);
fg_locate(17,0);
fg_text ("Do you wish to add a tilt? [y/n] ", 33);
//textcolor(8); // these textcolor routines draW the user's response in dark grey
yesno = getche();
yesno = toupper(yesno);
if (yesno == 'Y') {
fg_locate(18,0);
fg_text ("Please enter a tilt effect in millimetres. The terrain will be titled so", 73);
fg_locate(19,0);
fg_text ("that the top of the screen is raised while the bottom is lowered.", 65);
fg_locate(20,0);
fg_text ("Value (mm):", 11);
scanf ("%d", &tiltz);
// cprintf ("%d", tiltz);
tilt=1;
}
textcolor(15);
/*textcolor(15);
cprintf("You may load an Idrisi-generated Digital Elevation Model into the simulator.\n\r");
cprintf("This file must be of dimensions 127 x 127, and it must have been recorded in\n\r");
cprintf("integer format. Do you wish to load a 127 x 127 DEM for simulation?[y/n]:_\b");
textcolor(8); // these textcolor routines draW the user's response in dark grey
yesno = getche();
yesno = toupper(yesno);
textcolor(15);*/
/* The following if statement determines the course of the program as defined
* by the user's wish to read in a DEM created by the Idrisi GIS. If the
* result is that a DEM is desired, the lone function get_terrain() are
* called. Otherwise, the other set (run_bgi(), terrain_gen(),
* smooth_4stats(), and values_gen()) are run. These two courses reflect
* the two possible paths in establishing the terrain model and the
* evolution process: in the first case, values are read from the DEM, and
* the slope, aspect, and watershed values are calculated from that; in the
* second case, the fractal modelling network (see text) is first built, then
* the terrain model is built upon that, and lastly, the slope, aspect, and
* watershed values are derived. In either case the initialize() function
* is enacted first, and the values_gen() and erode() functions, which
* determine the slope, aspect, and watershed arrays in the first case, and
* perform the actual erosion simulation in thesecond. The final function
* calls (closegraph() and the farfree() functions) do the necessary wrap-up,
* first returning to text mode, then closing the dynamically allocated
* memory blocks used in the simulator.*/
initialize();
//if (yesno == 'Y')
// get_terrain();
//else
// {
setup_graphics();
terrain_gen();
smooth_4stats();
// }
values_gen(1);
erode();
fg_setmode(old_vmode);
erase();
farfree(f_model);
farfree(terrain);
farfree(river_a);
farfree(delta_z);
farfree(entrain);
farfree(lakes);
textcolor(8);
printf("\r\t\t");
cprintf("The TERRAIN EVOLUTION SIMULATOR by Michael Werneburg");
textcolor(7);
}
/////////////////////////////////////INITIALIZE()//FUNCTION///////////////////////////////////////////////////////////////////////
/* This function initializes the arrays f_model[] and terrain[]. These arrays are used to hold the fractal modelling network and
terrain elevation values, respectively. The first array is set to 15 for all points, and the second to 0. The reasoning for this
is as follows.
The function which follows initialize() - terrain_gen() - hunts the array f_model[] for zero values in its first pass. As there
are none but that at the point (63,63) (assigned in the second last line of the function), only that position is selected in the
first iteration of the function. The terrain[] array is set to 0 so that all points in the array have no elevation value at the
outset of the terrain_gen() function.
*/
int initialize() {
int array_place; // set memory block size to 127 by 127
f_model = (unsigned char huge*) farcalloc(16129, 1);
if (f_model == NULL) {
printf("error creating f_model array - insufficient memory!\n");
return(0);
}
terrain = (int huge*) farcalloc(16129, 2);
if (terrain == NULL) {
printf("error creating terrain array - insufficient memory!\n");
return(0);
}
lakes = (int huge*) farcalloc(16129, 2);
if (lakes == NULL) {
printf("error creating lakes array - insufficient memory!\n");
return(0);
}
for(array_place = 0; array_place < 16129; array_place++)
{
f_model[array_place] = 15; // flat fractal model framework
terrain[array_place] = 0; // no terrain
lakes[array_place] = 0; // no lakes
}
// 63 * 127 = 8001 + 63 = 8064
f_model [8064] = 0; // assigns a zero to (63,63)
terrain [8064] = 75; // no elevation at (63,63)
return(1); // This line indicates to main() that the function worked.
}
/////////////////////////////////SETUP_GRAPHICS()//FUNCTION///////////////////////////////////////////////////////////////////////
// This function sets up the video mode, depending on system capabilities.
void setup_graphics() {
unsigned char colour_NO;
fg_svgainit(0);
if (fg_testmode(25, 0)) { // tests for 640 x 480 resolution with 256 colours.
old_vmode = fg_getmode();
fg_setmode(25);
}
else if (fg_testmode(18,0)) { // tests for 640 x 480 resolution with 16 colours if 256 colour mode not possible
old_vmode = fg_getmode();
fg_setmode(18);
}
else {
printf ("\nThe Terrain Evolution Simulator requires VGA graphics capabilities");
printf ("\nto perform.");
exit(1);
}
//if (fg_getmode() == 25)
// {
setrgb(255, 63, 63, 63); // colour 255 is white
setrgb(254, 63, 0, 0); // colour 254 is red
// setrgb(253, 0, 63, 63); // colour 253 is cyan
setrgb(253, 0, 53, 53); // colour 253 is cyan
setrgb(252, 63, 63, 0); // colour 252 is yellow
setrgb(251, 63, 0, 63); // colour 251 is purple
setrgb(151, 0, 0, 0); // black in 236 colour mode
setrgb(152, 0, 0, 63);
setrgb(153, 63, 0, 0); // colour 153 is red
setrgb(154, 0, 63, 0);
setrgb(155, 20, 20, 20);
setcolor(255);
// }
//else
// {
// setrgb(15, 0, 63, 63); // colour 15 is cyan
// setcolor(15);
// }
fg_box (0,200,0,40);
fg_justify(0,0);
move_to(100, 20);
fg_print("Initializing video setup", 24);
//if (fg_getmode() == 25)
// {
for (colour_NO = 1; colour_NO < 55; colour_NO++) // 50 shades of grey
setrgb(colour_NO, (1.26*colour_NO), (1.26*colour_NO), (1.26*colour_NO));
for (colour_NO = 50; colour_NO < 100; colour_NO++) // 50 shades of yellow-green
setrgb(colour_NO, (0.63*(colour_NO-50)), (1.26*(colour_NO-50)), 0);
for (; colour_NO < 150; colour_NO++) // 50 shades of green
setrgb(colour_NO, 0, (1.26*(colour_NO-100)), (0.63*(colour_NO-100)));
// }
//else
// {
// for (colour_NO = 1; colour_NO < 15; colour_NO++)
// setrgb(colour_NO, colour_NO*4, colour_NO*4, colour_NO*4);
// }
maxx = fg_getmaxx(); // should be 640
maxy = fg_getmaxy(); // should be 480
lines = fg_getlines(); // should be 30 text lines on screen
move_to(100, 20);
setcolor(0);
fg_box (0,200,0,40);
fg_print("Initializing video setup", 24);
setcolor(155);
fg_box(127, 511, 47, 431);
setcolor(153);
}
//////////////////////////////////TERRAIN_GEN()//FUNCTION/////////////////////
/* The following block of nested for loops generates the fractal terrain.
* This is accomplished through searching the f_model[] array for a value
* equal to one less than the count of the iteration. That is, on the first
* pass (where the iteration counter "iter8") is equal to one, the array is
* searched for positions with the value 0. In this manner the position (63,
* 63) is found during the fisrt iteration. In the following iterations -
* indicated with an "iter8" value of 2 through 6 - values in the f_model[]
* array of 1 through 5 are hunted. With the discovery of the positions
* containing the proper iteration value(s), eight surrounding points in the
* array are assigned an iteration value equal to iter8. This indicates
* that these points are of the next lesser heirarchy (31, 31), (31, 63), (31, 94), (94, 63),
* (94, 94), (63, 94), (31, 94), and (31, 63) are found from the current
* central point of (63, 63). In the second iteration, the eight points
* surrounding each of these first eight positions are assigned with the
* D_xy value of 16, assigned in LOOP LINE 4. Since a number of these points
* overlap, only 40 positions are determined in iteration two, instead of the
* 64 expected (from 8 first-iteration points * 8 second iteration points).
* This goes on until the fifth generation, at which point the scale value is
* one, and all positions in the array are determined. As each of the new
* positions are found, an elevation value is calculated for that position.
* This elevation value is derived by adding or subtracting a number randomly
* generated for the new position to the elevation value stored in the point
* of higher heirarchy from which the new point's coordinates were found.
* In this fashion, all elevation values across the entire terrain[] array are
* ultimately derived from the value of 0 (which is equal to a real elevation
* of 25 m, being 25000 mm above the zero elevation at -25000 in TES units)
* seeded at (63, 63). For each random number generation, the range of
* possible values is determined by the current iteration, as an argument
* passed to the vari8() function listed below. For the first iteration,
* the range of values is -22400 to 22400; for the second it is -11200 to
* 11200, and so on, until in the final iteration the vertical range is -700
* to 700. This halving of vertical range with every iteration is in keeping
* with the halving of the horizontal distance involved in the derivation of
* the new points' array positions. Lastly, it was realized that the three
* positions with a lesser Y value, as well the two positions with a lesser
* X value than the central position in each iteration following the first
* may have already been determined by a previous interpolation of the same
* iteration.
* This is possible due to the overlap of points based on the scaling system.
* To demonstrate, the three positions (31,47), (47,47), and (47,63) found in
* the second iteration share as a common third-iteration "offspring"
* position the point (39, 55). An elevation value for this point will be
* first assigned as the value of (31, 47) +/- a random value between -11200
* and 11200. It will then be overwritten when the program re-evaluates it
* as the value of (47, 47) +/- ½ 11200, and again when the points around
* (47, 63) are determined. This problem gets worse in frequency with every
* iteration. To counteract this problem without creating as many as four
* more memory-swallowing arrays of 127 x 127 two-byte arrays, part of the
* code below checks first for a pre-existing value for each new position;
* it is for this reason that the terrain[] array was first set to zero for
* all points.
* If a non-zero value is found, it averages the elevation value it was going
* to assign (as found by the process described above) with the pre-existing
* value. This means that the value found from the (31, 47) position would
* be averaged with that found from (47, 47). This averaged value would then
* be averaged with that found from (63, 47), putting an unfortunate but
* inescapable weight on the last value.
* The procedure of averaging outlined above was the most realistic one found
* in creating a model by which each point could be traced in genesis to a
* nearby base point. */
void terrain_gen()
{
int rnd_value, // holder of randomly generated numbers
D_xy=64, // the scale of X and Y change with every iteration, defined as 32 divided by 2 to the power of iter8
tg_row,
tg_XY,
dx, dy,
iter8; // counter of iterations in generation of f_model[] array
unsigned int aplace;
randomize(); // this function initializes the random # generation routine for use.
//text_xy(253, maxx/2, maxy-30, "Generating fractal modelling framework");
for (iter8 = 1; iter8 < 7; iter8++) // LOOP LINE 1
{
D_xy/=2;
for (X = 0; X < 127; X++) // LOOP LINE 9
{
tg_row=X*127;
for (Y = 0; Y < 127; Y++) // LOOP LINE 11
{
tg_XY=tg_row+Y;
if (f_model[tg_XY] == (iter8 - 1)) // LOOP LINE 13
{
/* These lines assign four of the eight f_model[] heirarchy values found
* in each iteration, for those positions diagonally
* surrounding the central points found with each iteration. */
/* The following lines assign the f_model[] values to points in the array
* which are of a distance set by D_xy and on the same horizontal or
* vertical. */
/* These lines create the elevation values at each of the eight positions
* for which f_model[] values were found above. */
for (dx=-1;dx<2;dx++)
{
for (dy=-1;dy<2;dy++)
{
if ((dx || dy) && SANITY)
{
aplace=((X+(D_xy*dx))*127)+(Y+(D_xy*dy));
f_model[aplace] = iter8;
if (terrain[aplace])
terrain [aplace] = (vari8(iter8,X,Y) + terrain[aplace])/2;
else
terrain [aplace] = vari8(iter8,X,Y);
// if (f_model[tg_XY] == (iter8-1))
// draw(X, Y, (49 - (f_model[tg_XY]*8)));
}
}
}
}
} // end of y for loop
} // end of x for loop
} // end of for iter8 loop.
/* [Fix 9] Apply tilt uniformly to all cells after fractal generation.
The original applied tilt inside the fractal loop using a stale
'aplace' variable from the last iteration of the dx/dy loop,
so only one arbitrary neighbour per source cell was tilted. */
if (tiltz) {
int tx, ty;
for (tx = 0; tx < 127; tx++)
for (ty = 0; ty < 127; ty++)
terrain[tx * 127 + ty] += (int)((ty - 63) * (tiltz / 63.0));
}
/* These lines draw the legend of the f_model[] array to the screen beside the
* array. */
/*set_colour(155);
draw_rect(2, 118, 49, 178);
set_colour(254);
fg_box(0, 120, 47, 180);
set_colour(255);
fg_box(1, 119, 48, 179);
draw (-8, 7, 49);draw (-7, 7, 49);draw (-8, 8, 49);draw (-7, 8, 49);
text_xy(151, 50, 73,"Iteration 1");
draw (-8, 12, 41);draw (-7, 12, 41);draw (-8, 13, 41);draw (-7, 13, 41);
text_xy(151, 50, 88,"Iteration 2");
draw (-8, 17, 33);draw (-7, 17, 33);draw (-8, 18, 33);draw (-7, 18, 33);
text_xy(151, 50, 103,"Iteration 3");
draw (-8, 22, 25);draw (-7, 22, 25);draw (-8, 23, 25);draw (-7, 23, 25);
text_xy(151, 50, 118,"Iteration 4");
draw (-8, 27, 17);draw (-7, 27, 17);draw (-8, 28, 17);draw (-7, 28, 17);
text_xy(151, 50, 133,"Iteration 5");
draw (-8, 32, 9); draw (-7, 32, 9); draw (-8, 33, 9); draw (-7, 33, 9);
text_xy(151, 50, 148,"Iteration 6");
text_xy (151, 55, 170, "LEGEND");
text_xy (0, maxx/2, maxy-30, "Generating fractal modelling framework");
text_xy (253, maxx/2, maxy-30, "Fractal modelling framework complete");
text_xy (253, maxx/2, maxy-15, "press any key to continue");
noise = getche();
text_xy (0, maxx/2, maxy-15, "press any key to continue");
set_colour(0);
draw_rect(0, 120, 47, 180); */
text_xy (253, maxx/2, maxy-15, "Generating terrain model");
/* The following lines draw the terrain[] model to the screen, overwriting the
* previously drawn f_model[] array. */
for (X = 0; X < 127; X++)
{
tg_row=X*127;
for (Y = 0; Y < 127; Y++)
{
tg_XY=tg_row+Y;
if (terrain[tg_XY] > 0)
{
// draw (X, Y, (int)(24.5 + (terrain[tg_XY]/1000.0)));
draw (X, Y, (int)(50 + 3.5 * (terrain[tg_XY]/1000.0)));
if (lakes[tg_XY])
ldraw (X, Y, 152);
}
else
draw (X, Y, 253); // Ocean...
}
}
//legend(0); // This call causes the legend to be drawn.
text_xy (0, maxx/2, maxy-30, "Fractal modelling framework complete");
text_xy (0, maxx/2, maxy-15, "Generating terrain model");
text_xy (253, maxx/2, maxy-30, "Terrain model generated");
text_xy (253, maxx/2, maxy-15, "Press any key to continue");
noise = getche();
text_xy (0, maxx/2, maxy-15, "Press any key to continue");
text_xy (253, maxx/2, maxy-15, "Smoothening terrain model");
}
///////////////////////////////////VARI8()//FUNCTION/////////////////////////
/* This is the function - called in the terrain_gen() function above - which
* generates a new elevation value, which it then returns to the calling
* function. It takes as arguments the current scale factor (denoted in the
* parent function terrain_gen() as iter8), as well as the array position
* from which the currently active array position is derived; that is,
* (countX, countY) in terrain_gen(). The function begins by assigning the
* elevation value of the parent array position (midX, midY) to old_value.
* It then creates a new value randomly, defining the range through the
* "scale" variable, which is passed to the function as iter8() by
* terrain_gen(). If the random number generated was greater than half of
* the range, it (rnd_value) is added to old_value. If not, it is subtracted.
* The resulting value is then assigned to new_value, which is then
* returned to the calling function. */
// modified as of Nov 3, 1995 - I'm trying a new formula for the vertical
// difference with each pass. The basic formula is
// z(max)=44800/(1+(pow(0.5,(scale-1)))
int vari8(unsigned char scale,unsigned char midX,unsigned char midY)
{
int new_value, old_value, var_factor;
unsigned int rnd_value;
var_factor= (int) 10000 / (1+(pow(0.75, (scale -1))));
old_value = terrain[(midX*127) +midY]; // the elevation of the central calling position is copied to "old_value".
rnd_value = random(var_factor); // the random number (between zero and the value "var_factor") is assigned to rnd_value
new_value = old_value + ((var_factor/2) - rnd_value); // The new elevation value is assigned by adding a modified rnd_value to
// old_value. Since the modified rnd_value may be below zero, the new value
// may be below old_value.
return(new_value+500); // This line returns the new elevation value to the calling function.
}
//////////////////////////////////SMOOTH_4STATS()//FUNCTION//////////////////
/* The following code is designed to smoothen the terrain, to make its slopes
* more continuous, without jeopardizing the overall shape of the surface.
* This is accomplished through the use of a travelling window as outlined in
* the blurb above the values_gen() code. The window is used such that the
* current active position is set to the average of the four, six, or nine
* positions in the active window. This smoothing removes the unrealistic
* sharpness and obviously artificial terrain shapes which are a result of the
* terrain model generation process.
*/
int smooth_4stats()
{
long sum; // the sum of all terrain elevation values for positions in the window.
int divisor; // used to derive an average of all window elevation values as a divisor of "sum". [Fix 8: changed from
// unsigned int — see header.] This is assigned a value by determining the position of the active window
// in the river_a array.
char across, // the horizontal count across the active window.
along, // the vertical count across the active window
pass; // number of passes the smoothing algorithm makes
int s4_row,
s4_XY;
int huge *tsmooth; // memory block used to hold terrain values during smoothening
tsmooth = (int huge*) farcalloc(16129, 2);
if (terrain == NULL)
{
printf("error creating tsmooth array - insufficient memory!\n");
return(0);
}
//setfillstyle(1,0);
/* This code determines the position of the active window within the river_a
* array by finding the position of the active central point. In this manner,
* a window in the corner is assigned a 2 x 2 dimension for the calucation
* of sum and divisor values, an edge window is assigned a shape of 3 x 2
* positions, and all other windows 3 x 3. */
for (pass = 0; pass < 2; pass++)
{
for (X = 0; X < 127; X++)
{
s4_row=X*127;
for (Y = 0; Y < 127; Y++)
{
s4_XY=s4_row+Y;
sum = 0;
if (X == 0)
{ A = 0; B = 1; divisor = 2; }
else if (X == 126)
{ A = -1; B = 0; divisor = 2; }
else if ((X > 0) && (X < 126))
{ A = -1; B = 1; divisor = 3; }
else
{ A = 0; B = 0; divisor=1; }
if (Y == 0)
{ C = 0; D = 1; divisor *= 2; }
else if (Y == 126)
{ C = -1; D = 0; divisor *= 2; }
else if ((Y > 0) && (Y < 126))
{ C = -1; D = 1; divisor *= 3; }
else
{ A = 0; B = 0; divisor=1; }
/* Using the parameters determined above, the sum of all elevation values
* within the active window is found, and then divided by the appropriate
* value to come up with a sum. This is then assigned to the entrain model,
* and added to the running total avecount for later division by the size of
* the array (16129 positions). */
for (across = A; across <= B; across++)
for (along = C; along <= D; along++)
sum += terrain[((X+across)*127) + (Y+along)];
tsmooth[s4_row + Y] = (int)(sum / divisor);
// tsmooth[s4_row+Y] = terrain[s4_row+Y];
}
}
NOSMOOTH:
/* These lines draw the smoothed terrain model[] to the screen, and record
* the lowest and highest elevation values in the array. */
for (X = 0; X < 127; X++)
{
s4_row=X*127;
for (Y = 0; Y < 127; Y++)
{
s4_XY=s4_row+Y;
terrain[s4_XY] = tsmooth[s4_XY];
if (pass == 1)
{ // 24.5 is just a colour....
if (terrain[s4_XY] > 0)
{
// draw (X, Y, (int)(24.5 + (terrain[s4_XY]/1000.0)));
draw (X, Y, (int)(50 + 3.5 * (terrain[s4_XY]/1000.0)));
if (lakes[s4_XY])
ldraw (X,Y, 152);
}
else
draw (X,Y, 253); // Ocean
if (tsmooth[s4_XY] < lowcount)
lowcount = tsmooth[s4_XY];
else if (tsmooth[s4_XY] > highcount)
highcount = tsmooth[s4_XY];
avecount += tsmooth[s4_XY];
}
}
}
} // end of pass loop
farfree(tsmooth);
text_xy (0, maxx/2, maxy-30, "Terrain model generated");
text_xy (0, maxx/2, maxy-15, "Smoothening terrain model");
text_xy (253, maxx/2, maxy-30, "Terrain model smoothened");
text_xy (253, maxx/2, maxy-15, "press any key to continue");
// These lines display the highest, lowest, and average elevation values in the model.
hold_v = highcount/1000.0;
sprintf (output_h, "Highest elevation %.1f m", hold_v);
text_xy (253, maxx/2, 10, output_h);
hold_v = avecount / (16129000.0);
sprintf (output_a, "Average elevation %.1f m", hold_v);
text_xy (253, maxx/2, 22, output_a);
hold_v = lowcount/1000.0;
sprintf (output_l, "Lowest elevation %.1f m", hold_v);
text_xy (253, maxx/2, 35, output_l);
noise = getch();
text_xy (0, maxx/2, maxy-15, "press any key to continue");
text_xy (0, maxx/2, 10, output_h);
text_xy (0, maxx/2, 22, output_a);
text_xy (0, maxx/2, 35, output_l);
// These lines offer the choice to record the DEM as it currently is to a DEM.
if (!(ask))
{
text_xy (253, maxx/2, maxy-15, "Do you wish to save this model to an Idrisi DEM? [y/n]");
noise = getch();
text_xy (0, maxx/2, maxy-15, "Do you wish to save this model to an Idrisi DEM? [y/n]");
noise = toupper(noise);
if (noise == 'Y')
terrain_rec();
}
text_xy (253, maxx/2, maxy-15, "Now generating aspect and slope values");
return(1);
}
/////////////////////////////////VALUES_GEN()//FUNCTION///////////////////////
/* This function assigns slope and aspect values to each position in the
* terrain[] array. These are held in the delta_z[] and f_model[] arrays,
* respectively. These values are then used to determine the watershed which
* empties into each point in the array, which is equated linearly to the
* flow of water through that 10m x 10m position. At the same time, the
* amount of material which would be transported from each position - a
* function of the already derived slope and watershed values for that
* position. Note that f_model[] is reused here out of an unwillingness to
* pile up the variables.
*/
int values_gen(unsigned char isit_first)
{
int new_value; // variable used to identify locally lowest elevation & highest local slope.
int X1, // temporary X position as watersheds are calculated
Y1, // temporary Y position as watersheds are calculated
DX, // horizontal distance across window
DY; // vertical distance across window
double base, // holder of slope value for power function "pow()"
exponent = 1.4; // holder of exponent of slope for erosional equation
int min_dz = 1000,
vg_row,
tdx,
tdy,
A1,
B1,
C1,
D1,
lakec,
toggle,
max_dz = 0; // minimum and maximum slope value holders.
unsigned int vg_XY,
aplace;
/* During the first pass of this function only, these lines initiate the
* river_a[], delta_z[], and entrain[] arrays. */
if (isit_first)
{
/* [v1.30 FIX 1] Changed from unsigned char to unsigned int */
river_a = (unsigned int huge*) farcalloc(16129, 2);
if (river_a == NULL)
{
printf("error creating river_a array - insufficient memory!\n");
noise = getch();
return(0);
}
delta_z = (int huge*) farcalloc(16129, 2);
if (delta_z == NULL)
{
printf("error creating delta_z array - insufficient memory!\n");
noise = getch();
return(0);
}
entrain = (int huge*) farcalloc(16129, 2);
if (entrain == NULL)
{
printf("error creating entrain array - insufficient memory!\n");
noise = getch();
return(0);
}
}
for (X = 0; X < 127; X++)
{
vg_row=X*127;
for (Y = 0; Y < 127; Y++)
{
vg_XY=vg_row+Y;
river_a[vg_XY] = 1; // initializes the points in river_a[] to one.
delta_z[vg_XY] = 0; // initializes the points in delta_z[] to zero.
entrain[vg_XY] = 0; // initializes the points in entrain[] to zero.
/* The following lines determine the slope at each point in the array toward
* the aspect of greatest down-hill direction. These are assigned to the
* delta_z[] array for all points. The aspect of these recorded slopes are
* kept in the array f_model[] */
if (X == 0) { A = 0; B = 1; }
else if (X == 126) { A = -1; B = 0; }
else { A = -1; B = 1; }
if (Y == 0) { C = 0; D = 1; }
else if (Y == 126) { C = -1; D = 0; }
else { C = -1; D = 1; }
f_model[vg_XY]=0;
for (DX = A; DX <= B; DX++)
{
for (DY = C; DY <= D; DY++)
{
if (DX || DY)
{
aplace=((X+DX)*127)+Y+DY;
switch (DX*DY)
{
case (-1): // ...or
case (1): // We must fake a little trigonometry
new_value = (int)((((terrain[vg_XY]) - terrain[aplace])/1.414213562) + 0.5);
break;
case (0): // Simple math....
new_value = (terrain[vg_XY]) - terrain[aplace];
break;
}
if (new_value > delta_z[vg_XY])
{
delta_z[vg_XY] = new_value;
if ((DX == -1) && (DY == -1))
f_model[vg_XY] = 1; // the aspect value for each position is
else if ((DX == 0) && (DY == -1))
f_model[vg_XY] = 2; // found with these lines, which determine
else if ((DX == 1) && (DY == -1))
f_model[vg_XY] = 3; // the direction based on the X and Y offset,
else if ((DX == 1) && (DY == 0))
f_model[vg_XY] = 4; // which may range from -1 to 1 each. The
else if ((DX == 1) && (DY == 1))
f_model[vg_XY] = 5; // directions assigned range from 1 through 8
else if ((DX == 0) && (DY == 1))
f_model[vg_XY] = 6; // starting in the upper left.
else if ((DX == -1) && (DY == 1))
f_model[vg_XY] = 7;
else if ((DX == -1) && (DY == 0))
f_model[vg_XY] = 8;
}
}
}
} // end of for DY
// if (delta_z[vg_row+Y] == 0) // we're on a plain, so..
// f_model[vg_row+Y] = 0; // there is no aspect
/* [v1.31 FIX 7] Only detect lakes at true local minima
(at least one neighbour strictly higher), not flat plateaus. */
if (f_model[vg_XY] == 0)
{
int is_depression = 0;
for (DX = A; DX <= B && !is_depression; DX++)
for (DY = C; DY <= D && !is_depression; DY++)
if ((DX || DY) && ((Y+DY)<127)&&((X+DX)<127)&&((X+DX)>=0)&&((Y+DY)>=0))
if (terrain[((X+DX)*127)+Y+DY] > terrain[vg_XY])
is_depression = 1;
if (is_depression)
{
for (DX = A; DX <= B; DX++)
{
for (DY = C; DY <= D; DY++)
{
if ((DX || DY) && ((Y+DY)<127)&&((X+DX)<127)&&((X+DX)>=0)&&((Y+DY)>=0))
{
aplace=((X+DX)*127)+Y+DY;
if ((terrain[vg_XY] + river_a[vg_XY]) > terrain[aplace]
&& (terrain[vg_XY] <= terrain[aplace]))
{
lakes[aplace]=terrain[vg_XY] + river_a[vg_XY];
}
}
}
}
}
}
/* The following lines record the minimum and maximum slopes found, a process
* which is only enacted in the first call of the function, since doing so
* each time would slow the evolution process. */
if (isit_first)
{
hold_v += delta_z[vg_XY];
if (delta_z[vg_row+Y] > max_dz)
max_dz = delta_z[vg_XY];
else if (delta_z[vg_XY] < min_dz)
min_dz = delta_z[vg_XY];
}
}
}
/* These torturous lines add the ponding effect. This has a look at the
* entire terrain, and finds local lows. It then checks a ring of positions
* around the local low, and determines whether those spots are above the
* low in elevation, but below the elevation of the water filling the low.
* This is how pools form. */
for (X=0;X<127;X++)
{
for (Y=0;Y<127;Y++)
{
vg_XY=(X*127)+Y;
if (lakes[vg_XY])
{
A=0; B=0; C=0; D=0;
lakec=1;
while (lakec)
{
toggle=0; // Has at least one value changed?
if ((X+A)>0)
{
A--;
toggle++;
}
if ((X+B)<126)
{
B++;
toggle++;
}
if ((Y+C)>0)
{
C--;
toggle++;
}
if ((Y+D)<126)
{
D++;
toggle++;
}
lakec=0;
if (!toggle) break;
for (DX=A; DX<=B; DX++)
{
for (DY=C; DY<=D;DY++)
{
if ((DX==A || DX==B) || (DY==C || DY==D))
{
aplace=(X+DX)*127+(Y+DY);
if (terrain[aplace] < lakes[vg_XY]
&& terrain[aplace] >= terrain[vg_XY]
&& lakes[aplace] != lakes[vg_XY])
{
for (tdx=-1; tdx<=1;tdx++)
{
for (tdy=-1; tdy<=1; tdy++)
{
/* [v1.30 FIX 5] Added bounds check for tdx/tdy offsets */
if ((X+DX+tdx) >= 0 && (X+DX+tdx) < 127 &&
(Y+DY+tdy) >= 0 && (Y+DY+tdy) < 127)
{
if (lakes[(X+DX+tdx)*127+(Y+DY+tdy)] == lakes[vg_XY])
{
lakec++;
lakes[aplace] = lakes[vg_XY];
}
}
}
}
}
}
}
}
}
}
}
}
// These lines write the minimum, maximum, and mean slope values in the delta_z[] array to the screen.
if (isit_first)
{
base = tan((double)(max_dz/10000.0)) * 57.29577952;
sprintf(output_h, "Highest slope %.1f°", (float)base);
text_xy(253, maxx/2, 10, output_h);
hold_v /= 16129;
base = tan((double)(hold_v/10000.0)) * 57.29577952;
sprintf(output_a, "Average slope %.1f°", (float)base);
text_xy(253, maxx/2, 22, output_a);
base = tan((double)(min_dz/10000.0)) * 57.29577952;
sprintf(output_l, "Lowest slope %.1f°", (float)base);
text_xy(253, maxx/2, 35, output_l);
text_xy (0, maxx/2, maxy-30, "Terrain model smoothened");
text_xy (0, maxx/2, maxy-15, "Now generating aspect and slope values");
text_xy (253, maxx/2, maxy-30, "Slope and aspect model generated");
text_xy (253, maxx/2, maxy-15, "press any key to continue");
noise = getch();
text_xy (0, maxx/2, maxy-30, "Slope and aspect model generated");
text_xy (0, maxx/2, maxy-15, "press any key to continue");
text_xy (253, maxx/2, maxy-30, "Now generating watersheds & sediment values");
// text_xy (253, maxx/2, maxy-15, "This takes time - the dot traces progress");
}
/* This code determines the watershed values at each point by searching
* through the entire array for non-zero slopes. For each such value found,
* a subroutine follows the path downhill until a zero slope value is found.
* With every position found, the watershed value is incremented by one, to
* indicate that the first non-zero slope value position (at which the
* subroutine began) is contained within that position's watershed. This is
* accomplished through the determination of DX and DY from the aspect
* direction given in f_model[] for each position. These values, in turn,
* are added to the active X and Y coordinates so that the new, downhill
* (X,Y) are found, and then checked for slope value. In the case of
* non-zero slope, the value in f_model[] for that position is incremented by
* one, and the next downslope position is sought. For zero slope positions
* (local lows), no downslope position is sought.*/
for (X = 0; X < 127; X++)
{
for (Y = 0; Y < 127; Y++)
{
X1 = X;
Y1 = Y;
if (terrain[X1*127+Y1] > 0)
{
do {
vg_row=X1*127;
switch ((f_model[vg_row + Y1] % 10))
{
case 1: DX = -1; DY = -1; break; // These lines recreate the X and Y offset values used in the
case 2: DX = 0; DY = -1; break; // previous function for use below.
case 3: DX = 1; DY = -1; break;
case 4: DX = 1; DY = 0; break;
case 5: DX = 1; DY = 1; break;
case 6: DX = 0; DY = 1; break;
case 7: DX = -1; DY = 1; break;
case 8: DX = -1; DY = 0; break;
case 0: DX = 0; DY = 0; break;
}
// The following line adds one to the watershed value stored in river_a[] to the position downhill of the active position.
if (terrain[((X1+DX)*127)+(Y1+DY)] > 0)
river_a[((X1+DX)*127)+(Y1+DY)]++;
/* The following lines calculate the erosion due to wash and rilling which take
* place at each position. If several positions drain to the same position,
* that position is passed through by the algorithm several times, and the
* wash/rilling effect would be calculated each time, if the position were
* not flagged somehow after the first calculation. The flag used is the
* addition of ten to the f_model[] array for that position. While this has
* no effect on the aspect values held in that array, it allows for the
* determination of a previous calculation of that position's wash/rill
* erosion effects. The equation used is δ = (tan1.4θ * L0.6)/10, where
* θ is slope angle, and L is ground length of slope. */
if (f_model[(X1*127)+Y1] < 10)
{
base = (double) (delta_z[vg_row+Y1]/10000.0);
new_value = (int)(0.5 + (25.11886431 * pow(base, exponent)));
entrain[vg_row+Y1] += new_value;
/* [v1.31 FIX 6] Deposit to actual downhill neighbour
(X1+DX, Y1+DY). Original used entrain[vg_row+Y1+DY]
which only offset in Y, missing the DX row offset. */
entrain[((X1+DX)*127)+(Y1+DY)] -= new_value;
f_model[vg_row+Y1] += 10;
}
/* if (isit_first)
{
if (terrain[vg_row+Y] > 0)
{
//draw(X1,Y1,(int)(24.5 + (terrain[vg_row+Y1]/1000.0)));
draw (X, Y, (int)(50 + 3.5 * (terrain[vg_XY]/1000.0)));
if (lakes[vg_row+Y1])
ldraw(X1,Y1,152);
}
else
draw(X1,Y1,253); // Ocean
}*/
X1 += DX;
Y1 += DY;
/* [v1.30 FIX 3] Reordered while condition: bounds checks
BEFORE array access, and use X1*127+Y1 (not stale vg_row).
Original accessed f_model[vg_row+Y1] where vg_row was based
on the pre-increment X1, and checked arrays before bounds. */
} while ((X1>=0) && (Y1>=0) && (X1<=126) && (Y1<=126) && (f_model[X1*127+Y1]) && (delta_z[X1*127+Y1]));
}
}
}
/* At the end of every running of the values_gen() function, the following
* code draws the new terrain[] model to the screen, including the streams
* which might form upon the surface. No streams are assumed to form below
* a watershed area of 1000m², which yields a flow of 400 cubic metres per
* year if annual rainfall is assumed to approximate that of southern Ontario
* (800mm/yr), and that half of all precipitation is lost to evaporation and
* absorption by plants and soil. Assuming the annual precipitation value of
* 800 mm allows us to exclude rainsplash effects from the calculation as
* vegetation cover would be dense enough to fully shield any loose soil
* (see paper). */
for (X = 0; X < 129; X++)
{
vg_row=X*127;
for (Y = 0; Y < 127; Y++)
{
if (!isit_first && X < 127)
{
vg_XY=vg_row+Y;
if (terrain[vg_XY] > 0)
{
//draw (X, Y, (int)(24.5 + (terrain[vg_XY]/1000.0)));
draw (X, Y, (int)(50 + 3.5 * (terrain[vg_XY]/1000.0)));
if (lakes[vg_XY])
ldraw(X, Y, 152);
}
else
draw (X,Y,253); // Ocean
// if (((f_model[vg_XY] % 10) > 2) && ((f_model[vg_XY] % 10) < 7))
// draw(X, Y, (int)(124.5 + (terrain[vg_XY]/1000.0)));
// else
// draw(X, Y, (int)(74.5 + (terrain[vg_XY]/1000.0)));
}
if (river_a[((X-2)*127)+Y] > 10 && (X > 1) && (terrain[((X-2)*127)+Y] > 0))
{
putpixel(130 + ((X-2)*3), 50 + (Y*3), 152);
switch ((f_model[((X-2)*127) + Y]%10))
{
case 1:
putpixel(129 + (X-2)*3, 49 + Y*3, 152);
putpixel(128 + (X-2)*3, 48 + Y*3, 152);
break;
case 2:
putpixel(130 + (X-2)*3, 49 + Y*3, 152);
putpixel(130 + (X-2)*3, 48 + Y*3, 152);
break;
case 3:
putpixel(131 + (X-2)*3, 49 + Y*3, 152);
putpixel(132 + (X-2)*3, 48 + Y*3, 152);
break;
case 4:
putpixel(131 + (X-2)*3, 50 + Y*3, 152);
putpixel(132 + (X-2)*3, 50 + Y*3, 152);
break;
case 5:
putpixel(131 + (X-2)*3, 51 + Y*3, 152);
putpixel(132 + (X-2)*3, 52 + Y*3, 152);
break;
case 6:
putpixel(130 + (X-2)*3, 51 + Y*3, 152);
putpixel(130 + (X-2)*3, 52 + Y*3, 152);
break;
case 7:
putpixel(129 + (X-2)*3, 51 + Y*3, 152);
putpixel(128 + (X-2)*3, 52 + Y*3, 152);
break;
case 8:
putpixel(129 + (X-2)*3, 50 + Y*3, 152);
putpixel(128 + (X-2)*3, 50 + Y*3, 152);
break;
case 0:
ldraw((X-2), Y, 152);
break;
} // end of switch
}
}
}
if (isit_first)
{
text_xy (0, maxx/2, maxy-30, "Now generating watersheds & sediment values");
// text_xy (0, maxx/2, maxy-15, "This takes time - the dot traces progress");
text_xy (253, maxx/2, maxy-30, "Erosion model generated");
text_xy (253, maxx/2, maxy-15, "press any key to continue");
noise = getch();
text_xy(0, maxx/2, maxy-30, "Erosion model generated");
text_xy(0, maxx/2, maxy-15, "press any key to continue");
text_xy(0, maxx/2, 10, output_h);
text_xy(0, maxx/2, 22, output_a);
text_xy(0, maxx/2, 35, output_l);
text_xy(253, maxx/2, maxy-30, "Terrain evolution in progress");
// text_xy(253, maxx/2, maxy-15, "press any key to halt (then wait)");
text_xy(253, maxx/2, maxy-15, "press any key to halt");
legend(1);
}
return(1);
}
//////////////////////////////////////ERODE()//FUNCTION///////////////////////
/* This function is responsible for figuring the effects of erosion - as
* determined above - upon the model terrain. It also calculates the effect
* of erosion due to the presense of small streams, which form with a least
* watershed value of 1000m². The function runs until a key is pressed, and
* calls the values_gen() function with each iteration. */
// river_a = drainage area, f_model[] = aspect, delta_z[] = slope
void erode()
{
char keep_going = 1, // allows continuation of interrupted simulation
DX, // X offset in array
DY; // Y offset in array
int local_delta, // the amount which the local point changes.
e_row,
e_XY,
local_DX,
local_DY,
local_wouldbe,
local_set[3][3],
lset_sum=0,
local_lrec; // record of locally lowest point
/* [v1.30 FIX 1] eff_ra caps river_a to the original unsigned char range
for use in erosion formulas, preserving intended erosion rates while
allowing accurate watershed tracking with the wider type. */
unsigned char eff_ra;
unsigned int erosion_count = 1;
int clamp_i; /* [v1.30 FIX 2] loop counter for terrain clamping */
gotoxy(1,1);
while (keep_going)
{
while (!kbhit())
{
itoa(erosion_count, output_h, 10);
strcpy(output_a, "Year ");
strcat(output_a, output_h);
text_xy(255, maxx/2, 22, output_a);
lowcount = 31999;
highcount = 0;
for (X = 0; X < 127; X++)
{
e_row=X*127;
for (Y = 0; Y < 127; Y++)
{
e_XY=e_row+Y;
/* These lines calculate the erosion by channel formation, using the equation
* δ = a * W * A * bR * tanθ, where δ is denudation rate in mm/y, a is
* 1/40000mm/(mm*m²), W is watershed size in array units, A is size of array
* unit in m², b is runoff coefficient (0.5), R is runoff (800mm/yr), and θ
* is angle of slope. This equation boils down to δ = W * tanθ, and is
* coded as δ = river_a[(X*127)+Y] * (delta_z[(X*127)+Y]/10000), where
* delta_z[xy] is the change in elevation for the position (x,y), 10000 is
* the ground length of slope in mm, and river_a[xy] is the watershed at
* (x,y).
* Note that the rate of deposition of the material thus flushed out of the
* watershed is a function of the watershed size and the average slope of the
* watershed, which should, statistically speaking, be equal to the average
* slope of the entire model (and therefore is equal to 4°; tan(4°) is
* 0.069926812). This erosion takes place in the position within the
* watershed with a zero watershed; that is, the low to which the river feeds.
* It is therefore possible that considerably less material is being
* deposited to the mouth of system than is being eroded from the watershed,
* if the watershed is particularly steep.
* However, this was the only method of deposition based on channel entrainment
* which could be realistically modelled. */
/* [v1.30 FIX 1] Cap effective river_a to 255 for erosion math.
The original's unsigned char silently overflowed, accidentally
capping erosion rates for large watersheds. This explicit cap
preserves the same intended behavior. */
eff_ra = (river_a[e_XY] > 255) ? 255 : (unsigned char)river_a[e_XY];
// if ((river_a[e_row+Y] >= 10) && (delta_z[e_row+Y]))
if ((eff_ra >= 10) && delta_z[e_XY] && (terrain[e_XY] > 0))
{
entrain[e_XY] += (int) (0.5 + (eff_ra * (delta_z[e_XY])/10000.0));
}
else if (((lakes[e_XY] || eff_ra > 10) && (!delta_z[e_XY])) || terrain[e_XY] < 0) // an end of a river...
{
entrain[e_XY] -= (int) (0.5 + (eff_ra * 0.069926812));
}
// else
// {
// entrain[e_XY] += (int) (0.5 + (delta_z[e_XY])/10000.0);
// }
/* The following lines cause the "erosion" to occur on the terrain model.
* This is done by subtracting the amount of entrained material (in mm/yr)
* occuring at each position (and thus being removed) from the current
* elevation (in mm) for each year. If the amount of entrained material is
* negative - a feat which occurs due to deposition of material entrained
* from higher positions - the effect is to subtract a negative, thus
* increasing the elevation. */
local_delta=entrain[e_XY];
if (local_delta < 0 && (delta_z[e_XY] > 0))
{
lset_sum=0;
local_wouldbe=terrain[e_XY] - local_delta;
for (local_DX=-1;local_DX<=1;local_DX++)
{
for (local_DY=-1;local_DY<=1;local_DY++)
{
/* [v1.30 FIX 4] Reordered: bounds check (SANITY_TEST)
before array access, and corrected the condition
from (local_DX != 0 && local_DX != 0) which checked
DX twice, to (local_DX != 0 || local_DY != 0) which
correctly excludes only the center cell. */
if (SANITY_TEST && (local_DX != 0 || local_DY != 0) && (terrain[((X+local_DX)*127)+(Y+local_DY)] < local_wouldbe))
{
local_set[(local_DX+1)][(local_DY+1)]=1;
lset_sum++;
}
else
local_set[(local_DX+1)][(local_DY+1)]=0;
}
}
if (lset_sum)
{
lset_sum++;
local_delta /= lset_sum;
for (local_DX=-1;local_DX<=1;local_DX++)
{
for (local_DY=-1;local_DY<=1;local_DY++)
{
if ((local_set[(local_DX+1)][(local_DY+1)] || (!local_DX && !local_DY)) && SANITY_TEST)
terrain[((X+local_DX)*127)+(Y+local_DY)]-=local_delta;
}
}
}
else
terrain[e_XY] -= local_delta;
}
else if (local_delta) // erosion
terrain[e_XY] -= local_delta;
if (terrain[e_XY] < lowcount)
lowcount = terrain[e_XY];
else if (terrain[e_XY] > highcount)
highcount = terrain[e_XY];
}
}
/* [v1.30 FIX 2] Clamp terrain to +/-32000 after erosion.
The original's 16-bit int silently overflowed at +/-32767, which
accidentally prevented runaway erosion feedback (steep slope ->
more erosion -> steeper slope -> ...). With corrected arithmetic,
this explicit clamp serves the same stabilizing purpose. */
for (clamp_i = 0; clamp_i < 16129; clamp_i++)
{
if (terrain[clamp_i] > 32000) terrain[clamp_i] = 32000;
else if (terrain[clamp_i] < -32000) terrain[clamp_i] = -32000;
}
values_gen(0);
text_xy(0, maxx/2, 22, output_a);
erosion_count++;
}
noise = getch(); // This line reads in the key-press used to interrupt the program
text_xy(253, maxx/2, 18, "Do you wish to quit the simulation? [y/n]");
noise = getch();
text_xy(0, maxx/2, 18, "Do you wish to quit the simulation? [y/n]");
noise = toupper(noise);
if (noise == 'N') keep_going = 1;
else keep_going = 0;
}
printf ("\r");
text_xy (0, maxx/2, maxy-30, "Terrain evolution in progress");
//text_xy (0, maxx/2, maxy-15, "press any key to halt (then wait)");
text_xy (0, maxx/2, maxy-15, "press any key to halt");
text_xy (253, maxx/2, maxy-30, "Halted by user - this is the final terrain");
text_xy (253, maxx/2, maxy-15, "press any key to exit");
noise = getch();
text_xy (0, maxx/2, maxy-15, "press any key to exit");
if (!(ask)) {
text_xy (253, maxx/2, maxy-30, "Do you wish to save this model to an Idrisi DEM? [y/n]");
noise = getche();
text_xy (0, maxx/2, maxy-30, "Do you wish to save this model to an Idrisi DEM? [y/n]");
noise = toupper(noise);
if (noise == 'Y')
terrain_rec();
}
}
///////////////////////////////////TERRAIN_REC()//FUNCTION///////////////////
/* This function records the terrain model in use to disk, for use with the
* Idrisi GIS as a digital elevation model. The function works through
* opening a file in either text (ASCII) or binary format, as specified by
* the user. This file is then written-to for the 16129 positions in the
* terrain model. The numbers written are the Terrain Evolution Simulator
* vertical units, which are scaled so that each is equal to 1mm in elevation.
* The DEM created may be reimported to the TES with no change, through the
* use of a flag hidden in the .doc file which accompanies the DEM. The flag
* used is the insertion of the "¹" bullet, which is character number 249 in
* the ASCII set, and which looks like a raised period -> ¹ <-. This character
* is written to document file in place of the colon which should proceed the
* data format in the third line of the file. When a DEM is read into the
* TES, each string read before either "ascii" or "binary" is found is
* searched for the ¹ character
*/
void terrain_rec() {
FILE *ofp;
char f_outstring[9],
in_char = '¹',
D_type_string[6],
imgname[13],
docname[13];
int output,
tr_row,
input_pos;
int huge *terr_out;
lowcount = 32000;
highcount = -32000;
fg_locate(1,1);
//gotoxy(1,1);
do {
text_xy(254, maxx/2,15,"Please enter the output file name with no extension.:");
input_pos = 0;
in_char = getch();
while ((in_char != '\r') && (in_char != '\n')) {
f_outstring[input_pos] = in_char;
input_pos++;
in_char = getche();
}
f_outstring[input_pos] = '\0';
strcpy (imgname, f_outstring);
strcat (imgname, ".img");
strcpy (docname, f_outstring);
strcat (docname, ".doc");
text_xy(0, maxx/2,15,"Please enter the output file name with no extension.:");
} while (strchr(f_outstring, '.') != NULL);
do {
text_xy(254, maxx/2,15,"Do you wish to record to ASCII or binary format?[A/B]:");
D_type = getch();
text_xy(0, maxx/2,15,"Do you wish to record to ASCII or binary format?[A/B]:");
D_type = toupper(D_type);
if (!((D_type == 'A') || (D_type == 'B'))) {
beep(1000, 250);
text_xy(253, maxx/2, 25, "ONLY ENTER AN 'A' OR 'B'!");
}
} while (!((D_type == 'A') || (D_type == 'B')));
if (D_type == 'A') {
strcpy (D_type_string, "ascii");
ofp = fopen (imgname, "w"); // open file named tempname in ASCII mode...
for (X = 0 ; X < 127 ; X++) {
tr_row=X*127;
for(Y = 0 ; Y < 127 ; Y++) {
output = terrain[tr_row+Y];
fprintf (ofp, "%d\n", output); // values written to disk.
if (output < lowcount)
lowcount = output;
if (output > highcount)
highcount = output;
}
}
fclose(ofp); // .img file closed.
}
else {
strcpy (D_type_string, "binary");
ofp = fopen (imgname, "wb"); // open file named tempname in ASCII mode
terr_out = (int huge*) farcalloc(16129, 2);
if (terr_out == NULL) {
fg_setmode(old_vmode);
erase();
farfree(f_model);
farfree(terrain);
farfree(delta_z);
farfree(entrain);
farfree(river_a);
printf("error creating terr_out array - insufficient memory!\n");
getch();
abort();
}
for (X = 0 ; X < 127 ; X++) {
tr_row=X*127;
for(Y = 0 ; Y < 127 ; Y++) {
terr_out[tr_row+Y] = terrain[tr_row+Y];
if (terr_out[tr_row+Y] < lowcount)
lowcount = terr_out[tr_row+Y];
if (terr_out[tr_row+Y] > highcount)
highcount = terr_out[tr_row+Y];
}
}
for (output = 0; output < 16129; output++)
fwrite (&terr_out[output], 2, 1, ofp);
farfree(terr_out);
fclose(ofp); // .img file closed.
}
text_xy(0, maxx/2, 25, "ONLY ENTER AN 'A' OR 'B'!");
ofp = fopen(docname, "w"); // .doc file (
fprintf(ofp, "file title : Terrain generated by the Evolution Simulator\n");
fprintf(ofp, "data type : integer\n");
fprintf(ofp, "file type ¹ %s\n", D_type_string); // this line contains the hidden bullet....
fprintf(ofp, "columns : 127\n");
fprintf(ofp, "rows : 127\n");
fprintf(ofp, "ref. system : plane\n");
fprintf(ofp, "ref. units : m\n");
fprintf(ofp, "unit dist. : 10.0000000\n");
fprintf(ofp, "min. X : 0.0000000\n");
fprintf(ofp, "max. X : 1270.0000000\n");
fprintf(ofp, "min. Y : 0.0000000\n");
fprintf(ofp, "max. Y : 1270.0000000\n");
fprintf(ofp, "pos'n error : unknown\n");
fprintf(ofp, "resolution : 10.0\n");
fprintf(ofp, "min. value : %d\n", lowcount); // the lowest value recorded is inserted to the .doc file here
fprintf(ofp, "max. value : %d\n", highcount); // the highest value recorded is inserted to the .doc file here
fprintf(ofp, "value units : unspecified\n");
fprintf(ofp, "value error : unknown\n");
fprintf(ofp, "flag value : none\n");
fprintf(ofp, "flag def'n : none\n");
fclose(ofp);
lowcount = 32000;
highcount = -32000;
}
///////////////////////////////////////////LEGEND()//FUNCTION///////////////////////////////////////////
// The following function draws the legend to the screen (to go along with the terrain model).
void legend(int green_or_grey) {
float cohort;
int highlight=0,
count=0,
shadow=0;
char ostringl[10]; // elevation output string
//if (green_or_grey)
// {
highlight = 50;
shadow = 100;
// }
//for (shadow=0; shadow<254;shadow++)
// {
// set_colour(shadow);
// fg_rect(600+((shadow%2)*10),610+((shadow%2)*10),shadow*2,(shadow+1)*2);
// }
shadow=0;
set_colour(254);
fg_box(0,120,47,431);
fg_box(519,639,47,431);
move_to(0, 78);
fg_draw(120, 78);
fg_justify(1,0);
cohort=highcount/15000.0;
for (count=0;count<15;count++) {
draw (-8, 13+(count*5), (count*3.5)+highlight);
draw (-9, 13+(count*5), (count*3.5)+highlight);
draw (-8, 14+(count*5), (count*3.5)+highlight);
draw (-9, 14+(count*5), (count*3.5)+highlight);
sprintf (ostringl, "%.1f - %.1f", (count*cohort), ((count+1)*cohort));
text_xy(252,95, 91+(15*count), ostringl);
}
sprintf (ostringl, "%d: %.1f", highcount, cohort);
text_xy(252,95, 91+(15*25), ostringl);
fg_justify(0,0);
text_xy(253, 61, 58, "LEGEND");
text_xy(253, 61, 322, "elevation (m)");
set_colour(254);
move_to(0, 336);
fg_draw(120, 336);
//if (green_or_grey) // displays green-shade legend if green being displayed in 256 mode
// {
set_colour(152);
point(102, 362);
point(103, 363);
point(103, 364);
point(103, 365);
point(104, 366);
point(105, 367);
point(106, 368);
point(106, 369);
point(107, 369);
point(108, 370);
point(109, 371);
text_xy(252,61,367,"Stream");
set_colour(152);
draw (-6, 100, 152); draw(-6, 101, 152);
draw (-7, 100, 152); draw(-7, 101, 152);
draw (-8, 100, 152); draw(-8, 101, 152);
draw (-9, 100, 152); draw(-9, 101, 152);
text_xy(252,61,352, "Pond");
set_colour(152);
draw (-6, 110, 253); draw(-6, 111, 253);
draw (-7, 110, 253); draw(-7, 111, 253);
draw (-8, 110, 253); draw(-8, 111, 253);
draw (-9, 110, 253); draw(-9, 111, 253);
text_xy(252,61,382, "Sea");
if (fg_getmode() == 25) {
draw (-6, 112, 125); draw(-6, 111, 125);
draw (-7, 112, 125); draw(-7, 111, 125);
draw (-8, 112, 125); draw(-8, 111, 125);
draw (-9, 112, 125); draw(-9, 111, 125);
text_xy(252,61,385,"South-East ");
draw (-6, 117, 75); draw(-6, 116, 75);
draw (-7, 117, 75); draw(-7, 116, 75);
draw (-8, 117, 75); draw(-8, 116, 75);
draw (-9, 117, 75); draw(-9, 116, 75);
text_xy(252,61,400,"North-West ");
}
// }
draw (-7, 123, 49);
text_xy(252,61,417,"10m x 10m");
}
/////////////////////////////////////////INTRO()//FUNCTION////////////////////////////////////////////////////////////////////////
/* This function describes the workings of the program to the user, after the user runs the terrain evolution simulator with the
command "t_ev_sim :".
*/
void intro() {
textcolor(15);
cprintf("This program models the effects of erosion by surface water on a terrain com-\n\r");
cprintf("posed of unconsolidated material. The terrain model used may be derived in\n\r");
cprintf("one of two ways. The first method is to allow the program to create one using\n\r");
cprintf("a pseudo-fractal method. While this method results in realistic terrains, the\n\r");
cprintf("second method allows for the use of real elevation values, as imported from an\n\r");
cprintf("`Idrisi' Digital Elevation Model (DEM). In either case, the array size is\n\r");
cprintf("limited in size to 127 by 127 positions.\n\r");
cprintf("For the sake of relating the model to ground distances, the array positions\n\r");
cprintf("have been assigned the ground width of 10m on a side, giving an area for each\n\r");
cprintf("array position of 100m². The vertical scale of units has been set to 1 milli-\n\r");
cprintf("metre per unit. These values have been translated to metres for the legends\n\r");
cprintf("which explain the terrain images displayed during the program's course.\n\r");
cprintf("Erosion of the terrain surface is calculated for a number of processes which\n\r");
cprintf("occur under real conditions. These processes are those of wash, rill, and\n\r");
cprintf("channel erosion. Wash and rill erosion occur due to the runoff of the excess\n\r");
cprintf("rainfall which cannot be absorbed by the soil or lost due to evaporation.\n\r");
cprintf("For this model, the amount of runoff was set at 400mm/year, a rate reflecting\n\r");
cprintf("half of the yearly precipitation rate in Niagara. The 50% runoff value was\n\r");
cprintf("chosen as representative of a grass cover. The equation used to calculate\n\r");
cprintf("the combined effects of wash and rill erosion is δ = ((tanθ)^1.4 * L^0.6)/10,\n\r");
cprintf("where δ is removal rate in mm/y/100m² section, θ is angle of slope, and L is\n\r");
printf("\t\t\tpress any key to continue");
noise = getch();
cprintf("\rthe length of slope, which is 10m. \n\r");
cprintf("The second erosion process modelled is channel erosion, which is calculated\n\r");
cprintf("for positions in the terrain surface with a watershed value of 1000m² or\n\r");
cprintf("larger. This value was chosen as to reflect a degree of organization of\n\r");
cprintf("flow into gullies and channels. Unlike the wash/rill erosion, which occurs\n\r");
cprintf("primarily on slopes, this type of erosion occurs to the largest extent on\n\r");
cprintf("the places where flowing water gathers - such as valley and depressions.\n\r");
cprintf("The equation used to calculate this type of erosion is :\n\r");
cprintf("δ = a * W * A * bR * tanθ, where δ is erosion rate in mm/y/100m², a is an\n\r");
cprintf("erosion coefficient equal to 1/40000 mm/mm * m², W is number of array units in\n\r");
cprintf("the watershed, A is the area of each array unit (100m²), b is the runoff co-\n\r");
cprintf("efficient (0.5), R is annual rainfall (800mm/y), and θ is angle of slope.\n\r");
cprintf("The material removed due to erosion must go somewhere. The material removed\n\r");
cprintf("due to rill and wash is delivered during the same year's iteration to the\n\r");
cprintf("position downhill of each eroded position, while the material eroded due to\n\r");
cprintf("channel erosion is delivered to the low point into which the stream eventually\n\r");
cprintf("flows, also during the same year. Modifying both the rates of erosion and\n\r");
cprintf("subsequent deposition are the rules that a) no point may be eroded below the\n\r");
cprintf("elevation of its lowest neighbour, as running water will not appreciably carry\n\r");
cprintf("material uphill out of local lows; and b) deposition will not occur above the\n\r");
cprintf("elevation of the lowest neighbour, as material would be carried elsewhere upon\n\r");
printf("\t\t\tpress any key to continue");
noise = getch();
cprintf("\rsuch an occurrence.\n\r");
cprintf("During the evolution of the terrain, note the generalization of drainage as\n\r");
cprintf("small watersheds combine through stream capture and the sedimentation of local\n\r");
cprintf("declivities. This can result in a terrain with many small watersheds and short\r");
cprintf("streams evolving into one with only a handful of watersheds and streams which\n\r");
cprintf("cross the terrain model. Also note that the formation of ponds and lakes has\n\r");
cprintf("been left out of this simulation to show the effects of the surface. The\n\r");
cprintf("locations of ponds may be estimated by combining neighbouring points of depos-\n\r");
cprintf("ition (marked by a blue square on the terrain image) of similar elevation.\n\n\r");
textcolor(12);
cprintf(" -----USAGE OF THE TERRAIN EVOLUTION SIMULATOR COMMAND LINE-----\r\n\n");
textcolor(15);
cprintf("TO RUN THE TERRAIN SEVOLUTION SIMULATOR WITHOUT FRILLS, TYPE \"");
textcolor(11);
cprintf("t_ev_sim");
textcolor(15);
cprintf("\".\n\rTO SHOW THIS MESSAGE, TYPE \"");
textcolor(11);
cprintf ("t_ev_sim :");
textcolor(15);
cprintf ("\" AT THE PROMPT.\n\r");
cprintf("IF YOU'D LIKE THE OPTION OF EXPORTING TO DEM, TYPE \"");
textcolor(11);
cprintf ("t_ev_sim -o");
textcolor(15);
cprintf ("\" AT THE PROMPT.\n\r");
cprintf ("TO SUPRESS RIVERS DURING TERRAIN DRAWS, TYPE \"");
textcolor(11);
cprintf ("t_ev_sim -r");
textcolor(15);
cprintf ("\" AT THE PROMPT.\n\r");
cprintf ("Note: you may use any arguments together");
printf("\n\t\t\tpress any key to continue");
noise = getch();
clrscr();
}
////////////////////////////////////GET_TERRAIN()//FUNCTION/////////////////////////////////////////////////////////////////////
/* This function reads the terrain() array from a user-designated Idrisi DEM file. It works through prompting for DEM file
name (which may be entered with or without .img extension), and ensuring that the file exists. It then opens the .img file, and
reads in all of the values, assigning them to positions within the terrain[] array. THIS ONLY WORKS WITH 127 x 127 DEMs! Future
versions will hopefully be able to incorporate larger DEM files, and take into account different horizontal and vertical scales.
For now, all DEMs are assumed to be of 10m x 10m ground scale per position. Elevation values within the DEM are scaled over the
whole range of possible TES elevation values (approximately -25000 to +25000). This is accomplished through calculating the range
of elevation values in the DEM, and multiplying each elevation value in the DEM by a scaling value which extends this range over
the TES elevation range.
*/
int get_terrain() {
char was_gen_TES=0; // This varibale stores a flag determining whether the DEM on disk was generated by the TES!
int v_read,
gt_row,
gt_XY,
fileL = 32000, // This integer is used to hold the lowest value read in from the file
fileH = -32000; // This integer is used to hold the highest value read in from the file
float d_elev, // This variable is used to hold the range of elevation values read from the file
scale_elev; // This variable is used to hold the scaling value used in deriving TES elevation values
FILE *fp; // This is data type required for pointing to open files
int huge *in_terr;
/* The following do-while loop asks for a valid input file name until one is succesfully found. */
do {
cprintf ("\n\r\n\rPlease specify the DEM file name with no file extension\n\r");
cprintf (" :___________________________\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
scanf("%s%*c", &output_h); // This use of output_h is to hold the input file name
if (!strcmp(output_h, "QUIT") || !strcmp(output_h, "quit"))
abort();
strcpy(output_a, output_h);
if (!(strchr(output_a, '.'))) // this line checks for an extension in the file name entered
strcat(output_a, ".doc"); // if no extension was entered, the .doc extension is added to the filename
fp = fopen(output_a, "r");
if (fp == NULL)
printf ("\n\tThere is no such file! Type QUIT to abort program.\n");
else {
do {
fscanf (fp, "%s", &output_l);
if (strchr(output_l,'¹'))
was_gen_TES = 1;
} while (strcmp(output_l, "binary") && strcmp(output_l, "ascii"));
if (strchr(output_l, 'b') == NULL) D_type = 'A';
else D_type = 'B';
fclose(fp);
strcpy(output_a, output_h);
if (!(strchr(output_a, '.')))
strcat(output_a, ".img");
if (D_type == 'A')
fp = fopen(output_a, "r");
else
fp = fopen(output_a, "rb");
}
} while (fp == NULL);
setup_graphics();
/* The following lines read in the values from the specified file on disk. These are then converted to TES elevation values and
stored in the terrain[] array according to the data model indicated above. */
if (D_type == 'A')
{
for (X = 0; X < 127; X++) {
gt_row=X*127;
for (Y = 0; Y < 127; Y++) {
gt_XY=gt_row+Y;
fscanf (fp, "%d", &terrain[gt_XY]);
if (terrain[gt_XY] < fileL)
fileL = terrain[gt_XY];
if (terrain[gt_XY] > fileH)
fileH = terrain[gt_XY];
}
}
}
else {
in_terr = (int huge*) farcalloc(16129, 2);
if (in_terr == NULL) {
printf("error creating in_terr array - insufficient memory!\n");
return(0);
}
fread (in_terr, 2, 16129, fp);
for (X = 0; X < 127; X++) {
gt_row=X*127;
for (Y = 0; Y < 127; Y++) {
gt_XY=gt_row+Y;
terrain[gt_XY] = in_terr[gt_XY];
if (terrain[gt_XY] < fileL)
fileL = terrain[gt_XY];
if (terrain[gt_XY] > fileH)
fileH = terrain[gt_XY];
}
}
farfree(in_terr);
}
// The range of elevation values within the DEM is calculated here.
d_elev = fileH - fileL;
scale_elev = 48000.0 / d_elev; // The scaling value is calculated by dividing 48000 (just less than the full 50000) by the range
// of DEM elevations.
// The following lines convert the values read into the terrain[] array to TES values.
for (X = 0; X< 127; X++) {
gt_row=X*127;
for (Y = 0; Y < 127; Y++) {
gt_XY=gt_row+Y;
if (!was_gen_TES)
terrain[gt_XY] = -24990 + (scale_elev * (terrain[gt_XY] - fileL));
draw(X,Y,(int)((terrain[gt_XY]/1000.0) + 24.5));
if (terrain[gt_XY] < lowcount)
lowcount = terrain[gt_XY];
else if (terrain[gt_XY] > highcount)
highcount = terrain[gt_XY];
avecount += terrain[gt_XY];
}
}
text_xy (253, maxx/2, maxy-30, "DEM read");
text_xy (253, maxx/2, maxy-15, "press any key to continue");
legend(0); // This call to the legend() function draws the legend to the screen
// The following lines output the highest, average, and lowest elevation values.
hold_v = ((highcount/1000.0) + 25.0);
sprintf (output_h, "Highest elevation %.1f m", hold_v);
text_xy (253, maxx/2, 10, output_h);
hold_v = (avecount / 16129000.0) + 25.0;
sprintf (output_a, "Average elevation %.1f m", hold_v);
text_xy (253, maxx/2, 22, output_a);
hold_v = ((lowcount/1000.0) + 25.0);
sprintf (output_l, "Lowest elevation %.1f m", hold_v);
text_xy (253, maxx/2, 35, output_l);
noise = getch();
text_xy (0, maxx/2, 15, output_h);
text_xy (0, maxx/2, 25, output_a);
text_xy (0, maxx/2, 35, output_l);
text_xy (0, maxx/2, maxy-30, "DEM read");
text_xy (0, maxx/2, maxy-15, "press any key to continue");
text_xy (253, maxx/2, maxy-15, "Now generating aspect and slope values");
return(1);
}
/////////////////////////////////////DRAW()//FUNCTION/////////////////////////
/* This simple function draws a three pixel by three pixel block to the screen
* centred at the point determined by the arguments X_pos, Y_pos. It was
* made a function because it is used so often. */
void draw(int X_pos, int Y_pos, int col) {
int wha;
wha=fg_getcolor();
set_colour(col);
draw_rect (129 + X_pos*3, 131 + X_pos*3, 49 + Y_pos*3, 51+Y_pos*3);
set_colour(wha);
}
void ldraw (int X_pos, int Y_pos, int col) {
//setrgb(152,0,0,63);
col=152;
putpixel(130 + X_pos*3, 50 + Y_pos*3, col);
putpixel(129 + X_pos*3, 50 + Y_pos*3, col);
putpixel(131 + X_pos*3, 50 + Y_pos*3, col);
putpixel(130 + X_pos*3, 49 + Y_pos*3, col);
putpixel(130 + X_pos*3, 51 + Y_pos*3, col);
if (lakes[(X_pos-1)*127+Y_pos] || river_a[(X_pos-1)*127+Y_pos] > 10) {
putpixel(129 + X_pos*3, 49 + Y_pos*3, col);
putpixel(129 + X_pos*3, 51 + Y_pos*3, col);
}
if (lakes[X_pos*127+(Y_pos-1)] || river_a[X_pos*127+(Y_pos-1)] > 10) {
putpixel(129 + X_pos*3, 49 + Y_pos*3, col);
putpixel(131 + X_pos*3, 49 + Y_pos*3, col);
}
if (lakes[(X_pos+1)*127+Y_pos] || river_a[(X_pos+1)*127+Y_pos] > 10) {
putpixel(131 + X_pos*3, 49 + Y_pos*3, col);
putpixel(131 + X_pos*3, 51 + Y_pos*3, col);
}
if (lakes[X_pos*127+(Y_pos+1)] || river_a[X_pos*127+(Y_pos+1)] > 10) {
putpixel(129 + X_pos*3, 51 + Y_pos*3, col);
putpixel(131 + X_pos*3, 51 + Y_pos*3, col);
}
}
/* 4p27 array memory model (insufficient)
4p27_dma dynamically-allocated memory model
tmodel models surface & erosion
qmodel attempts at speeding up simulator
kmodel ditto
teevsim Terrain Evolution Simulator t_ev_sim.exe v 1.06
t_ev_sim Terrain Evolution Simulator t_ev_sim.exe v 1.07
tes1_08 Terrain Evolution Simulator t_ev_sim.exe v 1.08
tes1_09 " v 1.09 (larger z range, attempt at IBM8514 colour)
tes1_10 " v 1.10 (again, larger z range)
tes1_11 " v 1.11 (no reclassed DEMs - plateau effect,
fixed colours,
elev output in m,
fixed flipping of DEMs,
added legend for import
added stream erosion)
tes1_12 " v 1.12 (fixed vertical scale to 1 unit ≡ 1mm,
fiddled with read-in output, etc.,
stream drawing,
documentation)
t_ev_sim <-- end of 4p27...
tes1_14 TES v 1.13 (interruptable erosion simulation,
removal of multiple variable,
reads from binary)
tes1_15 reads, writes to binary, ascii
tes1_16a last of the grey types.
tes1_16b failed reattempt at IBM8514 colour....
tes1_16c green landscape, minimum overdeposition of 200mm allowed
tes1_17 fills potholes at depositional sites to lowest neighbour....
tes1_2 256 colours, legend(),
tes1_21.cpp fastgraf 16/256 colour capable
Reworked in 98/02 for improved river generation
Also attempts at lake formation.
tes_130.c v 1.30 corrected 2026: overflow fixes, bounds-check
fixes, deposition condition fix (see header).
tes_131.c v 1.31 corrected 2026: wash/rill deposition direction
fix, lake detection restricted to true depressions.
*/