/* 
	  Compile this program with:

		cc -o fracland fracland0.9.c -lm

-----------------------------------------------------------------------------

				Fracland


	A fractal landscape generator.

 	Copyright (C) 1992, Lawrence K. Coffin, John Stone, James Sullivan.
 	All rights reserved.
 
	Parts of this program are from "Numerical Recipes in C" by William H.
	Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling.

  This software may be freely copied, modified, and redistributed
  provided that this copyright notice is preserved on all copies.
 
  You may not distribute this software, in whole or in part, as part of
  any commercial product without the express consent of the authors.
  
  There is no warranty or other guarantee of fitness of this software
  for any purpose.  It is provided solely "as is".

 
  This is a program for generating fractal landscapes for Rayshade.  Output is
  in the Rayshade height field format as specified in the Rayshade 4.0 manual.
  The maximum and minimum possible altitudes range from a maximum of 1.0 to 
  a minimum of -1.0.

  This program was written by Larry Coffin, John Stone, and James Sullivan.
  Parts of this program are from "Numerical Recipies" by Someone Something.

 
  Useage:
 
  fracland [options] filename

   Scenefile:

	The height field file will be named "filename".  If no name
	is supplied, the output will be sent to "out.hf".

   Note: All errors and runtime comments are printed to stderr.
		(see the "-q" option)
 
   Options:

	-h 		Print a short summary of the commands
	-d divisions	Generate n subdivisions (default: 3)
	-s seed		Set the seed number for generating random noise in the
				altitudes of the vertces (default: time)
	-o rawoutput	Set the name for the raw output file (default: null)
	-q		Quiet.  Does not report statistics to stderr

  Seed value:

	The default for seed is the time at which the program is executed, this
	creates a new landscape each time.

  Raw output:

	The -o option creates a second file named "rawoutput".  This file is
	a listing of the verteces that make up the landscape a single vertex per
	line tab delimated.

	x1	y1	z1
	x2	y2	z2 ...

	This file is available for importation into a other programs such as
	contour plotting programs which may make it easier to plan scenes using
	the generated fractal landscape. (See also, the program cvhf which 
	displays height field files on Silicon Graphics machines.)

	Comments, suggestions, and bugs reports can be sent to:

	lcoffin@clciris.chem.umr.edu


*/


#include<stdio.h>
#include<math.h>
#include<float.h>
#include<stdlib.h>

float rv(void);
float **matrix(int,int,int,int);
void free_matrix(float**,int,int,int,int);

