Article 12756 of alt.sources:
Path: nntpd.lkg.dec.com!pa.dec.com!crl.dec.com!crl.dec.com!bloom-beacon.mit.edu!spool.mu.edu!newspump.wustl.edu!crcnews.unl.edu!backbone!backbone!wayne
From: wayne@backbone.uucp (Wayne Schlitt)
Newsgroups: alt.sources
Subject: XStar-2.0.0  Orbiting star system simulator    Part04/06
Date: Sun, 27 Aug 1995 23:14:06 GMT
Organization: The Backbone Cabal
Lines: 2727
Sender: wayne@backbone (Wayne Schlitt)
Message-ID: <WAYNE.95Aug27171406@backbone.uucp>
Reply-To: wayne@cse.unl.edu
NNTP-Posting-Host: cse.unl.edu
Originator: wayne@cse.unl.edu


Submitted-by: wayne@backbone.UUCP
Archive-name: XStar-2.0.0/part04


--- cut here ---

# Continuation of Shell Archive, created by backbone!wayne

# This is part 4

# To unpack the enclosed files, please use this file as input to the
# Bourne (sh) shell.  This can be most easily done by the command;
#     sh < thisfilename


# ---------- file animate.c ----------

filename="animate.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file animate.c...
fi

cat << 'END-OF-FILE' > $filename
#include "xstar.h"
#include "xstar_ext.h"
#include <sys/time.h>
#include <memory.h>


/*
** XStar -- a animated n-body solver for X windows
** Copyright (C) 1995  Wayne Schlitt (wayne@cse.unl.edu) and others
**
**  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 2 of the License, or
**  (at your option) any later version.
*/



/*
 * this is the main loop that animates the stars.  It controls the creation
 * of the star systems, the movement, and all user commands.
 */

