Article 12757 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    Part05/06
Date: Sun, 27 Aug 1995 23:14:50 GMT
Organization: The Backbone Cabal
Lines: 2816
Sender: wayne@backbone (Wayne Schlitt)
Message-ID: <WAYNE.95Aug27171450@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/part05


--- cut here ---

# Continuation of Shell Archive, created by backbone!wayne

# This is part 5

# 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 reset_hstep.c ----------

filename="reset_hstep.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file reset_hstep.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.
*/


/*
 * reset the hstep values to keep from overflowing the small hstep variable
 * in star_disp_type.  
 */

void reset_hstep( star_disp_type *sd )
{
    int		i;
    int		min_hstep;
    int		hstep_used;


    if( !num_disp_pt )
	return;
    
    /*
     * find the minimum hstep
     */
    min_hstep = sd->hsteps;

    for( i = 0; i < num_disp_pt; i++ )
    {
	if( sd->disp_pts[i].star == DISP_PT_UNUSED )
	    continue;
	
	if( sd->disp_pts[i].hstep < min_hstep )
	    min_hstep = sd->disp_pts[i].hstep;

    }

    if( sd->erase_hstep < min_hstep )
	min_hstep = sd->erase_hstep;


    /* make sure we won't have to call reset_hstep for at least 1k steps */
    hstep_used = sd->hsteps - sd->erase_hstep;
    if( hstep_used > MAX_HSTEP - 1024 )
	sd->erase_hstep += hstep_used - (MAX_HSTEP - 1024);
	

    /*
     * reset the hstep variables
     */
    sd->hsteps -= min_hstep;
    sd->raw_hsteps -= min_hstep / sd->hstep_scale;
    sd->erase_hstep -= min_hstep;

    for( i = 0; i < num_disp_pt; i++ )
    {
	if( sd->disp_pts[i].star == DISP_PT_UNUSED )
	    continue;
	
	sd->disp_pts[i].hstep -= min_hstep;
    }
}


END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file set_star_disp.c ----------

filename="set_star_disp.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file set_star_disp.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.
*/



/*
 * Set up the star display variables
 */

void set_star_disp( star_disp_type *sd, sys_param_type *sp )
{
    /* set up default erase_hstep adjust points */
    sd->num_disp_pt_32 = num_disp_pt / 32;
    sd->num_disp_pt_16 = num_disp_pt / 16;
    sd->num_disp_pt_8 = num_disp_pt / 8;
    sd->num_disp_pt_4 = num_disp_pt / 4;

    /* handle special cases */
    if( sp->do_bounce && sp->star_line )
	sd->hsteps = sd->num_steps = MAX_HSTEP/2;
	   
    else if( sp->do_bounce )
    {
	if( sp->no_speed )
	    if( sp->num_collapsar == 0 )
		sd->hsteps = sd->num_steps = 400 / sp->live_stars + 20;
	    else
		sd->hsteps = sd->num_steps = 20;
	else
	    sd->hsteps = sd->num_steps = 40 / sp->live_stars + 15;

	sd->num_disp_pt_32 = num_disp_pt * .875 + num_disp_pt / (32*4);
	sd->num_disp_pt_16 = num_disp_pt * .875 + num_disp_pt / (16*4);
	sd->num_disp_pt_8 = num_disp_pt * .875 + num_disp_pt / (8*4);
	sd->num_disp_pt_4 = num_disp_pt * .875 + num_disp_pt / (4*4);
    }
    else if( sp->few_stars )
    {
	sd->hsteps = sd->num_steps = 15;

	sd->num_disp_pt_32 = num_disp_pt * .75 + num_disp_pt / (32*4);
	sd->num_disp_pt_16 = num_disp_pt * .75 + num_disp_pt / (16*4);
	sd->num_disp_pt_8 = num_disp_pt * .75 + num_disp_pt / (8*4);
	sd->num_disp_pt_4 = num_disp_pt * .75 + num_disp_pt / (4*4);
    }
    else if( sp->star_circle )
	sd->hsteps = sd->num_steps = 400;

    else
	sd->hsteps = sd->num_steps = 200;


    /* if we don't have many display points, then let the tails grow slowly */
    if( sd->num_disp_pt_4 - sd->num_disp_pt_8 < 1024 )
	sd->num_disp_pt_4 += 1024 - (sd->num_disp_pt_4 - sd->num_disp_pt_8);
    
    
    /* set up other misc star_disp variables */
    sd->erase_disps = 0;
    sd->live_erase = sp->live_stars;

    sd->next_free = 1;
    sd->next_disp = 1;
    sd->next_erase = 0;
    sd->erase_hstep = 0;
    sd->updt_hstep = sd->hsteps;

    sd->hstep_scale = fv_inv * (1./4);

    sd->raw_num_steps = sd->num_steps / sd->hstep_scale;
    sd->raw_hsteps = sd->hsteps / sd->hstep_scale;
}



void set_buffer_factor( star_disp_type *sd, sys_param_type *sp )
{
    int		total_stars = sp->live_stars + sp->num_collapsar;
    
    sd->buffer_factor = max_disp_skip - total_stars*total_stars;
    if( sd->buffer_factor < 1 )
	sd->buffer_factor = 1;
}

END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file move_stars_taylor3.c ----------

filename="move_stars_taylor3.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file move_stars_taylor3.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.
*/



/*
 * taylor series expansion to the 2nd order with the 3rd order approximated
 * by using the derivative the the 2nd degree Lagrange polynomial.
 */


static int 	num_rk4;

void restart_taylor3( int n, double *m, point_2d *x, point_2d *v )
{
    num_rk4 = 3;
}



void init_taylor3( int n, double *m, point_2d *x, point_2d *v )
{
    init_rk4( n, m, x, v );
    
    restart_taylor3( n, m, x, v );
}



void move_stars_taylor3( point_2d *prev, point_2d *cur, double *m, 
		point_2d *vk[], point_2d *ak[],
		int max_stars, sys_param_type *sp, star_disp_type *sd )
{
    double		m_i, m_j;
    double		previ_x, previ_y;
    double		ai_x, ai_y;
    double		Dt_a;		/* derivative of the acceleration */

    register int        i, j, k;	/* star index */
    int			tx, ty;

    int			xpt, ypt;
    
    point_2d	*v_0 = vk[0];
    point_2d	*v_1 = vk[1];
    point_2d	*a_0 = ak[0];
    point_2d	*a_1 = ak[1];
    point_2d	*a_2 = ak[2];



    /*
     * we must use some other, self starting method for the first 3
     * calls to the taylor method
     */
    if( num_rk4 )
    {
	--num_rk4;
	move_stars_rk4( prev, cur, m, vk, ak, max_stars, sp, sd );
	return;
    }


    for( i = 0; i < max_stars; i++ )
    {
	m_i = m[i];
	if( m_i <= 0.0 )
	    continue;
	
	previ_x = prev[i].x;
	previ_y = prev[i].y;
	ai_x = a_0[i].x;
	ai_y = a_0[i].y;
	
	for( j = i+1; j < max_stars; j++ )
	{
	    double dx, dy, n2, f;
	    register double f_d;
	    
	    m_j = m[j];
	    if( m_j <= 0.0 )
		continue;
	    
	    dx = prev[j].x - previ_x;
	    dy = prev[j].y - previ_y;
	    n2 = dx*dx + dy*dy;
	    
	    if( n2 < collide_dist )
	    {
		a_0[i].x = ai_x;
		a_0[i].y = ai_y;
		collide( i, j, cur, prev, m, vk, ak, sp, sd );
		m_i = m[i];
		m_j = m[j];
		previ_x = prev[i].x;
		previ_y = prev[i].y;
		ai_x = a_0[i].x;
		ai_y = a_0[i].y;
	    }
	    else
	    {
		f = G / (sqrt(n2) * n2);
		
		f_d = dx * f;		
		ai_x += f_d * m_j;
		a_0[j].x -= f_d * m_i;
		
		f_d = dy * f;
		ai_y += f_d * m_j;
		a_0[j].y -= f_d * m_i;
	    }
	}
	
	/*
	 * Move the star
	 *
	 * You want each step to be as small as possible and this is done
	 * by making FV as large as possible.  If the steps are small
	 * enough, then we are approximating infinitesimal movements.
	 *
	 * The movement is based off of the formula:  (with t=1)
	 * x(t) = x0 + v t + 1/2 a t^2 + 1/(2*3) (Dt a) t^3
	 * v(t) = v0 + a t + 1/2 (Dt a) t^2
	 *
	 * In theory, these formulas should include an infinite number
	 * of derivatives of x(t), but the additional terms are not
	 * easy to calculate.  Even Dt a is just an estimate.
	 *
	 *
	 * note:
	 * I am currently using the derivative of the 2nd degree
	 * Lagrange polynomial to estimate Dt_a.  I tried using the
	 * 3rd degree Lagrange Polynomial, but it didn't work as well.
	 * Polynomials of high degree tend to not be very smooth so
	 * their derivatives can make things worse.
	 */
	
	if( m_i >= collapsar )
	{
	    /* I think it looks better if the collapsars don't move ;-) */
	    v_1[i].x = v_1[i].y = 0.0;
	    
	    a_0[i].x = a_0[i].y = 0.0;
	    ai_x = ai_y = 0.0;
	    a_1[i].x = a_1[i].y = 0.0;
	    
	    xpt = cur[i].x = prev[i].x;
	    ypt = cur[i].y = prev[i].y;
	}
	else
	{
	    a_0[i].x = ai_x;
/*	    Dt_a = (11./12)*ai_x - 1.5*a_1[i].x + .75*a_2[i].x
		- (1./6)*a_3[i].x; */
	    Dt_a = (3./4)*ai_x - a_1[i].x + (1./4)*a_2[i].x;
	    v_0[i].x = v_1[i].x + ai_x + Dt_a;
	    xpt = cur[i].x = previ_x + (v_1[i].x + ((1./2)*ai_x
						    + (1./3)*Dt_a));
	    
	    
	    a_0[i].y = ai_y;
/*	    Dt_a = (11./12)*ai_y - 1.5*a_1[i].y + .75*a_2[i].y
		- (1./6)*a_3[i].y; */
	    Dt_a = (3./4)*ai_y - a_1[i].y + (1./4)*a_2[i].y;
	    v_0[i].y = v_1[i].y + ai_y + Dt_a;
	    ypt = cur[i].y = previ_y + (v_1[i].y + ((1./2)*ai_y
						    + (1./3)*Dt_a));
	}
	ak[K-1][i].x = ak[K-1][i].y = 0;
	
#	include "plot_star.h"	
	
    }
}
END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file move_stars_rk4.c ----------

