[comp.sys.sgi] new chain hanging off cursor

architec@cutmcvax.cs.curtin.edu.au (Phil Dench ) (01/05/91)

For those that are interested, heres the latest version of that chain
hanging off cursor thing with new spider and wind options.

Phil

PS Ken, thanks for message. Lost address though so couldn't reply.

--------8<--------8<--------8<--------8<--------8<--------8<--------8<--------

/* $Header: /usr/people/envy/architec/src/chain/RCS/chain.c,v 1.4 91/01/05 13:24:42 architec Exp $ */

/****************************************************************************\
*                                                                            *
*                                                                            *
*                          Chain hanging off cursor                          *
*                                                                            *
*                                                                            *
*                          Original SDL code from                            *
*                                                                            *
*                           Alias/2 Version 2.4.1                            *
*            Scene Descriptiopn Language (SDL) Programming Manual            *
*                     Introduction to Dynamics Tutorial                      *
*                            Alias Research Inc.                             *
*                                                                            *
*                                                                            *
*                            Coerced into C/gl by                            *
*                                                                            *
*                                 Phil Dench                                 *
*                     architec@cutmcvax.cs.curtin.edu.au                     *
*                                                                            *
*                                                                            *
\****************************************************************************/

/*
 * $Log:	chain.c,v $
 * Revision 1.4  91/01/05  13:24:42  architec
 * added wind
 * 
 * Revision 1.3  91/01/04  18:20:32  architec
 * added crawling end (spider?) option
 * 
 * Revision 1.2  91/01/02  14:52:51  architec
 * clean up before release
 * 
 * Revision 1.1  91/01/02  14:06:18  architec
 * Initial revision
 * 
 */

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

#include <gl.h>
#include <device.h>

static char 	* prog_name;

static Boolean	end_crawls;
static Boolean	sleep_at_start;
static int		refresh_rate;
static int 		num_part;
static double	max_wind_force;

static Boolean	stuck_on_floor;
static Boolean	stuck_on_wall;
static Boolean	stuck_on_ceiling;

static double	wind_force;

#define 		RADIUS			1.0
#define 		SCALE 			20.0
#define			STICK_FORCE		0.5
#define			CRAWL_RATE		0.01

static void usage()
{
	(void) fprintf( stderr, "usage: %s [-c] [-n<num>] [-r<num>] [-s] [-w<num>]\n", prog_name);
}

static void help()
{
	(void) fprintf( stderr, "\n");
	usage();
	(void) fprintf( stderr, "\n");
	(void) fprintf( stderr, " where ...\n");
	(void) fprintf( stderr, "\n");
	(void) fprintf( stderr, "    -c         : end crawls along sides\n");
	(void) fprintf( stderr, "    -n<num>    : number of chain links [30]\n");
	(void) fprintf( stderr, "    -r<num>    : refresh rate          [6]\n");
	(void) fprintf( stderr, "    -s         : sleep at start\n");
	(void) fprintf( stderr, "    -w<num>    : max wind force        [1.0]\n");
	(void) fprintf( stderr, "\n");
}

static void parse_command_line( argc, argv)
	int argc;
	char ** argv;
{
	int arg;

	prog_name = (prog_name = strrchr( argv[ 0], '/')) ? (prog_name + 1) : argv[ 0];

	end_crawls = FALSE;
	sleep_at_start = FALSE;
	refresh_rate = 6;
	num_part = 30;
	max_wind_force = 1.0;

	for( arg = 1; arg < argc; arg++)
	{
		if( argv[ arg][ 0] == '-')
		{
			switch( argv[ arg][ 1])
			{
				case 'c':
					end_crawls = TRUE;
				break;

				case 'n':
					num_part = atoi( argv[ arg] + 2);
				break;

				case 'r':
					refresh_rate = atoi( argv[ arg] + 2);
					if( refresh_rate < 2)
					{
						refresh_rate = 2;
					}
				break;

				case 's':
					sleep_at_start = TRUE;
				break;

				case 'w':
					max_wind_force = atof( argv[ arg] + 2);
				break;

				case 'h':
					help();
					exit( 1);
				break;

				default:
					usage();
					exit( 1);
				break;
			}
		}
		else
		{
			usage();
			exit( 1);
		}
	}

	max_wind_force *= 0.01;
}

