/* rivers.c */
/* map generating program */
/* Copyright 2026 Torben Æ. Mogensen */

/* The program generates maps based on recursivesubdivision */
/* of triangles. The output is a colour BMP bitmap. */
/* It is assumed that pixels are square */
/* I have included a procedure to print the maps as .bmp */
/* on standard output or specified files. */

/* I have tried to avoid using machine specific features, so it should */
/* be easy to port the program to any machine. Beware, though that due */
/* to different precision on different machines, the same seed numbers */
/* can yield very different planets. */

/* The primitive user interface is primarily a result of portability concerns */

#ifdef THINK_C
#define macintosh 1
#endif

#ifdef macintosh
#include <console.h>
#include <unix.h>
#endif

#define O_BINARY 0x8000

#include <stdio.h>
#include <errno.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h> // Added boolean type

int BLACK = 0;
int WHITE = 1;
int BACK = 2;
int GRID = 3;
int OUTLINE1 = 4;
int OUTLINE2 = 5;
int LOWEST = 6;
int SEA = 7;
int LAND = 8;
int HIGHEST = 9;

//int debug = 0;

int nocols = 65536;

unsigned char rtable[65536], gtable[65536], btable[65536];

char cmdLine[] = "./planet -ps -s 0.31 -m 3 -L 30 -l 35 -G 10 -g 20 -w 320 -h 256 -o twodee.bmp";

#define PI 3.14159265358979323846
#define DEG2RAD 0.01745329251994329576923 /* pi/180 */

/* these three values can be changed to change world characteristica */

double M  = -.02;   /* initial altitude (slightly below sea level) */
double dd1 = 0.45;  /* weight for altitude difference */
double POWA = 1.0; /* power for altitude difference */
double dd2 = 0.035; /* weight for distance */
double POW = 0.47;  /* power for distance function */

int Depth; /* depth of subdivisions */
double r1,r2,r3,r4; /* seeds */
double longi,lat,scale;

int Width = 1024, Height = 1024; /* default map size */

int **col;  /* colour array */

int do_outline = 0;  /* if 1, draw coastal outline */
int do_bw = 0;       /* if 1, reduce map to black outline on white */
int contourstep = 0; /* if >0, # of colour steps between contour lines */
int *outx, *outy;


double cla, sla, clo, slo, rseed;

int matchMap = 0;
double matchSize = 0.1;

char colorsname[256] = "Olsson.col";

typedef struct Vert
{
  double x,y; /* coordinates */
  double h; /* altitude */
  double s; /* seed */
} vert;

typedef struct Edge
{
  double h, w; /* height and width of rivers. 0.0 means no river */
} edge;
edge None = {.h=0.0, .w=0.0};

vert inittriangle[3];
edge initedges[3];

int rivers = 1; // make rivers
double curves = 0.0; // make contour lines
double Wide = 150.0; /* when is a river wide */

double min3(double x, double y, double z) {
  return fmin(x,fmin(y,z));
}

double max3(double x, double y, double z) {
  return fmax(x,fmax(y,z));
}

double branchsize(double x, double y) { /* size of branching river */
  return 0.2*(x+y);
}

double magic = 100.0*PI;

double mix(double s1, double s2) { /* mix two seeds */
  return fmod((s1+magic)*(s2+magic), 2.0) - 1.0;
}

double nxt(double s) { /* next seed */
  return mix(s, s);
}

/* random value between a and b, tend towards middle */
double between(double a, double b, double seed) {
  return (a+b+seed*seed*seed*(a-b))/2.0;
}

/* distance between vertices */
double distance(vert a, vert b)
{
  double abx, aby;
  abx = a.x-b.x; aby = a.y-b.y;
  return sqrt(abx*abx+aby*aby);
}

double phi = (sqrt(5.0)-1.0)/2.0;
double phi1 = 1.0-(sqrt(5.0)-1.0)/2.0;

int islands = 0; // make islands in fjords

double **heights;
int **river;

// possibly extend river
edge extend(vert v, double e, double w,
	    vert vNear, vert vFar, double s, double dist) {
  edge e3;
  if  (v.h < fmin(0.0, e) && 0.0 < vFar.h) { // extend down
    e3.h = between(e, v.h, s); e3.w = 1.5*w;
  }
  else if (fabs(s) < 0.75 && e < min3(v.h, vNear.h, vFar.h)) { // extend up
    e3.h = between (e, min3(v.h, vNear.h, vFar.h), nxt(v.s));
    e3.w = 0.7*w;
  }
  else if (w > 2.0*Wide*dist) { // wide river
    e3.h = e; e3.w = w;
  }
  else e3 = None; // don't extend
  
  return e3;
}

