/****************************************************************************/
/*                                                                          */
/*  VolVis is a volume visualization system for investigating, manipulating */
/*  and rendering geometric and volumetric data.                            */
/*                                                                          */
/*  Copyright (C) 1993 by the Research Foundation of the State University   */
/*                            of New York                                   */
/*                                                                          */
/*  This program is free software; you can redistribute it and/or modify    */
/*  it under the terms of the GNU General Public License as published by    */
/*  the Free Software Foundation; either version 1, or (at your option)     */
/*  any later version.                                                      */
/*                                                                          */
/*  This program is distributed in the hope that it will be useful,         */
/*  but WITHOUT ANY WARRANTY; without even the implied warranty of          */
/*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           */
/*  GNU General Public License for more details.                            */
/*                                                                          */
/*  You should have received a copy of the GNU General Public License       */
/*  along with this program; if not, write to the Free Software             */
/*  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.               */
/*                                                                          */
/*  For information on VolVis, contact us at:                               */
/*                                                                          */
/*                volvis@cs.sunysb.edu                         (email)      */
/*                                                                          */
/*                Lisa Sobierajski & Ricardo Avila             (US Mail)    */
/*                Department of Computer Science                            */
/*                State University of New York at Stony Brook               */
/*                Stony Brook, New York  11794-4400                         */
/*                                                                          */
/****************************************************************************/



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

# include "C_volvis.h"

extern C_create_plane_equations();


C_identity3D(M)
C_Matrix		*M;
{
        int i, j;

        for (i=0; i < 4; i++)
                for (j=0; j < 4; j++)
                        M->matrix[i][j] = (float) (i == j);

}

C_scale3D(sx, sy, sz, M)
float 			sx, sy, sz;
C_Matrix		*M;
{
        C_identity3D( M );
        M->matrix[0][0] = sx;
        M->matrix[1][1] = sy;
        M->matrix[2][2] = sz;
}


C_translate3D(tx, ty, tz, M)
float 			tx, ty, tz;
C_Matrix		*M;
{
        C_identity3D( M );
        M->matrix[0][3] = tx;
        M->matrix[1][3] = ty;
        M->matrix[2][3] = tz;
}

C_rotate3D(axis, angle, M)
C_AxisType		axis;
float 			angle;
C_Matrix		*M;
{
        int top, bottom, left, right;

        C_identity3D( M );

        top     = (axis == C_X_AXIS ? 1 : 0);
        bottom  = (axis == C_Z_AXIS ? 1 : 2);
        left    = (axis == C_X_AXIS ? 1 : 0);
        right   = (axis == C_Z_AXIS ? 1 : 2);

        M->matrix[left][top]  = cos( (double) angle );
        M->matrix[right][top] = sin( (double) angle );

        if ( axis == C_Y_AXIS )
                M->matrix[right][top] = -(M->matrix[right][top]);

        M->matrix[left][bottom]  = -(M->matrix[right][top]);
        M->matrix[right][bottom] = M->matrix[left][top];

}

/****   PRE = M1 = [M2][M1],  POST = M1 = [M1][M2]   ****/

C_matrix_multiply( type, M1, M2 )

C_MatMultType		type;
C_Matrix		*M1;
C_Matrix		*M2;

{
	C_Matrix	temp;
	int		i, j, k;
	float		value;

	for (i=0; i<4; i++)
	  for (j=0; j<4; j++)
		temp.matrix[i][j] = M1->matrix[i][j];

	if ( type == C_POST_MULTIPLY )
	{
	  for (i=0; i<4; i++)
	    for (j=0; j<4; j++)
	    {
		value = 0.0;
		for (k=0; k<4; k++)
			value += M2->matrix[i][k] * temp.matrix[k][j];
		M1->matrix[i][j] = value;
	    }
	}

	else if ( type == C_PRE_MULTIPLY )
	{
	  for (i=0; i<4; i++)
	    for (j=0; j<4; j++)
	    {
		value = 0.0;
		for (k=0; k<4; k++)
			value += M2->matrix[k][j] * temp.matrix[i][k];
		M1->matrix[i][j] = value;
	    }
	}

	else
	{
		fprintf( stderr, "Invalid Matrix Multiplication!!\n");
		exit( C_ERROR );
	}

}

print_matrix( M )
C_Matrix	M;
{
	printf("\n");
	printf("%f  %f  %f %f\n",M.matrix[0][0], M.matrix[0][1], 
				 M.matrix[0][2], M.matrix[0][3]);
	printf("%f  %f  %f %f\n",M.matrix[1][0], M.matrix[1][1], 
				 M.matrix[1][2], M.matrix[1][3]);
	printf("%f  %f  %f %f\n",M.matrix[2][0], M.matrix[2][1], 
				 M.matrix[2][2], M.matrix[2][3]);
	printf("%f  %f  %f %f\n",M.matrix[3][0], M.matrix[3][1], 
				 M.matrix[3][2], M.matrix[3][3]);
	printf("\n");
}

print_fposition( P )
C_FPosition	P;
{
	printf("\n");
	printf(" %f %f %f\n", P.x, P.y, P.z );
}

C_transform_fposition( in_point, out_point, M )
C_FPosition		*in_point;
C_FPosition		*out_point;
C_Matrix		*M;

{

	out_point->x = in_point->x * M->matrix[0][0] +
		       in_point->y * M->matrix[0][1] +
		       in_point->z * M->matrix[0][2] +
		       		     M->matrix[0][3];
	
	out_point->y = in_point->x * M->matrix[1][0] +
		       in_point->y * M->matrix[1][1] +
		       in_point->z * M->matrix[1][2] +
		       		     M->matrix[1][3];
	
	out_point->z = in_point->x * M->matrix[2][0] +
		       in_point->y * M->matrix[2][1] +
		       in_point->z * M->matrix[2][2] +
		       		     M->matrix[2][3];

}
	