filename="move_stars_rk4.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file move_stars_rk4.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.
*/



/*
 * Runge-Kutta method of 4th order
 */

static point_2d	*ta_1, *ta_2, *ta_3;

static int	first_time = 1;

void restart_rk4( int n, double *m, point_2d *x, point_2d *v )
{
    return;
}


void init_rk4( int n, double *m, point_2d *x, point_2d *v )
{

    int		i;

    static int	last_n = -1;
    

    default_init( n, m, x, v );
    

    if( n != last_n )
    {
	last_n = n;
	
	/* allocate storage */
	if( !first_time )
	{
	    free( ta_1 );
	    free( ta_2 );
	    free( ta_3 );
	}
	
	ta_1 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	ta_2 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	ta_3 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );

	for( i = 0; i < n; i++ )
	{
	    ta_1[i].x = ta_1[i].y = 0;
	    ta_2[i].x = ta_2[i].y = 0;
	    ta_3[i].x = ta_3[i].y = 0;
	}
    }

    first_time = 0;
	
}


#define collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd ) do_collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd, __FILE__, __LINE__ )

static void do_collide_rk4( int i, int j,
			   point_2d *cur, point_2d *prev,
			   double *m,
			   point_2d *vk[], point_2d *ak[],
			   point_2d *ta_1, point_2d *ta_2, point_2d *ta_3,
			   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;
	}
	ta_1[i].x = (m[i]*ta_1[i].x + m[j]*ta_1[j].x)*tmass;
	ta_1[i].y = (m[i]*ta_1[i].y + m[j]*ta_1[j].y)*tmass;
	ta_2[i].x = (m[i]*ta_2[i].x + m[j]*ta_2[j].x)*tmass;
	ta_2[i].y = (m[i]*ta_2[i].y + m[j]*ta_2[j].y)*tmass;
	ta_3[i].x = (m[i]*ta_3[i].x + m[j]*ta_3[j].x)*tmass;
	ta_3[i].y = (m[i]*ta_3[i].y + m[j]*ta_3[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;
    }	    
    ta_1[j].x = ta_1[j].y = 0;
    ta_2[j].x = ta_2[j].y = 0;
    ta_3[j].x = ta_3[j].y = 0;

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



void move_stars_rk4( point_2d *prev, point_2d *cur, double *m, 
		point_2d *vk[], point_2d *ak[],
		int max_stars, sys_param_type *sp, star_disp_type *sd )
{
    double		m_i, m_j;
    double		previ_x, previ_y;
    double		vi_x, vi_y;
    double		a0i_x, a0i_y;
    double		a1i_x, a1i_y;
    double		ai_x, ai_y;
    double		a_sum1, a_sum2;

    register int        i, j, k;	/* star index */
    int			tx, ty;

    int			xpt, ypt;

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

    /*
     * initialize the ta_[1-3] arrays
     */
    for( i = 0; i < max_stars; i++ )
    {
	ta_1[i].x = ta_1[i].y = 0;
	ta_2[i].x = ta_2[i].y = 0;
	ta_3[i].x = ta_3[i].y = 0;
    }
    


    /*
     * move the stars
     */
    for( i = 0; i < max_stars; i++ )
    {
	m_i = m[i];
	if( m_i <= 0.0 )
	    continue;
	
	previ_x = prev[i].x;
	previ_y = prev[i].y;

	/*
	 * calculate the a_0 values
	 */
	ai_x = a_0[i].x;
	ai_y = a_0[i].y;
	for( j = i+1; j < max_stars; j++ )
	{
	    double dx, dy, n2, f;
	    register double f_d;
	    
	    m_j = m[j];
	    if( m_j <= 0.0 )
		continue;
	    
	    dx = prev[j].x - previ_x;
	    dy = prev[j].y - previ_y;
	    n2 = dx*dx + dy*dy;
	    
	    if( n2 < collide_dist )
	    {
		a_0[i].x = ai_x;
		a_0[i].y = ai_y;
		collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd );
		m_i = m[i];
		m_j = m[j];
		previ_x = prev[i].x;
		previ_y = prev[i].y;
		ai_x = a_0[i].x;
		ai_y = a_0[i].y;
	    }
	    else
	    {
		f = G / (sqrt(n2) * n2);
		
		f_d = dx * f;		
		ai_x += f_d * m_j;
		a_0[j].x -= f_d * m_i;
		
		f_d = dy * f;
		ai_y += f_d * m_j;
		a_0[j].y -= f_d * m_i;
	    }
	}
	a_0[i].x = ai_x;
	a_0[i].y = ai_y;
	
	
	/*
	 * calculate the ta_1 values
	 */
	ai_x = ta_1[i].x;
	ai_y = ta_1[i].y;
	vi_x = v_1[i].x;
	vi_y = v_1[i].y;

	for( j = i+1; j < max_stars; j++ )
	{
	    double dx, dy, n2, f;
	    register double f_d;
	    
	    m_j = m[j];
	    if( m_j <= 0.0 )
		continue;
	    
	    dx = prev[j].x - previ_x + .5*(v_1[j].x - vi_x);
	    dy = prev[j].y - previ_y + .5*(v_1[j].y - vi_y);

	    n2 = dx*dx + dy*dy;
	    
	    if( n2 < collide_dist )
	    {
		ta_1[i].x = ai_x;
		ta_1[i].y = ai_y;
		collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd );
		m_i = m[i];
		m_j = m[j];
		previ_x = prev[i].x;
		previ_y = prev[i].y;
		vi_x = v_1[i].x;
		vi_y = v_1[i].y;
		ai_x = ta_1[i].x;
		ai_y = ta_1[i].y;
	    }
	    else
	    {
		f = G / (sqrt(n2) * n2);
		
		f_d = dx * f;		
		ai_x += f_d * m_j;
		ta_1[j].x -= f_d * m_i;
		
		f_d = dy * f;
		ai_y += f_d * m_j;
		ta_1[j].y -= f_d * m_i;
	    }
	}
	ta_1[i].x = ai_x;
	ta_1[i].y = ai_y;
    }
    
    
    for( i = 0; i < max_stars; i++ )
    {
	m_i = m[i];
	if( m_i <= 0.0 )
	    continue;
	
	previ_x = prev[i].x;
	previ_y = prev[i].y;

	/*
	 * calculate the ta_2 values
	 */
	vi_x = v_1[i].x;
	vi_y = v_1[i].y;
	a0i_x = a_0[i].x;
	a0i_y = a_0[i].y;
	ai_x = ta_2[i].x;
	ai_y = ta_2[i].y;
	
	for( j = i+1; j < max_stars; j++ )
	{
	    double dx, dy, n2, f;
	    register double f_d;
	    
	    m_j = m[j];
	    if( m_j <= 0.0 )
		continue;
	    
	    dx = prev[j].x - previ_x + .5*(v_1[j].x - vi_x)
		+ .25*(a_0[j].x - a0i_x);
	    dy = prev[j].y - previ_y + .5*(v_1[j].y - vi_y)
		+ .25*(a_0[j].y - a0i_y);

	    n2 = dx*dx + dy*dy;
	    
	    if( n2 < collide_dist )
	    {
		ta_2[i].x = ai_x;
		ta_2[i].y = ai_y;
		collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd );
		m_i = m[i];
		m_j = m[j];
		previ_x = prev[i].x;
		previ_y = prev[i].y;
		vi_x = v_1[i].x;
		vi_y = v_1[i].y;
		a0i_x = a_0[i].x;
		a0i_y = a_0[i].y;
		ai_x = ta_2[i].x;
		ai_y = ta_2[i].y;
	    }
	    else
	    {
		f = G / (sqrt(n2) * n2);
		
		f_d = dx * f;		
		ai_x += f_d * m_j;
		ta_2[j].x -= f_d * m_i;
		
		f_d = dy * f;
		ai_y += f_d * m_j;
		ta_2[j].y -= f_d * m_i;
	    }
	}
	ta_2[i].x = ai_x;
	ta_2[i].y = ai_y;


	/*
	 * calculate the ta_3 values
	 */
	a1i_x = ta_1[i].x;
	a1i_y = ta_1[i].y;
	ai_x = ta_3[i].x;
	ai_y = ta_3[i].y;
	
	for( j = i+1; j < max_stars; j++ )
	{
	    double dx, dy, n2, f;
	    register double f_d;
	    
	    m_j = m[j];
	    if( m_j <= 0.0 )
		continue;
	    
	    dx = prev[j].x - previ_x + v_1[j].x - vi_x + .5*(ta_1[j].x - a1i_x);
	    dy = prev[j].y - previ_y + v_1[j].y - vi_y + .5*(ta_1[j].y - a1i_y);

	    n2 = dx*dx + dy*dy;
	    
	    if( n2 < collide_dist )
	    {
		ta_3[i].x = ai_x;
		ta_3[i].y = ai_y;
		collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd );
		m_i = m[i];
		m_j = m[j];
		previ_x = prev[i].x;
		previ_y = prev[i].y;
		vi_x = v_1[i].x;
		vi_y = v_1[i].y;
		a1i_x = ta_1[i].x;
		a1i_y = ta_1[i].y;
		ai_x = ta_3[i].x;
		ai_y = ta_3[i].y;
	    }
	    else
	    {
		f = G / (sqrt(n2) * n2);
		
		f_d = dx * f;		
		ai_x += f_d * m_j;
		ta_3[j].x -= f_d * m_i;
		
		f_d = dy * f;
		ai_y += f_d * m_j;
		ta_3[j].y -= f_d * m_i;
	    }
	}
	ta_3[i].x = ai_x;
	ta_3[i].y = ai_y;
	
	
	/*
	 * Move the star
	 */
	
	if( m_i >= collapsar )
	{
	    /* I think it looks better if the collapsars don't move ;-) */
	    v_1[i].x = v_1[i].y = 0.0;
	    
	    a_0[i].x = a_0[i].y = 0.0;
	    ai_x = ai_y = 0.0;
	    ta_1[i].x = ta_1[i].y = 0.0;
	    
	    xpt = cur[i].x = prev[i].x;
	    ypt = cur[i].y = prev[i].y;
	}
	else
	{
	    a_sum1 = ta_1[i].x + ta_2[i].x;
	    a_sum2 = (1./6)*(a_0[i].x + a_sum1);
	    xpt = cur[i].x = previ_x + (v_1[i].x + a_sum2);
	    v_0[i].x = v_1[i].x + a_sum2 + (1./6)*(a_sum1 + ai_x);
	    
	    a_sum1 = ta_1[i].y + ta_2[i].y;
	    a_sum2 = (1./6)*(a_0[i].y + a_sum1);
	    ypt = cur[i].y = previ_y + (v_1[i].y + a_sum2);
	    v_0[i].y = v_1[i].y + a_sum2 + (1./6)*(a_sum1 + ai_y);
	}
	ak[K-1][i].x = ak[K-1][i].y = 0;
	