// See https://tartarus.org/~simon/20110412-penrose/penrose.xhtml

// Triangle type 1: 36/72/72 degree angles, v2 at 36-degree angle
// Triangle type 2: 36/36/108 degree angles, v0 at 108-degree angle
// ei opposite vi

void tria1(vert v0, vert v1, vert v2,
	   edge e0, edge e1, edge e2,
	   double xmin, double xmax, double ymin, double ymax,
	   double epsilon, int type1) {
// v3 splits e0 into e4 and e5 at ratio phi (e4 longest)
// e3 between v0 and v3
// e4 between v2 and v3
// e5 between v1 and v3
  if (v0.x < xmin && v1.x < xmin && v2.x < xmin ||
      v0.x > xmax && v1.x > xmax && v2.x > xmax ||
      v0.y < ymin && v1.y < ymin && v2.y < ymin ||
      v0.y > ymax && v1.y > ymax && v2.y > ymax)
    return; // outside bounding box
  else {
    double dist = distance(v1, v2);

    // altitude variation depends on both horisontal and vertical distance
    double delta = 0.32*dist + 0.55*fabs(v1.h-v2.h);

    // attributes for v3
    double s3 = mix(v1.s, v2.s);
    double x3 = phi*v1.x+phi1*v2.x;
    double y3 = phi*v1.y+phi1*v2.y;

    edge e4, e5; double h3;

    if (e0.w == 0.0) {
      // no river on e0, v3 height by simple midpoint displacement
      e4 = e5 = None;
      h3 = phi*v1.h+phi1*v2.h + s3*delta;
    }
    else {
      double h = e0.h;
      double w = e0.w;
      // if islands=true, copy river (half width) to both edges
      // with low prob if dist relatively small and e0 under low sea
      if (islands && h < -0.005 && h > -0.05 && dist < 0.2
	  && fabs(nxt(s3)) > 0.9) {
	e4.h = h; e4.w = 0.5*w;
	e5.h = h; e5.w = 0.5*w;
	h3 = (h+v1.h+v2.h)/3.0 + s3*delta;
      } else {
	// copy river towards closest-height vertex
	// and v3 height from river height and opposite vertex height
	if (fabs(h-v1.h) > fabs(h-v2.h)) {
	  e4 = e0; e5 = None; h3 = (h+v1.h)/2.0 + s3*delta;
	} else {
	  e4 = None; e5 = e0; h3 = (h+v2.h)/2.0 + s3*delta;
	}
      }
    }
    vert v3;
    edge e3;
    v3.x=x3; v3.y=y3; v3.s=s3; v3.h=h3;

    if (dist < epsilon) { // plot point
      int i = (int) (round (((v0.x+v1.x+v2.x)/3.0-xmin) / epsilon));
      int j = (int) (round (((v0.y+v1.y+v2.y)/3.0-ymin) / epsilon));
      if (i < 0 || i >= Width || j < 0 || j >= Height)
	return; // outside bounding box
      else {
	if (heights[i][j] <= -2.0) heights[i][j] = h3;
	else heights[i][j] = (heights[i][j] + h3)/2.0;
	river[i][j] = river[i][j] || e0.w != 0.0 || e1.w != 0.0 || e2.w != 0.0;
	return;
      }
    } else { // subdivide
      // find attributes for e3
      e3 = None;
      double s03 = mix(v0.s, s3); // seed for e3
      
      if (e1.w == 0.0 && e2.w == 0.0 && e4.w == 0.0 && e5.w == 0.0) {
	if (rivers == 1 && fmax(v1.h, v2.h) > 0.1
	    && fmin(v1.h, v2.h) < min3(v0.h, v3.h, -0.05)) { // make new river
	  e3.h = between(fmin(v1.h, v2.h), fmin(v0.h, v3.h), s03);
	  e3.w = 0.5*dist;
	} else e3 = None;
      }

      else if (e1.w > 0.0 && e2.w == 0.0 && e4.w == 0.0 && e5.w == 0.0) 
	e3 = extend(v1, e1.h, e1.w, v0, v3, s03, dist);

      else if (e1.w == 0.0 && e2.w > 0.0 && e4.w == 0.0 && e5.w == 0.0) 
	e3 = extend(v2, e2.h, e2.w, v0, v3, s03, dist);

      else if (e1.w == 0.0 && e2.w == 0.0 && e4.w > 0.0 && e5.w == 0.0) 
	e3 = extend(v1, e4.h, e4.w, v3, v3, s03, dist);

      else if (e1.w > 0.0 && e2.w == 0.0 && e4.w == 0.0 && e5.w > 0.0) {
	e3.h = between (e1.h, e5.h, s03);
	e3.w = between(e1.w, e5.w, nxt(s03));
      }

      else if (e1.w == 0.0 && e2.w > 0.0 && e4.w > 0.0 && e5.w == 0.0) {
	e3.h = between (e2.h, e4.h, s03);
	e3.w = between(e2.w, e4.w, nxt(s03));
      }

      else if (e1.w > 0.0 && e2.w > 0.0 && e4.w == 0.0 && e5.w == 0.0) {
	e3.h = between (e1.h, e2.h, s03);
	e3.w = between(e1.w, e2.w, nxt(s03));
      }

      else if (e1.w > 0.0 && e2.w == 0.0 && e4.w > 0.0 && e5.w == 0.0) {
	if (fmax(e1.w, e4.w) > Wide*dist) { // wide river
	  // make wide branch
	  e3.h = between(fmin(e1.h, e4.h), min3(v1.h, v0.h, v3.h), nxt(s03));
	  e3.w = fmax(e1.w, e4.w);
	}
	else if (fabs(s03) < 2.75*dist &&
		 min3(v1.h, v0.h, v3.h) > fmin(e4.h, e1.h)) { // make branch
	  e3.h = between(fmin(e1.h, e4.h), min3(v1.h, v0.h, v3.h), nxt(s03));
	  e3.w = branchsize(e1.w, e4.w);
	}
	else e3 = None;
      }

      else if (e1.w == 0.0 && e2.w > 0.0 && e4.w == 0.0 && e5.w > 0.0) {
	if (fmax(e2.w, e5.w) > Wide*dist) { // wide river
	  // make wide branch
	  e3.h = between(fmin(e2.h, e5.h), min3(v2.h, v0.h, v3.h), nxt(s03));
	  e3.w = fmax(e2.w, e5.w);
	}
	else if (fabs(s03) < 2.75*dist &&
		 min3(v2.h, v0.h, v3.h) > fmin(e5.h, e2.h)) { // make branch
	  e3.h = between(fmin(e2.h, e5.h), min3(v2.h, v0.h, v3.h), nxt(s03));
	  e3.w = branchsize(e2.w, e5.w);
	}
	else e3 = None;
      }

      else if (e1.w > 0.0 && e2.w > 0.0 && e4.w > 0.0 && e5.w == 0.0) {
	e3.h = between(e2.h, fmin(e1.h, e4.h), s03);
	e3.w = between(e2.w, fmin(e1.w, e4.w), nxt(s03));
      }

      else if (e1.w > 0.0 && e2.w > 0.0 && e4.w == 0.0 && e5.w > 0.0) {
	e3.h = between(e1.h, fmin(e2.h, e5.h), s03);
	e3.w = between(e1.w, fmin(e2.w, e5.w), nxt(s03));
      }

      else if (e1.w == 0.0 && e2.w == 0.0 && e4.w > 0.0 && e5.w > 0.0)
	e3 = None;

      else if (e1.w > 0.0 && e2.w == 0.0 && e4.w > 0.0 && e5.w > 0.0) {
	e3.h = between(e1.h, e4.h, s03);
	e3.w = between(e1.w, e4.w, nxt(s03));
      }

      else if (e1.w == 0.0 && e2.w > 0.0 && e4.w > 0.0 && e5.w > 0.0) {
	e3.h = between(e2.h, e5.h, s03);
	e3.w = between(e2.w, e5.w, nxt(s03));
      }

      else if (e1.w > 0.0 && e2.w > 0.0 && e4.w > 0.0 && e5.w > 0.0)
	e3 = None;
    }
    if (type1) {
      tria1(v1, v3, v0, e3, e2, e5, xmin, xmax, ymin, ymax, epsilon, 1);
      tria1(v3, v2, v0, e1, e3, e4, xmin, xmax, ymin, ymax, epsilon, 0);
    } else {
      tria1(v3, v0, v2, e1, e4, e3, xmin, xmax, ymin, ymax, epsilon, 1);
      tria1(v3, v0, v1, e2, e5, e3, xmin, xmax, ymin, ymax, epsilon, 0);
    }
  }
}

