/*
 * mountain.c - creates a fractal mountain, using Carpenter's method with a
 *     different extension to square grids.  A pyramid of 4 glass spheres
 *     is added in front of the mountain.  None of the spheres are clipped.
 *     A few of the polygons are clipped.  One light source.
 *
 * Version:  2.2 (11/17/87)
 * Author:  Eric Haines, 3D/Eye, Inc.
 *
 * SIZE_FACTOR determines the number of objects output.
 *     Total triangular polygons = 2 * (4**SIZE_FACTOR)
 *
 *     SIZE_FACTOR	# triangles	# spheres
 *	    1		     8		     4
 *	    2		    32		     4
 *	    3		   128		     4
 *
 *	    6		  8192		     4
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "def.h"
#include "lib.h"

/* output format determines if polygons or true surfaces are used */
static int output_format = OUTPUT_CURVES;

/* Determine which raytracer we will use */
static int raytracer_format = OUTPUT_POVRAY;

/* size factor determines number of polygons */
static int size_factor = 8;

/* fractal dimension - affects variance of z.  Between 2 and 3 */
#define	FRACTAL_DIMENSION 2.2

/* change MOUNTAIN_NO to get a different mountain */
#define	MOUNTAIN_NO 21

/* lower left corner and width of mountain definitions */
#define	X_CORNER -1.0
#define	Y_CORNER -1.0
#define	WIDTH     2.0

/* hashing function to get a seed for the random number generator */
#define	hash_rand(A,B,C) ((((A)<<(23-(C))) + ((B)<<(15-(C))) + ((A)<<(7-(C))))\
                          & 0xffff)

static	long    num_pts ;
static  double	roughness ;

/* create a pyramid of crystal spheres */
void
create_spheres(center)
   COORD4 *center ;
{
    double  angle;
    COORD4  axis, new_pt, pt, sphere;
    long    i;
    MATRIX  mx;

    SET_COORD(axis, 1.0, 1.0, 0.0);
    lib_normalize_coord3(&axis);
    angle = acos((double)(-1.0/3.0));

    /* set center of pyramid and radius */
    SET_COORD4(pt, 0.0, 0.0, center->w * sqrt((double)(3.0/2.0)), 0.0);

    COPY_COORD(sphere, pt);
    ADD2_COORD(sphere, *center);
    sphere.w = center->w;
    lib_output_sphere(&sphere, output_format);

    lib_create_axis_rotate_matrix(mx, &axis, angle);
    lib_transform_coord(&new_pt, &pt, mx);

    for (i=0;i<3;++i) {
	lib_create_rotate_matrix(mx, Z_AXIS, (double)i * 2.0 * PI / 3.0);
	lib_transform_coord( &sphere, &new_pt, mx);
	ADD2_COORD(sphere, *center);
	sphere.w = center->w;
	lib_output_sphere(&sphere, output_format);
        }
}

/*
 * Build mountain section.  If at width > 1, split quadrilateral into four
 * parts.  Else if at width == 1, output quadrilateral as two triangles.
 */
