asthmatic_thematic

//================================================================================

// RELEASED UNDER GNU LICENSE
//   see www.gnu.org/licenses/gpl.html
// asthmatic_thematic

//================================================================================
// worldgen.c ====================================================================
//================================================================================

//include <cstdlib>
#include <time.h>
//include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include "math.h"

//================================================================================

// compiling and running
//   gcc worldgen.c ; ./a.out
// plotting data: use eg gnuplot
//   cd to correct directory:
//     cd 'dir'
//   set views:
//     set pm3d map
//     set zrange [0:10000]
//   map data
//     splot 't.txt' matrix
//     (or e.txt, or w.txt, ...)
// selecting types of data
//   either use gnuplot data ranges,
//   or go to "output" section below:
//     choose a structure to examine,
//     and set it to a high value
//       (for reproducibility, set
//       time now defn to a constant)
//

// this should be run at greater resolution: say 1028x1028
//   this would allow smoothing for more realistic thin rivers & lakes
//     help making lakes: topographic fill
//   tome coast would be smoother
// converting to tome format should be easy from this high resolution
//   use rectangles appropriately scaled to get 170x58 tome size
//   if shallow water in area, make whole tile shallow; same for towns etc
//   otherwise, use most common terrain as tile type

//================================================================================

#define TIMENOW (unsigned)time(0)

#define SIZE 512
//256
#define PIE 3.1415927

// RANDOM NUMBER FUNCTIONS
#define R_SZ (int)(rndmf(SIZE))
#define R_ONE rndmf(1.0)
#define R_1F1 (1.0-2.0*rndmf(1.0))
#define R_1I1 (1-(int)(rndmf(2.999999999999)))

// LENGTH SCALES
#define H_SCALE 2.0
// roughness goes up as scale below goes down
#define SCALE 0.55
#define DOWNSLOPE 0.005
#define DOWNVAR 1.0
#define WET_SCALE 1.8
#define BIAS_STRENGTH 0.0

// MOUNTAINS
#define N_MTN_RANGES 7
#define MTN_JIGGLE 0.7
#define MTN_WANDER 4.0
#define MTN_JUMP 0.07
#define MTN_WIDTH 0.04
#define MTN_ADD .006

// RIVERS AND SPRINGS
#define RIV_STEP 3
#define RIV_WIDTH 3
#define SM_WIDTH 5
#define N_SPRINGS 3
#define H_SPRING_MIN 0.4
#define H_SPRING_MAX 0.9
#define N_ANTI_SPRINGS 6
#define H_ANTI_MIN -0.3
#define H_ANTI_MAX -0.001
#define RIV_MAXL SIZE/4
#define LAKE_REPEAT 10
#define LAKE_WALKS 10

// FORESTS
#define N_FORESTS 20
#define N_FOR_JUMPS 30
#define N_FOR_WALKS 40
#define N_FOR_FAR 10
#define FOREST_WIDTH 2
#define H_FOREST_MIN 0.05
#define H_FOREST_MAX 0.7

// DESERTS
#define N_DESERTS 5
#define DESERT_WIDTH 8
#define N_DES_JUMPS 10
#define N_DES_WALKS 20
#define N_DES_FAR 4
#define SCORCH_J 0.1

// TUNDRA, GLACIERS
#define TUNDRA_J 0.8
#define N_GLACIER_WALKS 10
#define GLACIER_WIDTH 8

// TOPOGRAPHY
#define N_CONE 9
#define H_CONE 1.8
#define SCALE_CONE 0.01
#define H_HILLS 0.2
#define H_MTN   0.5
#define H_PEAK  0.9

// TERRAINS
//   including shallow and deep water
#define PLAIN     1
#define DESERT    2
#define WOODS     3
#define HILLS     4
#define TUNDRA    5
#define SHALLOW   6

#define SCORCH    12
#define FOREST    13
#define MOUNTAIN  14
#define GLACIER   15
#define DEEP      16

#define PEAK      24

#define DIFFIMPASS 10

// TOWNS AND SITES