static void open_graphics()
{
	/*
	noborder();
	prefposition( 0, 1, 0, 1);
	*/

	/*
	foreground();
	*/

	(void) winopen( prog_name);

	color( BLACK);
	clear();

	fullscrn();

	ortho2( 0.0, 1280.0 / SCALE, 0.0, 1024.0 / SCALE);

	drawmode( PUPDRAW);
}

static void draw_chain( num_part, x, y, old_x, old_y)
	int		num_part;
	double 	* x,
			* y,
			* old_x,
			* old_y;
{
	int part;

	/* Loop for each object outputting the geometry */

	color( PUP_CLEAR);

	move( old_x[ 0], old_y[ 0], 0.0);
	for( part = 1; part < num_part; part++)
	{
		draw( old_x[ part], old_y[ part], 0.0);
	}

	color( PUP_WHITE);

	move( old_x[ 0] = x[ 0], old_y[ 0] = y[ 0], 0.0);
	for( part = 1; part < num_part; part++)
	{
		draw( old_x[ part] = x[ part], old_y[ part] = y[ part], 0.0);
	}
}

static void calc_chain( num_part, x, y, vx, vy, mass)
	int 	num_part;
	double	* x,
			* y,
			* vx,
			* vy,
			* mass;
{
	double forcex[ 64], forcey[ 64];

	/* Define the control parameters for the particles and the springs */

	static int numint = 10;

	static double radius = RADIUS;

	static double timescale = 2.1;
	static double gravity = 0.125 / 30.0;
	static double ymax = 1022.0 / SCALE;
	static double ymin =    2.0 / SCALE;
	static double xmax = 1278.0 / SCALE;
	static double xmin =    2.0 / SCALE;
	static double elasticity = 0.95;
	static double floor_damping = 0.1;
	static double air_damping = 0.999;
	static double bond = 4.0;
	static double damping = 4.0;
	static double bend = 0.15;
	static double bend_damping = 0.15;

	double rad2 = radius * radius;
	double floor_thickness = radius / 2.0;

	double a = bond;
	double b = bond / (radius * radius);

	double dt = timescale * 1.0 / numint;
	double dt2 = 0.5 * dt * dt;

	int sub;
	int part;

	/* Loop for each subinterval of time */

	for( sub = 0; sub < numint; sub++)
	{
		for( part = 0; part < num_part; part++)
		{
			forcex[ part] = 0.0;
			forcey[ part] = 0.0;
		}

		/* Loop for each part to compute the effect of the neighbour forces */

		for( part = 0; part < num_part - 1; part++)
		{
			double dx, dy, r2;
			int other;

			other = part + 1;

			/* Compute the interparticle attraction */
				
			dx = x[ part] - x[ other];
			dy = y[ part] - y[ other];
				
			r2 = dx * dx + dy * dy;

			/* original code  - causes chain to stretch - pd
			if( r2 < 4.0 * rad2)
			*/
			{
				double r, r3, r5;
				double f, scf;
				double dvx, dvy;

				r = sqrt( r2);
				if( r < 0.0001) r = 0.0001;
								
				/* Compute the attraction-repulsion term */
					
				r3 = r * r2;
				r5 = r3 * r2;					

				/* the original calc - causes chain to stretch - pd
				scf = (radius - 0.5 * r) / radius;
				*/

				scf = 0.5;

				dvx = vx[ part] - vx[ other];
				dvy = vy[ part] - vy[ other];
				
				/* the original calc - causes chain to stretch - pd
				f = ((a / r5 - b / r3) - damping *
						 (dvx * dx + dvy * dy) / r2) * scf;
				*/

				f = ((a / r2 - b / r) - damping *
						 (dvx * dx + dvy * dy) / r2) * scf;
				
				forcex[ part] = forcex[ part] + f * dx + wind_force * fabs( dy);
				forcey[ part] = forcey[ part] + f * dy;

				forcex[ other] = forcex[ other] - f * dx + wind_force * fabs( dy);
				forcey[ other] = forcey[ other] - f * dy;

				/* Compute the floor drag-factor term */

				if( y[ part] < ymin + floor_thickness)
				{
					double atn;

					atn = floor_damping * (ymin + floor_thickness -
										 y[ part]) / (floor_thickness);
					forcex[ part] = forcex[ part] - atn * vx[ part];
					forcey[ part] = forcey[ part] - atn * vy[ part];
				}
			}
		}

		/* Loop for each particle computing the torques */
		
		for( part = 1; part < num_part - 1; part++)
		{
			int last, next;
			double xmean, ymean;
			double dx, dy;
			double dxb, dyb;
			double dvx, dvy;
			double vxmean, vymean;
			double d;

			last = part - 1;
			next = part + 1;
			xmean = (x[ last] + x[ part] + x[ next]) / 3.0;
			ymean = (y[ last] + y[ part] + y[ next]) / 3.0;
			dx = x[ part] - xmean;
			dy = y[ part] - ymean;
			d = sqrt( dx * dx + dy * dy);
			dxb = dx * bend;
			dyb = dy * bend;
			
			vxmean = (vx[ last] + vx[ part] + vx[ next]) / 3.0;
			vymean = (vy[ last] + vy[ part] + vy[ next]) / 3.0;
			dvx = vx[ part] - vxmean;
			dvy = vy[ part] - vymean;

			dxb = dxb + bend_damping * dvx;
			dyb = dyb + bend_damping * dvy;
			
			forcex[ part] = forcex[ part] - dxb;
			forcey[ part] = forcey[ part] - dyb;

			dxb = dxb * 0.5;
			dyb = dyb * 0.5;

			forcex[ last] = forcex[ last] + dxb;
			forcey[ last] = forcey[ last] + dyb;

			forcex[ next] = forcex[ next] + dxb;
			forcey[ next] = forcey[ next] + dyb;
		}

		/* Loop for each particle applying the forces to the position */

		for( part = 0; part < num_part; part++)
		{
			double ax, ay;
			double mag;

			ax = forcex[ part] / mass[ part];
			ay = - gravity + forcey[ part] / mass[ part];

			/*
			if( part > 0 || frame > 900)
			*/

			if( (part == (num_part - 1)) && (stuck_on_floor || stuck_on_wall || stuck_on_ceiling))
			{
				mag = sqrt( (forcex[ part] * forcex[ part]) + (forcey[ part] * forcey[ part]));

				if( mag > STICK_FORCE)
				{
					stuck_on_floor = FALSE;
					stuck_on_wall = FALSE;
					stuck_on_ceiling = FALSE;
				}
			}

			if( (part == (num_part - 1)) && stuck_on_floor)
			{
				if( x[ part] < ((xmax + xmin) / 2.0))
				{
					x[ part] = x[ part] - CRAWL_RATE * (1.2 - mag / STICK_FORCE);
				}
				else
				{
					x[ part] = x[ part] + CRAWL_RATE * (1.2 - mag / STICK_FORCE);
				}

				vx[ part] = 0.0;
				vy[ part] = 0.0;
			}
			else if( (part == (num_part - 1)) && stuck_on_wall)
			{
				y[ part] = y[ part] + CRAWL_RATE * (1.2 - mag / STICK_FORCE);

				vx[ part] = 0.0;
				vy[ part] = 0.0;
			}
			else if( (part == (num_part - 1)) && stuck_on_ceiling)
			{
				if( x[ part] < ((xmax + xmin) / 2.0))
				{
					x[ part] = x[ part] + CRAWL_RATE * (1.2 - mag / STICK_FORCE);
				}
				else
				{
					x[ part] = x[ part] - CRAWL_RATE * (1.2 - mag / STICK_FORCE);
				}

				vx[ part] = 0.0;
				vy[ part] = 0.0;
			}
			else if( part > 0)
			{			
				x[ part] = x[ part] + vx[ part] * dt + ax * dt2;
				y[ part] = y[ part] + vy[ part] * dt + ay * dt2;
				
				vx[ part] = vx[ part] + ax * dt;
				vy[ part] = vy[ part] + ay * dt;

				vx[ part] *= air_damping;
				vy[ part] *= air_damping;
			}
			else
			{
				x[ part] = getvaluator( MOUSEX) / SCALE;
				y[ part] = getvaluator( MOUSEY) / SCALE;

				vx[ part] = 0.0;
				vy[ part] = 0.0;
			}

			/* Bounce off the floor */

			if( end_crawls && (part == (num_part - 1)))
			{
				if( y[ part] > ymax)
				{
					stuck_on_ceiling = TRUE;
					stuck_on_wall = FALSE;
				}
				else if( y[ part] < ymin)
				{
					stuck_on_floor = TRUE;
				}

				if( ((x[ part] > xmax) || (x[ part] < xmin)) && !stuck_on_ceiling)
				{
					stuck_on_wall = TRUE;
					stuck_on_floor = FALSE;
				}
			}
			else
			{
				if( y[ part] > ymax)
				{
					y[ part] = ymax + (ymax - y[ part]) * elasticity;
					vy[ part] = -vy[ part] * elasticity;
				}
				else if( y[ part] < ymin)
				{
					y[ part] = ymin + (ymin - y[ part]) * elasticity;
					vy[ part] = -vy[ part] * elasticity;
				}

				if( x[ part] > xmax)
				{
					x[ part] = xmax + (xmax - x[ part]) * elasticity;
					vx[ part] = -vx[ part] * elasticity;
				}
				else if( x[ part] < xmin)
				{
					x[ part] = xmin + (xmin - x[ part]) * elasticity;
					vx[ part] = -vx[ part] * elasticity;
				}
			}
		}
	}
}