void readcolors(FILE *colfile, char* colorsname)
{
  int crow, cNum = 0, oldcNum, i;

  if (NULL == (colfile = fopen(colorsname, "r")))
    {
      fprintf(stderr,
              "Cannot open %s\n",
              colorsname);
      exit(1);
    }


  /* Format of colour file is a sequence of lines       */
  /* each consisting of four integers:                  */
  /*   colour_number red green blue                     */
  /* where 0 <= colour_number <= 65535                  */
  /* and 0<= red, green, blue <= 255                    */
  /* The colour numbers must be increasing              */
  /* The first colours have special uses:               */
  /* 0 is usually black (0,0,0)                         */
  /* 1 is usually white (255,255,255)                   */
  /* 2 is the background colour                         */
  /* 3 is used for latitude/longitude grid lines        */
  /* 4 and 5 are used for outlines and contour lines    */
  /* 6 upwards are used for altitudes                   */
  /* Halfway between 6 and the max colour is sea level  */
  /* Shallowest sea is (max+6)/2 and land is above this */
  /* With 65536 colours, (max+6)/2 = 32770              */
  /* Colours between specified are interpolated         */

  for (crow = 0; !feof(colfile); crow++)
    {
      int       rValue, 
                gValue, 
                bValue,
                result = 0;

      oldcNum = cNum;  /* remember last colour number */
      result = fscanf(colfile, " %d %d %d %d",
                      &cNum, &rValue, &gValue, &bValue);

      if (result > 0)
        {
          if (cNum < oldcNum) cNum = oldcNum;
          if (cNum > 65535) cNum = 65535;
          rtable[cNum] = rValue;
          gtable[cNum] = gValue;
          btable[cNum] = bValue;
          /* interpolate colours between oldcNum and cNum */
          for (i = oldcNum+1; i<cNum; i++) {
            rtable[i] = (rtable[oldcNum]*(cNum-i)+rtable[cNum]*(i-oldcNum))
                        / (cNum-oldcNum+1);
            gtable[i] = (gtable[oldcNum]*(cNum-i)+gtable[cNum]*(i-oldcNum))
                        / (cNum-oldcNum+1);
            btable[i] = (btable[oldcNum]*(cNum-i)+btable[cNum]*(i-oldcNum))
                        / (cNum-oldcNum+1);
          }
        }
    }

  fclose(colfile);

  nocols = cNum+1;
  if (nocols < 10) nocols = 10;

  HIGHEST = nocols - 1;
  SEA = (HIGHEST+LOWEST)/2;
  LAND = SEA+1;

  for (i = cNum+1; i<nocols; i++) {
    /* fill up rest of colour table with last read colour */
    rtable[i] = rtable[cNum];
    gtable[i] = gtable[cNum];
    btable[i] = btable[cNum];
  }
}