C_transform_hg_fposition( in_point, out_point, M )
C_FPosition		*in_point;
C_FPosition		*out_point;
C_Matrix		*M;

{
	float		w;

	out_point->x = in_point->x * M->matrix[0][0] +
		       in_point->y * M->matrix[0][1] +
		       in_point->z * M->matrix[0][2] +
		       		     M->matrix[0][3];
	
	out_point->y = in_point->x * M->matrix[1][0] +
		       in_point->y * M->matrix[1][1] +
		       in_point->z * M->matrix[1][2] +
		       		     M->matrix[1][3];
	
	out_point->z = in_point->x * M->matrix[2][0] +
		       in_point->y * M->matrix[2][1] +
		       in_point->z * M->matrix[2][2] +
		       		     M->matrix[2][3];

	 w 	     = in_point->x * M->matrix[3][0] +
		       in_point->y * M->matrix[3][1] +
		       in_point->z * M->matrix[3][2] +
		       		     M->matrix[3][3];

	if ( w )
	{
		out_point->x /= w;
		out_point->y /= w;
		out_point->z /= w;
	}
}
	
C_transform_iposition( in_point, out_point, M )
C_IPosition		*in_point;
C_IPosition		*out_point;
C_Matrix		*M;

{

	out_point->x = (int) ( (float)( in_point->x ) * M->matrix[0][0] +
		               (float)( in_point->y ) * M->matrix[0][1] +
		               (float)( in_point->z ) * M->matrix[0][2] +
		       		                        M->matrix[0][3] );
	
	out_point->y = (int) ( (float)( in_point->x ) * M->matrix[1][0] +
		               (float)( in_point->y ) * M->matrix[1][1] +
		               (float)( in_point->z ) * M->matrix[1][2] +
		       		                        M->matrix[1][3] );

	out_point->z = (int) ( (float)( in_point->x ) * M->matrix[2][0] +
		               (float)( in_point->y ) * M->matrix[2][1] +
		               (float)( in_point->z ) * M->matrix[2][2] +
		       		                        M->matrix[2][3] );
}
	
C_transform_fvector( in_vector, in_origin, out_vector, out_origin, M )
C_FVector		*in_vector;
C_FPosition		*in_origin;
C_FVector		*out_vector;
C_FPosition		*out_origin;
C_Matrix		*M;

{
	float		temp;


	out_origin->x = in_origin->x * M->matrix[0][0] +
		        in_origin->y * M->matrix[0][1] +
		        in_origin->z * M->matrix[0][2] +
				       M->matrix[0][3];
	
	out_origin->y = in_origin->x * M->matrix[1][0] +
		        in_origin->y * M->matrix[1][1] +
		        in_origin->z * M->matrix[1][2] +
				       M->matrix[1][3];
	
	out_origin->z = in_origin->x * M->matrix[2][0] +
		        in_origin->y * M->matrix[2][1] +
		        in_origin->z * M->matrix[2][2] +
				       M->matrix[2][3];



	out_vector->x = (in_origin->x + in_vector->x) * M->matrix[0][0] +
		        (in_origin->y + in_vector->y) * M->matrix[0][1] +
		        (in_origin->z + in_vector->z) * M->matrix[0][2] +
				       			M->matrix[0][3] -
				       			out_origin->x;
	
	out_vector->y = (in_origin->x + in_vector->x) * M->matrix[1][0] +
		        (in_origin->y + in_vector->y) * M->matrix[1][1] +
		        (in_origin->z + in_vector->z) * M->matrix[1][2] +
				    		   	M->matrix[1][3] -
				       			out_origin->y;
	
	out_vector->z = (in_origin->x + in_vector->x) * M->matrix[2][0] +
		        (in_origin->y + in_vector->y) * M->matrix[2][1] +
		        (in_origin->z + in_vector->z) * M->matrix[2][2] +
				       			M->matrix[2][3] -
				       			out_origin->z;

	C_Normalize( out_vector->x, out_vector->y, out_vector->z, temp );
}


C_transform_c_sys( in_c_sys, out_c_sys, M )
C_CoordSys		*in_c_sys;
C_CoordSys		*out_c_sys;
C_Matrix		*M;

{

	C_transform_fvector( &(in_c_sys->x_axis),
			     &(in_c_sys->origin),
			     &(out_c_sys->x_axis),
			     &(out_c_sys->origin), M );

	C_transform_fvector( &(in_c_sys->y_axis),
			     &(in_c_sys->origin),
			     &(out_c_sys->y_axis),
			     &(out_c_sys->origin), M );

	C_transform_fvector( &(in_c_sys->z_axis),
			     &(in_c_sys->origin),
			     &(out_c_sys->z_axis),
			     &(out_c_sys->origin), M );

}

	


void C_copy_matrix( M1, M2 )
C_Matrix		*M1, *M2;

{
	int i, j;

	for (i = 0; i < 4; i++)
		for (j = 0; j < 4; j++)
			M2->matrix[i][j] = M1->matrix[i][j];
}


void C_copy_c_sys( c_sys1, c_sys2 )
C_CoordSys	*c_sys1, *c_sys2;
{
	c_sys2->origin.x = c_sys1->origin.x;
	c_sys2->origin.y = c_sys1->origin.y;
	c_sys2->origin.z = c_sys1->origin.z;
	c_sys2->x_axis.x = c_sys1->x_axis.x;
	c_sys2->x_axis.y = c_sys1->x_axis.y;
	c_sys2->x_axis.z = c_sys1->x_axis.z;
	c_sys2->y_axis.x = c_sys1->y_axis.x;
	c_sys2->y_axis.y = c_sys1->y_axis.y;
	c_sys2->y_axis.z = c_sys1->y_axis.z;
	c_sys2->z_axis.x = c_sys1->z_axis.x;
	c_sys2->z_axis.y = c_sys1->z_axis.y;
	c_sys2->z_axis.z = c_sys1->z_axis.z;
}

	  