main(int argc, char *argv[])
{
	int x, y, i, in, i2;
	int c, size;

	int help = 0;
	int quiet = 0;
	int subdiv = 3;
	int seed = time(0);
	int output = 0;

	long int n;
	FILE *out, *rawout;
	float rc, **al;
	char *rawname;
	char *outfile;
	extern char *optarg;
	extern int optind;

/*
 * Check program arguments and reset variables accordingly
 */

	while((c = getopt(argc, argv, "hqd:s:o:")) != EOF)
		switch(c) {
		case 'h': help = 1;
			break;
		case 'q': quiet = 1;
			break;
		case 'd': subdiv = atoi(optarg);
			break;
		case 's': seed = atoi(optarg);
			break;
		case 'o': rawname = optarg;
			output = 1;
			break;
		case '?': help = 1;
	}
	if((outfile = argv[optind]) == NULL) {
		outfile = "out.hf";
	}

/*
 * Print out the useage and some other stuff about fracland
 */	
	if(help) {
		fprintf(stderr,"\n");
		fprintf(stderr,"  Fracland v0.9\n");
		fprintf(stderr,"    Copyright 1992 by Lawrence Coffin, John Stone and James Sullivan.\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"  Useage:\n");
		fprintf(stderr,"\n");
 		fprintf(stderr,"  fracland [options] filename\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"   Scenefile:\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"	The height field file will be named \"filename\".  If no name\n");
		fprintf(stderr,"	is supplied, the output will be sent to \"out.hf\".\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"   Note: All errors and runtime comments are printed to stderr.\n");
		fprintf(stderr,"		(see the \"-q\" option)\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"   Options:\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"	-h 		Print a short summary of the commands\n");
		fprintf(stderr,"	-d divisions	Generate n subdivisions (default: 3)\n");
		fprintf(stderr,"	-s seed		Set the seed number for generating random noise in the\n");
		fprintf(stderr,"				altitudes of the vertces (default: time)\n");
		fprintf(stderr,"	-o rawoutput	Set the name for the raw output file (default: null)\n");
		fprintf(stderr,"	-q		Quiet.  Does not report statistics to stderr\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"  Seed value:\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"	The default for seed is the time at which the program is executed, this\n");
		fprintf(stderr,"	creates a new landscape each time.\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"  Raw output:\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"	The -o option creates a second file named \"rawoutput\".  This file is\n");
		fprintf(stderr,"	a listing of the verteces that make up the landscape a single vertex per\n");
		fprintf(stderr,"	line, tab delimated.\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"	x1	y1	z1\n");
		fprintf(stderr,"	x2	y2	z2 ...\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"	This file is available for importation into a other programs such as\n");
		fprintf(stderr,"	contour plotting programs which may make it easier to plan scenes using\n");
		fprintf(stderr,"	the generated fractal landscape.\n");
		fprintf(stderr,"\n");
	exit(1);
	}

	if((out=fopen(outfile,"w")) == NULL) {
		fprintf(stderr,"Error opening %s.\n",outfile);
		exit(1);
	}
	
	if(output){
		if((rawout=fopen(rawname,"w")) == NULL) {
			fprintf(stderr,"Error opening %s.\n",rawname);
			exit(1);
		}
	}

	srand(seed);
	n=pow(2,subdiv);
/*
 * Allocate memory for the altitude matrix
 */
	al = matrix(0,n,0,n);

/*
 * Generate the altitudes
 */

	al[0][0] = 0;
	al[0][n] = 0;
	al[n][0] = 0;
	al[n][n] = 0;
	for (i=0; i<subdiv; i++){
		in = n/pow(2,i);
		i2 = in/2;
		rc = 1/pow(2,i);
		for (x=0; x<n; x+=in){
			for (y=0; y<n; y+=in){
				al[x+i2][y   ] = (al[x   ][y   ] + al[x+in][y   ])/2 + rv()*rc;
				al[x   ][y+i2] = (al[x   ][y   ] + al[x   ][y+in])/2 + rv()*rc;
				al[x+in][y+i2] = (al[x+in][y   ] + al[x+in][y+in])/2 + rv()*rc;
				al[x+i2][y+in] = (al[x   ][y+in] + al[x+in][y+in])/2 + rv()*rc;
				al[x+i2][y+i2] = (al[x][y] + al[x+in][y] + al[x][y+in] + al[x+in][y+in])/4 + rv()*rc;
			}
		}
	}
	size = n+1;

	if(!(quiet)) {
		fprintf(stderr,"\n   Number of subdivisions is:\t%i\n",subdiv);
		fprintf(stderr,"   Number of points per side is:\t%i\n", size);
		fprintf(stderr,"   Seed is:\t\t\t%i\n",seed);
		if(output){
			fprintf(stderr,"   Raw output file is:\t%s\n",rawname);
		}
		fprintf(stderr,"\n");
	}

/*
 * Print out the points into the heightfield file
 */

	fwrite(&size,sizeof(size),1,out);
	i = 0;
        for (y=0; y<=n; y++){
		for (x=0; x<=n; x++){
			fwrite(&al[x][y],sizeof(al[x][y]),1,out);
			if(output){
				fprintf(rawout,"%6.3f\t%6.3f\t%6.3f\n",((float)(i%size)/(float)(size-1)),((float)(i/size)/(float)(size-1)),al[x][y]);
				i++;
			}
		}
	}
	
/*
 * Free memory used for the altitude matrix
 */

	free_matrix(al,0,n,0,n);

	if(output){
		fclose(rawout);
	}


	fclose(out);
	
}

/*
 * Generate a random value from -0.5 to 0.5
 */

float rv(void){
	float y;
	y = (float)rand();
	y /= RAND_MAX;
	y -= 0.5;
	return (y);
}

/*
 * Generate an error message if error in allocating matrix
 */

void nrerror(error_text)
char error_text[];
{
	void exit();

	fprintf(stderr,"Numerical Recipes run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}

/*
 * Allocate memory for a matrix
 */

float **matrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
	int i;
	float **m;

	m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float*));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float));
		if (!m[i]) nrerror("allocation failure 2 in matrix()");
		m[i] -= ncl;
	}
	return m;
}

/*
 * Free memory used by a matrix
 */

void free_matrix(m,nrl,nrh,ncl,nch)
float **m;
int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}