/****************************************************************************/
/*                                                                          */
/*  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"
# include "C_raytrace.h"

# define debug(x) 


int C_rt_geo_sphere( world, rt_info, iptr, closest_distance,
		     shade_info )
C_World				*world;
C_RayTraceInfo			*rt_info;
C_IntersectionList		*iptr;
float				*closest_distance;
C_ShadeInfo			*shade_info;
{
	extern C_transform_fposition(); 
	extern C_transform_fvector(); 

	C_Volume		*volume;
	float			dot_product;	
	float			t, distance;
	C_FPosition		new_ray_origin;
	C_FPosition		point_one;
	C_FPosition		intersection_point;
	C_FPosition		junk_point;
	C_FVector		new_ray_direction;
	C_FVector		normal;
	C_FVector		world_normal;
	float			current_intersection;
        float                   diffx, diffy, diffz;

	float			quad_eq_a,
				quad_eq_b,
				quad_eq_c,
				t1, t2,
				radical_term,
				local_distance;

	C_GeoData		*table_entry;


	volume = world->volume[iptr->volume];

	C_transform_fvector( &(rt_info->ray_direction), 
			       &(rt_info->ray_origin), 
			       &(new_ray_direction), &(new_ray_origin),
	 		       &(volume->wtol_voxels) );

	table_entry = volume->data.geo_data;

	quad_eq_b = 2.0*( 
		new_ray_direction.x * ( new_ray_origin.x - 0.5 ) + 
		new_ray_direction.y * ( new_ray_origin.y - 0.5 ) + 
		new_ray_direction.z * ( new_ray_origin.z - 0.5 ) ); 

	quad_eq_c = C_Square( new_ray_origin.x ) + 
		    C_Square( new_ray_origin.y ) + 
		    C_Square( new_ray_origin.z ) + 
		    0.5  - ( new_ray_origin.x  + 
		    new_ray_origin.y  + new_ray_origin.z  );

	radical_term = C_Square( quad_eq_b ) - 
			       4.0 * quad_eq_c;

	local_distance = sqrt( (double) (
		C_Square( iptr->position.x - new_ray_origin.x ) +
		C_Square( iptr->position.y - new_ray_origin.y ) +
		C_Square( iptr->position.z - new_ray_origin.z ) ) ); 

	if ( radical_term == 0.0 )
	{
	    	t = quad_eq_b / -2.0;
	}
	else if ( radical_term > 0.0 )
	{
		t1 = (-quad_eq_b - 
		      sqrt( (double) radical_term ) ) / 2.0;

		t2 = (-quad_eq_b + 
		      sqrt( (double) radical_term ) ) / 2.0;

		if ( t1 > local_distance && (t1 < t2  || t2 < local_distance ) )
		{
			t = t1;
		}
		else if (t2 > local_distance)
		{
			t = t2;
		}
		else
		{
			t = -1.0;
		}
	}
	else t = -1.0;

	if (t > local_distance && ( rt_info->ray_type == C_PRIMARY_RAY ||
		         t > 0.0001 ))
	{
		point_one.x = new_ray_origin.x + new_ray_direction.x * t;
		point_one.y = new_ray_origin.y + new_ray_direction.y * t;
		point_one.z = new_ray_origin.z + new_ray_direction.z * t;

		C_transform_fposition( &point_one,
			&intersection_point, 
			&(volume->ltow_voxels) );

		diffx = intersection_point.x - rt_info->ray_origin.x;
		diffy = intersection_point.y - rt_info->ray_origin.y;
		diffz = intersection_point.z - rt_info->ray_origin.z;

	       	current_intersection = (float) sqrt ( 
	    		C_Square( diffx ) +
	    		C_Square( diffy ) +
			C_Square( diffz ) );

		if ( current_intersection < *closest_distance &&
		     current_intersection > iptr->distance &&
		     current_intersection < iptr->max_distance )
		
		{
		    *closest_distance = current_intersection;

		    normal.x = point_one.x - 0.5;
		    normal.y = point_one.y - 0.5;
		    normal.z = point_one.z - 0.5;

		    C_transform_fvector( 
			    &(normal), &(point_one),
			    &(world_normal), 
			    &(junk_point),
			    &(volume->ltow_voxels));

	  	    shade_info->world_distance = current_intersection;
	  	    shade_info->volume_num = iptr->volume;
	  	    shade_info->world_point.x = intersection_point.x;
	  	    shade_info->world_point.y = intersection_point.y;
	  	    shade_info->world_point.z = intersection_point.z;
	  	    shade_info->local_point.x = point_one.x;
	  	    shade_info->local_point.y = point_one.y;
	  	    shade_info->local_point.z = point_one.z;
		    shade_info->world_normal.x = world_normal.x;
		    shade_info->world_normal.y = world_normal.y;
		    shade_info->world_normal.z = world_normal.z;
		}
	}

}