C_create_volume_matrices( volume )
C_Volume		*volume;

{

	float		theta_one, theta_two, theta_three;
	float		trans_x, trans_y, trans_z;
	float		scale_x, scale_y, scale_z;
	C_Matrix	temp;
	C_FPosition	temp_pos;
	C_FVector	new_y_axis;


	/****    Create the wtol_units matrix    ****/


	C_identity3D( &(volume->wtol_units) );

	C_translate3D(  -(volume->c_sys.origin.x), 
			-(volume->c_sys.origin.y), 
			-(volume->c_sys.origin.z), &temp );

	C_matrix_multiply( C_POST_MULTIPLY, 
			&(volume->wtol_units), &temp );


	theta_one = -1.0*(float)(atan2( (double)volume->c_sys.x_axis.y, 
			         (double)volume->c_sys.x_axis.x ));



	C_rotate3D( C_Z_AXIS, theta_one, &temp );

	C_matrix_multiply( C_POST_MULTIPLY, 
			&(volume->wtol_units), &temp );


	theta_two = (float) asin( (double)(volume->c_sys.x_axis.z) );


	C_rotate3D( C_Y_AXIS, theta_two, &temp );

	C_matrix_multiply( C_POST_MULTIPLY, 
			&(volume->wtol_units), &temp );

	C_transform_fvector( &(volume->c_sys.y_axis), 
			     &(volume->c_sys.origin),
			     &new_y_axis, 
			     &temp_pos,
			     &(volume->wtol_units) );

	theta_three = -1.0*(float)(atan2( (double)(new_y_axis.z), 
					  (double)(new_y_axis.y) ));


	C_rotate3D( C_X_AXIS, theta_three, &temp );

	C_matrix_multiply( C_POST_MULTIPLY, 
			&(volume->wtol_units), &temp );

	

	/****    Create the ltow_units matrix    ****/

	C_identity3D( &(volume->ltow_units) );

	C_translate3D(  volume->c_sys.origin.x, 
			volume->c_sys.origin.y, 
			volume->c_sys.origin.z, &temp );

	C_matrix_multiply( C_PRE_MULTIPLY, 
			&(volume->ltow_units), &temp );

	C_rotate3D( C_Z_AXIS, -theta_one, &temp );

	C_matrix_multiply( C_PRE_MULTIPLY, 
			&(volume->ltow_units), &temp );

	C_rotate3D( C_Y_AXIS, -theta_two, &temp );

	C_matrix_multiply( C_PRE_MULTIPLY, 
			&(volume->ltow_units), &temp );

	C_rotate3D( C_X_AXIS, -theta_three, &temp );

	C_matrix_multiply( C_PRE_MULTIPLY, 
			&(volume->ltow_units), &temp );


	/****    Create the ltow_voxels matrix    ****/


	C_copy_matrix(  &(volume->ltow_units),
			&(volume->ltow_voxels) );

        scale_x = (float)volume->x_size_units/(float)(volume->x_size_voxels-1);
        scale_y = (float)volume->y_size_units/(float)(volume->y_size_voxels-1);
        scale_z = (float)volume->z_size_units/(float)(volume->z_size_voxels-1);

        C_scale3D( scale_x, scale_y, scale_z, &temp );

        C_matrix_multiply( C_PRE_MULTIPLY, &(volume->ltow_voxels), &temp );


	/****    Create the wtol_voxels matrix    ****/

	C_copy_matrix(  &(volume->wtol_units),
			&(volume->wtol_voxels) );

	scale_x = (float)(volume->x_size_voxels-1)/(float)volume->x_size_units;
	scale_y = (float)(volume->y_size_voxels-1)/(float)volume->y_size_units;
	scale_z = (float)(volume->z_size_voxels-1)/(float)volume->z_size_units;

	C_scale3D( scale_x, scale_y, scale_z, &temp );

	C_matrix_multiply( C_POST_MULTIPLY, &(volume->wtol_voxels), &temp );


}



C_create_view_matrices( view )
C_View		*view;