#	include "plot_star.h"	
	
    }
}
END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file move_stars_gpemce8.c ----------

filename="move_stars_gpemce8.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file move_stars_gpemce8.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.
*/



/*
 * Gragg's method of Polynomial-Extrapolation with Modifications for
 * Conservation of Energy.  (i.e., this is a discrete mechanics method.)
 *
 * This method is from book _Numerical_Solutions_of_the_N-body_Problem_ by
 * Andrzej Marciniak.  This appears to be his favorite method.
 */

#define H	(1.)			/* step size			*/
#define P	(4)			/* half the order of the approx */
#define MAX_SKIP (32)			/* max # of correction skips	*/


/*
 * Notes:
 *
 * H should not be modified.  To change the step size, modify fv instead.
 *
 * P is one half the order of the approximation.  With P=4, FV=.25 you get
 * a very accurate approximation.  With P=2, FV=.25, you get an approximation
 * that is slightly worse than the taylor approximation with FV=2, but the
 * taylor method is much faster.  According to Andrzej's book, going beyond
 * P=4 won't get you much additional precision due to rounding errors.
 */


static int	ibeta[P+1];
static int	ibeta_2[P+1];
static double	beta[P+1];
static double	beta_2_inv[P+1];
static double	beta_inv[P+1];

static double	e0_2, a1, orig_a1, last_a1;
static double	alpha1[P+1], alpha2[P+1], bk[P+1], f[P], f1[P];
static point_2d	*x1, *y, *z, *txk[P+1];
static point_2d	*v1, *u, *w, *tvk[P+1];
static point_2d	*sum2, *sum3;
static double	*sum, *sum1;

static int	first_time = 1;

static int	eps_skip;
static int	max_eps_skip;


void restart_gpemce8( int n, double *m, point_2d *x, point_2d *v )
{
    return;
}

void init_gpemce8( int n, double *m, point_2d *x, point_2d *v )
{

    int		i, j, k, l;
    double	a, c;
    double	sum4, sum5;

    static int	last_n = -1;
    

    default_init( n, m, x, v );
    

    if( n != last_n )
    {
	last_n = n;
	
	/* allocate storage */
	if( !first_time )
	{
	    free( x1 );
	    free( y );
	    free( z );
	    free( v1 );
	    free( u );
	    free( w );
	    free( sum );
	    free( sum1 );
	    free( sum2 );
	    free( sum3 );
	    
	    for( k = 0; k < P+1; k++ )
	    {
		free( txk[k] );
		free( tvk[k] );
	    }
	}
	
	x1 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	y = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	z = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	v1 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	u = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	w = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	sum = (double *)Safe_Malloc( sizeof( double ) * n );
	sum1 = (double *)Safe_Malloc( sizeof( double ) * n );
	sum2 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	sum3 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	
	for( k = 0; k < P+1; k++ )
	{
	    txk[k] = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	    tvk[k] = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	}
    }
    

    /* calculate the total energy of the system (times 2) */
    e0_2 = 0;
    for( i = 0; i < n; i++ )
    {
	if( m[i] <= 0. )
	    continue;
	
	sum4 = 0;
	
	for( j = 0; j < n; j++ )
	{
	    if( m[j] <= 0. )
		continue;
	
	    if( i == j )
		continue;

	    
	    sum4 += m[j] / sqrt( (x[i].x - x[j].x)*(x[i].x - x[j].x)
				+ ( x[i].y - x[j].y )*( x[i].y - x[j].y ));
	}

	e0_2 += m[i]*( v[i].x*v[i].x + v[i].y*v[i].y - G*sum4 );
    }


    if( first_time )
    {
	/* step 1 */
	
	ibeta[0] = 1;
	ibeta[1] = 2;
	ibeta[2] = 3;
	for( k = 3; k < P+1; k++ )
	    ibeta[k] = 2*ibeta[k-2];

	for( k = 0; k < P+1; k++ )
	{
	    ibeta_2[k] = 2*ibeta[k];
	    beta[k] = ibeta[k];
	    beta_inv[k] = 1. / beta[k];
	    beta_2_inv[k] = .5 * beta_inv[k];
	}
	

	/* step 2 */
	for( k = 0; k < P+1; k++)
	{
	    bk[k] = 1;
	    
	    for( l = 0; l < P+1; l++ )
	    {
		if( k == l )
		    continue;
		
		bk[k] *= ibeta[l]*ibeta[l];
	    }
	}
	
	for( i = 0; i < P; i++ )
	    f[i] = f1[i] = 1;
	
	for( l = 0; l < P-1; l++ )
	{
	    for( i = 0; i < P-1; i++ )
	    {
		if( l == i )
		    continue;
		
		a = bk[l] - bk[i];
		f[l] *= (bk[P-1] - bk[i]) / a;
		f1[l] *= (bk[P] - bk[i]) / a;
	    }
	    
	    f[l] *= bk[P-1]/bk[l];
	    f1[l] *= bk[P]/bk[l];
	}
	
	a = c = 0;
	
	for( l = 0; l < P; l++ )
	{
	    a += f1[l];
	    c += f[l];
	}
	
	a = 1 - a;
	c = 1 - c;
	alpha1[P] = 0;
	alpha1[P-1] = 1/c;
	alpha2[P] = 1;
	alpha2[P-1] = -a/c;
	
	for( l = 0; l < P-1; l++ )
	{
	    alpha1[l] = -f[l]/c;
	    alpha2[l] = a*f[l]/c - f1[l];
	}
	
	sum4 = sum5 = 0;
	
	for( k = 0; k < P+1; k++ )
	{
	    a = pow( beta[k], 2*P );
	    sum4 += alpha1[k]/a;
	    sum5 += alpha2[k]/a;
	}
	
	orig_a1 = -sum4/sum5;
    }
    

    a1 = orig_a1;
    last_a1 = a1;
    eps_skip = -1;
    max_eps_skip = 4;
    
    first_time = 0;
	
}