#define T_START   50
#define T_HUMAN   51
#define T_SANDM   52
#define T_ELF     53
#define T_DWARF   54
#define T_SNOWM   55
#define T_RIVM    56

#define T_BURNM   62
#define T_TREEM   63
#define T_MTNM    64
#define T_ICEM    65
#define T_OCEANM  66

/*
later: ease of moving for races

hum  1  2   3  4  5  6 12  13  14  15
elf  3  1  13  2  4  5  6  12  14  15
dwf  4  1  14  2  3  5  6  12  13  15
sno  5 15  14 ...
*/

void fill1DFractArray (float *, int, int, float, float);
void fill2DFractArray (float *, int, int, float, float);

#include "dmnd_square.c"

int main() {


  float * ffaa;
  int hsize=SIZE/2;
  float h1 = 1.0;
  float h2 = 0.4;
  float h_main = 0.5;
  float h_polar = 1.5;

  float dist, angle;
  int i, j;
  int t;
  int i1,j1,i2,j2;
  float divs;
  int ii,jj;
  int k, kk;
  int irndrad, jrndrad;

  float elev[SIZE][SIZE];
  int terr[SIZE][SIZE];
  float wet[SIZE][SIZE];
  unsigned nowtime;

  nowtime = TIMENOW;
  srand(nowtime);

  //================================================================
  // set basic topography
  ffaa = alloc2DFractArray(SIZE);
  fill2DFractArray (ffaa, SIZE, TIMENOW, H_SCALE, SCALE);

  for (i=0; i<SIZE; i++) {
    for (j=0; j<SIZE; j++) {
      elev[i][j] = ffaa[(i*SIZE)+j] + h_main;
    }
  }

  //================================================================
  // set fault/tectonic bias a la continents
  float bias_angle;
  bias_angle = R_ONE*2*PIE;

  //================================================================
  // scale down radially from center
  int igeo, jgeo;
  float set_dist, set_slop;
  float * fdist, * fslop;
  float angle_norm, angle_rem;
  float leave_pole;
  float bias_near;

  fdist = alloc1DFractArray(16); //17 //distance to geo center
  fslop = alloc1DFractArray(16); //17 //slope to use

  for (k=0; k<2; k++) {
    fill1DFractArray(fdist, 16, TIMENOW, SIZE/3.0, .1);
    fill1DFractArray(fslop, 16, TIMENOW, DOWNSLOPE, .2);

    igeo = R_SZ/3. + SIZE/3.;
    jgeo = R_SZ/3. + SIZE/3.;

    for (i=0; i<SIZE; i++) {
      for (j=0; j<SIZE; j++) {
        dist = distance(i,igeo,j,jgeo);
        angle = atan2(i-igeo, j-jgeo) + PIE;
        angle_norm = angle * (16./(2.*PIE));

        // interpolate
        angle_rem = angle_norm-(int)angle_norm;
        set_dist = (1-angle_rem)*fdist[(int)angle_norm] +
          angle_rem*fdist[(int)angle_norm+1];
        set_dist += SIZE/3.;

        // fault bias...
        bias_near = 1 + BIAS_STRENGTH*cos(angle-bias_angle)*(1+0.1*R_1F1);

        set_slop =
          -DOWNSLOPE * (1+DOWNVAR*fslop[(int)(angle*16./360.)]);

        // ARCTIC: adjust scale down to get more land up north
        if (j>3.*SIZE/4.) set_slop *= (1 - h_polar*0.5*4.*(j-3.*SIZE/4.)/SIZE);

        if (dist > set_dist) elev[i][j] += set_slop*(dist-set_dist);
      }
    }
  }

  //================================================================
  // raise and lower topographic random cones
  float loci, locj, cone_height, cone_rad;

  for (k=0; k<N_CONE; k++) {
    loci = R_SZ;
    locj = R_SZ;
    cone_height = H_CONE    * (1.25 - 2*R_ONE);
    cone_rad =    SIZE/6.0  * (1.5 -    R_ONE);
    fill1DFractArray(fdist, 16, TIMENOW, cone_rad, SCALE_CONE);
    for (i=0; i<SIZE; i++) {
      for (j=0; j<SIZE; j++) {
        //  if (distance(loci,i,locj,j) < radcone) elev[i][j] -= h_cone;
        dist = distance(i,loci,j,locj);
        angle = 360/(2*PIE)*atan2(i-loci, j-locj)+180.;
        angle_norm = angle*16./360.;
        angle_rem = angle_norm-(int)angle_norm;

        // interpolate
        set_dist = (1-angle_rem)*fdist[(int)angle_norm] +
          angle_rem*fdist[(int)angle_norm+1];
        set_dist += cone_rad;

        if (dist < set_dist)
          elev[i][j] -=
            cone_height * pow((set_dist-dist)/set_dist,2);
      }
    }
  }

  //================================================================
  // could raise arctic arctic a bit more
  //  for (i=0; i<SIZE; i++)
  //    for (j=3.*SIZE/4.; j<SIZE; j++)
  //      elev[i][j] += h_polar * R_ONE * 4.*(j-3.*SIZE/4.)/SIZE;

  //================================================================
  // add mountain ranges
  float * fi;
  float * fj;
  float mtn_wid;
  int kkk, l1, l2;
  int mtn_width;
  int n, kkkmax;
  fi = ((float *) malloc (sizeof(float) * 7)); //alloc1DFractArray(6); //7
  fj = ((float *) malloc (sizeof(float) * 7)); //alloc1DFractArray(6); //7

  for (n=0; n<N_MTN_RANGES; n++) {
  i1 = R_SZ;
  j1 = R_SZ;
  i2 = R_SZ;
  j2 = R_SZ;

  segment(fi, i1, i2, 6); //fi[6] = i2;
  segment(fj, j1, j2, 6); //fj[6] = j2;

  for (kkk=0; kkk<6; kkk++) {
    divs = distance(fi[kkk],fi[kkk+1],fj[kkk],fj[kkk+1]);
    i1 = fi[kkk]   + R_1I1 * MTN_JIGGLE * (R_SZ/80. + 1);
    i2 = fi[kkk+1] + R_1I1 * MTN_JIGGLE * (R_SZ/80. + 1);
    j1 = fj[kkk]   + R_1I1 * MTN_JIGGLE * (R_SZ/80. + 1);
    j2 = fj[kkk+1] + R_1I1 * MTN_JIGGLE * (R_SZ/80. + 1);

    divs = distance(i1,i2,j1,j2);
    for (k=0; k <= divs; k++) {
      ii = k*(i2-i1)/divs + i1;
      jj = k*(j2-j1)/divs + j1;

      for (kk=0; kk<MTN_WANDER*SIZE/15.; kk++) {
        irndrad = MTN_JUMP  * 3*R_1F1 * SIZE/20.;
        ii+=irndrad;
        jrndrad = MTN_JUMP  * 3*R_1F1 * SIZE/20.;
        jj+=jrndrad;
        mtn_wid = MTN_WIDTH * 5*R_1F1 * SIZE/10.;

        for (l1=-mtn_wid; l1<mtn_wid; l1++) {
          for (l2=-mtn_wid; l2<mtn_wid; l2++) {
            if ((ii+l1)>=0 && (ii+l1)<SIZE && (jj+l2)>=0 && (jj+l2)<SIZE) {
                elev[ii+l1][jj+l2] += MTN_ADD * R_ONE;
            }
          }
        }
      }
    }
  }
  }

  //================================================================
  // raise map so there is 70% land: quick and dirty version
  //   (should use quicksort; choose 70% item; add constant so it's zero)
  int num_above = 0;
  float correction = 0.0;
  float h_corr = 0.3;

  for (i=0; i<SIZE; i++)
    for (j=0; j<SIZE; j++)
      if (elev[i][j] > 0) num_above++;
  if (num_above > 0.85*SIZE*SIZE) correction=-h_corr;
  if (num_above < 0.60*SIZE*SIZE) correction=+h_corr;
  if (correction != 0.0)
    for (i=0; i<SIZE; i++)
      for (j=0; j<SIZE; j++)
        elev[i][j] += correction;

  //================================================================
  // ERODE SHORELINES (to make a bit smoother edges)

  //================================================================
  // SMALL DECREASE TO EDGES (to minimize need for borders, below)

  //================================================================
  // SCALE ABOVEGROUND TO ONE
  float maxval = -100.0;
  for (i=0; i<SIZE; i++) for (j=0; j<SIZE; j++)
    if (elev[i][j]>maxval) maxval=elev[i][j];
  for (i=0; i<SIZE; i++) for (j=0; j<SIZE; j++)
    if (elev[i][j]>0) elev[i][j] /= maxval;

  //================================================================================
  // END OF TOPOGRAPHY==============================================================
  // START TERRAIN==================================================================
  //================================================================================

  //================================================================
  // FORM TERRAIN MAP

  // convert to initial terrain
  for (i=0; i<SIZE; i++) {
    for (j=0; j<SIZE; j++) {
      t = DEEP;
      if (elev[i][j] > 0 && elev[i][j] < H_HILLS) t = PLAIN; //plain
      if (elev[i][j] > H_HILLS && elev[i][j] < H_MTN) t = HILLS; //hill
      if (elev[i][j] > H_MTN && elev[i][j] < H_PEAK) t = MOUNTAIN;  //mtn
      if (elev[i][j] > H_PEAK) t = PEAK; // peak
      terr[i][j] = t;
    }
  }

  //================================================================
  // ADD SPRINGS

  int notmid;
  int spr_x[N_SPRINGS];
  int spr_y[N_SPRINGS];
  int x, y;
  int antispr_x[N_ANTI_SPRINGS];
  int antispr_y[N_ANTI_SPRINGS];

  // i guess this could loop infinitely in tiny maps: careful
  for (i=0; i<N_SPRINGS; i++) {
    notmid = 1;
    while (notmid) {
      x = R_SZ;
      y = R_SZ;
      if (elev[x][y]>H_SPRING_MIN && elev[x][y]<H_SPRING_MAX && y>SIZE/5.) notmid=0;
    }
    spr_x[i] = x;
    spr_y[i] = y;
  }
  for (i=0; i<N_ANTI_SPRINGS; i++) {
    notmid = 1;
    while (notmid) {
      x = R_SZ;
      y = R_SZ;
      if (elev[x][y]>H_ANTI_MIN && elev[x][y]<H_ANTI_MAX && y>SIZE/5.) notmid=0;
    }
    antispr_x[i] = x;
    antispr_y[i] = y;
  }

  //================================================================
  // RUN RIVERS & LAKES
  // could use nelder mead or something better here

  float minval = 1000000.0;
  int mindex=0;
  int notsea, notlake, nsteps;
  // these are dangerous definitions: floats vs ints
  float spring_elev, lastelev, angle_down, new_x, new_y, xx, yy;
  float i_grad, j_grad, old_i, old_j;
  int orig_x, orig_y;

  for (kkk=0; kkk<N_SPRINGS; kkk++) {
    xx = spr_x[kkk]; yy = spr_y[kkk];
    notsea = notlake = 1;
    nsteps = 0;
    while (1) {
      nsteps++;

      i_grad = 0.0;
      for (i=-SM_WIDTH; i<=SM_WIDTH; i++)
        for (j=-SM_WIDTH; j<=SM_WIDTH; j++)
          if ((xx+i)<SIZE && (xx+i)>=0 && (yy+j)<SIZE && (yy+j)>=0)
            i_grad += i*elev[(int)xx+i][(int)yy+j];
      j_grad = 0.0;
      for (i=-SM_WIDTH; i<=SM_WIDTH; i++)
        for (j=-SM_WIDTH; j<=SM_WIDTH; j++)
          if ((xx+i)<SIZE && (xx+i)>=0 && (yy+j)<SIZE && (yy+j)>=0)
            j_grad += j*elev[(int)xx+i][(int)yy+j];
      angle = atan2(j_grad,i_grad);

      new_x = xx - cos(angle)*RIV_STEP;
      new_y = yy - sin(angle)*RIV_STEP;

      if ((new_x>=SIZE-1) || (new_x<1) || (new_y>=SIZE-1) || (new_y<1))
        {printf(" -escape- "); break;}

      // fill in river
      divs = distance(xx,new_x,yy,new_y);
      for (k=0; k<divs; k++) {
        ii = k*(new_x-xx)/divs + xx;
        jj = k*(new_y-yy)/divs + yy;
        for (i=0; i<RIV_WIDTH; i++) for (j=0; j<RIV_WIDTH; j++) terr[ii+i][jj+j] = SHALLOW;
      }
      xx = new_x; yy = new_y;

      if (nsteps > 100) {printf(" -deadlooplake- "); notlake=0; break;}

      if (distance(spr_x[kkk],xx,spr_y[kkk],yy)>2*RIV_MAXL) {printf(" -long- "); break;}

      // if ocean, quit
      if (elev[(int)new_x][(int)new_y] < 0) {notsea=0; printf(" -ocean- "); break;}
    }
    // lake fill (very very simple right now)
    // future improvements:
    // do we want to avoid lakes right near shore?
    // any easy way to fill topographically?
    //   take endpoint x,y of river
    //   grab radius of 5 around;
    //     if min_elev is not on border, fill from min_elev to half min-max
    //   connect river endpoint x,y to min_elev to make sure river hits lake
    // flow from lakes:
    //   half the time,
    //     use segments to jump to random near point which is lower in elev
    //   restart river from there
    // bigger issue: dont apply shallow immed; make linked list;
    //   if river doesn't make it far, dont bother changing to shallow
    //   if succeeds, apply linked list to make shallow
    // also make random lakes
    if (notlake==0 && elev[(int)new_x][(int)new_y]>0.15) {
      orig_x = xx; orig_y = yy;
      for (i1=0; i1<LAKE_REPEAT; i1++) {
        i = orig_x; j = orig_y;
        for (i2=0; i2<LAKE_WALKS; i2++) {
          terr[i][j] = terr[i+1][j] = terr[i][j+1] = terr[i+1][j+1] = SHALLOW;
          i += R_1I1; j += R_1I1;
        }
      }
    }
  }

  //anti-springs
  for (kkk=0; kkk<N_ANTI_SPRINGS; kkk++) {
    xx = antispr_x[kkk]; yy = antispr_y[kkk];
    notsea = notlake = 1;
    nsteps = 0;
    while (1) {
      nsteps++;

      i_grad = 0.0;
      for (i=-SM_WIDTH; i<=SM_WIDTH; i++)
        for (j=-SM_WIDTH; j<=SM_WIDTH; j++)
          if ((xx+i)<SIZE && (xx+i)>=0 && (yy+j)<SIZE && (yy+j)>=0)
            i_grad += i*elev[(int)xx+i][(int)yy+j];
      j_grad = 0.0;
      for (i=-SM_WIDTH; i<=SM_WIDTH; i++)
        for (j=-SM_WIDTH; j<=SM_WIDTH; j++)
          if ((xx+i)<SIZE && (xx+i)>=0 && (yy+j)<SIZE && (yy+j)>=0)
            j_grad += j*elev[(int)xx+i][(int)yy+j];
      angle = atan2(j_grad,i_grad);

      new_x = xx + cos(angle)*RIV_STEP;
      new_y = yy + sin(angle)*RIV_STEP;

      if ((new_x>=SIZE-1) || (new_x<1) || (new_y>=SIZE-1) || (new_y<1))
        {printf(" -escape- "); break;}

      // fill in river
      divs = distance(xx,new_x,yy,new_y);
      for (k=0; k<divs; k++) {
        ii = k*(new_x-xx)/divs + xx;
        jj = k*(new_y-yy)/divs + yy;
        if (elev[ii][jj]>0)
          for (i=0; i<RIV_WIDTH; i++) for (j=0; j<RIV_WIDTH; j++) terr[ii+i][jj+j] = SHALLOW;
      }
      xx = new_x; yy = new_y;

      if (nsteps > 100) {printf(" -deadlooplake- "); notlake=0; break;}
      if (distance(antispr_x[kkk],xx,antispr_y[kkk],yy)>RIV_MAXL) {printf(" -long- "); break;}

    }
    if (notlake==0 && elev[(int)new_x][(int)new_y]>0.15) {
      orig_x = xx; orig_y = yy;
      for (i1=0; i1<LAKE_REPEAT; i1++) {
        i = orig_x; j = orig_y;
        for (i2=0; i2<LAKE_WALKS; i2++) {
          terr[i][j] = terr[i+1][j] = terr[i][j+1] = terr[i+1][j+1] = SHALLOW;
          i += R_1I1; j += R_1I1;
        }
      }
    }
  }

  //================================================================
  // DETERMINE RANDOM WETNESS MAP
  fill2DFractArray (ffaa, SIZE, TIMENOW, 1, WET_SCALE);

  minval =  100000.0;
  maxval = -100000.0;
  for (i=0; i<SIZE; i++) {
    for (j=0; j<SIZE; j++) {
      wet[i][j] = ffaa[(i*SIZE)+j];
      if (wet[i][j]>maxval) maxval=wet[i][j];
      if (wet[i][j]<minval) minval=wet[i][j];
    }
  }

  // normalize to 0-1
  for (i=0; i<SIZE; i++) {
    for (j=0; j<SIZE; j++) {
      wet[i][j] -= minval;
      wet[i][j] /= (maxval-minval);
    }
  }

  // make river areas humid
  for (i=0; i<SIZE; i++)
    for (j=0; j<SIZE; j++)
      if (terr[i][j] == SHALLOW)
        for (ii=i-8; ii<i+8; ii++)
          for (jj=j-8; jj<j+8; jj++) wet[ii][jj] += 0.1;

  // make south more arid
  for (i=0; i<SIZE; i++)
    for (j=0; j<SIZE/4.; j++)
      wet[i][j] *= j/(SIZE/4.);

  //================================================================
  // RUN FORESTS

  int for_x[N_FORESTS], for_y[N_FORESTS];

  // locations of forest seeds
  for (i=0; i<N_FORESTS; i++) {
    notmid = 1;
    while (notmid) {
      x = R_SZ;
      y = R_SZ;
      //      printf("(%i %i)\n", x, y);
      if (elev[x][y]>H_FOREST_MIN && elev[x][y]<H_FOREST_MAX && wet[x][y]>0.6)
        notmid=0;
    }
    for_x[i] = x;
    for_y[i] = y;
    //    printf("F(%i %i);[%.2f]. ",x,y,wet[x][y]);
  }

  int nfj;
  //multiple random walks from same point; jump to new point, repeat
  //check for success by humidity
  for (kkk=0; kkk<N_FORESTS; kkk++) {
    x = for_x[kkk];
    y = for_y[kkk];
    nfj = N_FOR_JUMPS * R_ONE;
    for (i=0; i<nfj; i++) {
      for (j=0; j<N_FOR_WALKS; j++) {
        for (ii=-FOREST_WIDTH; ii<=FOREST_WIDTH; ii++) {
          for (jj=-FOREST_WIDTH; jj<=FOREST_WIDTH; jj++)
            if (terr[x+ii][y+jj] <= TUNDRA && wet[x+ii][y+jj]>(R_ONE+0.1)
                && elev[x+ii][y+jj]>0 && elev[x+ii][y+jj]<H_MTN)
              terr[x+ii][y+jj] = WOODS;
        }
        x += R_1I1; y += R_1I1;
      }
      // jump to new point; forests spread along latitudes
      x += 1.5 * R_1I1 * R_ONE * N_FOR_FAR;
      y += .75 * R_1I1 * R_ONE * N_FOR_FAR;
      // more forests if humid
      if (wet[x][y]>0.75 && R_ONE>0.1) i--;
    }
  }

  //================================================================
  // PLACE DESERTS

  int des_x[N_DESERTS], des_y[N_DESERTS];

  // desert seeds
  for (i=0; i<N_DESERTS; i++) {
    notmid = 1;
    while (notmid) {
      x = R_SZ;
      y = R_SZ;
      //      printf("(%i %i)\n", x, y);
      if (wet[x][y]<0.5 && elev[x][y]>0 && y/SIZE<0.9) notmid=0;
    }
    des_x[i] = x;
    des_y[i] = y;
    //    printf("D(%i %i);[%.2f]. ",x,y,wet[x][y]);
  }

  //multiple random walks from same point; jump to new point, repeat
  for (kkk=0; kkk<N_DESERTS; kkk++) {
    x = des_x[kkk];
    y = des_y[kkk];
    nfj = N_DES_JUMPS * R_ONE;
    for (i=0; i<nfj; i++) {
      for (j=0; j<N_DES_WALKS; j++) {
        for (ii=-DESERT_WIDTH; ii<=DESERT_WIDTH; ii++)
          for (jj=-DESERT_WIDTH; jj<=DESERT_WIDTH; jj++)
            if (terr[x+ii][y+jj] <= TUNDRA &&
                wet[x+ii][y+jj]<.5 || wet[x+ii][y+jj]<(R_ONE*R_ONE) &&
                (abs(ii)+abs(jj)) <= DESERT_WIDTH && elev[x+ii][y+jj]>0 &&
                elev[x+ii][y+jj] < H_MTN)
              if (wet[x+ii][y+jj] < 0.15 && R_ONE<0.3 && y<3.*SIZE/4.)
                terr[x+ii][y+jj] = SCORCH; else
                terr[x+ii][y+jj] = DESERT;
        x += 2*R_1I1; y += 2*R_1I1;
      }
      x += R_1I1 * R_ONE * N_DES_FAR;
      y += R_1I1 * R_ONE * N_DES_FAR;
      // deserts in dry areas and along shores and low elevations
      if ((wet[x][y]<0.25 || elev[x][y]<0.04) && R_ONE>0.1) i--;
    }
  }

  //================================================================
  // PLACE TUNDRA
  //  note: this turns some mountains (but not peaks) into glaciers
  for (j=SIZE-1; j>TUNDRA_J*SIZE; j--)
    for (i=0; i<SIZE; i++)
      if ((j-TUNDRA_J*SIZE)/(SIZE*(1-TUNDRA_J))*pow(elev[i][j]+0.6,3.0)>R_ONE
           && elev[i][j]>0 && terr[i][j]!=PEAK)
        if (terr[i][j]==MOUNTAIN) terr[i][j]=GLACIER; else terr[i][j] = TUNDRA;

  //================================================================
  // SET BOUNDARIES
  // glacier in north;
  // scorching sands in south
  // future boundaries: mountains in east or west; forests in west or east;
  //   (place forest on wetter edge, mountain on higher)

  // glacier
  // multiple random walks from top
  for (i=0; i<SIZE; i++) {
    j = SIZE-1;
    x = i; y = j;
    if (elev[x][y]>0 || R_ONE>0.9) {
      for (k=0; k<N_GLACIER_WALKS; k++) {
        for (ii=-GLACIER_WIDTH; ii<=GLACIER_WIDTH; ii++) {
          for (jj=-GLACIER_WIDTH; jj<=GLACIER_WIDTH; jj++) {
            if ((abs(ii)+abs(jj))<=GLACIER_WIDTH
                && (y+jj)<SIZE && (x+ii)>=0 && (x+ii)<SIZE)
              terr[x+ii][y+jj] = GLACIER;
          }
        }
        if (pow(elev[i][j]+0.6,3.0)>R_ONE
            && (y)<SIZE && (x)>=0 && (x)<SIZE && R_ONE>0.1)
          k--;
        x += R_1I1; y-= (int)(rndmf(2)); // y += R_1I1;
      }
    }
  }

  // scorching sands
  float * jdist;
  int jmax;
  float jpow;
  jdist = alloc1DFractArray(SIZE);
  fill1DFractArray(jdist, SIZE, TIMENOW, 1, .05);

  for (i=0; i<SIZE; i++) {
    jmax = SIZE*SCORCH_J*(jdist[i]+1.0)*(1+0.05*(R_1F1));
    for (j=0; j<jmax; j++) {
      jpow = pow((jmax-j)/jmax+0.7,5.0);
      if (elev[i][j]>0 && terr[i][j]!=PEAK &&
          (jpow>R_ONE || (elev[i][j]<0.075 && R_ONE>0.1) ||
           (j<jmax*0.1 && R_ONE>0.3))) {
        if (j>0.8*jmax && R_ONE>0.5) terr[i][j] = DESERT; else terr[i][j] = SCORCH;
      }
    }
  }

  //================================================================================
  // END OF TERRAIN=================================================================
  // START CIVILIZATIONS============================================================
  //================================================================================

  //================================================================
  // place home towns

  int town_x[SHALLOW+1], town_y[SHALLOW+1];

  // starting town
  while (1) {
    x = R_SZ/2. + SIZE/4.;
    y = R_SZ/2. + SIZE/4.;
    if (elev[x][y] < 0.5) break;
  }
  town_x[0] = x; town_y[0] = y;
  terr[x][y] = T_START;

  // normal towns
  for (i=PLAIN; i <= SHALLOW; i++) {
    while (1) {
      x = R_SZ;
      y = R_SZ;
      if (terr[x][y]==i) break;
    }
    town_x[i] = x; town_y[i] = y;
    terr[x][y] = T_START+i;
  }

  // extreme towns

  // report
  printf("\nTOWNS\n");
  for (i=0; i <= SHALLOW; i++)
    printf("%i: (%i %i)[%i].\n", i,town_x[i],town_y[i], terr[town_x[i]][town_y[i]]);

  /*
  // check accessibility to normal towns
  for (i=1; i <= SHALLOW; i++) {
  if (!exist_path(town_x[0],town_y[0], town_x[i],town_y[i]))
  printf("!!!NO PATH!!! (%i)", i);
  */

  // place other dungeons, etc., and check accessibility

  //================================================================================
  //================================================================================
  //================================================================================
  // FINISHED: JUST OUTPUT

  // J, THEN I, to make map look correct
  FILE *fp_elev;
  FILE *fp_terr;
  FILE *fp_wet;
  fp_elev = fopen("e.txt", "wb");
  fp_terr = fopen("t.txt", "wb");
  fp_wet  = fopen("w.txt", "wb");

  //rare odd data in terrain: mash to deep
  //where could this possibly come from...?
  for (i=0; i<SIZE; i++) {
    for (j=0; j<SIZE; j++) {
      if (terr[i][j]<0 || terr[i][j]>1000)
        {terr[i][j]=DEEP; printf("mt..");}
    }
  }
  //havent seen anything wrong in elev, but mash just in case
  for (i=0; i<SIZE; i++) {
    for (j=0; j<SIZE; j++) {
      if (elev[i][j]<-1000 || elev[i][j]>1.1)
        {elev[i][j]=0.5; printf("me..");}
    }
  }


  //================================================================
  //debugging only: to make certain features pop on terrain map
  for (i=0; i<SIZE; i++) {
    for (j=0; j<SIZE; j++) {
      if (terr[i][j]==WOODS)  terr[i][j]=70;
      if (terr[i][j]==DESERT)  terr[i][j]=140;
      if (terr[i][j]==SCORCH)  terr[i][j]=180;
      if (terr[i][j]==DEEP)  terr[i][j]=-1;
    }
  }
  //debugging only: to make certain features pop on elevation map
  /*
  for (i=0; i<SIZE; i++) {
    for (j=0; j<SIZE; j++) {
      //      if (terr[i][j]==70) elev[i][j]=1.4;
    }
  }
  */

  //================================================================
  // OUTPUT
  for (j=0; j<SIZE; j++) {
    for (i=0; i<SIZE; i++) {
      fprintf(fp_elev, "%.2f  ", elev[i][j]);
      fprintf(fp_terr, "%2i   ", terr[i][j]);
      fprintf(fp_wet,  "%.2f  ",  wet[i][j]);
    }
    fprintf(fp_elev, "\n");
    fprintf(fp_terr, "\n");
    fprintf(fp_wet,  "\n");
  }
  fclose(fp_elev); fclose(fp_terr);

  printf("\nnowtime: %u. DONE.\n", nowtime);
}