void printbmp(outfile) /* prints picture in BMP format */
FILE *outfile;
{
  int i,j,c,s0, s, W1;

  fprintf(outfile,"BM");

  W1 = (3*Width+3);
  W1 -= W1 % 4;
  s0 = (strlen(cmdLine)+strlen("Command line:\n\n")+3) & 0xffc;
  s = s0 + 54+W1*Height; /* file size */
  putc(s&255,outfile);
  putc((s>>8)&255,outfile);
  putc((s>>16)&255,outfile);
  putc(s>>24,outfile);

  putc(0,outfile);
  putc(0,outfile);
  putc(0,outfile);
  putc(0,outfile);

  putc(54,outfile); /* offset to data */
  putc(0,outfile);
  putc(0,outfile);
  putc(0,outfile);

  putc(40,outfile); /* size of infoheader */
  putc(0,outfile);
  putc(0,outfile);
  putc(0,outfile);

  putc(Width&255,outfile);
  putc((Width>>8)&255,outfile);
  putc((Width>>16)&255,outfile);
  putc(Width>>24,outfile);

  putc(Height&255,outfile);
  putc((Height>>8)&255,outfile);
  putc((Height>>16)&255,outfile);
  putc(Height>>24,outfile);

  putc(1,outfile);  /* no. of planes = 1 */
  putc(0,outfile);

  putc(24,outfile);  /* bpp */
  putc(0,outfile);

  putc(0,outfile); /* no compression */
  putc(0,outfile);
  putc(0,outfile);
  putc(0,outfile);

  putc(0,outfile); /* image size (unspecified) */
  putc(0,outfile);
  putc(0,outfile);
  putc(0,outfile);

  putc(0,outfile); /* h. pixels/m */
  putc(32,outfile);
  putc(0,outfile);
  putc(0,outfile);

  putc(0,outfile); /* v. pixels/m */
  putc(32,outfile);
  putc(0,outfile);
  putc(0,outfile);

  putc(0,outfile); /* colours used (unspecified) */
  putc(0,outfile);
  putc(0,outfile);
  putc(0,outfile);


  putc(0,outfile); /* important colours (all) */
  putc(0,outfile);
  putc(0,outfile);
  putc(0,outfile);

  for (j=Height-1; j>=0; j--) {
    for (i=0; i<Width; i++) {
      putc(btable[col[i][j]],outfile);
      putc(gtable[col[i][j]],outfile);
      putc(rtable[col[i][j]],outfile);
    }
    for (i=3*Width; i<W1; i++) putc(0,outfile);
  }

  fprintf(outfile,"Command line:\n%s\n",cmdLine);
  fclose(outfile);
}