extern double drand48();

static void calc_wind()
{
	static double dwind_force = 0.0;

	dwind_force += ((drand48() - 0.5) * max_wind_force / 1000.0);
	wind_force += dwind_force;

	if( ((wind_force > max_wind_force) && (dwind_force > 0.0)) ||
		((wind_force < -max_wind_force) && (dwind_force < 0.0)))
	{
		dwind_force = 0.0;
	}
}

int main( argc, argv)
	int 	argc;
	char 	** argv;
{
	/* Simulate interparticle attraction for several bodies */

	/* Define the arrays to store the particle descriptions */

	double x[ 64], y[ 64], vx[ 64], vy[ 64], mass[ 64];
	double old_x[ 64], old_y[ 64];

	int part;

	Boolean paused = FALSE;
	Boolean quit = FALSE;

	stuck_on_floor = FALSE;
	stuck_on_wall = FALSE;
	stuck_on_ceiling = FALSE;

	wind_force = 0.0;

	parse_command_line( argc, argv);

	/* sleep for a bit */
	if( sleep_at_start)
	{
		(void) sleep( 120);
	}

	open_graphics();

	old_x[ 0] = x[ 0] = getvaluator( MOUSEX) / SCALE;
	old_y[ 0] = y[ 0] = getvaluator( MOUSEY) / SCALE;

	for( part = 1; part < num_part; part++)
	{
		old_x[ part] = x[ part] = x[ part - 1] + RADIUS;
		old_y[ part] = y[ part] = y[ part - 1];

		vx[ part] = 0.0;
		vy[ part] = 0.0;
		mass[ part] = 1.0;
	}

	qdevice( TIMER0);
	noise( TIMER0, refresh_rate);

	qdevice( WINFREEZE);
	qdevice( WINTHAW);
	qdevice( WINQUIT);

	while( !quit)
	{
		long dev;
		short val;

		switch( dev = qread( &val))
		{
			case TIMER0:
				if( !paused && !( getbutton( MENUBUTTON)))
				{
					calc_wind();
					calc_chain( num_part, x, y, vx, vy, mass);

					if( !( getbutton( MENUBUTTON)))
					{
						draw_chain( num_part, x, y, old_x, old_y);
					}
				}
			break;

			case WINTHAW:
				paused = FALSE;
			break;

			case WINFREEZE:
				paused = TRUE;
				color( PUP_CLEAR);
				clear();
			break;

			case WINQUIT:
				quit = TRUE;
				color( PUP_CLEAR);
				clear();
			break;

			case REDRAW:
				endfullscrn();
				drawmode( NORMALDRAW);
				color( BLACK);
				clear();

				fullscrn();
				ortho2( 0.0, 1280.0 / SCALE, 0.0, 1024.0 / SCALE);
				drawmode( PUPDRAW);
			break;
		}
	}

	return( 0);
}