{

	float		theta_one, theta_two, theta_three;
	float		trans_x, trans_y, trans_z;
	float		scale_x, scale_y, scale_z;
	C_Matrix	temp;
	C_FPosition	temp_pos;
	C_FVector	new_y_axis;

	/****	 Set width and height in units if perspective    ****/

	if (view->fov != 0.0)
        {
                view->height_units = 2.0*(float)(tan((double)(view->fov/2.0)));
                view->width_units = view->height_units *
                    (((float)view->width_pixels)/((float)view->height_pixels));
        }

	/****    Create the wtol_units matrix    ****/
	C_identity3D( &(view->wtol_units) );

	C_translate3D(  -(view->c_sys.origin.x), 
			-(view->c_sys.origin.y), 
			-(view->c_sys.origin.z), &temp );

	C_matrix_multiply( C_POST_MULTIPLY, 
			&(view->wtol_units), &temp );

	theta_one = -1.0*(float)(atan2( (double)(view->c_sys.x_axis.y), 
			         	(double)(view->c_sys.x_axis.x)) );

	C_rotate3D( C_Z_AXIS, theta_one, &temp );

	C_matrix_multiply( C_POST_MULTIPLY, 
			&(view->wtol_units), &temp );

	theta_two = (float)(asin( (double)(view->c_sys.x_axis.z) ));

	C_rotate3D( C_Y_AXIS, theta_two, &temp );

	C_matrix_multiply( C_POST_MULTIPLY, 
			&(view->wtol_units), &temp );

	C_transform_fvector( &(view->c_sys.y_axis), 
			     &(view->c_sys.origin),
			     &new_y_axis, 
			     &temp_pos,
			     &(view->wtol_units) );

	theta_three = -1.0*(float)(atan2( (double)(new_y_axis.z), 
					  (double)(new_y_axis.y)) );

	C_rotate3D( C_X_AXIS, theta_three, &temp );

	C_matrix_multiply( C_POST_MULTIPLY, 
			&(view->wtol_units), &temp );

	/****    Create the ltow_units matrix    ****/

	C_identity3D( &(view->ltow_units) );

	C_translate3D(  view->c_sys.origin.x, 
			view->c_sys.origin.y, 
			view->c_sys.origin.z, &temp );

	C_matrix_multiply( C_PRE_MULTIPLY, 
			&(view->ltow_units), &temp );

	C_rotate3D( C_Z_AXIS, -theta_one, &temp );

	C_matrix_multiply( C_PRE_MULTIPLY, 
			&(view->ltow_units), &temp );

	C_rotate3D( C_Y_AXIS, -theta_two, &temp );

	C_matrix_multiply( C_PRE_MULTIPLY, 
			&(view->ltow_units), &temp );

	C_rotate3D( C_X_AXIS, -theta_three, &temp );

	C_matrix_multiply( C_PRE_MULTIPLY, 
			&(view->ltow_units), &temp );

	/****    Create the ltow_pixels matrix    ****/

	C_copy_matrix(  &(view->ltow_units),
			&(view->ltow_pixels) );

        scale_x = (float)view->width_units / view->width_pixels;
        scale_y = (float)view->height_units / view->height_pixels;
        scale_z = 1.0; 

        C_scale3D( scale_x, scale_y, scale_z, &temp );

        C_matrix_multiply( C_PRE_MULTIPLY, &(view->ltow_pixels), &temp );

	/****    Create the wtol_pixels matrix    ****/

	C_copy_matrix(  &(view->wtol_units),
			&(view->wtol_pixels) );

	scale_x = (float)view->width_pixels / view->width_units;
	scale_y = (float)view->height_pixels / view->height_units;
	scale_z = 1.0; 

	C_scale3D( scale_x, scale_y, scale_z, &temp );

	C_matrix_multiply( C_POST_MULTIPLY, &(view->wtol_pixels), &temp );
}

C_rotate_volume_world( volume, axis, angle, reference )
C_Volume		*volume;
C_AxisType		axis;
float			angle;
C_FPosition		*reference;
{
	C_Matrix	mat;

	C_translate3D( -(reference->x), -(reference->y), 
		       -(reference->z), &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	C_rotate3D( axis, angle, &mat );

	
	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );
			
	C_translate3D( reference->x, reference->y, reference->z, &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	C_transform_c_sys( &(volume->orig_c_sys),
			   &(volume->c_sys), &(volume->transform) );

	C_create_volume_matrices( volume );

	C_create_plane_equations( volume );
}


C_rotate_volume_local( volume, axis, angle, reference )
C_Volume		*volume;
C_AxisType		axis;
float			angle;
C_FPosition		*reference;
{
	C_Matrix	mat;

	/* mat = T */
	C_translate3D( (volume->c_sys.origin.x - reference->x),
		       (volume->c_sys.origin.y - reference->y),
		       (volume->c_sys.origin.z - reference->z), &mat );	

	/* transform = orig_transform * T */
	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	/* mat = R */
	C_rotate3D( axis, angle, &mat );

	/* mat = wtol * R */
	C_matrix_multiply( C_PRE_MULTIPLY, &mat, &(volume->wtol_units) ); 	

	/* mat = wtol * R * ltow */
	C_matrix_multiply( C_POST_MULTIPLY, &mat, &(volume->ltow_units) );

	/* transform = orig_transform * T * wtol * R * ltow */
	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	/* mat = -T */
	C_translate3D( (reference->x - volume->c_sys.origin.x),
		       (reference->y - volume->c_sys.origin.y),
		       (reference->z - volume->c_sys.origin.z), &mat );

	/* transform = orig_transform * T * wtol * R * ltow * -T */
	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	C_transform_c_sys( &(volume->orig_c_sys),
			   &(volume->c_sys), &(volume->transform) );

	C_create_volume_matrices( volume );

	C_create_plane_equations( volume );

}

C_rotate_volume_view( volume, view, axis, angle, reference )
C_Volume		*volume;
C_View			*view;
C_AxisType		axis;
float			angle;
C_FPosition		*reference;
{
	C_Matrix	mat;

	/* mat = T */
	C_translate3D( (view->c_sys.origin.x - reference->x),
		       (view->c_sys.origin.y - reference->y),
		       (view->c_sys.origin.z - reference->z), &mat );	

	/* transform = orig_transform * T */
	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	/* mat = R */
	C_rotate3D( axis, angle, &mat );

	/* mat = wtol * R */
	C_matrix_multiply( C_PRE_MULTIPLY, &mat, &(view->wtol_units) ); 	

	/* mat = wtol * R * ltow */
	C_matrix_multiply( C_POST_MULTIPLY, &mat, &(view->ltow_units) );

	/* transform = orig_transform * T * wtol * R * ltow */
	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	/* mat = -T */
	C_translate3D( (reference->x - view->c_sys.origin.x),
		       (reference->y - view->c_sys.origin.y),
		       (reference->z - view->c_sys.origin.z), &mat );

	/* transform = orig_transform * T * wtol * R * ltow * -T */
	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	C_transform_c_sys( &(volume->orig_c_sys),
			   &(volume->c_sys), &(volume->transform) );

	C_create_volume_matrices( volume );

	C_create_plane_equations( volume );
}

C_rotate_volume_volume( volume, csys_volume, axis, angle, reference )
C_Volume		*volume;	/* Rotate This Volume */
C_Volume		*csys_volume;	/* Use This Volumes Coord. Sys. */
C_AxisType		axis;		/* Rotate About This Axis */
float			angle;		/* Rotate By This Angle */
C_FPosition		*reference;	/* Rotate Using This Ref. Point */
{
	C_Matrix	mat;

	/* mat = T */
	C_translate3D( (csys_volume->c_sys.origin.x - reference->x),
		       (csys_volume->c_sys.origin.y - reference->y),
		       (csys_volume->c_sys.origin.z - reference->z), &mat );	

	/* transform = orig_transform * T */
	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	/* mat = R */
	C_rotate3D( axis, angle, &mat );

	/* mat = wtol * R */
	C_matrix_multiply( C_PRE_MULTIPLY, &mat, &(csys_volume->wtol_units) );

	/* mat = wtol * R * ltow */
	C_matrix_multiply( C_POST_MULTIPLY, &mat, &(csys_volume->ltow_units) );

	/* transform = orig_transform * T * wtol * R * ltow */
	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	/* mat = -T */
	C_translate3D( (reference->x - csys_volume->c_sys.origin.x),
		       (reference->y - csys_volume->c_sys.origin.y),
		       (reference->z - csys_volume->c_sys.origin.z), &mat );

	/* transform = orig_transform * T * wtol * R * ltow * -T */
	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );


	C_transform_c_sys( &(volume->orig_c_sys),
			   &(volume->c_sys), &(volume->transform) );

	C_create_volume_matrices( volume );

	C_create_plane_equations( volume );
}