void array2bmp() {
  int ncols = nocols-6;
  for (int i=0; i<Width; i++)
    for (int j=0; j<Height; j++) {
      if (heights[i][j] <= -2.0) col[i][j] = 0;
      else if (river[i][j] && heights[i][j] >= -0.003)
	col[i][j] = 6 + ncols/2;
      else if (curves>0.0 && heights[i][j] > 0.0 &&
            i>0 && i<Width-1 && j>0 && j<Height-1 &&
            ((int) (heights[i][j]*curves) > (int) (heights[i-1][j]*curves) ||
             (int) (heights[i][j]*curves) > (int) (heights[i+1][j]*curves) ||
             (int) (heights[i][j]*curves) > (int) (heights[i][j-1]*curves) ||
             (int) (heights[i][j]*curves) > (int) (heights[i][j+1]*curves)))
	col[i][j] = 4 + ((int) (heights[i][j]*curves)) % 2;
      else {
	int c = 6 + ncols*(heights[i][j]+1.0)/2.0;
	if (c < 6) c = 0;
	else if (c >= ncols) c = 1;
	col[i][j] = c;
      }
    }
}

int main(int argc, const char* argv[] )
{
  double seed = 0.5;
  double xmid = 0.5;
  double ymid = 0.5;
  double scale = 1.0;
  int i = 0;
  char outfilename[100] = "land";
  char outfilename1[105];
  char colorsname[100] = "OlssonW.col";

  double h00 = -3.0;
  double h01 = -3.0;
  double h10 = -3.0;
  double h11 = -3.0;
  
  // read command line parameters:
  // -s <float> : set seed
  // -x <float> : set x of centre (between 0.0 and 1.0)
  // -y <float> : set y of centre (between 0.0 and 1.0)
  // -m <float> : set magnification
  // -c <float> : set contour lines
  // -o <file name> : set output file name (without .bmp extension)
  // -C <file name> : set colour file name (without .col extension)
  // -r : flip river generation (default on)
  // -i : flip islands in fjords (default off)
  // -N <float> : set altitude of north point
  // -E <float> : set altitude of east point
  // -S <float> : set altitude of south point
  // -W <float> : set altitude of west point

  while (i<argc) {
    if (argv[i][0] == '-')
      switch (argv[i][1]) {
      case 's' : if (++i < argc &&  sscanf(argv[i],"%lf",&seed)) break;
      case 'x' : if (++i < argc &&  sscanf(argv[i],"%lf",&xmid)) break;
      case 'y' : if (++i < argc &&  sscanf(argv[i],"%lf",&ymid)) break;
      case 'm' : if (++i < argc &&  sscanf(argv[i],"%lf",&scale)) break;
      case 'c' : if (++i < argc &&  sscanf(argv[i],"%lf",&curves)) break;
      case 'N' : if (++i < argc &&  sscanf(argv[i],"%lf",&h01)) break;
      case 'E' : if (++i < argc &&  sscanf(argv[i],"%lf",&h00)) break;
      case 'W' : if (++i < argc &&  sscanf(argv[i],"%lf",&h10)) break;
      case 'S' : if (++i < argc &&  sscanf(argv[i],"%lf",&h11)) break;
      case 'w' : if (++i < argc &&  sscanf(argv[i],"%d",&Width)) break;
      case 'h' : if (++i < argc && sscanf(argv[i],"%d",&Height)) break;
      case 'o' : if (++i < argc && sscanf(argv[i],"%255[^\n]",outfilename)) break;
      case 'C' : if (++i < argc && sscanf(argv[i],"%255[^\n]",colorsname)) break;
      case 'r' : { rivers = 1-rivers; break; }
      case 'i' : { islands = 1-islands; break; }
      }
    i++;
  }
  if (xmid < 0.0) xmid = 0.0;
  if (xmid > 1.0) xmid = 1.0;
  if (ymid < 0.0) ymid = 0.0;
  if (ymid > 1.0) ymid = 1.0;
  ymid = 1.0 - ymid;
  if (scale < 1.0) scale = 1.0;
  if (rivers == 0) islands = 0; // no islands if no rivers 

  double seed00 = mix(seed, PI);
  double seed01 = mix(seed, 2.234551);
  double seed10 = mix(seed, 5.629339);
  double seed11 = mix(seed, 9.075109);
  double xmin = fmax(0.0, xmid - 0.5 / scale);
  double ymin = fmax(0.0, ymid - 0.5 / scale);
  double xmax = fmin(1.0, xmid + 0.5 / scale);
  double ymax = fmin(1.0, ymid + 0.5 / scale);
  double cs = cos(36.0*PI/180.0);
  double sn = sin(36.0*PI/180.0);
  double cs1 = 0.5 * cs / sn;
  double w = 0.5 / sn;

  if (h00 == -3.0) h00 = between(-1.0,1.0,seed00);
  if (h01 == -3.0) h01 = between(-1.0,1.0,seed01);
  if (h10 == -3.0) h10 = between(-1.0,1.0,seed10);
  if (h11 == -3.0) h11 = between(-1.0,1.0,seed11);

  edge river0;

  if (rivers == 1 && fmax(h00, h11) > 0.1
      && fmin(h00, h11) < min3(h01, h10, -0.1)) { // make new river
    river0.h = between(fmin(h00, h11), fmin(h01, h10), nxt(seed11));
    river0.w = 0.5;
  } else river0 = None;

  col = (int**)calloc(Width,sizeof(int*));
  if (col == 0) {
    fprintf(stderr, "Memory allocation failed.");
    exit(1);
  }
  for (int i=0; i<Width; i++) {
    col[i] = (int*)calloc(Height,sizeof(int));
    if (col[i] == 0) {
      fprintf(stderr,
              "Memory allocation failed at %d out of %d cols\n",
              i+1,Width);
      exit(1);
    }
  }

  heights = (double**)calloc(Width,sizeof(double*));
  if (heights == 0) {
    fprintf(stderr, "Memory allocation failed.");
    exit(1);
  }
  for (int i=0; i<Width; i++) {
    heights[i] = (double*)calloc(Height,sizeof(double));
    if (heights[i] == 0) {
      fprintf(stderr,
              "Memory allocation failed at %d out of %d cols\n",
              i+1,Width);
      exit(1);
    }
  }

  river = (int**)calloc(Width,sizeof(int*));
  if (river == 0) {
    fprintf(stderr, "Memory allocation failed.");
    exit(1);
  }
  for (int i=0; i<Width; i++) {
    river[i] = (int*)calloc(Height,sizeof(int));
    if (river[i] == 0) {
      fprintf(stderr,
              "Memory allocation failed at %d out of %d cols\n",
              i+1,Width);
      exit(1);
    }
  }

  for (int i = 0; i<Width; i++)
    for (int j = 0; j<Height; j++) {
      heights[i][j] = -3.0;
      river[i][j] = 0;
    }
  
  vert v1, v2, v3;

  v1.x = cs1; v1.y = 1.0; v1.s = seed11; v1.h = h11;
  v2.x = w; v2.y = 0.5; v2.s = seed10; v2.h = h10;
  v3.x = 0.0; v3.y = 0.5; v3.s = seed00; v3.h = h00;

  tria1(v1, v2, v3, river0, None, None,
        xmin, xmax, ymin, ymax,
        (1.0 / (double) Width / scale),
        1);

  v1.x = cs1; v1.y = 0.0; v1.s = seed01; v1.h = h01;
  v2.x = w; v2.y = 0.5; v2.s = seed10; v2.h = h10;
  v3.x = 0.0; v3.y = 0.5; v3.s = seed00; v3.h = h00;

  tria1(v1, v2, v3, river0, None, None,
        xmin, xmax, ymin, ymax,
        (1.0 / (double) Width / scale),
        1);


  FILE *outfile, *colfile = NULL;
  readcolors(colfile, colorsname);

  array2bmp();

  sprintf(outfilename1,"%s.bmp",outfilename);
  
  outfile = fopen(outfilename1,"wp");
  
  printbmp(outfile);
}