void
Animate()
{
    int			num_erase;
    int			bind_check;
    point_2d		*cur, *prev;	/* current and previous positions */
    point_2d		*temp_points;
    point_2d		*vk[K];		/* velocities		*/
    double              *m;             /* masses */
    point_2d		*ak[K];		/* accelerations	*/
    XEvent              xev;
    int			found;
    int			i, k;
    double		found_mass, tmass, dist;

    sys_param_type	sp;
    star_disp_type	sd;

    int			color_skip = 0;

    int			paused = FALSE;

    memset( &sp, 0, sizeof( sp ) );
    memset( &sd, 0, sizeof( sd ) );
    sd.next_free = 1;


    /* Allocate memory. */
    if( num_disp_pt )
    {
	sd.max_points = 512;		/* replotting screen needs large buf */

	sd.disp_pts = (disp_point_type *)Safe_Malloc( sizeof( disp_point_type )
						  * num_disp_pt );

	sd.hash_index = (ushort_t *)Safe_Malloc( sizeof( ushort_t )
					     * HASH_TABLE_SIZE );

    }
    else
    {
	sd.max_points = max_stars*10;	/* will rarely be more than 1/4 full */
	sd.tmp_pts = (XPoint *) Safe_Malloc(sizeof(XPoint) * sd.max_points);
	sd.pixels = (int *) Safe_Malloc(sizeof(int) * sd.max_points);
    }
    sd.points = (XPoint *) Safe_Malloc(sizeof(XPoint) * sd.max_points);
    cur = (point_2d *) Safe_Malloc( sizeof( point_2d ) * max_stars );
    prev = (point_2d *) Safe_Malloc( sizeof( point_2d ) * max_stars );
    sd.last_disp = (XPoint *) Safe_Malloc( sizeof( XPoint ) * max_stars );
    m = (double *) Safe_Malloc(sizeof(double) * max_stars);
    sd.star_color = (int *)Safe_Malloc( sizeof( int ) * max_stars );

    for( k = 0; k < K; k++ )
    {
	vk[k] = (point_2d *)Safe_Malloc( sizeof( point_2d ) * max_stars );
	ak[k] = (point_2d *)Safe_Malloc( sizeof( point_2d ) * max_stars );
    }


    /* set up colors */
    sd.color_number = lrand48() % NUM_COLORS;
    if( rotate_colors || multi_colors)
	init_colors( colors, color_gcs, &sd, NUM_COLORS );


    /* initialize the system */
    set_redraw_full( &sd );
    clear_disp_pt( &sd );
    bind_check = 0;
    num_erase = 0;
    
    init_stars( prev, cur, m, vk, ak, max_stars, &sp, &sd );


    /*
     * Seemingly endless loop.
     */
    while( TRUE )
    {
	sd.raw_num_steps++;
	sd.raw_hsteps++;
	sd.num_steps = sd.raw_num_steps * sd.hstep_scale;
	sd.hsteps = sd.raw_hsteps * sd.hstep_scale;
	if( sd.hsteps >= MAX_HSTEP )
	    reset_hstep( &sd );
	
	sd.num_visible = 0;
	
	/* Age the the points. */
	temp_points = prev;
	prev = cur;
	cur = temp_points;
	
	temp_points = vk[K-1];
	for( k = K-1; k > 0; k-- )
	    vk[k] = vk[k-1];
	vk[0] = temp_points;
	
	temp_points = ak[K-1];
	for( k = K-1; k > 0; k-- )
	    ak[k] = ak[k-1];
	ak[0] = temp_points;
	

	/* calculate and display the new locations of each star */
	if( sp.do_bounce )
	    check_bounce( prev, cur, m, vk, ak, max_stars, &sp, &sd );

	move_fna( prev, cur, m, vk, ak, max_stars, &sp, &sd );

	
	/* 
	 * if there are not enough stars to be interesting, then add
	 * a new star or start over with a new system.
	 */
	if( sp.live_stars <= sp.min_stars
	   || ( sp.num_add > 1 && sp.live_stars <= sp.min_stars + 1 )
	   )
	{
	    new_star( prev, cur, m, vk, ak, max_stars, &sp, &sd, 0 );
	    continue;
	}
	

        /*
	 * Display stars
	 *
	 * some of the terms in these formulas should really be adjustable
	 * depending on the speed of the computer and the quality of the
	 * X server and compiler.  max_disp_skip in particular.
	 */
	sd.num_disp_skipped++;
	if( num_disp_pt )
	{
	    int		pts_to_plot = DIFF( sd.next_free, sd.next_disp );
	    
	    if( pts_to_plot > sd.num_visible
	       && ( sd.num_disp_skipped > sd.buffer_factor /* long time between disps */
		   || pts_to_plot >= 4*sd.num_visible		/* jerky lines */
		   || sp.live_stars >= DIFF( sd.next_erase, sd.next_free )	/* buf full */
		   )
	       )
		update_screen( &sp, &sd, &display );
	}
	else
	{
	    if( sd.points_used > sd.num_visible
	       && ( sd.num_disp_skipped > sd.buffer_factor /* long time between disps */
		   || sd.points_used >= 5*sd.num_visible	/* jerky lines */
		   || sd.points_used + stars + MAX_COLLAPSAR >= sd.max_points	/* buf full */
		   )
	       )
		update_screen( &sp, &sd, &display );
	}	
	
        /*
	 * Check for events.
	 */
	sd.num_poll_skipped++;
        if( sd.num_poll_skipped >= POLL_SKIP )
        {
	    int		check_event;
	    
	    sd.num_poll_skipped = 0;
	
	    check_event = XEventsQueued( display.dpy, QueuedAfterFlush );
	    if( !check_event )
	    {
		/* Delay so we don't use all of the cpu time. */
#ifdef USE_USLEEP
		if( delay != 0 ) usleep(delay);
		if( paused ) usleep(3000*1000);
#else
		if( delay != 0 ) nap(0,delay);
		if( paused ) nap(3,0);
#endif
	    }
	
	
	    /* check to see if the star system is still bound */
	    if( verbose > 1 && (bind_check++)*fv_inv > 128 )
	    {
		sys_const	sc;
		
		bind_check = 0;
		
		calc_sys_const( &sc, max_stars, cur, m, vk[0] );
		if( (sc.k + sc.u)*fv*fv/fm > -1E-7 )
		{
		    fprintf( stderr, "%s:%d  Error:  too much energy, system is not bound.\n", __FILE__, __LINE__ );
		    prt_sys_const( max_stars, cur, m, vk[0] );
		    prt_sys_const_delta( orig_sc, max_stars, cur, m, vk[0] );
		}
		else if( verbose > 3 )
		{
/*		    prt_sys_const( max_stars, cur, m, vk[0] ); */
		    prt_sys_const_delta( orig_sc, max_stars, cur, m, vk[0] );
		}
	    }


	    /* change the color of the stars */
 	    if( (color_skip++)*fv_inv > 80000./(POLL_SKIP*NUM_COLORS) )
	    {
		color_skip = 0;
		sd.color_number = (sd.color_number + 1) % NUM_COLORS;

		plot_collapsars( cur, m, &sp, &sd, &display );
	    }
	
	
	    /* Screen saver/background:  restart after 7 minutes
	     * with out a star dieing */
	    if( timeout || root )
	    {
		double	lost_stars;
		
		time_t	cur_time = time(NULL);
		
		if( cur_time - sd.tstamp > 420*( 1 + sp.do_bounce ) )
		    init_stars( prev, cur, m, vk, ak, max_stars, &sp, &sd );
		else if( cur_time - sd.tstamp > 210*( 1 + sp.do_bounce )
			&& sp.num_add > 1 )
		    new_star( prev, cur, m, vk, ak, max_stars, &sp, &sd, 0 );

		
		if( sd.live_erase < sp.live_stars )
		    sd.live_erase = sp.live_stars;
		
		/*
		 * should we clear the screen?
		 */
		if( !num_disp_pt && !sp.star_line )
		    if( sp.do_bounce )
		    {
			if( sd.erase_disps*fv_inv
			   > (rotate_colors ? 190.0 : 100.0) )
			    erase = TRUE;
		    }
		    else
		    {
			lost_stars = sd.live_erase - sp.live_stars;
			
			if( sd.erase_disps >
			   fv*( (rotate_colors ? 3800.0 : 2000.0)
			       - (1200*lost_stars*lost_stars)
			       / (sp.live_stars+2.0)
			       )
			   )
			    erase = TRUE;
		    }
	    }
	
	
	    /* if we can't see any stars, then we should start over */
	    if( sd.num_seen == 0 )
	    {
		if( verbose > 1 )
		    fprintf( stderr, "%s:%d  no stars are visible...\n", __FILE__, __LINE__ );
		
		new_star( prev, cur, m, vk, ak, max_stars, &sp, &sd, 0 );
	    }
	    sd.num_seen = 0;
	
	
	    /*
	     * check for X events and user commands
	     */

	    do
	    {
		if( check_event )
		{
		    XNextEvent(display.dpy, &xev);
		    HandleEvent(&xev, &sd);
		}
		
		/* Add a collapsar or a star */
		if( add || add_star )
		{
		    found = -1;
		    for( i = 0; i < max_stars; i++ )
			if( m[i] <= 0 )
			{
			    found = i;
			    break;
			}
		    
		    if( (found == -1 || stars + MAX_COLLAPSAR >= max_stars )
		       && stars + MAX_COLLAPSAR <= MAX_NUM_STARS )
		    {
			max_stars++;
			if( !num_disp_pt && sd.max_points < max_stars*10 )
			    sd.max_points = max_stars*10;
			
			if( rotate_colors || multi_colors)
			    init_colors( colors, color_gcs, &sd, NUM_COLORS );
			
			
			sd.points = (XPoint *) Safe_Realloc(sd.points, sizeof(XPoint) * sd.max_points );
			sd.star_color = (int *)Safe_Realloc( sd.star_color, sizeof( int ) * max_stars );
			if( !num_disp_pt )
			{
			    sd.tmp_pts = (XPoint *) Safe_Realloc(sd.tmp_pts, sizeof(XPoint) * sd.max_points );
			    sd.pixels = (int *) Safe_Realloc(sd.pixels, sizeof(int) * sd.max_points);
			}
			cur = (point_2d *) Safe_Realloc(cur, sizeof(point_2d) * max_stars);
			prev = (point_2d *) Safe_Realloc(prev, sizeof(point_2d) * max_stars);
			sd.last_disp = (XPoint *) Safe_Realloc( sd.last_disp, sizeof(XPoint) * max_stars );
			m = (double *) Safe_Realloc(m, sizeof(double) * max_stars);
			for( k = 0; k < K; k++ )
			{
			    vk[k] = (point_2d *)Safe_Realloc( vk[k], sizeof( point_2d ) * max_stars );
			    ak[k] = (point_2d *)Safe_Realloc( ak[k], sizeof( point_2d ) * max_stars );
			}
			
			found = max_stars - 1;
			m[found] = 0;
		    }
		    
		    if( found != -1 )
		    {
			sp.num_add++;
			stars++;
			
			new_star( prev, cur, m, vk, ak, max_stars, &sp, &sd,
				 add );
		    }
		    
		    
		    add = FALSE;
		    add_star = FALSE;
		    
		}
		
		/* delete a star */
		if( del_star )
		{
		    del_star = FALSE;
		    
		    /* kill the least massive one */
		    found_mass = collapsar*10;
		    for( i = 0, found = -1; i < max_stars; i++ )
		    {
			dist = sqrt(cur[i].x*cur[i].x + cur[i].y*cur[i].y);
			if( dist < 50 )
			    dist = 50;
			
			tmass = m[i] / dist;
			
			if( m[i] > 0 && tmass < found_mass )
			{
			    found = i;
			    found_mass = tmass;
			}
		    }
		    
		    if( found != -1 && stars > 2 )
		    {
			stars--;
			sp.live_stars--;
			sd.tstamp = time(NULL);
			set_buffer_factor( &sd, &sp );
			if( verbose > 1 )
			    fprintf( stderr, "%s:%d  live_stars=%d  m[%d]=%g\n", __FILE__, __LINE__, sp.live_stars, found, m[found] );
			
			m[found] = 0.0;
			cur[found].x = cur[found].y = prev[found].x = prev[found].y = -1.0;
			vk[0][found].x = vk[0][found].y = 0.0;
			
			if( sp.live_stars <= sp.min_stars )
			    init_stars( prev, cur, m, vk, ak, max_stars, &sp, &sd );
			
			init_fna( max_stars, m, cur, vk[0] );
		    }
		}
		
		
		/* Erase trails */
		if( erase )
		{
		    erase = FALSE;
		    
		    update_screen( &sp, &sd, &display );
		    XFlush( display.dpy );
		    sd.points_plotted += sd.points_used;
		    sd.total_points += sd.points_plotted;
		    num_erase++;
		    sd.points_plotted = 0;
		    sd.points_used = 0;
		    
		    XFillRectangle(display.dpy, display.win, display.erase_gc,
				   0,0, winW, winH);
		    
		    sd.erase_disps = 0;
		    sd.live_erase = sp.live_stars;
		    
		    clear_disp_pt( &sd );
		    plot_collapsars( cur, m, &sp, &sd, &display );
		}
		
		
		/* new start trails */
		if( new_start )
		{
		    new_start = FALSE;
		    init_stars( prev, cur, m, vk, ak, max_stars, &sp, &sd );
		}
		
		
		/* stop updating */
		if( pause_updt )
		{
		    pause_updt = FALSE;
		    paused = !paused;
		}
		
		
		/* Clean up and shut down. */
		if( stop )
		{
		    stop = FALSE; /* reset the "stop" variable */
		    
		    
		    XFillRectangle(display.dpy, display.win, display.erase_gc,
				   0,0, winW, winH);
		    XFlush( display.dpy );
		    
		    if( !root )
			XDestroyWindow(display.dpy, display.win);
		    
		    /* Free Memory */
		    free(sd.points);
		    free(sd.star_color);
		    if( !num_disp_pt )
		    {
			free(sd.tmp_pts);
			free(sd.pixels);
		    }
		    free(cur);
		    free(prev);
		    free(sd.last_disp);
		    free(m);
		    for( k = 0; k < K; k++ )
		    {
			free( vk[k] );
			free( ak[k] );
		    }
		    
		    /* Free the graphics contexts. */
		    XFreeGC(display.dpy, display.star_gc);
		    XFreeGC(display.dpy, display.erase_gc);
		    
		    if( verbose > 0 && last_init != 0 )
		    {
			sd.tstamp = time(NULL);
			if( verbose == 1 )
			    fprintf( stderr,
				    "%5d/%2d/%2d %4d %8.2f %2d %3d %3d %5d   | %4d %5d %7d %4d%9d\n",
				    orig_sp.live_stars, orig_sp.live_stars-sp.mono_stars,
				    orig_sp.num_add, sp.num_collapsar,
				    sp.size, sp.star_circle, sp.star_line, sp.do_bounce,
				    sp.no_speed,
				    sd.tstamp - last_init, num_collide, num_edge,
				    sp.num_add, sd.num_steps );
			else
			    fprintf( stderr, "%s:%d  total life span of the star system: %d sec (%d steps)\n", __FILE__, __LINE__, sd.tstamp - last_init, sd.num_steps );
			
			if( verbose > 3 )
			{
			    prt_sys_const( max_stars, cur, m, vk[0] );
			    prt_sys_const_delta( orig_sc, max_stars, cur, m, vk[0] );
			}
		    }
		    
		    return;
		}
		
		if( screen_exposed )
		{
		    screen_exposed = FALSE;
		    
		    redraw_screen( cur, m, &sp, &sd, &display );

		    set_redraw_full( &sd );
		}
		
		if( do_redraw )
		{
		    do_redraw = FALSE;
		    
		    XFillRectangle(display.dpy, display.win, display.erase_gc,
				   0,0, winW, winH);
		    
		    redraw_screen( cur, m, &sp, &sd, &display );
		}
		
		if( check_event )
		    check_event = XEventsQueued( display.dpy, QueuedAfterReading );
	    }
	    while( check_event );
	}
    }
}
END-OF-FILE

if [ "$filename" != "/dev/null" ]
then
  size=`wc -c < $filename`

  if [ $size != 13223 ]
  then
    echo $filename changed - should be 13223 bytes, not $size bytes
  fi

  chmod 660 $filename
fi

# ---------- file collide.c ----------

filename="collide.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file collide.c...
fi

cat << 'END-OF-FILE' > $filename
#include "xstar.h"
#include "xstar_ext.h"
#include <sys/time.h>



/*
** XStar -- a animated n-body solver for X windows
** Copyright (C) 1995  Wayne Schlitt (wayne@cse.unl.edu)
**
**  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 2 of the License, or
**  (at your option) any later version.
*/