void move_stars_gpemce8( point_2d *x, point_2d *x_new, double *m, 
		point_2d *vk[], point_2d *ak[],
		int max_stars, sys_param_type *sp, star_disp_type *sd )
{
    double	eps = 1./(1024. * 1024. * 1024. * 1024. );
    
    int		ibetak_2;
    double	betak_2_inv;
    double	betak_inv;
    double	betak;

    int		p1;
    int		i, j, k;
    
    double	ap, b, c, df, ff, r, r1, r2;

    int		eps_calc;

    double	dx, dy;
    double	r_inv, r3_inv;

    double	m_i;

    int		tx, ty;

    int		xpt, ypt;
    
    point_2d	*v_new = vk[0];
    point_2d	*v = vk[1];
    
    

    /* step 3 */
    for( k = 0; k < P+1; k++ )
    {
	ibetak_2 = ibeta_2[k];
	betak = beta[k];
	betak_inv = beta_inv[k];
	betak_2_inv = beta_2_inv[k];
	
	for( i = 0; i < max_stars; i++ )
	{
	    if( m[i] <= 0. )
		continue;
	    
	    /* step 4 */
	    z[i] = x[i];
	    u[i] = v[i];

	    /* step 5 */
	    y[i].x = z[i].x + H*betak_2_inv*u[i].x;
	    y[i].y = z[i].y + H*betak_2_inv*u[i].y;

	    sum2[i].x = sum2[i].y = 0;
	}

	/* step 5 (continued) */
	for( i = 0; i < max_stars; i++ )
	{
	    if( m[i] <= 0. )
		continue;
	
	    for( j = i+1; j < max_stars; j++ )
	    {
		if( m[j] <= 0. )
		    continue;
		
		dx = z[i].x - z[j].x;
		dy = z[i].y - z[j].y;

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

	    w[i].x = u[i].x - H*G*betak_2_inv*sum2[i].x;
	    w[i].y = u[i].y - H*G*betak_2_inv*sum2[i].y;
	}
    
	
	for( p1 = 2; p1 < ibetak_2; p1++ )
	{
	    for( i = 0; i < max_stars; i++ )
	    {
		if( m[i] <= 0. )
		    continue;
		
		x1[i] = y[i];
		v1[i] = w[i];

		y[i].x = z[i].x + H*betak_inv*v1[i].x;
		y[i].y = z[i].y + H*betak_inv*v1[i].y;

		sum2[i].x = sum2[i].y = 0;
	    }
	    
	    for( i = 0; i < max_stars; i++ )
	    {
		if( m[i] <= 0. )
		    continue;
		
		for( j = i+1; j < max_stars; j++ )
		{
		    if( m[j] <= 0. )
			continue;
		
		    dx = x1[i].x - x1[j].x;
		    dy = x1[i].y - x1[j].y;
		    r3_inv = 1 / (dx*dx + dy*dy);
		    r3_inv = sqrt( r3_inv ) * r3_inv;
		
		    sum2[i].x += m[j] * dx * r3_inv;
		    sum2[j].x -= m[i] * dx * r3_inv;
		    sum2[i].y += m[j] * dy * r3_inv;
		    sum2[j].y -= m[i] * dy * r3_inv;
		}
		
		w[i].x = u[i].x - H*G*betak_inv*sum2[i].x;
		w[i].y = u[i].y - H*G*betak_inv*sum2[i].y;
	    }

	    for( i = 0; i < max_stars; i++ )
	    {
		z[i] = x1[i];
		u[i] = v1[i];
	    }
	}

	
	/* step 6 */
	for( i = 0; i < max_stars; i++ )
	{
	    sum2[i].x = sum2[i].y = 0;
	    sum3[i].x = sum3[i].y = 0;
	}

	for( i = 0; i < max_stars; i++ )
	{
	    if( m[i] <= 0. )
		continue;
		
	    for( j = i+1; j < max_stars; j++ )
	    {
		if( m[j] <= 0. )
		    continue;
		
		dx = y[i].x - y[j].x;
		dy = y[i].y - y[j].y;
		r3_inv = 1 / (dx*dx + dy*dy);
		r3_inv = sqrt( r3_inv ) * r3_inv;
		
		sum2[i].x += m[j] * dx * r3_inv;
		sum2[j].x -= m[i] * dx * r3_inv;
		sum2[i].y += m[j] * dy * r3_inv;
		sum2[j].y -= m[i] * dy * r3_inv;
		
		dx = z[i].x - z[j].x + H*betak_inv*(w[i].x - w[j].x);
		dy = z[i].y - z[j].y + H*betak_inv*(w[i].y - w[j].y);
		r3_inv = 1 / (dx*dx + dy*dy);
		r3_inv = .5 * sqrt( r3_inv ) * r3_inv;
		
		sum3[i].x += m[j] * dx * r3_inv;
		sum3[j].x -= m[i] * dx * r3_inv;
		sum3[i].y += m[j] * dy * r3_inv;
		sum3[j].y -= m[i] * dy * r3_inv;
	    }

	    sum3[i].x += sum2[i].x;
	    sum3[i].y += sum2[i].y;

	    x1[i].x = (y[i].x + z[i].x)/2
		+ H*betak_2_inv*(w[i].x + u[i].x/2
				 - H*G*betak_2_inv*sum2[i].x );
	    x1[i].y = (y[i].y + z[i].y)/2
		+ H*betak_2_inv*(w[i].y + u[i].y/2
				 - H*G*betak_2_inv*sum2[i].y );

	    v1[i].x = (w[i].x + u[i].x)/2 - H*G*betak_2_inv*sum3[i].x;
	    v1[i].y = (w[i].y + u[i].y)/2 - H*G*betak_2_inv*sum3[i].y;
	}

	for( i = 0; i < max_stars; i++ )
	{
	    if( m[i] <= 0. )
		continue;
	    
	    txk[k][i] = x1[i];
	    tvk[k][i] = v1[i];
	}
    }
    
    /* step 8 */
    for( eps_calc = 0; eps_skip <= 1 && eps_calc < 16; eps_calc++ )
    {
	ap = a1;

	ff = df = 0;
	
	for( i = 0; i < max_stars; i++ )
	    sum[i] = sum1[i] = 0;

	for( i = 0; i < max_stars; i++ )
	{
	    if( m[i] <= 0. )
		continue;
	    
	    
	    for( j = i+1; j < max_stars; j++ )
	    {
		if( m[j] <= 0. )
		    continue;
		
		
		r1 = r2 = 0;
		for( k = 0; k < P+1; k++ )
		{
		    b = txk[k][i].x - txk[k][j].x;
		    r1 += (alpha1[k] + alpha2[k] * a1) * b;
		    r2 += alpha2[k] * b;
		}
		
		r = r1*r1;
		c = r1 * r2;
		
		r1 = r2 = 0;
		for( k = 0; k < P+1; k++ )
		{
		    b = txk[k][i].y - txk[k][j].y;
		    r1 += (alpha1[k] + alpha2[k] * a1) * b;
		    r2 += alpha2[k] * b;
		}
		
		r += r1*r1;
		c += r1 * r2;
		
		r3_inv = 1 / r;
		r_inv = sqrt( r3_inv );
		r3_inv *= c*r_inv;
		
		sum[i] += m[j] * r_inv;
		sum1[i] += m[j] * r3_inv;
		sum[j] += m[i] * r_inv;
		sum1[j] += m[i] * r3_inv;
	    }
	    
	    
	    r = c = 0;
	    
	    r1 = r2 = 0;
	    
	    for( k = 0; k < P+1; k++ )
	    {
		b = tvk[k][i].x;
		r1 += (alpha1[k] + alpha2[k] * a1) * b;
		r2 += alpha2[k] * b;
	    }
	    
	    r += r1*r1;
	    c += r1 * r2;
	    
	    r1 = r2 = 0;
	    
	    for( k = 0; k < P+1; k++ )
	    {
		b = tvk[k][i].y;
		r1 += (alpha1[k] + alpha2[k] * a1) * b;
		r2 += alpha2[k] * b;
	    }
	    
	    r += r1*r1;
	    c += r1 * r2;
	    
	    ff += m[i]*(r - G*sum[i]);
	    df += m[i]*(2*c + G*sum1[i]);
	}
	
	a1 -= (ff - e0_2) / df;

	if( fabs( a1 - ap ) < eps )
	    break;
    }
/*    if( eps_calc > 0 || fabs( a1 - last_a1 ) > .125*eps ) */
    if( eps_calc > 0 )
    {
	last_a1 = a1;
	if( eps_skip >= 0 )
	    max_eps_skip = (max_eps_skip + 1) >> 1;

	eps_skip = -2;
    }
    /*
     * we don't need to update a1 _every_ time, so skip a few.  updating a1
     * twice in a row is important.  It keeps the oscillations down.
     */
    eps_skip++;
    if( eps_skip > max_eps_skip )
    {
	eps_skip = 0;
	if( max_eps_skip < MAX_SKIP )
	    max_eps_skip += 1 + (max_eps_skip >> 3);
    }


    /* step 9 */
    for( i = 0; i < max_stars; i++ )
    {
	if( m[i] <= 0. || m[i] >= collapsar )
	{
	    x_new[i] = x[i];
	    v_new[i] = v[i];
	}
	else
	{
	    x_new[i].x = x_new[i].y = 0;
	    v_new[i].x = v_new[i].y = 0;
	}
    }

    for( k = 0; k < P+1; k++ )
    {
	b = alpha1[k] + alpha2[k] * a1;
	
	for( i = 0; i < max_stars; i++ )
	{
	    if( m[i] <= 0. || m[i] >= collapsar )
		continue;
		
	    x_new[i].x += b * txk[k][i].x;
	    x_new[i].y += b * txk[k][i].y;

	    v_new[i].x += b * tvk[k][i].x;
	    v_new[i].y += b * tvk[k][i].y;
	}
    }

    /*
     * plot the points
     */
    for( i = 0; i < max_stars; i++ )
    {
	m_i = m[i];
	if( m_i <= 0.0 )
	    continue;

	/* check for collisions */
	for( j = i+1; j < max_stars; j++ )
	{
	    if( m[j] <= 0.0 )
		continue;
	    
	    dx = x_new[i].x - x_new[j].x;
	    dy = x_new[i].y - x_new[j].y;
	    
	    if( dx*dx + dy*dy < collide_dist )
		collide( i, j, x, x_new, m, vk, ak, sp, sd );
	}
	
		
	xpt = x_new[i].x;
	ypt = x_new[i].y;
	
#	define cur	x_new
#	define prev	x
#	define v_0	v
#	include "plot_star.h"	
	
    }
}
END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file move_stars_ab7.c ----------

filename="move_stars_ab7.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file move_stars_ab7.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.
*/



/*
 * Adam-Bashford method of 7th order
 */

static int 	num_rk4;

void restart_ab7( int n, double *m, point_2d *x, point_2d *v )
{
    num_rk4 = 7;
}



void init_ab7( int n, double *m, point_2d *x, point_2d *v )
{
    init_rk4( n, m, x, v );
    
    restart_ab7( n, m, x, v );
}



void move_stars_ab7( point_2d *prev, point_2d *cur, double *m, 
		point_2d *vk[], point_2d *ak[],
		int max_stars, sys_param_type *sp, star_disp_type *sd )
{
    double		m_i, m_j;
    double		previ_x, previ_y;
    double		ai_x, ai_y;

    register int        i, j, k;	/* star index */
    int			tx, ty;

    int			xpt, ypt;
    
    point_2d	*v_0 = vk[0];
    point_2d	*v_1 = vk[1];
    point_2d	*v_2 = vk[2];
    point_2d	*v_3 = vk[3];
    point_2d	*v_4 = vk[4];
    point_2d	*v_5 = vk[5];
    point_2d	*v_6 = vk[6];
    point_2d	*v_7 = vk[7];
    point_2d	*a_0 = ak[0];
    point_2d	*a_1 = ak[1];
    point_2d	*a_2 = ak[2];
    point_2d	*a_3 = ak[3];
    point_2d	*a_4 = ak[4];
    point_2d	*a_5 = ak[5];
    point_2d	*a_6 = ak[6];


    /*
     * we must use some other, self starting method for the first 7
     * calls to the Adam-Bashford 7th order method
     */
    if( num_rk4 )
    {
	--num_rk4;
	move_stars_rk4( prev, cur, m, vk, ak, max_stars, sp, sd );
	return;
    }


    for( i = 0; i < max_stars; i++ )
    {
	m_i = m[i];
	if( m_i <= 0.0 )
	    continue;
	
	previ_x = prev[i].x;
	previ_y = prev[i].y;
	ai_x = a_0[i].x;
	ai_y = a_0[i].y;
	
	for( j = i+1; j < max_stars; j++ )
	{
	    double dx, dy, n2, f;
	    register double f_d;
	    
	    m_j = m[j];
	    if( m_j <= 0.0 )
		continue;
	    
	    dx = prev[j].x - previ_x;
	    dy = prev[j].y - previ_y;
	    n2 = dx*dx + dy*dy;
	    
	    if( n2 < collide_dist )
	    {
		a_0[i].x = ai_x;
		a_0[i].y = ai_y;
		collide( i, j, cur, prev, m, vk, ak, sp, sd );
		m_i = m[i];
		m_j = m[j];
		previ_x = prev[i].x;
		previ_y = prev[i].y;
		ai_x = a_0[i].x;
		ai_y = a_0[i].y;
	    }
	    else
	    {
		f = G / (sqrt(n2) * n2);
		
		f_d = dx * f;		
		ai_x += f_d * m_j;
		a_0[j].x -= f_d * m_i;
		
		f_d = dy * f;
		ai_y += f_d * m_j;
		a_0[j].y -= f_d * m_i;
	    }
	}
	
	/*
	 * Move the star
	 */
	
	if( m_i >= collapsar )
	{
	    /* I think it looks better if the collapsars don't move ;-) */
	    v_1[i].x = v_1[i].y = 0.0;
	    
	    a_0[i].x = a_0[i].y = 0.0;
	    ai_x = ai_y = 0.0;
	    a_1[i].x = a_1[i].y = 0.0;
	    
	    xpt = cur[i].x = prev[i].x;
	    ypt = cur[i].y = prev[i].y;
	}
	else
	{
	    a_0[i].x = ai_x;
	    v_0[i].x = v_1[i].x + (198721/60480.)*a_0[i].x - (447288/60480.)*a_1[i].x + (705549/60480.)*a_2[i].x - (688256/60480.)*a_3[i].x + (407139/60480.)*a_4[i].x - (134472/60480.)*a_5[i].x + (19087/60480.)*a_6[i].x;
	    xpt = cur[i].x = prev[i].x + (198721/60480.)*v_1[i].x - (447288/60480.)*v_2[i].x + (705549/60480.)*v_3[i].x - (688256/60480.)*v_4[i].x + (407139/60480.)*v_5[i].x - (134472/60480.)*v_6[i].x + (19087/60480.)*v_7[i].x;

	    a_0[i].y = ai_y;
	    v_0[i].y = v_1[i].y + (198721/60480.)*a_0[i].y - (447288/60480.)*a_1[i].y + (705549/60480.)*a_2[i].y - (688256/60480.)*a_3[i].y + (407139/60480.)*a_4[i].y - (134472/60480.)*a_5[i].y + (19087/60480.)*a_6[i].y;
	    ypt = cur[i].y = prev[i].y + (198721/60480.)*v_1[i].y - (447288/60480.)*v_2[i].y + (705549/60480.)*v_3[i].y - (688256/60480.)*v_4[i].y + (407139/60480.)*v_5[i].y - (134472/60480.)*v_6[i].y + (19087/60480.)*v_7[i].y;
	    

	    vk[K-1][i].x = vk[K-1][i].y = 0;
	}
	ak[K-1][i].x = ak[K-1][i].y = 0;
	
#	include "plot_star.h"	
	
    }
}
END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file move_stars_am7.c ----------

filename="move_stars_am7.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file move_stars_am7.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.
*/



/*
 * 6th order Adam-Bashford predictor with a 7th order Adam-Moulton corrector, 
 */


static point_2d	*tloc, *tv, *ta;

static int 	num_rk4;

static int	first_time = 1;

void restart_am7( int n, double *m, point_2d *x, point_2d *v )
{
    num_rk4 = 5;
}


void init_am7( int n, double *m, point_2d *x, point_2d *v )
{

    int		i;

    static int	last_n = -1;
    

    init_rk4( n, m, x, v );
    
    restart_am7( n, m, x, v );


    if( n != last_n )
    {
	last_n = n;
	
	/* allocate storage */
	if( !first_time )
	{
	    free( tloc );
	    free( tv );
	    free( ta );
	}
	
	tloc = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	tv = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );
	ta = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n );

	for( i = 0; i < n; i++ )
	{
	    tloc[i].x = tloc[i].y = 0;
	    tv[i].x = tv[i].y = 0;
	    ta[i].x = ta[i].y = 0;
	}
    }

    first_time = 0;
	
}