C_translate_volume_local( volume, tx, ty, tz )
C_Volume		*volume;
float			tx;
float			ty;
float			tz;
{
	C_FPosition	wpoint;
	C_FPosition	lpoint;
	C_Matrix	mat;

	lpoint.x = tx;
	lpoint.y = ty;
	lpoint.z = tz;

	C_transform_fposition( &lpoint, &wpoint, &(volume->ltow_units) ); 

	C_translate3D( ( wpoint.x - volume->c_sys.origin.x ),
		       ( wpoint.y - volume->c_sys.origin.y ),
		       ( wpoint.z - volume->c_sys.origin.z ), &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	C_transform_c_sys( &(volume->orig_c_sys),
			   &(volume->c_sys), &(volume->transform) );

	C_create_volume_matrices( volume );

	C_create_plane_equations( volume );
}

C_translate_volume_world( volume, tx, ty, tz )
C_Volume		*volume;
float			tx;
float			ty;
float			tz;
{
	C_Matrix	mat;

	C_translate3D( tx, ty, tz, &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	C_transform_c_sys( &(volume->orig_c_sys),
			   &(volume->c_sys), &(volume->transform) );

	C_create_volume_matrices( volume );

	C_create_plane_equations( volume );

}

C_translate_volume_view( volume, view, tx, ty, tz )
C_Volume		*volume;
C_View			*view;
float			tx;
float			ty;
float			tz;
{
	C_FPosition	wpoint;
	C_FPosition	lpoint;
	C_Matrix	mat;

	lpoint.x = tx;
	lpoint.y = ty;
	lpoint.z = tz;

	C_transform_fposition( &lpoint, &wpoint, &(view->ltow_units) ); 

	C_translate3D( ( wpoint.x - view->c_sys.origin.x ),
		       ( wpoint.y - view->c_sys.origin.y ),
		       ( wpoint.z - view->c_sys.origin.z ), &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	C_transform_c_sys( &(volume->orig_c_sys),
			   &(volume->c_sys), &(volume->transform) );

	C_create_volume_matrices( volume );

	C_create_plane_equations( volume );
}

C_translate_volume_volume( volume, csys_volume, tx, ty, tz )
C_Volume		*volume;
C_Volume		*csys_volume;
float			tx;
float			ty;
float			tz;
{
	C_FPosition	wpoint;
	C_FPosition	lpoint;
	C_Matrix	mat;

	lpoint.x = tx;
	lpoint.y = ty;
	lpoint.z = tz;

	C_transform_fposition( &lpoint, &wpoint, &(csys_volume->ltow_units) ); 

	C_translate3D( ( wpoint.x - csys_volume->c_sys.origin.x ),
		       ( wpoint.y - csys_volume->c_sys.origin.y ),
		       ( wpoint.z - csys_volume->c_sys.origin.z ), &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );

	C_transform_c_sys( &(volume->orig_c_sys),
			   &(volume->c_sys), &(volume->transform) );

	C_create_volume_matrices( volume );

	C_create_plane_equations( volume );
}

C_scale_volume_local( volume, sx, sy, sz, reference )
C_Volume		*volume;
float			sx;
float			sy;
float			sz;
C_FPosition		*reference;
{
	C_Matrix	mat;
	C_FPosition	wref;
	C_FPosition	new_reference;
	C_FPosition	loc_ref;

	if ( sx <= 0.0 || sy <= 0.0 || sz <= 0.0 )
	{
		C_error_message("Invalid Scale Factor!!!\n");
		return;
	}

	volume->x_size_units *= sx;
	volume->y_size_units *= sy;
	volume->z_size_units *= sz;

	C_transform_fposition( reference, &loc_ref, &(volume->wtol_units) );

	loc_ref.x *= sx;
	loc_ref.y *= sy;
	loc_ref.z *= sz;
	

	C_transform_fposition(&loc_ref, &new_reference, &(volume->ltow_units));

	C_translate3D( (reference->x - new_reference.x), 
		       (reference->y - new_reference.y),
		       (reference->z - new_reference.z), &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(volume->transform), &mat );
	
	C_transform_c_sys( &(volume->orig_c_sys),
			   &(volume->c_sys), &(volume->transform) );


	C_create_volume_matrices( volume );

	C_create_plane_equations( volume );

}


C_rotate_view_world( view, axis, angle, reference )
C_View			*view;
C_AxisType		axis;
float			angle;
C_FPosition		*reference;
{
	C_Matrix	mat;


	C_translate3D( -(reference->x), -(reference->y), 
		       -(reference->z), &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );

	C_rotate3D( axis, angle, &mat );

	
	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );
			
	C_translate3D( reference->x, reference->y, 
		       reference->z, &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );

	C_transform_c_sys( &(view->orig_c_sys),
			   &(view->c_sys), &(view->transform) );

	C_create_view_matrices( view );

}

C_rotate_view_local( view, axis, angle, reference )
C_View			*view;
C_AxisType		axis;
float			angle;
C_FPosition		*reference;
{
	C_Matrix	mat;

	C_translate3D( (view->c_sys.origin.x - reference->x),
		       (view->c_sys.origin.y - reference->y),
		       (view->c_sys.origin.z - reference->z), &mat );	

	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );

	C_rotate3D( axis, angle, &mat );

	C_matrix_multiply( C_PRE_MULTIPLY, &mat, &(view->wtol_units) ); 	

	C_matrix_multiply( C_POST_MULTIPLY, &mat, &(view->ltow_units) );

	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );

	C_translate3D( (reference->x - view->c_sys.origin.x),
		       (reference->y - view->c_sys.origin.y),
		       (reference->z - view->c_sys.origin.z), &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );

	C_transform_c_sys( &(view->orig_c_sys),
			   &(view->c_sys), &(view->transform) );

	C_create_view_matrices( view );
}