/*
 * This routine is called (via the macro collide()) when two stars get
 * close enough that they should merge.  The minimal separation
 * distance is important.  If you don't try to collapse two nearby
 * stars together soon enough, they reach very high spin velocities.
 * Then, the slightest disturbance can cause them to step over and
 * past each other in a single movement.  In reality they would have
 * had to have passed through each other to do this last step so they
 * should have collided.  After this last step, they are now far
 * enough apart and have such high velocities, that their attraction
 * can no longer keep them together and they fly off at very high
 * speeds.
 *
 * The smaller the FV value, the larger this value needs to be,
 * because movement is coarser and each step can be larger.  The collide
 * distance can be changed with the -l command line option.
 *
 * Be aware that besides this overstep artifact caused by discrete
 * movement, there is also the sling shot effect which is quite real.
 * The sling shot effect is actually used for interplanetary probes
 * and such.  It works by transferring some of the energy from a
 * larger body (planet, star) to the smaller body (probe), allowing it
 * to gain a fair amount of speed.  You can usually tell the
 * difference between the results of the two phenomenon, because the
 * overstep artifact shoots stars off at a very high speed, and the
 * sling shot effect only shoots stars off at a high speed.
 *
 * Now that I can calculate the star systems constants of motion, you
 * can always distinguish these two phenomenons by checking the total
 * energy of the system.  If the energy has increased, then it is the
 * overstep artifact.  If the energy has remained constant, then it is
 * the sling shot effect.
 *
 *
 * Since collapsars don't move, you will never see the sling shot
 * effect from a collapsar.  Because collapsars are heavy, you will
 * often see the overstep artifact from them.
 */

void do_collide( int i, int j, point_2d *cur, point_2d *prev, double *m,
		point_2d *vk[], point_2d *ak[],
		sys_param_type *sp, star_disp_type *sd,
		char *file, int line 
		)
{
    /* collision:  totally inelastic -> combine */
    int		k;
    double tmass;
    double orig_m_i = m[i];
    
    num_collide++;
    if( verbose > 3 )
    {
	prt_sys_const( max_stars, prev, m, vk[1] );
	prt_sys_const_delta( orig_sc, max_stars, prev, m, vk[1] );
    }


    tmass = 1. / (m[i] + m[j]);
    
    for( k = 0; k < K; k++ )
    {
	vk[k][i].x = (m[i]*vk[k][i].x + m[j]*vk[k][j].x)*tmass;
	vk[k][i].y = (m[i]*vk[k][i].y + m[j]*vk[k][j].y)*tmass;
    }
    
    
    /* don't create a new collapsar */
    if( m[i] >= collapsar )
	m[i] += m[j];
    else if ( m[j] >= collapsar )
    {
	m[i] += m[j];

	cur[i] = cur[j];
	prev[i] = prev[j];
    }
    else
    {
	/* the new location of the combined stars */
	prev[i].x = (m[i]*prev[i].x + m[j]*prev[j].x)*tmass;
	prev[i].y = (m[i]*prev[i].y + m[j]*prev[j].y)*tmass;
	
	
	for( k = 0; k < K; k++ )
	{
	    ak[k][i].x = (m[i]*ak[k][i].x + m[j]*ak[k][j].x)*tmass;
	    ak[k][i].y = (m[i]*ak[k][i].y + m[j]*ak[k][j].y)*tmass;
	}
	
	m[i] += m[j];
	if( m[i] >= collapsar )
	{
	    /* This _should_, by construction, never happen. */
	    if( verbose > 0 )
		fprintf( stderr, "%s:%d  combined mass too large: m[%d]=%g\n",
			__FILE__, __LINE__, i, m[i]/fm );
	    
	    m[i] = collapsar * .999;
	}
    }
    
    --sp->live_stars;
    sd->tstamp = time(NULL);
    set_buffer_factor( sd, sp );
    if( verbose > 1 )
	fprintf( stderr, "%s:%d  live_stars=%d  stars: %d,%d  m: %g+%g=%g\n", file, line, sp->live_stars, i, j, orig_m_i/fm, m[j]/fm, m[i]/fm );
    
    m[j] = 0.0;
    cur[j].x = cur[j].y = prev[j].x = prev[j].y = -1.0;
    for( k = 0; k < K; k++ )
    {
	vk[k][j].x = vk[k][j].y = 0;
	ak[k][j].x = ak[k][j].y = 0;
    }	    

    init_fna( max_stars, m, prev, vk[1] );
}


END-OF-FILE

if [ "$filename" != "/dev/null" ]
then
  size=`wc -c < $filename`

  if [ $size != 4353 ]
  then
    echo $filename changed - should be 4353 bytes, not $size bytes
  fi

  chmod 660 $filename
fi

# ---------- file sys_const.c ----------

filename="sys_const.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file sys_const.c...
fi

cat << 'END-OF-FILE' > $filename
#include "xstar.h"
#include "xstar_ext.h"


/*
** XStar -- a animated n-body solver for X windows
** Copyright (C) 1995  Wayne Schlitt (wayne@cse.unl.edu)
**
**  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 2 of the License, or
**  (at your option) any later version.
*/



/*
 * this file contains a bunch of routines that deal with the constants
 * of motion that all star systems have.  (total energy, angular momentum,
 * linear momentum, etc)
 */


double kinetic_energy( int i, int n, point_2d *cur, double *m, point_2d *v )
{
    return .5 * m[i] * (v[i].x*v[i].x + v[i].y*v[i].y);
}


double binding_energy( int i, int n, point_2d *cur, double *m, point_2d *v )
{
    double	dx, dy;
    double	sum;
    int		j;

    sum = 0;
    for( j = 0; j < n; j++ )
    {
	if( m[j] <= 0. )
	    continue;

	if( i == j )
	    continue;

	dx = cur[i].x - cur[j].x;
	dy = cur[i].y - cur[j].y;
	sum += m[j]/sqrt( dx*dx + dy*dy );
    }

    return -.5 * G * m[i] * sum;
}

    


void calc_sys_const( sys_const *sc,
		    int n, point_2d *cur, double *m,  point_2d *v )
{
    int		i;
    double	m_i;
    double	k, u;
    
    sc->k = 0;
    sc->u = 0;
    sc->e = 0;
    sc->cm.x = 0;
    sc->cm.y = 0;
    sc->p.x = 0;
    sc->p.y = 0;
    sc->l = 0;
    sc->m = 0;

    for( i = 0; i < n; i++ )
    {
	m_i = m[i];
	if( m_i <= 0. )
	    continue;

	k = kinetic_energy( i, n, cur, m, v );
	u = binding_energy( i, n, cur, m, v );

	sc->k += k;

	sc->u += u;
	sc->e += k + u;

	sc->cm.x += m_i * cur[i].x;
	sc->cm.y += m_i * cur[i].y;

	sc->p.x += m_i * v[i].x;
	sc->p.y += m_i * v[i].y;

	sc->l += m_i * (cur[i].x * v[i].y - cur[i].y * v[i].x);

	sc->m += m_i;
    }

    sc->cm.x /= sc->m;
    sc->cm.y /= sc->m;

}


void do_prt_sys_const( int n, point_2d *cur, double *m,  point_2d *v, char *file, int line )
{
    sys_const	sc;

    calc_sys_const( &sc, n, cur, m, v );

    fprintf( stderr, "%s:%d  k=%9.3e + u=%9.3e = %9.3e  l=%8.5g  m=%7.4g\n", file, line, sc.k*fv*fv/fm, sc.u*fv*fv/fm, (sc.k + sc.u)*fv*fv/fm, sc.l*fv/fm, sc.m/fm );
    fprintf( stderr, "%s:%d  cm=(%6g,%6g)  p=(%6g,%6g)\n", file, line, sc.cm.x, sc.cm.y, sc.p.x*fv/fm, sc.p.y*fv/fm );
}



#define per_diff(a,b) ( (a) ? 100.*((a)-(b))/(a) : (a)-(b) )

void do_prt_sys_const_delta( sys_const sc0, int n, point_2d *cur, double *m,  point_2d *v, char *file, int line )
{
    sys_const	sc1;

    calc_sys_const( &sc1, n, cur, m, v );

    fprintf( stderr, "%s:%d  k + u = %9.3e%%  l=%8.5g%%\n", file, line, per_diff(sc0.k + sc0.u,sc1.k + sc1.u), per_diff(sc0.l,sc1.l), per_diff(sc0.m,sc1.m) );
    fprintf( stderr, "%s:%d  cm=(%6g,%6g)  p=(%6g,%6g)\n", file, line, sc0.cm.x-sc1.cm.x, sc0.cm.y-sc1.cm.y, sc0.p.x-sc1.p.x, sc0.p.y-sc1.p.y );
}

END-OF-FILE

if [ "$filename" != "/dev/null" ]
then
  size=`wc -c < $filename`

  if [ $size != 2850 ]
  then
    echo $filename changed - should be 2850 bytes, not $size bytes
  fi

  chmod 660 $filename
fi

# ---------- file default_init.c ----------

filename="default_init.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file default_init.c...
fi

cat << 'END-OF-FILE' > $filename
#include "xstar.h"
#include "xstar_ext.h"

/*
** XStar -- a animated n-body solver for X windows
** Copyright (C) 1995  Wayne Schlitt (wayne@cse.unl.edu)
**
**  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 2 of the License, or
**  (at your option) any later version.
*/



/*
 * what every movement routine needs to do upon initialization.
 */

void default_init( int n, double *m, point_2d *x, point_2d *v )
{
    if( verbose > 3 )
    {
	prt_sys_const( n, x, m, v );
	prt_sys_const_delta( orig_sc, n, x, m, v );
    }
    if( verbose > 0 )
	calc_sys_const( &orig_sc, n, x, m, v );
}
END-OF-FILE

if [ "$filename" != "/dev/null" ]
then
  size=`wc -c < $filename`

  if [ $size != 731 ]
  then
    echo $filename changed - should be 731 bytes, not $size bytes
  fi

  chmod 660 $filename
fi

# ---------- file set_sys_param.c ----------

filename="set_sys_param.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file set_sys_param.c...
fi

cat << 'END-OF-FILE' > $filename
#include "xstar.h"
#include "xstar_ext.h"

    
/*
** XStar -- a animated n-body solver for X windows
** Copyright (C) 1995  Wayne Schlitt (wayne@cse.unl.edu)
**
**  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 2 of the License, or
**  (at your option) any later version.
*/