void
grow_mountain(width, ll_x, ll_y, ll_fz, lr_fz, ur_fz, ul_fz)
   long width, ll_x, ll_y;
   double ll_fz, lr_fz, ur_fz, ul_fz;
{
   long    half_width, iz;
   double  l_fx, r_fx, l_fy, u_fy;
   double  lower_fz, right_fz, upper_fz, left_fz, middle_fz;
   long    num_tri, num_tri_vert, num_vert;
   double  rise_height, hside_length;
   COORD4  tri_vert[3];


    if (width == 1) {
	/* calculate x and y coordinates of corners */
	l_fx = X_CORNER + (double)ll_x * WIDTH / (double)num_pts;
	r_fx = X_CORNER + (double)(ll_x+1) * WIDTH / (double)num_pts;
	l_fy = Y_CORNER + (double)ll_y * WIDTH / (double)num_pts;
	u_fy = Y_CORNER + (double)(ll_y+1) * WIDTH / (double)num_pts;

	/* output two triangles for section */
	for (num_tri = 0; num_tri < 2; ++num_tri ) {
	    for (num_vert = 0; num_vert < 3; ++num_vert ) {
		num_tri_vert = (num_vert + num_tri * 2 ) % 4;
		if (num_tri_vert == 0) {
		    SET_COORD(tri_vert[num_vert], l_fx, l_fy, ll_fz);
		    }
		else if ( num_tri_vert == 1 ) {
		    SET_COORD(tri_vert[num_vert], r_fx, l_fy, lr_fz);
		    }
		else if ( num_tri_vert == 2 ) {
		    SET_COORD(tri_vert[num_vert], r_fx, u_fy, ur_fz);
		    }
		else {
		    SET_COORD(tri_vert[num_vert], l_fx, u_fy, ul_fz);
		    }
	        }
	    lib_output_polygon(3, tri_vert);
	    }
        }
    else {
	/* subdivide edges and move in z direction */
	half_width = width >> 1;
	hside_length = (double)half_width * WIDTH / (double)num_pts;
	rise_height = hside_length * roughness;

	/* for each midpoint, find z */
	iz = MOUNTAIN_NO + hash_rand(ll_x + half_width, ll_y, size_factor);
	lower_fz = (ll_fz + lr_fz) / 2.0 + rise_height * lib_gauss_rand(iz);
	iz = MOUNTAIN_NO + hash_rand(ll_x+width, ll_y+half_width, size_factor);
	right_fz = (lr_fz + ur_fz) / 2.0 + rise_height * lib_gauss_rand(iz);
	iz = MOUNTAIN_NO + hash_rand(ll_x+half_width, ll_y+width, size_factor);
	upper_fz = (ur_fz + ul_fz) / 2.0 + rise_height * lib_gauss_rand(iz);
	iz = MOUNTAIN_NO + hash_rand( ll_x, ll_y + half_width, size_factor);
	left_fz  = (ul_fz + ll_fz) / 2.0 + rise_height * lib_gauss_rand(iz);
	iz = MOUNTAIN_NO + hash_rand(ll_x+half_width, ll_y+half_width,
                                     size_factor);
	middle_fz = (ll_fz+lr_fz+ur_fz+ul_fz) / 4.0 +
                    1.4142136 * rise_height * lib_gauss_rand(iz);

        /* For raytracers without autobounding of primitives, we can
           often create a heirarchy of bounding boxes */
        switch (raytracer_format) {
           case OUTPUT_POVRAY: fprintf(stdout, "composite {\n"); break;
           case OUTPUT_QRT:    fprintf(stdout, "BEGIN_BBOX\n"); break;
           }

	/* check subsections for subdivision or output */
	grow_mountain(half_width, ll_x, ll_y,
		      ll_fz, lower_fz, middle_fz, left_fz ) ;
	grow_mountain(half_width, ll_x+half_width, ll_y,
		      lower_fz, lr_fz, right_fz, middle_fz);
	grow_mountain(half_width, ll_x+half_width, ll_y+half_width,
		      middle_fz, right_fz, ur_fz, upper_fz);
	grow_mountain(half_width, ll_x, ll_y+half_width,
		      left_fz, middle_fz, upper_fz, ul_fz);

        switch (raytracer_format) {
           case OUTPUT_POVRAY:
              fprintf(stdout, "bounded_by { ");
              l_fx = X_CORNER+(double)ll_x*WIDTH / (double)num_pts;
              r_fx = X_CORNER+(double)(ll_x+width+1)*WIDTH/(double)num_pts;
              l_fy = Y_CORNER+(double)ll_y*WIDTH/(double)num_pts;
              u_fy = Y_CORNER+(double)(ll_y+width+1)*WIDTH/(double)num_pts;
              fprintf(stdout, "box { <%g %g %g> <%g %g %g> } } }\n",
                      l_fx, l_fy, MIN(ll_fz, lr_fz) - rise_height * width,
                      r_fx, u_fy, MAX(ur_fz, ul_fz) + rise_height * width);
              break;
           case OUTPUT_QRT:    fprintf(stdout, "END_BBOX\n"); break;
           }
        }
}

void
main(argc, argv)
   int argc;
   char *argv[];
{
    COORD4  back_color, obj_color;
    COORD4  center, light;
    COORD4  from, at, up;
    double  ratio;

   /* Start by defining which raytracer we will be using */
   switch (argc) {
      case 3:
	 raytracer_format = atoi(argv[2]);
      case 2:
	 size_factor = atoi(argv[1]);
	 break;
      defaut:
	 break;
      }
   lib_set_output_file(stdout);
   if (raytracer_format == OUTPUT_RTRACE ||
       raytracer_format == OUTPUT_PLG)
      lib_set_raytracer(OUTPUT_DELAYED);
   else
      lib_set_raytracer(raytracer_format);
   lib_set_polygonalization(3, 3);

    /* output viewpoint */
    SET_COORD( from, -1.6, 1.6, 1.7);
    SET_COORD( at, 0.0, 0.0, 0.0 );
    SET_COORD( up, 0.0, 0.0, 1.0 );
    lib_output_viewpoint(&from, &at, &up, 45.0, 1.0, 0.01, 512, 512);

    /* output background color - UNC sky blue */
    SET_COORD(back_color, 0.078, 0.361, 0.753);
    lib_output_background_color(&back_color);

    /* output light source */
    SET_COORD4(light, -100.0, -100.0, 100.0, 1.0);
    lib_output_light(&light);

    /* set up crystal sphere color - clear white */
    SET_COORD(obj_color, 1.0, 1.0, 1.0 );
    lib_output_color(NULL, &obj_color, 0.0, 0.3, 0.2, 0.7, 5.0, 0.9, 1.5);

    /* output crystal spheres */
    SET_COORD4(center, -0.8, 0.8, 1.00, 0.17);
    create_spheres(&center);

    /* set up mountain color - grey */
    SET_COORD(obj_color, 0.5, 0.45, 0.35);
    lib_output_color(NULL, &obj_color, 0.2, 0.8, 0.0, 0.0, 0.0, 0.0, 1.0);

    /* grow mountain */
    num_pts = 1 << size_factor;
    ratio = 2.0 / exp((double)(log((double)2.0) / (FRACTAL_DIMENSION-1.0)));
    roughness = sqrt((double)(SQR(ratio) - 1.0));
    grow_mountain(num_pts, 0, 0, 0.0, 0.0, 0.0, 0.0);

   /* Make sure everything is cleaned up */
   if (raytracer_format == OUTPUT_RTRACE ||
       raytracer_format == OUTPUT_PLG) {
      lib_set_raytracer(raytracer_format);
      lib_flush_definitions();
      }
}