static void collide_am7( int i, int j,
			point_2d *cur, point_2d *prev, point_2d *tloc,
			double *m,
			point_2d *vk[], point_2d *ak[],
			point_2d *tv, point_2d *ta,
			sys_param_type *sp, star_disp_type *sd
			)
{
    /* 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;
    }
    tv[i].x = (m[i]*tv[i].x + m[j]*tv[j].x)*tmass;
    tv[i].y = (m[i]*tv[i].y + m[j]*tv[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;
	tloc[i].x = (m[i]*tloc[i].x + m[j]*tloc[j].x)*tmass;
	tloc[i].y = (m[i]*tloc[i].y + m[j]*tloc[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;
	}
	ta[i].x = (m[i]*ta[i].x + m[j]*ta[j].x)*tmass;
	ta[i].y = (m[i]*ta[i].y + m[j]*ta[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 = tloc[j].x = tloc[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;
    }	    
    tv[j].x = tv[j].y = 0;
    ta[j].x = ta[j].y = 0;

    init_fna( max_stars, m, tloc, tv );
}


void move_stars_am7( point_2d *prev, point_2d *cur, double *m, 
		point_2d *vk[], point_2d *ak[],
		int max_stars, sys_param_type *sp, star_disp_type *sd )
{
    double		m_i, m_j;
    double		previ_x, previ_y;
    double		tloci_x, tloci_y;
    double		ai_x, ai_y;

    register int        i, j, k;	/* star index */
    int			tx, ty;

    int			xpt, ypt;
    
    point_2d	*v_0 = vk[0];
    point_2d	*v_1 = vk[1];
    point_2d	*v_2 = vk[2];
    point_2d	*v_3 = vk[3];
    point_2d	*v_4 = vk[4];
    point_2d	*v_5 = vk[5];
    point_2d	*v_6 = vk[6];
    point_2d	*a_0 = ak[0];
    point_2d	*a_1 = ak[1];
    point_2d	*a_2 = ak[2];
    point_2d	*a_3 = ak[3];
    point_2d	*a_4 = ak[4];
    point_2d	*a_5 = ak[5];


    /*
     * we must use some other, self starting method for the first 7
     * calls to the Adam-Bashford predictor, Adam-Moulton corrector method
     */
    if( num_rk4 )
    {
	--num_rk4;
	move_stars_rk4( prev, cur, m, vk, ak, max_stars, sp, sd );
	return;
    }


    /*
     * predictor step
     */
    for( i = 0; i < max_stars; i++ )
    {
	m_i = m[i];
	if( m_i <= 0.0 )
	    continue;
	
	previ_x = prev[i].x;
	previ_y = prev[i].y;
	ai_x = a_0[i].x;
	ai_y = a_0[i].y;
	
	for( j = i+1; j < max_stars; j++ )
	{
	    double dx, dy, n2, f;
	    register double f_d;
	    
	    m_j = m[j];
	    if( m_j <= 0.0 )
		continue;
	    
	    dx = prev[j].x - previ_x;
	    dy = prev[j].y - previ_y;
	    n2 = dx*dx + dy*dy;
	    
	    if( n2 < collide_dist )
	    {
		a_0[i].x = ai_x;
		a_0[i].y = ai_y;
		collide( i, j, cur, prev, m, vk, ak, sp, sd );
		m_i = m[i];
		m_j = m[j];
		previ_x = prev[i].x;
		previ_y = prev[i].y;
		ai_x = a_0[i].x;
		ai_y = a_0[i].y;
	    }
	    else
	    {
		f = G / (sqrt(n2) * n2);
		
		f_d = dx * f;		
		ai_x += f_d * m_j;
		a_0[j].x -= f_d * m_i;
		
		f_d = dy * f;
		ai_y += f_d * m_j;
		a_0[j].y -= f_d * m_i;
	    }
	}
	
	/*
	 * Move the star
	 */
	
	if( m_i >= collapsar )
	{
	    /* I think it looks better if the collapsars don't move ;-) */
	    v_1[i].x = v_1[i].y = 0.0;
	    
	    a_0[i].x = a_0[i].y = 0.0;
	    ai_x = ai_y = 0.0;
	    a_1[i].x = a_1[i].y = 0.0;

	    tloc[i].x = prev[i].x;
	    tloc[i].y = prev[i].y;
	}
	else
	{
	    a_0[i].x = ai_x;
	    a_0[i].y = ai_y;

	    tv[i].x = v_1[i].x + (4227/1440.)*a_0[i].x - (7673/1440.)*a_1[i].x + (9482/1440.)*a_2[i].x - (6798/1440.)*a_3[i].x + (2627/1440.)*a_4[i].x - (425/1440.)*a_5[i].x;
	    tv[i].y = v_1[i].y + (4227/1440.)*a_0[i].y - (7673/1440.)*a_1[i].y + (9482/1440.)*a_2[i].y - (6798/1440.)*a_3[i].y + (2627/1440.)*a_4[i].y - (425/1440.)*a_5[i].y;
	    
	    tloc[i].x = prev[i].x + (4227/1440.)*v_1[i].x - (7673/1440.)*v_2[i].x + (9482/1440.)*v_3[i].x - (6798/1440.)*v_4[i].x + (2627/1440.)*v_5[i].x - (425/1440.)*v_6[i].x;
	    tloc[i].y = prev[i].y + (4227/1440.)*v_1[i].y - (7673/1440.)*v_2[i].y + (9482/1440.)*v_3[i].y - (6798/1440.)*v_4[i].y + (2627/1440.)*v_5[i].y - (425/1440.)*v_6[i].y;
	    
	}
    }


    /*
     * corrector step
     */
    for( i = 0; i < max_stars; i++ )
    {
	m_i = m[i];
	if( m_i <= 0.0 )
	    continue;
	
	tloci_x = tloc[i].x;
	tloci_y = tloc[i].y;
	ai_x = ta[i].x;
	ai_y = ta[i].y;
	
	for( j = i+1; j < max_stars; j++ )
	{
	    double dx, dy, n2, f;
	    register double f_d;
	    
	    m_j = m[j];
	    if( m_j <= 0.0 )
		continue;
	    
	    dx = tloc[j].x - tloci_x;
	    dy = tloc[j].y - tloci_y;
	    n2 = dx*dx + dy*dy;
	    
	    if( n2 < collide_dist )
	    {
		ta[i].x = ai_x;
		ta[i].y = ai_y;
		collide_am7( i, j, cur, prev, tloc, m, vk, ak, tv, ta, sp, sd );
		m_i = m[i];
		m_j = m[j];
		tloci_x = tloc[i].x;
		tloci_y = tloc[i].y;
		ai_x = ta[i].x;
		ai_y = ta[i].y;
	    }
	    else
	    {
		f = G / (sqrt(n2) * n2);
		
		f_d = dx * f;		
		ai_x += f_d * m_j;
		ta[j].x -= f_d * m_i;
		
		f_d = dy * f;
		ai_y += f_d * m_j;
		ta[j].y -= f_d * m_i;
	    }
	}
	
	/*
	 * Move the star
	 */
	
	if( m_i >= collapsar )
	{
	    /* I think it looks better if the collapsars don't move ;-) */
	    v_1[i].x = v_1[i].y = 0.0;
	    
	    ta[i].x = ta[i].y = 0.0;
	    ai_x = ai_y = 0.0;
	    a_1[i].x = a_1[i].y = 0.0;
	    
	    xpt = cur[i].x = tloc[i].x;
	    ypt = cur[i].y = tloc[i].y;
	}
	else
	{
	    ta[i].x = ai_x;
	    ta[i].y = ai_y;

	    v_0[i].x = v_1[i].x + (19087/60480.)*ta[i].x + (65112/60480.)*a_0[i].x - (46461/60480.)*a_1[i].x + (37504/60480.)*a_2[i].x - (20211/60480.)*a_3[i].x + (6312/60480.)*a_4[i].x - (863/60480.)*a_5[i].x;
	    v_0[i].y = v_1[i].y + (19087/60480.)*ta[i].y + (65112/60480.)*a_0[i].y - (46461/60480.)*a_1[i].y + (37504/60480.)*a_2[i].y - (20211/60480.)*a_3[i].y + (6312/60480.)*a_4[i].y - (863/60480.)*a_5[i].y;

	    xpt = cur[i].x = prev[i].x + (19087/60480.)*tv[i].x + (65112/60480.)*v_1[i].x - (46461/60480.)*v_2[i].x + (37504/60480.)*v_3[i].x - (20211/60480.)*v_4[i].x + (6312/60480.)*v_5[i].x - (863/60480.)*v_6[i].x;
	    ypt = cur[i].y = prev[i].y + (19087/60480.)*tv[i].y + (65112/60480.)*v_1[i].y - (46461/60480.)*v_2[i].y + (37504/60480.)*v_3[i].y - (20211/60480.)*v_4[i].y + (6312/60480.)*v_5[i].y - (863/60480.)*v_6[i].y;

	    vk[K-1][i].x = vk[K-1][i].y = 0;
	}
	ta[i].x = ta[i].y = 0;
	ak[K-1][i].x = ak[K-1][i].y = 0;
	
#	include "plot_star.h"	
	
    }
}
END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file dump_sys.c ----------