/*
 * This routine decides what kind of star system to create.  It
 * selects from either a normal star system, or one of many special
 * kinds of star systems.  It then sets all the appropriate parameters
 * and lets set_xmva() do the real work.  (Just deciding what do
 * do seems like a fairly good sized routine...)
 */


/* #define DEBUG */


void set_sys_param( sys_param_type *sp, int max_stars )
{
    int		tstars;
#if defined( DEBUG )
    static int  seed = 1;
    long	junk;
    int		i;
#endif
    
#if defined( DEBUG )
    /*
     * set up the random number generator to a predictable state...
     */
    srand48( seed++ );

    /* seeds won't change much, so use the rng a few times to let the
     * sequences diverge a little bit... */
    for( i = 0; i < 48; i++ )
	junk = lrand48();
#endif
    
#if 0
/* ** debug ** */
    /* for easier debugging */
    sp->star_distrib = DIST_UNIFORM;
    sp->num_collapsar = 0;
    sp->star_circle = 1;
    sp->star_line = 0;
    sp->do_bounce = 1;
    sp->no_speed = 1;
    sp->few_stars = 0;
    sp->live_stars = 3.5 + drand48() * drand48() * (stars-3);
    sp->cir_dist = 50;
    sp->min_stars = 1;
    sp->num_add = 0;
    sp->drift = 0;
#else

    sp->live_stars = 2*stars;
    sp->num_collapsar = 0;
    sp->star_circle = 0;
    sp->star_line = 0;
    sp->do_bounce = 0;
    sp->no_speed = 0;
    sp->cir_dist = 0;
    sp->few_stars = 0;
    sp->drift = -1;

    /* try creating a star system, but don't violate these restrictions */
    while( sp->live_stars > stars
	  || sp->num_collapsar + sp->live_stars > max_stars
	  || (sp->few_stars && sp->live_stars*2 > max_stars)
	  )
    {
	if( (lrand48() % 5) <= 1 )
	{
	    /* create a normal star system */
	    sp->star_distrib = lrand48() % 3;
	    sp->num_collapsar = 0;
	    sp->star_circle = 0;
	    sp->star_line = 0;
	    sp->do_bounce = 0;
	    sp->no_speed = 0;
	    sp->cir_dist = 0;
	    sp->few_stars = 0;
	    sp->live_stars = stars;
	    sp->min_stars = 1;
	    sp->num_add = 0;
	}
	else
	{
	    switch( lrand48() % 8 )
	    {
	    case 0:
		/* 10% of the time, create a small number of stars */
		sp->star_distrib = lrand48() % 3;
		sp->num_collapsar = 0;
		sp->star_circle = 0;
		sp->star_line = 0;
		sp->do_bounce = 0;
                sp->no_speed = 0;
		sp->cir_dist = 0;
		sp->few_stars = 1;
		sp->live_stars = 2 + lrand48() % 4;
		sp->min_stars = 1;
		sp->num_add = 0;
		break;
		
	    case 1:
		/* 10% of the time, create a circular star system */
		tstars = ( stars > 25 ) ? 25 : stars-3;
		sp->star_distrib = DIST_UNIFORM;
		sp->num_collapsar = 0;
		sp->star_circle = 1;
		sp->star_line = 0;
		sp->do_bounce = 0;
                sp->no_speed = 0;
		sp->few_stars = 0;
		sp->live_stars = 3.5 + drand48() * drand48() * tstars;
		sp->cir_dist = 50 + RAND(80);
		sp->min_stars = 1;
		sp->num_add = 0;
		sp->drift = (lrand48() % 5) == 0;
		break;
		
	    case 2:
		/* 10% of the time, create a linear star system */
		sp->star_distrib = DIST_UNIFORM;
		sp->num_collapsar = 0;
		sp->star_circle = 0;
		sp->star_line = 1;
		sp->do_bounce = 1;
		sp->no_speed = 1;
		sp->few_stars = 0;
		switch( lrand48() % 6 )
		{
		case 0:
		    sp->live_stars = 2;
		    break;

		case 1:
		case 2:
		    sp->live_stars = 4;
		    break;

		default:
		    tstars = ( stars > 25 ) ? 25 : stars-2;
		    sp->live_stars = 5 + drand48() * drand48() * drand48() * tstars;
		    break;
		}
		sp->cir_dist = 40 + drand48() * 60;
		sp->min_stars = 1;
		sp->num_add = 0;
		sp->drift = 1;

		/* due to bugs in the bounce code, this should be even */
		sp->live_stars += sp->live_stars & 1;

		/* if we only have 2 stars, put a spin on it */
		if( sp->live_stars == 2 )
		{
		    if( (lrand48() % 4) )
		    {
			sp->no_speed = 0;
			sp->drift = (lrand48() % 3) > 0;
		    }
		}
		else
		    sp->rnd_spacing = (lrand48() % 3) == 0;

		break;
		
	    case 3:
	    case 4:
		/* 20% of the time, create bouncing stars */
		sp->star_distrib = lrand48() % 3;
		sp->num_collapsar = 0;
		sp->star_circle = 0;
		sp->star_line = 0;
		sp->do_bounce = 1;
                sp->no_speed = (lrand48() % 3) == 0;
		sp->cir_dist = 0;
		sp->few_stars = 0;
		sp->live_stars = 3 + lrand48() % 6;
		sp->min_stars = 1;
		sp->num_add = 0;
		sp->drift = 1;

		if( sp->no_speed && sp->star_distrib == DIST_CENTER )
		    sp->star_distrib = DIST_UNIFORM;

		/* every once and a while, throw in some collapsars */
		if( lrand48() % 3 == 0 )
		    sp->num_collapsar = lrand48() % (MAX_COLLAPSAR-1);
		break;


	    default:
		/* 50% of the time, create collapsars */
		
		sp->star_distrib = 1 + lrand48() % 2;

		switch( lrand48() % 10 )
		{
		case 0:
		    sp->num_collapsar = 1;
		    sp->min_stars = 1;
		    break;
		    
		case 1:
		case 2:
		    sp->num_collapsar = 2;
		    sp->min_stars = 0;
		    break;
		    
		case 3:
		case 4:
		case 5:
		    sp->num_collapsar = 3;
		    sp->min_stars = 0;
		    break;
		    
		case 6:
		case 7:
		    sp->num_collapsar = 4;
		    sp->min_stars = 0;
		    break;

		case 8:
		    sp->num_collapsar = 5;
		    sp->min_stars = 0;
		    break;
		    
		case 9:
		default:
		    sp->num_collapsar = 6;
		    sp->min_stars = 0;
		    break;
		}
		
		sp->star_circle = 0;
		sp->star_line = 0;
		sp->do_bounce = 0;
                sp->no_speed = 0;
		sp->cir_dist = 0;
		sp->few_stars = 0;
		sp->live_stars = stars;
		
		
		sp->num_add = stars / 3;
		sp->live_stars -= sp->num_add;
		sp->drift = 0;
	    }
	}	
    }	
#endif

    /* 15-25% binary stars */
    if( sp->live_stars > 9 && !sp->star_circle && !sp->star_line )
	sp->mono_stars = sp->live_stars * (.925 - .05 * drand48())
	    + drand48()*2;
    else
	sp->mono_stars = sp->live_stars;
    if( sp->mono_stars > sp->live_stars )
	sp->mono_stars = sp->live_stars;


    /* should stars drift? */
    if( sp->drift == -1 )
	if( sp->live_stars <= 5 )
	{
	    if( sp->live_stars == 2 )
		sp->drift = (lrand48() % 3) < 2;
	    else
		sp->drift = 1;
	}
	else
	    sp->drift = 0;


    /*
     * calculate the size and amount of energy the system should have
     */
    sp->size = (sp->live_stars + sp->num_collapsar) * STARAREA;
    switch( sp->star_distrib )
    {
    case DIST_CENTER:
	sp->size = sqrt( sp->size ) * .6666;
	break;

    case DIST_UNIFORM:
	sp->size = sqrt( sp->size ) * .64;
	break;

    case DIST_RING:
	sp->size = sqrt( sp->size ) * .58;
	break;

    default:
	sp->size = 0;
	fprintf( stderr, "%s:%d  Error!  star distribution=%d\n",
		sp->star_distrib );
	exit( 1 );
	break;
    }

    if( sp->do_bounce )
	sp->size *= .5;


    /*
     * calculate the energy factor for some types of systems
     */
    sp->energy_fact = -1;
    if( sp->star_circle )
    {
	if( lrand48() % 16 )
	    sp->energy_fact = drand48() * .31 + .54;
	else
	    sp->energy_fact = .5;	/* perfect circle		*/
    }

    else if( sp->live_stars <= 2 )
    {
	if( sp->star_line )
	    sp->energy_fact = drand48() * .04 + .0005;
	else if( (lrand48() % 8) == 0 )
	    sp->energy_fact = .5;
	else if( lrand48() % 2 )
	    sp->energy_fact = drand48() * .9 + .06;
	else
	    sp->energy_fact = drand48() * .4 + .06;
    }
    sp->calc_energy_fact = ( sp->energy_fact == -1 );
    
	
    
    if( verbose > 1 )
	fprintf( stderr, "%s:%d  num stars=%d  num collapsars=%d  size=%f\n",
		__FILE__, __LINE__, sp->live_stars, sp->num_collapsar,
		sp->size );
}
END-OF-FILE

if [ "$filename" != "/dev/null" ]
then
  size=`wc -c < $filename`

  if [ $size != 7729 ]
  then
    echo $filename changed - should be 7729 bytes, not $size bytes
  fi

  chmod 660 $filename
fi

# ---------- file set_xmva.c ----------

filename="set_xmva.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file set_xmva.c...
fi