C_rotate_view_volume( view, csys_volume, axis, angle, reference )
C_View			*view;
C_Volume		*csys_volume;
C_AxisType		axis;
float			angle;
C_FPosition		*reference;
{
	C_Matrix	mat;

	C_translate3D( (csys_volume->c_sys.origin.x - reference->x),
		       (csys_volume->c_sys.origin.y - reference->y),
		       (csys_volume->c_sys.origin.z - reference->z), &mat );	

	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );

	C_rotate3D( axis, angle, &mat );

	C_matrix_multiply(C_PRE_MULTIPLY, &mat, &(csys_volume->wtol_units) ); 	

	C_matrix_multiply(C_POST_MULTIPLY, &mat,&(csys_volume->ltow_units));

	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );

	C_translate3D( (reference->x - csys_volume->c_sys.origin.x),
		       (reference->y - csys_volume->c_sys.origin.y),
		       (reference->z - csys_volume->c_sys.origin.z), &mat );

	C_transform_c_sys( &(view->orig_c_sys),
			   &(view->c_sys), &(view->transform) );

	C_create_view_matrices( view );
}

C_translate_view_local( view, tx, ty, tz )
C_View			*view;
float			tx;
float			ty;
float			tz;
{
	C_FPosition	wpoint;
	C_FPosition	lpoint;
	C_Matrix	mat;


	lpoint.x = tx;
	lpoint.y = ty;
	lpoint.z = tz;

	C_transform_fposition( &lpoint, &wpoint, &(view->ltow_units) ); 


	C_translate3D( ( wpoint.x - view->c_sys.origin.x ),
		       ( wpoint.y - view->c_sys.origin.y ),
		       ( wpoint.z - view->c_sys.origin.z ), &mat );


	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );


	C_transform_c_sys( &(view->orig_c_sys),
			   &(view->c_sys), &(view->transform) );

	C_create_view_matrices( view );
}
	
C_translate_view_world( view, tx, ty, tz )
C_Volume		*view;
float			tx;
float			ty;
float			tz;
{
	C_Matrix	mat;


	C_translate3D( tx, ty, tz, &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );

	C_transform_c_sys( &(view->orig_c_sys),
			   &(view->c_sys), &(view->transform) );

	C_create_view_matrices( view );
}

C_translate_view_volume( view, csys_volume, tx, ty, tz )
C_View			*view;
C_Volume		*csys_volume;
float			tx;
float			ty;
float			tz;
{
	C_FPosition	wpoint;
	C_FPosition	lpoint;
	C_Matrix	mat;


	lpoint.x = tx;
	lpoint.y = ty;
	lpoint.z = tz;

	C_transform_fposition( &lpoint, &wpoint, &(csys_volume->ltow_units) ); 


	C_translate3D( ( wpoint.x - csys_volume->c_sys.origin.x ),
		       ( wpoint.y - csys_volume->c_sys.origin.y ),
		       ( wpoint.z - csys_volume->c_sys.origin.z ), &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );


	C_transform_c_sys( &(view->orig_c_sys),
			   &(view->c_sys), &(view->transform) );

	C_create_view_matrices( view );
}

C_scale_view_local( view, sx, sy, sz, reference )
C_View			*view;
float			sx;
float			sy;
float			sz;
C_FPosition		*reference;

{
	C_Matrix	mat;
	C_FPosition	wref;
	C_FPosition	new_reference;
	C_FPosition	loc_ref;


	view->width_units *= sx;
	view->height_units *= sy;

	C_transform_fposition( reference, &loc_ref, &(view->wtol_units) );

	loc_ref.x *= sx;
	loc_ref.y *= sy;
	loc_ref.z *= sz;
	

	C_transform_fposition( &loc_ref, &new_reference, &(view->ltow_units) );

	C_translate3D( (reference->x - new_reference.x), 
		       (reference->y - new_reference.y),
		       (reference->z - new_reference.z), &mat );

	C_matrix_multiply( C_POST_MULTIPLY, &(view->transform), &mat );
	
	C_transform_c_sys( &(view->orig_c_sys),
			   &(view->c_sys), &(view->transform) );


	C_create_view_matrices( view );

}