filename="dump_sys.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file dump_sys.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 dumps a star system into a format suitable for including
 * in another C program.  This routine was used to create the init_sys_*()
 * routines.
 */

void dump_sys( point_2d *cur, point_2d *prev, double *m,
	      point_2d *vk[], point_2d *ak[],
	      int max_stars, sys_param_type *sp )
{
    int		i, k;
    
    /*
     * dump the variables so that we can recreate this system later.
     */

    fprintf( stderr, "    sp->star_distrib     = %d;\n", sp->star_distrib );
    fprintf( stderr, "    sp->size             = %.16e;\n", sp->size );
    fprintf( stderr, "    sp->num_collapsar    = %d;\n", sp->num_collapsar );
    fprintf( stderr, "    sp->mono_stars       = %d;\n", sp->mono_stars );
    fprintf( stderr, "\n" );
    fprintf( stderr, "    sp->star_circle      = %d;\n", sp->star_circle );
    fprintf( stderr, "    sp->cir_dist         = %d;\n", sp->cir_dist );
    fprintf( stderr, "\n" );
    fprintf( stderr, "    sp->star_line        = %d;\n", sp->star_line );
    fprintf( stderr, "    sp->rnd_spacing      = %d;\n", sp->rnd_spacing );
    fprintf( stderr, "\n" );
    fprintf( stderr, "    sp->do_bounce        = %d;\n", sp->do_bounce );
    fprintf( stderr, "    sp->no_speed         = %d;\n", sp->no_speed );
    fprintf( stderr, "    sp->few_stars        = %d;\n", sp->few_stars );
    fprintf( stderr, "    sp->drift            = %d;\n", sp->drift );
    fprintf( stderr, "\n" );
    fprintf( stderr, "    sp->energy_fact      = %d;\n", sp->energy_fact );
    fprintf( stderr, "    sp->calc_energy_fact = %d;\n", sp->calc_energy_fact );
    fprintf( stderr, "\n" );
    fprintf( stderr, "    sp->min_stars        = %d;\n", sp->min_stars );
    fprintf( stderr, "    sp->num_add          = %d;\n", sp->num_add );
    fprintf( stderr, "    sp->live_stars       = %d;\n", sp->live_stars );
    
    for( i = 0; i < max_stars; i++ )
    {
	if( m[i] <= 0 )
	    continue;
	
	fprintf( stderr, "\n/* star %d */\n", i );
	fprintf( stderr, "\n    m[%d]        = %.16e;\n", i, m[i]/fm );
	fprintf( stderr, "    cur[%d].x    = %.16e;\n", i, cur[i].x );
	fprintf( stderr, "    cur[%d].y    = %.16e;\n", i, cur[i].y );
	fprintf( stderr, "    prev[%d].x   = %.16e;\n", i, prev[i].x );
	fprintf( stderr, "    prev[%d].y   = %.16e;\n", i, prev[i].y );

	for( k = 0; k < K; k++ )
	{
	    if( vk[k][i].x != 0 || vk[k][i].y != 0 )
	    {
		fprintf( stderr, "    vk[%d][%d].x  = %.16e;\n",
			k, i, vk[k][i].x*fv );
		fprintf( stderr, "    vk[%d][%d].y  = %.16e;\n",
			k, i, vk[k][i].y*fv );
	    }
	    
	    if( ak[k][i].x != 0 || ak[k][i].y != 0 )
	    {
		fprintf( stderr, "    ak[%d][%d].x  = %.16e;\n",
			k, i, ak[k][i].x*fv*fv );
		fprintf( stderr, "    ak[%d][%d].y  = %.16e;\n",
			k, i, ak[k][i].y*fv*fv );
	    }
	    
	}
    }
}