cat << 'END-OF-FILE' > $filename
#include "xstar.h"
#include "xstar_ext.h"


/*
** XStar -- a animated n-body solver for X windows
** Copyright (C) 1995  Wayne Schlitt (wayne@cse.unl.edu)
**
**  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 2 of the License, or
**  (at your option) any later version.
*/

/* note:  new_star() has similar functionality... change both places */

/*
 * The routine set_xmva sets up the star system's location (x), mass
 * (m), velocity (v) and acceleration (a).  The star system is
 * carefully crafted so that the stars orbit each other and the star
 * system doesn't drift off the screen.
 *
 *
 * The algorithm that sets up the star system is somewhat crude.  It
 * picks a direction that the star will move by combining the vector
 * needed to orbit the center and a vector needed to orbit the nearest
 * star.  The magnitude of the velocity is set by assuming that having
 * the kinetic energy to be no more than half the binding energy is a
 * good value.  This seems to work fairly well, but it still has too
 * much voodoo magic in it for my liking.
 *
 * Over the course of developing xstar, I tried many different methods
 * of setting up the initial star system.  Xgrav just set random
 * velocities and it worked well enough to be interesting.  My first
 * attempt at an improvement was to use all sorts of magic numbers and
 * weird ad-hoc formulas to try and make these random settings more
 * useful.  It worked better than xgrav's method, but it didn't scale
 * well and it clearly had room for improvement.  I found that just
 * making every star orbit the center was a major improvement, but
 * many other tweaks failed.
 *
 * I tried setting up the stars so that the velocity of each star was
 * the sum of the velocities needed to orbit each star in the star
 * system.  I had thought that would make a "perfect" orbit, but I was
 * mistaken.  The resulting velocities were too large.  The velocities
 * could be fudged to ok values, but that didn't solve all the the
 * problems.  The resulting star system tended to not orbit the
 * center, and there were still a lot of early collisions.  This was
 * before I had the binding energy calculation working, so maybe with
 * this correction factor this method might work better.
 *
 * My next though was that for each star, the remaining stars set up
 * a force field, and the gradient of that field would tell me which
 * direction the change in acceleration would be the greatest.  By picking
 * a velocity perpendicular to the gradient, the star would be moving
 * in a direction of constant acceleration and you could use the formula
 * a = v^2 / r to determine the speed.  I still don't think this would
 * create the perfect orbit for the star since the force field changes
 * with time.
 *
 * I tried to think of some sort of iterative approach to the problem,
 * where I would pick an initial approximate velocity vector for each
 * star and then let the velocities converge toward a perfect a perfect
 * orbit, but I couldn't think of any way of doing this.
 *
 * So, while the current star system is a hack, it is the best hack
 * that I know how to make.
 *
 * Heck, it is possible that if you just used random velocities like xgrav
 * had, and use the binding energy correction, that the star system
 * would work fairly well.  The binding energy correction was clearly
 * a very powerful tool to making interesting star systems.
 *
 *
 * Note:  I was just thinking...  Right now we create systems that rotate,
 * yielding systems with a large angular momentum.   Since angular momentum
 * is conserved, that means after a bunch of collisions, we must still have
 * a large angular momentum therefore the star system must be very spread
 * out.  If we want to keep things on the screen, we should try to keep the
 * angular momentum close to zero.  Or, am I missing something...
 */

double sqr( double a )
{
    return a*a;
}