C_rotate_light_world( light, axis, angle, reference )
C_Light			*light;
C_AxisType		axis;
float			angle;
C_FPosition		*reference;
{
	C_Matrix	mat;


	switch ( light->light_type )
	{

	  case C_POINT_LIGHT:


		C_translate3D( -(reference->x), -(reference->y), 
		       		-(reference->z), &mat );

		C_matrix_multiply( C_POST_MULTIPLY, 
			   	&(light->light.light_point->transform), &mat );

		C_rotate3D( axis, angle, &mat );
		
		C_matrix_multiply( C_POST_MULTIPLY, 
		  		&(light->light.light_point->transform), &mat );
			
		C_translate3D( reference->x, reference->y, 
		       	       reference->z, &mat );

		C_matrix_multiply( C_POST_MULTIPLY, 
			   	&(light->light.light_point->transform), &mat );

		C_transform_fposition(
				&(light->light.light_point->orig_light_pos),
			      	&(light->light.light_point->light_pos), 
			      	&(light->light.light_point->transform) );

		break;

	  default:

		C_error_message("Internal Error: MAT_ROTLIGHT1");
	}

}

C_rotate_light_local( light, axis, angle, reference )
C_Light			*light;
C_AxisType		axis;
float			angle;
C_FPosition		*reference;
{
	C_Matrix	mat;


	switch ( light->light_type )
	{

	  case C_POINT_LIGHT:

		C_translate3D( -(reference->x), -(reference->y), 
		       		-(reference->z), &mat );

		C_matrix_multiply( C_POST_MULTIPLY, 
			   	&(light->light.light_point->transform), &mat );

		C_rotate3D( axis, angle, &mat );
		
		C_matrix_multiply( C_POST_MULTIPLY, 
		  		&(light->light.light_point->transform), &mat );
			
		C_translate3D( reference->x, reference->y, reference->z, &mat );

		C_matrix_multiply( C_POST_MULTIPLY, 
			   	&(light->light.light_point->transform), &mat );

		C_transform_fposition(
				&(light->light.light_point->orig_light_pos),
			      	&(light->light.light_point->light_pos), 
			      	&(light->light.light_point->transform) );
		break;

	  default:

		C_error_message("Internal Error: MAT_ROTLIGHT2");
	}
}

C_rotate_light_view( light, view, axis, angle, reference )
C_Light			*light;
C_View			*view;
C_AxisType		axis;
float			angle;
C_FPosition		*reference;
{
	C_Matrix	mat;

	switch ( light->light_type )
	{
	  case C_POINT_LIGHT:

		C_translate3D(  (view->c_sys.origin.x - reference->x),
				(view->c_sys.origin.y - reference->y),
				(view->c_sys.origin.z - reference->z), &mat );	

		C_matrix_multiply( C_POST_MULTIPLY, 
			   	&(light->light.light_point->transform), &mat );

		C_rotate3D( axis, angle, &mat );

		C_matrix_multiply( C_PRE_MULTIPLY, &mat, &(view->wtol_units) );

		C_matrix_multiply( C_POST_MULTIPLY, &mat, &(view->ltow_units) );
		
		C_matrix_multiply( C_POST_MULTIPLY, 
		  		&(light->light.light_point->transform), &mat );
			
		C_translate3D(  (reference->x - view->c_sys.origin.x),
				(reference->y - view->c_sys.origin.y),
				(reference->z - view->c_sys.origin.z), &mat );	

		C_matrix_multiply( C_POST_MULTIPLY, 
			   	&(light->light.light_point->transform), &mat );

		C_transform_fposition(
				&(light->light.light_point->orig_light_pos),
			      	&(light->light.light_point->light_pos), 
			      	&(light->light.light_point->transform) );
		break;

	  default:

		C_error_message("Internal Error: MAT_ROTLIGHT3");
	}

}

C_rotate_light_volume( light, csys_volume, axis, angle, reference )
C_Light			*light;
C_Volume		*csys_volume;
C_AxisType		axis;
float			angle;
C_FPosition		*reference;
{
	C_Matrix	mat;

	switch ( light->light_type )
	{

	  case C_POINT_LIGHT:

		C_translate3D(  (csys_volume->c_sys.origin.x - reference->x),
				(csys_volume->c_sys.origin.y - reference->y),
				(csys_volume->c_sys.origin.z - reference->z), 
				&mat );

		C_matrix_multiply( C_POST_MULTIPLY, 
			   	&(light->light.light_point->transform), &mat );

		C_rotate3D( axis, angle, &mat );

		C_matrix_multiply( C_PRE_MULTIPLY, &mat, 
					&(csys_volume->wtol_units) );

		C_matrix_multiply( C_POST_MULTIPLY, &mat, 
					&(csys_volume->ltow_units) );
		
		C_matrix_multiply( C_POST_MULTIPLY, 
		  		&(light->light.light_point->transform), &mat );
			
		C_translate3D(  (reference->x - csys_volume->c_sys.origin.x),
				(reference->y - csys_volume->c_sys.origin.y),
				(reference->z - csys_volume->c_sys.origin.z), 
				&mat );

		C_matrix_multiply( C_POST_MULTIPLY, 
			   	&(light->light.light_point->transform), &mat );

		C_transform_fposition(
				&(light->light.light_point->orig_light_pos),
			      	&(light->light.light_point->light_pos), 
			      	&(light->light.light_point->transform) );
		break;

	  default:

		C_error_message("Internal Error: MAT_ROTLIGHT4");
	}
}