END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file init_sys_4.c ----------

filename="init_sys_4.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file init_sys_4.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.
*/



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

    sp->star_distrib  = 0;
    sp->size          = 2.1400000000000000e+02;
    sp->num_collapsar = 0;
    sp->mono_stars    = 4;

    sp->star_circle   = 0;
    sp->cir_dist      = 0;

    sp->star_line     = 0;

    sp->do_bounce     = 0;
    sp->no_speed      = 0;
    sp->few_stars     = 0;
    sp->drift         = 0;

    sp->min_stars     = 1;
    sp->num_add       = 0;
    sp->live_stars    = 4;

/* star 0 */

    m[0]        = 2.7865168539325840e+00;
    cur[0].x    = 1.5835647927008775e+01;
    cur[0].y    = 1.8923118463283608e+01;
    prev[0].x   = 1.5835647927008775e+01;
    prev[0].y   = 1.8923118463283608e+01;
    vk[0][0].x  = 9.0941882306722933e-03;
    vk[0][0].y  = 1.6257009059518187e-03;
    ak[0][0].x  = 6.9931581889377990e-06;
    ak[0][0].y  = -8.0395052303351502e-06;

/* star 1 */

    m[1]        = 5.4831460674157304e+00;
    cur[1].x    = 6.0521453125534116e+01;
    cur[1].y    = 9.4644832984345086e+00;
    prev[1].x   = 6.0521453125534116e+01;
    prev[1].y   = 9.4644832984345086e+00;
    vk[0][1].x  = 2.5826514208378845e-03;
    vk[0][1].y  = -1.4348317268012971e-02;
    ak[0][1].x  = -8.3617935918065052e-06;
    ak[0][1].y  = -4.3824681501481231e-07;

/* star 2 */

    m[2]        = 2.5168539325842696e+00;
    cur[2].x    = 8.3340088264855723e+00;
    cur[2].y    = -2.3248024535545774e+01;
    prev[2].x   = 8.3340088264855723e+00;
    prev[2].y   = -2.3248024535545774e+01;
    vk[0][2].x  = -1.4385862343579154e-02;
    vk[0][2].y  = -1.3381252605483343e-03;
    ak[0][2].x  = 3.0889420909050185e-06;
    ak[0][2].y  = 9.3058808330760419e-06;

/* star 3 */

    m[3]        = 5.2134831460674160e+00;
    cur[3].x    = -7.6139068509249114e+01;
    cur[3].y    = -8.8449218374312917e+00;
    prev[3].x   = -7.6139068509249114e+01;
    prev[3].y   = -8.8449218374312917e+00;
    vk[0][3].x  = -6.3202458968509645e-04;
    vk[0][3].y  = 1.4867554009648933e-02;
    ak[0][3].x  = 3.5653642534100792e-06;
    ak[0][3].y  = 2.6539749189938063e-07;


    /* scale things to the current accuracy */

    for( i = 0; i < max_stars; i++ )
    {
	if( m[i] <= 0 )
	    continue;
	
	m[i] *= fm;

	for( k = 0; k < K; k++ )
	{
	    vk[k][i].x *= fv_inv;
	    vk[k][i].y *= fv_inv;

	    ak[k][i].x *= fv_inv*fv_inv;
	    ak[k][i].y *= fv_inv*fv_inv;
	}
    }
}
END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file init_sys_8.c ----------

filename="init_sys_8.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file init_sys_8.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.
*/



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

    sp->star_distrib  = 0;
    sp->size          = 3.0300000000000000e+02;
    sp->num_collapsar = 0;
    sp->mono_stars    = 8;

    sp->star_circle   = 0;
    sp->cir_dist      = 0;

    sp->star_line     = 0;

    sp->do_bounce     = 0;
    sp->no_speed      = 0;
    sp->few_stars     = 0;
    sp->drift         = 0;

    sp->min_stars     = 1;
    sp->num_add       = 0;
    sp->live_stars    = 8;

/* star 0 */

    m[0]        = 1.1468208092485550e+00;
    cur[0].x    = -2.7038883602489072e+00;
    cur[0].y    = 6.4089177018335889e+01;
    prev[0].x   = -2.7038883602489072e+00;
    prev[0].y   = 6.4089177018335889e+01;
    vk[0][0].x  = 1.3226751051071546e-02;
    vk[0][0].y  = -2.2643948440514610e-03;
    ak[0][0].x  = -2.1461637019132762e-06;
    ak[0][0].y  = -3.8329229061867005e-06;

/* star 1 */

    m[1]        = 2.2566473988439308e+00;
    cur[1].x    = 6.0566200308691172e+01;
    cur[1].y    = 5.0696810406423374e+01;
    prev[1].x   = 6.0566200308691172e+01;
    prev[1].y   = 5.0696810406423374e+01;
    vk[0][1].x  = 1.1326416168617813e-02;
    vk[0][1].y  = -1.2983888267850539e-02;
    ak[0][1].x  = -3.7978097624060328e-06;
    ak[0][1].y  = -4.9601977809216730e-06;

/* star 2 */

    m[2]        = 1.0358381502890175e+00;
    cur[2].x    = -1.3325368021270076e+01;
    cur[2].y    = 4.3795680059746633e+00;
    prev[2].x   = -1.3325368021270076e+01;
    prev[2].y   = 4.3795680059746633e+00;
    vk[0][2].x  = -1.4657477777562395e-03;
    vk[0][2].y  = 8.3098025452551811e-03;
    ak[0][2].x  = 1.4837373960419059e-06;
    ak[0][2].y  = 3.2846997147130252e-06;