void set_xmva( point_2d *cur, point_2d *prev, double *m,
	      point_2d *vk[], point_2d *ak[],
	      int max_stars, sys_param_type *sp )
{
    int		i, j;

    point_2d	tp;			/* total momentum		*/
    point_2d	cm;			/* center of mass		*/
    double	tm;			/* total mass			*/

    double	dx, dy;

    double	dist;
    double	max_dist, cur_dist;

    int		cur_dist_idx;
    int		tot_stars;


    double	theta, theta0;

    point_2d	*v_0 = vk[0];

    point_2d	drift_amount;

    cm.x = cm.y = 0;
    tm = 0;
    
    
    /*
     * calculate the amount to drift
     */
    if( sp->drift )
    {
	if( sp->star_line )
	{
	    drift_amount.x = 0;
	    do
		drift_amount.y = (drand48() - .5 ) * fv_inv * .0025;
	    while( fabs( drift_amount.y ) < fv_inv * .0005 );
	}
	else
	{
	    drift_amount.x = (drand48() - .5 ) * fv_inv * .001;
	    drift_amount.y = (drand48() - .5 ) * fv_inv * .001;
	}
    }
    

    /*
     * create the star system
     */
    if( sp->star_line )
    {
	/*
	 * create stars in a straight line
	 */
	
	cm.x = cm.y = 0;
	tm = 0;

	for( i = 0; i < sp->live_stars; i++ )
	{
	    if( sp->rnd_spacing )
		prev[i].x = cur[i].x = sp->cir_dist * ( i + drand48() - .5 );
	    else
		prev[i].x = cur[i].x =
		    sp->cir_dist * (i - sp->live_stars*.5 + .5);
	    prev[i].y = cur[i].y = 0;

	    m[i] = fm*1;
	

	    /* keep track of what we have done for later adjustments */
	    tm += m[i];

	    cm.x += m[i]*prev[i].x;
	    cm.y += m[i]*prev[i].y;
	}
    }
	
    else if( sp->star_circle )
    {
	/*
	 * create stars in a circle
	 */

	cm.x = cm.y = 0;
	tm = 0;

	theta = 2 * M_PI/ sp->live_stars;
	theta0 = M_PI * drand48();

	for( i = 0; i < sp->live_stars; i++ )
	{
	    prev[i].x = cur[i].x = sp->cir_dist * cos( i * theta + theta0 );
	    prev[i].y = cur[i].y = sp->cir_dist * sin( i * theta + theta0 );

	    m[i] = fm*1;

	    /* keep track of what we have done for later adjustments */
	    tm += m[i];

	    cm.x += m[i]*prev[i].x;
	    cm.y += m[i]*prev[i].y;
	}
    }

    else
    {
	/*
	 * create a bunch of stars
	 */
	
	max_dist = far_dist * far_dist;
	for( i = 0; i < sp->live_stars; i++ )
	{
	    if( i < sp->mono_stars )
	    {
		do
		{
		    /*
		     * use appropriate distribution of stars
		     */
		    theta = 2*M_PI * drand48();
		    switch( sp->star_distrib )
		    {
		    case DIST_CENTER:
			dist = drand48() * sp->size;
			break;
			
		    case DIST_UNIFORM:
			dist = sqrt( drand48() ) * sp->size;
			break;
			
		    case DIST_RING:
			dist = sqrt(sqrt( drand48() )) * sp->size;
			break;
			
		    default:
			/* this shouldn't happen... */
			dist = 100;
		    }
		    prev[i].x = cur[i].x = dist * cos( theta );
		    prev[i].y = cur[i].y = dist * sin( theta );

		    
		    /*
		     * don't put stars too close together
		     */
		    cur_dist = far_dist * far_dist;
		    for( j = 0; j < i; j++ )
		    {
			dx = prev[i].x - prev[j].x;
			dy = prev[i].y - prev[j].y;
			dist = dx*dx + dy*dy;
			
			if( dist < cur_dist )
			    cur_dist = dist;
		    }
		}
		while( cur_dist < 30*30
		      && ( cur_dist < 20*20 || drand48() < .5 )
		      && ( cur_dist < 10*10 || drand48() < .75 )
		      );
		
		m[i] = fm*90.+fm*RAND(120);
	    }
	    else
	    {
		/*
		 * create binary stars on the outside
		 */
		cur_dist = -1;
		cur_dist_idx = -1;
		for( j = 0; j < sp->mono_stars; j++ )
		{
		    dist = prev[j].x*prev[j].x + prev[j].y*prev[j].y;
		    if( dist < max_dist && dist > cur_dist )
		    {
			cur_dist = dist;
			cur_dist_idx = j;
		    }
		}
		max_dist = cur_dist * .9999;
		
		theta = drand48() * M_PI * 2;
		dist = 6 + drand48() * 8;
		
		prev[i].x = cur[i].x = cur[cur_dist_idx].x + dist * cos(theta);
		prev[i].y = cur[i].y = cur[cur_dist_idx].y + dist * sin(theta);
		
		
		/* cut the mass down on binary stars */
		m[i] = (fm*90.+fm*RAND(120))*.5;
		
		m[cur_dist_idx] *= .5;
		
		tm -= m[cur_dist_idx];
		cm.x -= m[cur_dist_idx]*prev[cur_dist_idx].x;
		cm.y -= m[cur_dist_idx]*prev[cur_dist_idx].y;
		
		
		/* make only half of the outer most stars binary stars */
		cur_dist = -1;
		for( j = 0; j < sp->mono_stars; j++ )
		{
		    dist = prev[j].x*prev[j].x + prev[j].y*prev[j].y;
		    if( dist < max_dist && dist > cur_dist )
			cur_dist = dist;
		}
		max_dist = cur_dist * .9999;
	    }
	    
	    
	    /* keep track of what we have done for later adjustments */
	    tm += m[i];
	    
	    cm.x += m[i]*prev[i].x;
	    cm.y += m[i]*prev[i].y;
	}
    }


    /*
     * relocate the center of mass to the origin
     */
    cm.x /= tm;
    cm.y /= tm;

    for( i = 0; i < sp->live_stars; i++ )
    {
	prev[i].x = cur[i].x = prev[i].x - cm.x;
	prev[i].y = cur[i].y = prev[i].y - cm.y;
    }


    /*
     * scale down the total mass by the number of stars
     *
     * I am really not sure what should be done here.
     *
     * On one hand, if you don't take into account how many stars there
     * are, and you make the average star a certain mass, then after
     * several (many?) collisions, you can end up with a very massive
     * star and you will run into a lot of overstep problems.  You
     * will also be increasing the total mass of the system as the
     * number of stars increases and that will make stars move faster
     * and thus change the accuracy of the system.
     *
     * On the other hand, if you make the star system have a constant
     * value and adjust the mass of each star accordingly, then for
     * very small star systems, the individual stars will be quite
     * massive and you will run into the overstep problem.  Also, even
     * though the total mass of the star system is constant, with lots
     * of stars the mass will be very spread out and the stars will
     * end up moving slower, thus changing the accuracy of the system.
     *
     * So, the mass of the star system should scale with the number of
     * stars, but sub-linearly.  what this formula is, I don't know
     * and I don't know how to figure it out.
     *
     * For right now, it appears that keeping a constant star system
     * mass is better than a simple linear scale, so that's what I am
     * doing.
     */
    if( sp->num_collapsar || sp->do_bounce )
	tm = (fm * 4) / tm;
    else if( sp->few_stars && sp->live_stars > 2 )
	tm = (fm * 6) / tm;
    else
	tm = (fm * 16) / tm;
    
	
    if( tm <= 0 )
    {
	fprintf( stderr, "%s:%d  Error!  tm=%g  should be >0\n",
		__FILE__, __LINE__, tm );
	tm = fm*.1;
    }

    for( i = 0; i < sp->live_stars; i++ )
    {
	m[i] *= tm;
	if( m[i] >= collapsar )
	{
	    /* This _should_, by construction, never happen. */
	    if( verbose > 0 )
		fprintf( stderr, "%s:%d  scaled mass too large: m[%d]=%g\n",
			__FILE__, __LINE__, i, m[i]/fm );

	    m[i] = collapsar * .999;
	}
    }


    /*
     * create the collapsars
     */
    tot_stars = sp->live_stars + sp->num_collapsar;
    if( sp->num_collapsar )
    {
	cm.x = cm.y = 0;
	tm = 0;
	
	for( i = sp->live_stars; i < tot_stars; i++ )
	{
	    /*
	     * don't put stars too close together and use the appropriate
	     * distribution of stars
	     */
	    do
	    {
		theta = 2*M_PI * drand48();
		switch( sp->star_distrib )
		{
		case DIST_CENTER:
		    dist = drand48() * sp->size;
		    break;
			
		case DIST_UNIFORM:
		    dist = sqrt( drand48() ) * sp->size;
		    break;
			
		case DIST_RING:
		    dist = sqrt(sqrt( drand48() )) * sp->size;
		    break;
			
		default:
		    dist = 100;
		}
		if( !sp->do_bounce )
		{
		    prev[i].x = cur[i].x = dist * .5 * cos( theta );
		    prev[i].y = cur[i].y = dist * .5 * sin( theta );
		}
		else
		{
		    prev[i].x = cur[i].x = dist * cos( theta );
		    prev[i].y = cur[i].y = dist * sin( theta );
		}

		cur_dist = far_dist * far_dist;
		for( j = sp->live_stars; j < i; j++ )
		{
		    dx = prev[i].x - prev[j].x;
		    dy = prev[i].y - prev[j].y;
		    dist = dx*dx + dy*dy;

		    if( dist < cur_dist )
			cur_dist = dist;
		}
	    }
	    while( cur_dist < 30*30
		  && ( cur_dist < 20*20 || drand48() < .5 )
		  && ( cur_dist < 10*10 || drand48() < .75 )
		  );

	    v_0[i].x = 0.0;
	    v_0[i].y = 0.0;
	    m[i] = collapsar * 1.0001;
	    
	    /* keep track of what we have done for later adjustments */
	    tm += m[i];
	    
	    cm.x += m[i]*prev[i].x;
	    cm.y += m[i]*prev[i].y;
	}
	
	
	/*
	 * relocate the collapsars' center of mass to the origin
	 */
	cm.x /= tm;
	cm.y /= tm;
	
	for( i = sp->live_stars; i < tot_stars; i++ )
	{
	    prev[i].x = cur[i].x = prev[i].x - cm.x;
	    prev[i].y = cur[i].y = prev[i].y - cm.y;
	}
    }


    /*
     * set the velocities
     */
    for( i = 0; i < sp->live_stars; i++ )
	set_va( cur, prev, m, vk, ak, sp, i );
    


    /*
     * make sure the system has no overall momentum
     */
    tm = 0;
    tp.x = tp.y = 0;
    for( i = 0; i < sp->live_stars; i++ )
    {
	tp.x += v_0[i].x*m[i];
	tp.y += v_0[i].y*m[i];

	tm += m[i];
    }

    tp.x /= tm;
    tp.y /= tm;

    /* ok, if there aren't many stars, let them drift around */
    if( sp->drift )
    {
	tp.x += drift_amount.x;
	tp.y += drift_amount.y;
    }

    for( i = 0; i < sp->live_stars; i++ )
    {
	v_0[i].x -= tp.x;
	v_0[i].y -= tp.y;
    }


    /*
     * if we are drifting, move the starting locations of the stars
     * in the opposite direction
     */

    if( sp->drift && ( !sp->do_bounce || sp->star_line ) )
    {
	double		gap;
	point_2d	edge;
	point_2d	offset;

	if( sp->star_line )
	    gap = 40;
	else
	    gap = sp->size * 1.5;
	    
	edge.x = ( max_x < -min_x ) ? max_x : -min_x;
	edge.y = ( max_y < -min_y ) ? max_y : -min_y;

	if( gap > edge.x )
	    gap = edge.x;
	if( gap > edge.y )
	    gap = edge.y;

	if( fabs( drift_amount.x / edge.x ) > fabs( drift_amount.y / edge.y ) )
	{
	    /* the x direction is the limiting factor */
	    offset.x = ( drift_amount.x > 0 ) ? -edge.x + gap : edge.x - gap;
	    offset.y = offset.x * drift_amount.y / drift_amount.x;
	}
	else
	{
	    /* the y direction is the limiting factor */
	    offset.y = ( drift_amount.y > 0 ) ? -edge.y + gap : edge.y - gap;
	    offset.x = offset.y * drift_amount.x / drift_amount.y;
	}


	for( i = 0; i < tot_stars; i++ )
	{
	    prev[i].x = cur[i].x = prev[i].x - offset.x;
	    prev[i].y = cur[i].y = prev[i].y - offset.y;
	}

    }
}


void set_va( point_2d *cur, point_2d *prev, double *m,
	    point_2d *vk[], point_2d *ak[],
	    sys_param_type *sp,
	    int i )
{
    int		j;

    double	tm;			/* total mass			*/

    double	dx, dy;
    double	tdx, tdy;

    double	dist, tdist;

    double	k, u;
    double	v_fac;

    int		tot_stars;


    point_2d	*v_0 = vk[0];
    point_2d	*a_0 = ak[0];

    double	max_energy;

    double	a_mag, q_a_mag;
    double	f;


    tot_stars = sp->live_stars + sp->num_collapsar;

    
    /*
     * initialize the acceleration array
     */
    for( j = 0; j < tot_stars; j++ )
    {
	double f, f_d;

	if( j == i )
	    continue;

	if( m[j] <= 0.0 )
	    continue;
	    
	dx = prev[j].x - prev[i].x;
	dy = prev[j].y - prev[i].y;
	dist = dx*dx + dy*dy;
	    
	if( dist >= collide_dist )
	{
	    f = G / (sqrt(dist) * dist);
		
	    f_d = dx * f;		
	    a_0[i].x += f_d * m[j];
		
	    f_d = dy * f;
	    a_0[i].y += f_d * m[j];
	}
    }
    
    
    /*
     * set the velocities
     */
    dx = cur[i].x;
    dy = cur[i].y;
    dist = sqrt( dx*dx + dy*dy );
    
    if( dist > 4 && !sp->star_circle )
    {
	/*
	 * assume that the stars form a uniformly dense sphere and
	 * calculate the gravitational force on this particular star.
	 *
	 * we assume that all stars closer to the center than this one
	 * are uniformly distributed and thus can be treated as a
	 * point mass at the center.  we also assume that the stars
	 * further out than this one are uniformly distributed and thus
	 * have no affect on this star.  Of course, none of this is true.
	 */
	tm = 0;
	dist *= .99999999;
	for( j = 0; j < tot_stars; j++ )
	{
	    if( m[j] <= 0 )
		continue;
	    
	    tdx = cur[j].x;
	    tdy = cur[j].y;
	    tdist = sqrt( tdx*tdx + tdy*tdy );
	    
	    if( tdist < dist )
		tm += m[j];
	}
	if( tm < .5*fm )
	    tm = .5*fm;
	
	if( sp->num_collapsar )
	    tm *= .30;
	else
	    tm *= .45;
	
	
	/* vector tangent to a perfectly circular orbit */
	f = sqrt( G * tm / dist);
	
	v_0[i].x = (dy / dist) * f;
	v_0[i].y = (-dx / dist) * f;
	
    }
    
    
    /*
     * add a twist to the velocity to avoid any close by stars
     */
    a_mag = sqrt( a_0[i].x * a_0[i].x + a_0[i].y * a_0[i].y );
    if( a_mag*fv*fv > 1E-9 )
    {
	q_a_mag = sqrt( a_mag );		
	
	if( sp->star_circle )
	{
	    /* orbit center of the stars system */
	    v_0[i].x += (-a_0[i].y * sqrt(sp->cir_dist)) / q_a_mag;
	    v_0[i].y += (a_0[i].x * sqrt(sp->cir_dist)) / q_a_mag;
	}
	else
	{
	    dist = far_dist*far_dist;
	    for( j = 0; j < tot_stars; j++ )
	    {
		if( i == j )
		    continue;

		if( m[j] <= 0 )
		    continue;
		
		tdx = (prev[i].x*m[i] + prev[j].x*m[j]) / (m[i] + m[j])
		    - prev[i].x;
		tdy = (prev[i].y*m[i] + prev[j].y*m[j]) / (m[i] + m[j])
		    - prev[i].y;
		
		tdist = tdx*tdx + tdy*tdy;
		
		if( tdist < dist )
		    dist = tdist;
	    }
	    dist = sqrt( sqrt( dist ) );


	    /* orbit closest stars with radius of 'dist' to the C of M */
	    /* derived from a = v^2/rcm and  a = r12/|r12| * G * m1 / r12^2 */
	    v_0[i].x += (-a_0[i].y * dist) / q_a_mag;
	    v_0[i].y += (a_0[i].x * dist) / q_a_mag;


	    /* give systems with collapsars a random twist */
	    if( sp->num_collapsar )
	    {
		v_0[i].x *= 1 + (drand48() - .4)/2;
		v_0[i].y *= 1 + (drand48() - .4)/2;
	    }
	}
    }
    
    
    /*
     * make sure this star is still bound
     */
    if( !sp->no_speed )
    {
	k = kinetic_energy( i, tot_stars, cur, m, v_0 );
	u = binding_energy( i, tot_stars, cur, m, v_0 );
	
	if( sp->calc_energy_fact )
	{
	    if( sp->num_collapsar )
	    {
		sp->energy_fact = drand48() * .25 + .30;
	    }
	    else
		switch( sp->star_distrib )
		{
		case DIST_CENTER:
		    sp->energy_fact = drand48() * .1 + .50;
		    break;
		    
		case DIST_UNIFORM:
		    sp->energy_fact = drand48() * .1 + .45;
		    break;
		    
		case DIST_RING:
		    sp->energy_fact = drand48() * .1 + .40;
		    break;
		    
		default:
		    sp->energy_fact = 0;
		    fprintf( stderr, "%s:%d  Error!  star distribution=%d\n",
			    sp->star_distrib );
		    exit( 1 );
		    break;
		}
	}
	max_energy = -u * sp->energy_fact;
	
	
	if( k > max_energy || sp->star_circle )
	{
	    v_fac = sqrt( max_energy / k );
	    
	    v_0[i].x *= v_fac;
	    v_0[i].y *= v_fac;
	}
#if 0
	else
	    v_fac = -1;
	    
	    k = kinetic_energy( i, tot_stars, cur, m, v_0 );
	    u = binding_energy( i, tot_stars, cur, m, v_0 );
	    
	    fprintf( stderr, "%s:%d %d  k=%9.6f u=%9.6f  -k/u=%9.6f   |%7.4f,%7.4f|  v_fac=%g\n", __FILE__, __LINE__, i, k*fv*fv/fm, u*fv*fv/fm, -k/u, fv*v_0[i].x, fv*v_0[i].y, v_fac );
#endif
	    
    }
    
    
    /* sometimes it's nice to star with no velocity */
    if( sp->no_speed )
	v_0[i].x = v_0[i].y = 0;
}
END-OF-FILE