C_translate_light_world( light, tx, ty, tz )
C_Light			*light;
float			tx;
float			ty;
float			tz;
{
	C_Matrix	mat;

	switch ( light->light_type )
	{

	  case C_POINT_LIGHT:

		C_translate3D( tx, ty, tz, &mat );

		C_matrix_multiply( C_POST_MULTIPLY, 
				   &(light->light.light_point->transform), 
				   &mat );

		C_transform_fposition( 
				&(light->light.light_point->orig_light_pos),
			   	&(light->light.light_point->light_pos), 
				&(light->light.light_point->transform) );
		break;

	  default:

		C_error_message("Internal Error: MAT_TRANSLIGHT1");
	}
}

C_translate_light_local( light, tx, ty, tz )
C_Light			*light;
float			tx;
float			ty;
float			tz;
{
	C_Matrix	mat;

	switch ( light->light_type )
	{

	  case C_POINT_LIGHT:

		C_translate3D( tx, ty, tz, &mat );

		C_matrix_multiply( C_POST_MULTIPLY, 
				   &(light->light.light_point->transform), 
				   &mat );

		C_transform_fposition( 
				&(light->light.light_point->orig_light_pos),
			   	&(light->light.light_point->light_pos), 
				&(light->light.light_point->transform) );
		break;

	  default:

		C_error_message("Internal Error: MAT_TRANSLIGHT2");
	}
}

C_translate_light_view( light, view, tx, ty, tz )
C_Light			*light;
C_View			*view;
float			tx;
float			ty;
float			tz;
{
	C_FPosition	wpoint;
	C_FPosition	lpoint;
	C_Matrix	mat;

	switch ( light->light_type )
	{

	  case C_POINT_LIGHT:

		lpoint.x = tx;
		lpoint.y = ty;
		lpoint.z = tz;

		C_transform_fposition( &lpoint, &wpoint, &(view->ltow_units) ); 

		C_translate3D(  ( wpoint.x - view->c_sys.origin.x ),
		       		( wpoint.y - view->c_sys.origin.y ),
		       		( wpoint.z - view->c_sys.origin.z ), &mat );

		C_matrix_multiply( C_POST_MULTIPLY, 
				   &(light->light.light_point->transform), 
				   &mat );

		C_transform_fposition( 
				&(light->light.light_point->orig_light_pos),
			   	&(light->light.light_point->light_pos), 
				&(light->light.light_point->transform) );
		break;

	  default:

		C_error_message("Internal Error: MAT_TRANSLIGHT3");
	}
}

C_translate_light_volume( light, csys_volume, tx, ty, tz )
C_Light			*light;
C_Volume		*csys_volume;
float			tx;
float			ty;
float			tz;
{
	C_FPosition	wpoint;
	C_FPosition	lpoint;
	C_Matrix	mat;

	switch ( light->light_type )
	{

	  case C_POINT_LIGHT:

		lpoint.x = tx;
		lpoint.y = ty;
		lpoint.z = tz;

		C_transform_fposition(  &lpoint, &wpoint, 
					&(csys_volume->ltow_units) ); 

		C_translate3D(  ( wpoint.x - csys_volume->c_sys.origin.x ),
		       		( wpoint.y - csys_volume->c_sys.origin.y ),
		       		( wpoint.z - csys_volume->c_sys.origin.z ), 
				&mat );

		C_matrix_multiply( C_POST_MULTIPLY, 
				   &(light->light.light_point->transform), 
				   &mat );

		C_transform_fposition( 
				&(light->light.light_point->orig_light_pos),
			   	&(light->light.light_point->light_pos), 
				&(light->light.light_point->transform) );
		break;

	  default:

		C_error_message("Internal Error: MAT_TRANSLIGHT4");
	}
}

void	C_null3D( M )
C_Matrix	*M;
{
	int	i, j;

	for ( i = 0; i < 4; i++)
	{
		for ( j = 0; j < 4; j++)
		{
			M->matrix[i][j] = 0.0;
		}
	}

	M->matrix[3][3] = 1.0;
}

void	C_make_center_display_matrix( dx, dy, M)
float	dx,
	dy;
C_Matrix	*M;
{
	C_identity3D( M);

	M->matrix[0][3] += dx;
	M->matrix[1][3] += dy;
}

void	C_make_perspective_matrix( view, M, near, far )
C_View		*view;
C_Matrix	*M;
float		far,
		near;
{
	float	theta,
		f_height_pixels,
		f_width_pixels;

	f_height_pixels = (float) (view->height_pixels);
	f_width_pixels  = (float) (view->width_pixels);

	C_identity3D( M);

	theta = view->fov /2.00;	

	theta = 1.00/ tan( theta);		/* gets the cotangent */

	M->matrix[0][0] = theta/( f_width_pixels/ f_height_pixels);

	M->matrix[1][1] = theta;

	M->matrix[2][2] = ( far + near )/(far - near );

	M->matrix[2][3] = -1.0 * (2.0 * far * near)/(far - near);

	M->matrix[3][2] = 1.00;

	M->matrix[3][3] = 0.00;

}	

void	C_transform_point( M1, P)
C_Matrix	*M1,
		*P;
{
	C_Matrix	temp;
	
	C_copy_matrix( P, &temp);

	P->matrix[0][0] = temp.matrix[0][0] * M1->matrix[0][0] +
				temp.matrix[1][1] * M1->matrix[0][1] +
					temp.matrix[2][2] * M1->matrix[0][2];

	P->matrix[1][1] = temp.matrix[0][0] * M1->matrix[1][0] +
				temp.matrix[1][1] * M1->matrix[1][1] +
					temp.matrix[2][2] * M1->matrix[1][2];

	P->matrix[2][2] = temp.matrix[0][0] * M1->matrix[2][0] +
				temp.matrix[1][1] * M1->matrix[2][1] +
					temp.matrix[2][2] * M1->matrix[2][2];

}
