ivp_buoyancy_solver.cxx

hl2 source code. Do not use it illegal.

		    ball_projected_surface_content_under = surface_content_down;
		    goto calculate_impulse; //already finished here
				
		} else {
		    center_of_segment_os.set(geom_center_os);
		    center_of_segment_os.add_multiple(&direction_down, length_p1p2*length_p1p2*length_p1p2 / (12 * surface_content_down));
		}
	    }
			
			
	    /*** compute part of ellipse above surface ***/
			
			
	    //compute center of whole ellipse
	    IVP_U_Float_Point projected_center_on_surface_os;
	    surface_os->proj_on_plane(geom_center_os, &projected_center_on_surface_os);
	    IVP_U_Float_Point center_of_ellipse_os;
	    speed_plane_os.proj_on_plane(&projected_center_on_surface_os, &center_of_ellipse_os);
			
	    IVP_U_Float_Point edge_O_center_os;
	    edge_O_center_os.subtract(&center_of_ellipse_os, &O_os);
	    IVP_DOUBLE b = edge_O_center_os.real_length();
	    IVP_U_Float_Point edge_Q_center_os;
	    edge_Q_center_os.subtract(&center_of_ellipse_os, &Q_os);
	    IVP_DOUBLE part_b = edge_Q_center_os.real_length();	
	    IVP_DOUBLE radius_of_water_disection = IVP_Inline_Math::ivp_sqrtf(IVP_Inline_Math::fabsd(radius*radius - distance*distance));

	    IVP_DOUBLE surface_content_of_cut_off_ellipse_segment;
	    IVP_U_Float_Point center_of_part_ellipse_os;
	    IVP_DOUBLE surface_content_up;
	    if (b < buoyancy_eps) {  //safety check
		//if this is the case, then speed_plane_os is nearly orthogonal to surface_os, so this part of the dampened surface can be assumed 0
		surface_content_of_cut_off_ellipse_segment = 0;
		surface_content_up = 0.0f; // OS @@@@ ???
		center_of_part_ellipse_os.set(0.0f, 0.0f, 0.0f);
	    } else {

		surface_content_of_cut_off_ellipse_segment = IVP_Inline_Math::fabsd((radius_of_water_disection*b*IVP_Inline_Math::save_acosf(part_b / b) - part_b*(p1p2_os.real_length()*0.5f)));  //segment of ellipse under surface_os
	    
		IVP_DOUBLE surface_content_of_whole_ellipse = (IVP_PI*radius_of_water_disection*b);
		surface_content_up = IVP_Inline_Math::fabsd(surface_content_of_whole_ellipse - surface_content_of_cut_off_ellipse_segment);
			
		//compute the center of the part of the ellipse above surface
		IVP_DOUBLE h = b - part_b;
		IVP_DOUBLE s = 2.0f * IVP_Inline_Math::ivp_sqrtf(IVP_Inline_Math::fabsd(2*h*b - h*h));
		IVP_DOUBLE arc = IVP_Inline_Math::fabsd(2*IVP_Inline_Math::fast_asin( s/(2.0f * b) )*b);
		IVP_DOUBLE A = IVP_Inline_Math::fabsd(0.5f * (arc * b - s*(b - h)));

		IVP_U_Float_Point center_of_cut_off_part_of_ellipse;  //center of small cut-off part of ellipse under surface
		if (A < buoyancy_eps) { //safety check
		    center_of_cut_off_part_of_ellipse.set(0.0f, 0.0f, 0.0f);
		    surface_content_of_cut_off_ellipse_segment = 0.0f;
		} else {
		    edge_O_center_os.normize();  //savety check above (for b)
		    center_of_cut_off_part_of_ellipse.set(&center_of_ellipse_os);
		    center_of_cut_off_part_of_ellipse.add_multiple(&edge_O_center_os, s*s*s / (12 * A));
		}

		if (surface_content_up < buoyancy_eps) { //safety check
		    center_of_part_ellipse_os.set(0.0f, 0.0f, 0.0f);
		    surface_content_up = 0.0f;
		} else {
		    center_of_part_ellipse_os.set_multiple(&center_of_cut_off_part_of_ellipse, -surface_content_of_cut_off_ellipse_segment);
		    center_of_part_ellipse_os.add_multiple(&center_of_ellipse_os, surface_content_of_whole_ellipse);
		    center_of_part_ellipse_os.mult( 1.0f / surface_content_up );
		}
	    }
			
	    /*** compute results ***/
	    ball_projected_surface_content_under = surface_content_up + surface_content_down;

	    if (ball_projected_surface_content_under < buoyancy_eps) {
		ball_projected_surface_content_under = 0.0f;
		//ball_point_of_impulse is still set to (0,0,0)
	    } else {
		ball_point_of_impulse.set_multiple(&center_of_segment_os, surface_content_down);
		ball_point_of_impulse.add_multiple(&center_of_part_ellipse_os, surface_content_up);
		ball_point_of_impulse.mult(1.0f / ball_projected_surface_content_under);
	    }
			
	} else { //ball flies out of the water
			
	    /*** calc ball_projected_surface_content_under ***/
			
	    //compute center of whole ellipse
	    IVP_U_Float_Point projected_center_on_surface_os;
	    surface_os->proj_on_plane(geom_center_os, &projected_center_on_surface_os);
	    IVP_U_Float_Point center_of_ellipse_os;
	    speed_plane_os.proj_on_plane(&projected_center_on_surface_os, &center_of_ellipse_os);
			
	    //calc various ellipse related values
	    IVP_U_Float_Point edge_O_center_os;
	    edge_O_center_os.subtract(&center_of_ellipse_os, &O_os);
	    IVP_DOUBLE b = edge_O_center_os.real_length();
	    IVP_U_Float_Point edge_Q_center_os;
	    edge_Q_center_os.subtract(&center_of_ellipse_os, &Q_os);
	    IVP_DOUBLE radius_of_water_disection = IVP_Inline_Math::ivp_sqrtf(IVP_Inline_Math::fabsd(radius*radius - distance*distance));
	    IVP_DOUBLE height = IVP_Inline_Math::fabsd(b - edge_Q_center_os.real_length());
	    //IVP_DOUBLE length_p1p2 = p1p2_os.real_length();
			
	    //calc circle segment under water
	    IVP_U_Float_Point MQ_os;
	    MQ_os.subtract(&Q_os, geom_center_os);
	    IVP_DOUBLE h = radius - MQ_os.real_length();  //height of disected segment of circle
	    IVP_DOUBLE arc = 2*IVP_Inline_Math::fast_asin(length_p1p2 / (2.0f * radius)) * radius;  //length of arc of circle segment
	    IVP_DOUBLE surface_content_of_circle_segment = IVP_Inline_Math::fabsd(0.5f * (arc*radius - length_p1p2*(radius - h)));

	    if (surface_content_of_circle_segment < buoyancy_eps) {  //safety check
		//if this is the case, then there is near to nothing to damp
		ball_projected_surface_content_under = 0.0f;
		goto calculate_impulse;
	    }
			
	    //calc segment of ellipse under water
	    IVP_DOUBLE surface_content_of_ellipse_segment;
	    if (b < buoyancy_eps) { //safety check
		//if this is the case then surface_os and speed_plane_os are nearly orthogonal, so the projected ellipse has nearly vanished
		surface_content_of_ellipse_segment = 0.0f;
	    } else {
		surface_content_of_ellipse_segment = b*radius_of_water_disection *
		    IVP_Inline_Math::save_acosf(edge_Q_center_os.real_length() / b) - 
		    edge_Q_center_os.real_length()*(length_p1p2 * 0.5f);
	    }
			
	    ball_projected_surface_content_under = surface_content_of_circle_segment - surface_content_of_ellipse_segment;
	    if (ball_projected_surface_content_under < buoyancy_eps) {  //safety check
		ball_projected_surface_content_under = 0.0f;
		goto calculate_impulse;
	    }
			
	    /*** calc ball_point_of_impulse ***/
			
	    //IVP_U_Float_Point direction;
	    //direction.set_multiple(&edge_O_center_os, -1.0f);  //turn direction
	    //direction.normize();

	    //calc mass center of circle segment
	    IVP_U_Float_Point center_of_circle_segment;
	    center_of_circle_segment.set(geom_center_os);
	    center_of_circle_segment.add_multiple(&direction_down, length_p1p2*length_p1p2*length_p1p2 / (12*surface_content_of_circle_segment)); //safety check for surface_content_of_circle_segment already above
			
	    //calc mass center of ellipse segment
	    IVP_U_Float_Point center_of_ellipse_segment;
	    IVP_DOUBLE s = 2.0f * IVP_Inline_Math::ivp_sqrtf(IVP_Inline_Math::fabsd(2.0f * height * b - height*height));
	    IVP_DOUBLE ellipse_arc = 2*IVP_Inline_Math::fast_asin(s / (2.0f*b))*b;
	    IVP_DOUBLE A = 0.5f * (ellipse_arc * b - s*(b - height));
	    if ( (b<buoyancy_eps) || (A<buoyancy_eps) ) { //safety check
		//if this is the case then surface_os and speed_plane_os are nearly orthogonal, so the projected ellipse has nearly vanished
		center_of_ellipse_segment.set(0.0f, 0.0f, 0.0f);
	    } else {
		center_of_ellipse_segment.set(&center_of_ellipse_os);
		center_of_ellipse_segment.add_multiple(&direction_down, s*s*s / (12.0f * A));
	    }

	    ball_point_of_impulse.set_multiple(&center_of_ellipse_segment, -surface_content_of_ellipse_segment);
	    ball_point_of_impulse.add_multiple(&center_of_circle_segment, surface_content_of_circle_segment);
	    ball_point_of_impulse.mult(1.0f / ball_projected_surface_content_under);
			
	}
		
    }
    break;
    case 2: { //center lies above the surface and speed_plane_os disects ball only above the surface
    case_2:
	//check in which direction the ball is heading
	IVP_DOUBLE factor = surface_os->dot_product(&speed_plane_os);
	if (factor < 0.0f) {
			
	    //compute ball_point_of_impulse
	    IVP_U_Float_Point projected_center_on_surface_os;
	    surface_os->proj_on_plane(geom_center_os, &projected_center_on_surface_os);
	    speed_plane_os.proj_on_plane(&projected_center_on_surface_os, &ball_point_of_impulse);
			
	    //compute ball_projected_surface_content_under
	    ball_projected_surface_content_under = IVP_Inline_Math::fabsd(factor*IVP_PI*(radius*radius - distance*distance));
			
	} else {
			
	    //no dampening needed
	    ball_projected_surface_content_under = 0.0f;
	    //ball_point_of_impulse.set(0.0f, 0.0f, 0.0f);
			
	}
    }
    break;
    case 3: { //center lies under the surface and speed_plane_os disects ball under and above the surface

	//check if surface_os and speed_plane_os are nearly parallel
	//if so, the we can switch to case 4
	IVP_U_Float_Point direction_up;
	speed_plane_os.proj_on_plane(surface_os, &direction_up);
	direction_up.mult(-1.0f);
	if (direction_up.quad_length() < buoyancy_eps) {
	    //if this is the case, then surface_os and speed_plane_os are nearly parallel, so we can switch to case 4
	    goto case_4;
	}
	direction_up.normize();
		
	IVP_U_Float_Point p1p2_os;
	p1p2_os.subtract(p2_os, p1_os);
		
	//calc point in the middle of the straight that results from the intersection between surface_os and speed_plane_os
	IVP_U_Float_Point Q_os;  //point between p1_os and p2_os
	Q_os.set(p1_os);
	Q_os.add_multiple(&p1p2_os, 0.5f);

	IVP_DOUBLE length_p1p2 = p1p2_os.real_length();
	if (length_p1p2 < buoyancy_eps) {  //check for safety reasons
	    //in this case Q lies on the surface of the ball
	    ball_projected_surface_content_under = IVP_PI * radius * radius;
	    ball_point_of_impulse.set(geom_center_os);
	    break;
	}
	//calc projected highest point O_os on arc of ellipse
	IVP_U_Float_Hesse p1p2_hesse_os;
	p1p2_hesse_os.set(&p1p2_os);
	p1p2_hesse_os.calc_hesse_val(&Q_os);
	p1p2_hesse_os.normize();
		
	IVP_U_Float_Point tmp_O1_os;
	IVP_U_Float_Point tmp_O2_os;
	//int erg =
	compute_disection_points_with_ball(surface_os, &p1p2_hesse_os, geom_center_os, radius, &tmp_O1_os, &tmp_O2_os);
	//	    IVP_ASSERT( erg == 1 );
	//decide which point is the right one and project it onto speed_plane_os => point O_os
	IVP_U_Float_Point tmp_O1O2_os;
	tmp_O1O2_os.subtract(&tmp_O2_os, &tmp_O1_os);
	IVP_U_Float_Point O_os;
	if (speed_plane_os.dot_product(&tmp_O1O2_os) >= 0.0f) {
	    speed_plane_os.proj_on_plane(&tmp_O1_os, &O_os);
	} else {
	    speed_plane_os.proj_on_plane(&tmp_O2_os, &O_os);
	}
		
	//compute center of whole ellipse
	IVP_U_Float_Point projected_center_on_surface_os;
	surface_os->proj_on_plane(geom_center_os, &projected_center_on_surface_os);
	IVP_U_Float_Point center_of_ellipse_os;
	speed_plane_os.proj_on_plane(&projected_center_on_surface_os, &center_of_ellipse_os);
		
		
	//check in which direction the ball is heading
	IVP_DOUBLE factor = surface_os->dot_product(&speed_plane_os);
	if (factor <= buoyancy_eps) {  //ball falls into the medium
			
	    /*** compute circle segment under surface ***/
			
	    IVP_U_Float_Point MQ_os;
	    MQ_os.subtract(&Q_os, geom_center_os);
	    IVP_DOUBLE h = radius - MQ_os.real_length();  //height of disected part of circle
	    IVP_DOUBLE arc = 2*IVP_Inline_Math::fast_asin(length_p1p2 / (2.0f * radius)) * radius;  //length of arc which is cut off speed_plane_os by surface_os
			
	    IVP_DOUBLE cut_off_part_above_surface = IVP_Inline_Math::fabsd(0.5f * (arc*radius - length_p1p2*(radius - h)));
	    IVP_DOUBLE surface_content_down = IVP_Inline_Math::fabsd(IVP_PI*radius*radius - cut_off_part_above_surface);

	    IVP_ASSERT(surface_content_down > buoyancy_eps);
			
	    IVP_U_Float_Point center_of_segment_os;  //center of part under surface

	    if (cut_off_part_above_surface < buoyancy_eps) { //safety check
		//nearly the whole surface to damp is in the medium
		ball_projected_surface_content_under = IVP_PI * radius * radius;
		ball_point_of_impulse.set(geom_center_os);
	    } else {
			
		if (factor > -buoyancy_eps) { //speed_plane_os is orthogonal on surface_os
		    IVP_U_Float_Point center_of_cut_off_part_above_surface;
		    center_of_cut_off_part_above_surface.set(geom_center_os);
		    center_of_cut_off_part_above_surface.add_multiple(surface_os, -(length_p1p2*length_p1p2*length_p1p2 / (12 * cut_off_part_above_surface)));
		    center_of_segment_os.set_multiple(&center_of_cut_off_part_above_surface, -cut_off_part_above_surface);
		    center_of_segment_os.add_multiple(geom_center_os, IVP_PI*radius*radius);
		    center_of_segment_os.mult( 1.0f / surface_content_down );
				
		    //if both surfaces are orthogonal, then nothing else has to be computed
		    ball_projected_surface_content_under = surface_content_down;
		    ball_point_of_impulse = center_of_segment_os;
		    goto calculate_impulse;
		} else {
		    //IVP_U_Float_Point direction_up;
		    //direction_up.subtract(&O_os, geom_center_os);
		    //direction_up.normize();
				
		    IVP_U_Float_Point center_of_cut_off_part_above_surface;
		    center_of_cut_off_part_above_surface.set(geom_center_os);
		    center_of_cut_off_part_above_surface.add_multiple(&direction_up, (length_p1p2*length_p1p2*length_p1p2) / (12 * cut_off_part_above_surface));
		    //compute center_of_segment_os by weighed surface contents
		    center_of_segment_os.set_multiple(&center_of_cut_off_part_above_surface, -cut_off_part_above_surface);
		    center_of_segment_os.add_multiple(geom_center_os, IVP_PI*radius*radius);
		    center_of_segment_os.mult( 1.0f / surface_content_down );
		}
	    }
			
			
	    /*** compute part of ellipse above surface ***/
			
	    //compute surface content of ellipse segment above surface
	    IVP_U_Float_Point edge_O_center_os;