if [ "$filename" != "/dev/null" ]
then
  size=`wc -c < $filename`

  if [ $size != 18996 ]
  then
    echo $filename changed - should be 18996 bytes, not $size bytes
  fi

  chmod 660 $filename
fi

# ---------- file update_screen.c ----------

filename="update_screen.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file update_screen.c...
fi

cat << 'END-OF-FILE' > $filename
#include "xstar.h"
#include "xstar_ext.h"

/*
** XStar -- a animated n-body solver for X windows
** Copyright (C) 1995  Wayne Schlitt (wayne@cse.unl.edu)
**
**  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 2 of the License, or
**  (at your option) any later version.
*/



/* update the screen with anything that hasn't been plotted */
void update_screen( sys_param_type *sp, star_disp_type *sd, disp *display )
{

    if( num_disp_pt )
    {
	int		i;

	int		num_pts, pp, next_color;
	int		pts_to_plot;
	int		free_pts;
	int		stop_erase, stop_err;
	int		new_steps;
	

	/*
	 * erase the old points
	 */

	/* need to have at least a minimal amount of slots free */
	stop_erase = sd->next_erase;
	
	free_pts = DIFF( sd->next_erase, sd->next_free );
	if( free_pts < sp->live_stars*2+2 )
	{
	    stop_erase = SUM( stop_erase, (sp->live_stars*2+2) - free_pts );

	    if( sd->erase_hstep < sd->disp_pts[ stop_erase ].hstep )
		sd->erase_hstep = sd->disp_pts[ stop_erase ].hstep;
	}

	new_steps = sd->num_steps - sd->updt_hstep;
	sd->erase_hstep += new_steps;

	/* modify erase_hstep based on how full/empty the buffer is */
	if( free_pts < sp->live_stars*8+2 )
	    sd->erase_hstep += new_steps << 2;

	else if( free_pts < sd->num_disp_pt_32 )
	    sd->erase_hstep += new_steps << 1;

	else if( free_pts < sd->num_disp_pt_16 )
	    sd->erase_hstep += new_steps;

	else if( free_pts < sd->num_disp_pt_8 )
	    sd->erase_hstep += new_steps >> 2;

	else if( free_pts < sd->num_disp_pt_4 )
	    sd->erase_hstep -= (new_steps >> 5) + 1;

	else
	    sd->erase_hstep -= new_steps >> 1;


	    
	/* find where to stop erasing stars at */
	stop_err = PREV( sd->next_erase );
	for( ;; stop_erase = NEXT( stop_erase ) )
	{
	    if( stop_erase == stop_err )
	    {
		if( verbose > 0 )
		    fprintf( stderr, "%s:%d step%7d  Error:  buffer overflow.\n", __FILE__, __LINE__, sd->num_steps );
		stop_erase = PREV( stop_erase );
		sd->erase_hstep = sd->disp_pts[ stop_erase ].hstep;
		break;
	    }
	    
	    if( sd->disp_pts[ stop_erase ].star == DISP_PT_UNUSED )
		continue;

	    if( sd->disp_pts[ stop_erase ].hstep <= sd->erase_hstep )
		continue;

	    break;
	}
	

	/* erase the stars */
	if( stop_erase != sd->next_erase )
	{
	    num_pts = 0;
	    for (pp = sd->next_erase; pp != stop_erase; pp = NEXT(pp) )
	    {
		if( sd->disp_pts[pp].star == DISP_PT_UNUSED )
		    continue;
	    
		sd->points[num_pts++] = sd->disp_pts[pp].pt;

		/* the following is not needed and screws up the redraw code */
/*		sd->disp_pts[pp].star = DISP_PT_UNUSED; */
		
		/*
		 * delete this point out of the hash table
		 */
	    {
		int		hash_loc;
		int		hashed_idx;
		int		found = FALSE;
		int		xpt = sd->disp_pts[pp].pt.x;
		int		ypt = sd->disp_pts[pp].pt.y;
		
		
		for( hash_loc = PT_HASH( xpt, ypt ),
		    hashed_idx = sd->hash_index[ hash_loc ];
		    
		    hashed_idx != HASH_UNUSED;
		    
		    hash_loc = NEXT_H( hash_loc ),
		    hashed_idx = sd->hash_index[ hash_loc ] )
		{
		    if( hashed_idx == HASH_SEARCH )
			continue;
		    
		    if( sd->disp_pts[ hashed_idx ].pt.x != xpt
		       || sd->disp_pts[ hashed_idx ].pt.y != ypt
		       )
			continue;
		    
		    /* we have found the entry, now delete it */
		    found = TRUE;
		    if( sd->hash_index[ NEXT_H( hash_loc ) ] == HASH_UNUSED )
		    {
			do
			{
			    sd->hash_index[ hash_loc ] = HASH_UNUSED;
			    hash_loc = PREV_H( hash_loc );
			}
			while( sd->hash_index[ hash_loc ] == HASH_SEARCH );
		    }
		    else
			sd->hash_index[ hash_loc ] = HASH_SEARCH;
			
			
		    break;
		}

		if( !found && verbose > 0 )
		    fprintf( stderr, "%s:%d step%7d  Error: (%d,%d) was not found in the hash table!  hash_loc=%d\n", __FILE__, __LINE__, sd->num_steps, xpt, ypt, PT_HASH( xpt, ypt ) );
		
	    }
		
		/* see if the buffer is too full */
		if( num_pts >= sd->max_points )
		{
		    XDrawPoints( display->dpy, display->win, display->erase_gc,
				sd->points, num_pts, CoordModeOrigin );

		    num_pts = 0;
		}
	    }
	    
	    if( num_pts )
	    {
		XDrawPoints( display->dpy, display->win, display->erase_gc,
			    sd->points, num_pts, CoordModeOrigin );
		num_pts = 0;
	    }

	    if( sd->next_erase == sd->next_disp )
		sd->next_disp = pp;
	    sd->next_erase = pp;

	    free_pts = DIFF( sd->next_erase, sd->next_free );
	    if( verbose > 0 &&  free_pts < sp->live_stars*2+2 )
		fprintf( stderr, "%s:%d step%7d  free_pts=%d  min=%d  live_stars=%d\n", __FILE__, __LINE__, sd->num_steps, free_pts, sp->live_stars*2+2, sp->live_stars );
	}


	
	/*
	 * display the new points
	 */
	pts_to_plot = DIFF( sd->next_free, sd->next_disp );
	if( pts_to_plot > 0 )
	{
	    
	    for( i = 0; i < NUM_COLORS; i = next_color )
	    {
		next_color = NUM_COLORS;
		
		num_pts = 0;
		for (pp = sd->next_disp; pp != sd->next_free; pp = NEXT(pp) )
		{
		    if( sd->disp_pts[pp].star == DISP_PT_UNUSED )
			continue;
		    
		    if( rotate_colors )
		    {
			int	color = sd->disp_pts[pp].color;

			if( i == color )
			    sd->points[num_pts++] = sd->disp_pts[pp].pt;
			else if( i < color && color < next_color )
			    next_color = color;
		    }
		    else if( multi_colors )
		    {
			int	color = sd->star_color[ sd->disp_pts[pp].star ];
			if( i == color )
			    sd->points[num_pts++] = sd->disp_pts[pp].pt;
			else if( i < color && color < next_color )
			    next_color = color;
		    }
		    else
		    {
			sd->points[num_pts++] = sd->disp_pts[pp].pt;
		    }
		    
		    
		    if( num_pts >= sd->max_points )
		    {
			if( rotate_colors || multi_colors )
			    XDrawPoints( display->dpy, display->win,
					color_gcs[i], sd->points, num_pts,
					CoordModeOrigin );
			else 
			    XDrawPoints( display->dpy, display->win,
					display->star_gc, sd->points, num_pts,
					CoordModeOrigin );
			
			num_pts = 0;
		    }
		}
		
		if( num_pts )
		{
		    if( rotate_colors || multi_colors )
			XDrawPoints( display->dpy, display->win,
				    color_gcs[i], sd->points, num_pts,
				    CoordModeOrigin );
		    else 
			XDrawPoints( display->dpy, display->win,
				    display->star_gc, sd->points, num_pts,
				    CoordModeOrigin );
		    num_pts = 0;
		}
	    }
	    sd->next_disp = sd->next_free;
	    
	    
	    /* if we're not generating much data, then force it out quickly */
	    if( sd->num_disp_skipped > 5 + sd->buffer_factor
	       || (pts_to_plot >= 4*sd->num_visible
		   && sd->num_disp_skipped*4 > sd->buffer_factor
		   )
	       )
		XFlush( display->dpy );
	    
	    sd->points_plotted += pts_to_plot;
	    
	}

	sd->num_disp_skipped = 0;
	sd->updt_hstep = sd->num_steps;
    }	


    else
    {
	sd->erase_disps++;
	
	if( rotate_colors )
	    XDrawPoints(display->dpy, display->win,
			color_gcs[ sd->color_number ],
			sd->points, sd->points_used, CoordModeOrigin );
	else if( multi_colors )
	{
	    int		num_pts, pp, next_color;
	    int		i;
	    
	    
	    for( i = 0; i < NUM_COLORS; i = next_color )
	    {
		next_color = NUM_COLORS;
		
		num_pts = 0;
		for (pp = 0; pp < sd->points_used; pp++)
		{
		    if( i == sd->pixels[pp] )
			sd->tmp_pts[num_pts++] = sd->points[pp];
		    else if( i < sd->pixels[pp]
			    && sd->pixels[pp] < next_color )
			next_color = sd->pixels[pp];
		}
		
		if( num_pts )
		{
		    XDrawPoints(display->dpy, display->win,
				color_gcs[ i ],
				sd->tmp_pts, num_pts, CoordModeOrigin );
		    num_pts = 0;
		}
	    }
	}
	else
	    XDrawPoints(display->dpy, display->win,
			display->star_gc,
			sd->points, sd->points_used, CoordModeOrigin );
	    
	
	/* if we're not generating much data, then force it out quickly */
	if( sd->num_disp_skipped > 5 + sd->buffer_factor
	   || (sd->points_used >= 5*sd->num_visible
	       && sd->num_disp_skipped*4 > sd->buffer_factor
	       )
	   )
	    XFlush( display->dpy );
	
	sd->points_plotted += sd->points_used;
	sd->points_used = 0;
	sd->num_disp_skipped = 0;
    }
}