/* star 3 */

    m[3]        = 2.7560693641618497e+00;
    cur[3].x    = 5.0296086642545227e+01;
    cur[3].y    = 3.3756236820588086e+00;
    prev[3].x   = 5.0296086642545227e+01;
    prev[3].y   = 3.3756236820588086e+00;
    vk[0][3].x  = -3.7597216177549879e-03;
    vk[0][3].y  = -8.3720019727285024e-03;
    ak[0][3].x  = -1.9157531920650712e-06;
    ak[0][3].y  = 3.8658856402276140e-06;

/* star 4 */

    m[4]        = 2.1271676300578037e+00;
    cur[4].x    = 6.8270243736874789e+01;
    cur[4].y    = -1.1595121023145049e+02;
    prev[4].x   = 6.8270243736874789e+01;
    prev[4].y   = -1.1595121023145049e+02;
    vk[0][4].x  = -1.4273624489641538e-02;
    vk[0][4].y  = -8.7326800300004132e-03;
    ak[0][4].x  = -9.7523640637911302e-07;
    ak[0][4].y  = 1.7221688018067123e-06;

/* star 5 */

    m[5]        = 2.0531791907514454e+00;
    cur[5].x    = -6.1922907126650045e+01;
    cur[5].y    = 8.2905813493962910e+01;
    prev[5].x   = -6.1922907126650045e+01;
    prev[5].y   = 8.2905813493962910e+01;
    vk[0][5].x  = 1.4057051300353198e-02;
    vk[0][5].y  = 9.7559665347273592e-03;
    ak[0][5].x  = 4.3248908536138410e-06;
    ak[0][5].y  = -6.5951618951828798e-06;

/* star 6 */

    m[6]        = 1.9236994219653181e+00;
    cur[6].x    = -4.8033093344416187e+01;
    cur[6].y    = 4.7578502331255180e+01;
    prev[6].x   = -4.8033093344416187e+01;
    prev[6].y   = 4.7578502331255180e+01;
    vk[0][6].x  = 3.8568338976656485e-03;
    vk[0][6].y  = 9.4050873297924203e-03;
    ak[0][6].x  = 2.2103324433326556e-06;
    ak[0][6].y  = 3.6631219712243281e-06;

/* star 7 */

    m[7]        = 2.7005780346820809e+00;
    cur[7].x    = -6.8161132128329925e+01;
    cur[7].y    = -8.0295095516064634e+01;
    prev[7].x   = -6.8161132128329925e+01;
    prev[7].y   = -8.0295095516064634e+01;
    vk[0][7].x  = -1.2873818798120700e-02;
    vk[0][7].y  = 9.9296098370418636e-03;
    ak[0][7].x  = 1.3764904099135169e-06;
    ak[0][7].y  = 1.6155758780025582e-06;


    /* scale things to the current accuracy */

    for( i = 0; i < max_stars; i++ )
    {
	if( m[i] <= 0 )
	    continue;
	
	m[i] *= fm;

	for( k = 0; k < K; k++ )
	{
	    vk[k][i].x *= fv_inv;
	    vk[k][i].y *= fv_inv;

	    ak[k][i].x *= fv_inv*fv_inv;
	    ak[k][i].y *= fv_inv*fv_inv;
	}
    }
}
END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file init_sys_8b.c ----------

filename="init_sys_8b.c"

if [ -f $filename ]
then
  echo File \"$filename\" already exists\!  Skipping...
  filename=/dev/null		# throw it away
else
  echo extracting file init_sys_8b.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.
*/



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

    sp->star_distrib  = 0;
    sp->size          = 3.0300000000000000e+02;
    sp->num_collapsar = 0;
    sp->mono_stars    = 8;

    sp->star_circle   = 0;
    sp->cir_dist      = 0;

    sp->star_line     = 0;

    sp->do_bounce     = 0;
    sp->no_speed      = 0;
    sp->few_stars     = 0;
    sp->drift         = 0;

    sp->min_stars     = 1;
    sp->num_add       = 0;
    sp->live_stars    = 8;

/* star 0 */

    m[0]        = 3.2820512820512824e+00;
    cur[0].x    = -5.9426018201326272e+01;
    cur[0].y    = 1.7647796986427664e+01;
    prev[0].x   = -5.9426018201326272e+01;
    prev[0].y   = 1.7647796986427664e+01;
    vk[0][0].x  = 4.8740542740093326e-03;
    vk[0][0].y  = 1.0900511132255125e-02;
    ak[0][0].x  = 2.5173820214687024e-06;
    ak[0][0].y  = -2.4452551992655210e-06;

/* star 1 */

    m[1]        = 2.0968660968660968e+00;
    cur[1].x    = 9.3599198133270978e-01;
    cur[1].y    = -4.0668403188234521e+01;
    prev[1].x   = 9.3599198133270978e-01;
    prev[1].y   = -4.0668403188234521e+01;
    vk[0][1].x  = -7.7175472328441120e-03;
    vk[0][1].y  = -3.7383560620019456e-03;
    ak[0][1].x  = -3.2319296217630353e-06;
    ak[0][1].y  = 1.7991971056596420e-06;

/* star 2 */

    m[2]        = 2.0512820512820515e+00;
    cur[2].x    = -5.9384119374432700e+01;
    cur[2].y    = 9.8236339712784101e+01;
    prev[2].x   = -5.9384119374432700e+01;
    prev[2].y   = 9.8236339712784101e+01;
    vk[0][2].x  = 1.5656891177842101e-02;
    vk[0][2].y  = 8.0608082916574265e-03;
    ak[0][2].x  = 8.2457296929936326e-07;
    ak[0][2].y  = -3.3782755206993198e-06;

/* star 3 */

    m[3]        = 9.3447293447293456e-01;
    cur[3].x    = -7.0405575139421671e+01;
    cur[3].y    = -2.7001090804229296e+01;
    prev[3].x   = -7.0405575139421671e+01;
    prev[3].y   = -2.7001090804229296e+01;
    vk[0][3].x  = -7.6525422837526603e-03;
    vk[0][3].y  = 1.6937754974762298e-02;
    ak[0][3].x  = 6.5073079181275713e-06;
    ak[0][3].y  = 3.6029649981805933e-06;

/* star 4 */

    m[4]        = 2.2336182336182335e+00;
    cur[4].x    = 9.2302766801034323e+01;
    cur[4].y    = -6.0862441826543247e+01;
    prev[4].x   = 9.2302766801034323e+01;
    prev[4].y   = -6.0862441826543247e+01;
    vk[0][4].x  = -1.0083930782485945e-02;
    vk[0][4].y  = -1.3905687652693736e-02;
    ak[0][4].x  = -2.1482127913348859e-06;
    ak[0][4].y  = 2.1385211861947785e-06;

/* star 5 */

    m[5]        = 7.5213675213675213e-01;
    cur[5].x    = -7.2473741939590468e+00;
    cur[5].y    = 1.6015571178554193e+01;
    prev[5].x   = -7.2473741939590468e+00;
    prev[5].y   = 1.6015571178554193e+01;
    vk[0][5].x  = 8.7707293778696745e-03;
    vk[0][5].y  = -9.8131165679121637e-04;
    ak[0][5].x  = -4.0382422553598351e-06;
    ak[0][5].y  = -3.0253523635815457e-06;

/* star 6 */

    m[6]        = 2.9173789173789175e+00;
    cur[6].x    = 8.7255538788646888e+01;
    cur[6].y    = 2.6925933985080562e+01;
    prev[6].x   = 8.7255538788646888e+01;
    prev[6].y   = 2.6925933985080562e+01;
    vk[0][6].x  = 5.5994912063467676e-03;
    vk[0][6].y  = -1.4944457539263970e-02;
    ak[0][6].x  = -2.0118695310450290e-06;
    ak[0][6].y  = -1.6850985387665771e-06;

/* star 7 */

    m[7]        = 1.7321937321937322e+00;
    cur[7].x    = -4.3062922104262000e+01;
    cur[7].y    = -5.9796442419559682e+01;
    prev[7].x   = -4.3062922104262000e+01;
    prev[7].y   = -5.9796442419559682e+01;
    vk[0][7].x  = -1.4541568021132618e-02;
    vk[0][7].y  = 8.7153380803802905e-03;
    ak[0][7].x  = 2.5674930287103935e-06;
    ak[0][7].y  = 5.9061580476882484e-06;


    /* scale things to the current accuracy */

    for( i = 0; i < max_stars; i++ )
    {
	if( m[i] <= 0 )
	    continue;
	
	m[i] *= fm;

	for( k = 0; k < K; k++ )
	{
	    vk[k][i].x *= fv_inv;
	    vk[k][i].y *= fv_inv;

	    ak[k][i].x *= fv_inv*fv_inv;
	    ak[k][i].y *= fv_inv*fv_inv;
	}
    }
}
END-OF-FILE

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

  if [ $size != 4465 ]
  then
    echo $filename changed - should be 4465 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.