/* clear out the old data */
void clear_disp_pt( star_disp_type *sd )
{
    int		i;

    if( num_disp_pt )
    {
	/* mark the entries as unused */
	for( i = 0; i < num_disp_pt; i++ )
	    sd->disp_pts[i].star = DISP_PT_UNUSED;

	for( i = 0; i < HASH_TABLE_SIZE; i++ )
	    sd->hash_index[i] = HASH_UNUSED;
    }
    else
	sd->points_used = 0;

    sd->num_disp_skipped = 0;
    sd->num_poll_skipped = 0;
    sd->num_seen = 1;

    sd->points_plotted = 0;
    sd->total_points = 0;
}




/* redraw the screen from scratch */
void redraw_screen( point_2d *cur, double *m, sys_param_type *sp, star_disp_type *sd, disp *display )
{
    int		i;
    
    int		num_pts, pp, next_color;
    

    /* redraw the collapsars */
    plot_collapsars( cur, m, sp, sd, display );
    
	    
    /* redraw the stars */
    if( !num_disp_pt )
	return;
    
    for( i = 0; i < NUM_COLORS; i = next_color )
    {
	next_color = NUM_COLORS;
	
	num_pts = 0;
	for (pp = sd->next_erase; pp != sd->next_free; pp = NEXT(pp) )
	{
	    XPoint	*pt = &sd->disp_pts[pp].pt;

	    if( sd->disp_pts[pp].star == DISP_PT_UNUSED )
		continue;
		    
	    if( pt->x < sd->redraw_min_x || pt->x > sd->redraw_max_x
	       || pt->y < sd->redraw_min_y || pt->y > sd->redraw_max_y
	       )
		continue;

	    
	    if( rotate_colors )
	    {
		int	color = sd->disp_pts[pp].color;
		
		if( i == color )
		    sd->points[num_pts++] = sd->disp_pts[pp].pt;
		else if( i < color && color < next_color )
		    next_color = color;
	    }
	    else if( multi_colors )
	    {
		int	color = sd->star_color[ sd->disp_pts[pp].star ];

		if( i == color )
		    sd->points[num_pts++] = sd->disp_pts[pp].pt;
		else if( i < color && color < next_color )
		    next_color = color;
	    }
	    else
		sd->points[num_pts++] = sd->disp_pts[pp].pt;
	    
	    
	    if( num_pts >= sd->max_points )
	    {
		if( rotate_colors || multi_colors )
		    XDrawPoints( display->dpy, display->win,
				color_gcs[i], sd->points, num_pts,
				CoordModeOrigin );
		else 
		    XDrawPoints( display->dpy, display->win,
				display->star_gc, sd->points, num_pts,
				CoordModeOrigin );
		
		num_pts = 0;
	    }
	}
	
	if( num_pts )
	{
	    if( rotate_colors || multi_colors )
		XDrawPoints( display->dpy, display->win,
			    color_gcs[i], sd->points, num_pts,
			    CoordModeOrigin );
	    else 
		XDrawPoints( display->dpy, display->win,
			    display->star_gc, sd->points, num_pts,
			    CoordModeOrigin );
	    num_pts = 0;
	}
    }
    
    XFlush( display->dpy );
}


void set_redraw_full( star_disp_type *sd )
{
    sd->redraw_min_x = min_x + center_x;
    sd->redraw_max_x = max_x + center_x;
    sd->redraw_min_y = min_y + center_y;
    sd->redraw_max_y = max_y + center_y;
}

END-OF-FILE

if [ "$filename" != "/dev/null" ]
then
  size=`wc -c < $filename`

  if [ $size != 10691 ]
  then
    echo $filename changed - should be 10691 bytes, not $size bytes
  fi

  chmod 660 $filename
fi

# ---------- file plot_collapsars.c ----------

filename="plot_collapsars.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file plot_collapsars.c...
fi

cat << 'END-OF-FILE' > $filename
#include "xstar.h"
#include "xstar_ext.h"

/*
** XStar -- a animated n-body solver for X windows
** Copyright (C) 1995  Wayne Schlitt (wayne@cse.unl.edu)
**
**  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 2 of the License, or
**  (at your option) any later version.
*/



/* plot the collapsars on the screen */
void plot_collapsars( point_2d *cur, double *m, sys_param_type *sp, star_disp_type *sd, disp *display )
{

    int		i;

    int		x, y;

    int		num_pts;

    XPoint	*disp_buf;


    /* find a good scratch array */
    if( num_disp_pt )
	disp_buf = sd->points;
    else
	disp_buf = sd->tmp_pts;
    
    

    /* draw the collapsars */
    
    for( i = 0; i < max_stars; i++ )
    {
	if( m[i] < collapsar )
	    continue;

	
	num_pts = 0;
	for( x = -1; x < 2; x++ )
	    for( y = -1; y < 2; y++ )
	    {
		disp_buf[num_pts].x = cur[i].x + center_x + x;
		disp_buf[num_pts].y = cur[i].y + center_y + y;
		num_pts++;
		
	    }
	
	if( rotate_colors )
	    XDrawPoints( display->dpy, display->win,
			color_gcs[ sd->color_number ],
			disp_buf, num_pts,
			CoordModeOrigin );
	else if( multi_colors )
	    XDrawPoints( display->dpy, display->win,
			color_gcs[ sd->star_color[ i ]],
			disp_buf, num_pts,
			CoordModeOrigin );
	else 
	    XDrawPoints( display->dpy, display->win,
			display->star_gc, disp_buf, num_pts,
			CoordModeOrigin );
	
    }
}

END-OF-FILE

if [ "$filename" != "/dev/null" ]
then
  size=`wc -c < $filename`

  if [ $size != 1509 ]
  then
    echo $filename changed - should be 1509 bytes, not $size bytes
  fi

  chmod 660 $filename
fi

echo end of this archive file....
exit 0
--- cut here ---
-- 
Wayne Schlitt can not assert the truth of all statements in this
article and still be consistent.