Article 12755 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    Part03/06
Date: Sun, 27 Aug 1995 23:13:17 GMT
Organization: The Backbone Cabal
Lines: 2000
Sender: wayne@backbone (Wayne Schlitt)
Message-ID: <WAYNE.95Aug27171317@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/part03


--- cut here ---

# Continuation of Shell Archive, created by backbone!wayne

# This is part 3

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

filename="xstar.c"

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

cat << 'END-OF-FILE' > $filename
/*
** XStar -- a animated n-body solver for X windows
** Copyright (C) 1995  Wayne Schlitt (wayne@cse.unl.edu) and others
** Version 2.0.0  8/27/95
**
** Insofar as I have any claim to this program (and I haven't much since it's
** a hack of somebody elses) the GNU General Public License would apply.
**
**  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 program is distributed in the hope that it will be useful,
**  but WITHOUT ANY WARRANTY; without even the implied warranty of
**  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
**  GNU General Public License for more details.
**
**  You should have received a copy of the GNU General Public License
**  along with this program; if not, write to the Free Software
**  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
**
**
** XStar was based on XGrav 1.0.1 (ca 1/20/95)
** written by David Flater (dave@case50.ncsl.nist.gov)
** The README file for XGrav is in README.xgrav.
**
** XGrav was based on XSwarm 2.3 (ca 1/4/91)
** written by Jeff Butterworth (butterwo@cs.unc.edu)
**
** XSwarm was based on psychoII (date and author unknown)
** and on meltdown (date and author also unknown)
**
** psychoII was possibly based on ico (date and author unknown)
**
** ico was probably a blatant rip-off of some other program :->
*/

#include "xstar.h"
#include "xstar_ext.h"

/* These are needed for the nap function. */
#include <sys/time.h>
#include <signal.h>

/* This makes it work with virtual root window managers. */
#include "vroot.h"

#include <time.h>
#include <assert.h>
#include <string.h>
#include "patchlevel.h"

/* Who needs the this header file?? */
#if defined(NEEDS_GETOPT_H)
#include <getopt.h>
#endif



int
main(argc, argv)
int argc;
char **argv;
{
    char        *geometry = NULL;


    Initialize();

    Parse_Arguments(argc, argv, &geometry);

    /* The program acts as a screen saver if timeout is non-zero. */
    if( !timeout )
    {
        Create_Window(geometry);
        Animate();
    }
    else
    {
        Change_Screen_Saver(False);
        while(TRUE)
        {
            Wait_For_Idleness();
            Create_Big_Window();
            Animate();
        }
    }

    return 0;
}


void
Quit( int signal )
{
    if( timeout )
    {
        fprintf(stderr, "Restoring screen saver state.\n");
        Change_Screen_Saver(True);
    }

    fprintf(stderr, "Terminating because of signal %d. ", signal);
    switch( signal )
    {
        case SIGHUP:    fprintf(stderr, "(Hangup)\n"); break;
        case SIGINT:    fprintf(stderr, "(Interrupt)\n"); break;
        case SIGQUIT:   fprintf(stderr, "(Quit)\n"); break;
        case SIGILL:    fprintf(stderr, "(Illegal Instruction)\n"); break;
        case SIGFPE:    fprintf(stderr, "(Floating Point Exception)\n"); break;
        case SIGKILL:   fprintf(stderr, "(Kill)\n"); break;
        case SIGBUS:    fprintf(stderr, "(Bus Error)\n"); break;
        case SIGSEGV:   fprintf(stderr, "(Segmentation Violation)\n"); break;
        case SIGTERM:   fprintf(stderr, "(Termination from Kill)\n"); break;
        case SIGSTOP:   fprintf(stderr, "(Stop not from tty)\n"); break;
        case SIGTSTP:   fprintf(stderr, "(Stop from tty)\n"); break;
#ifdef SIGXCPU
        case SIGXCPU:   fprintf(stderr, "(Exceeded CPU Time Limit)\n"); break;
#endif
        default:        fprintf(stderr, "(Unexpected Signal)\n");
    }

    exit(-signal);
}


void
Initialize()
{
    double	rand_test;
    
    /* Get the random number generator ready. */
    srand48(time(NULL));
    sleep(1);				/* let other xtstars seed w/same time*/
    rand_test = drand48();
    if( rand_test > 1 || rand_test < 0 )
    {
	fprintf( stderr, "Error:  drand48 returned %g, which is <0 or >1\n" );
	fprintf( stderr, "Most likely, drand48() was not given a prototype when compiled." );
	fprintf( stderr, "Try defining RAND48_PROTO and recompling." );
	exit( 1 );
    }
    

    signal(SIGHUP, Quit);
    signal(SIGINT, Quit);
    signal(SIGQUIT, Quit);
    signal(SIGILL, Quit);
    signal(SIGFPE, Quit);
    signal(SIGKILL, Quit);
    signal(SIGBUS, Quit);
    signal(SIGSEGV, Quit);
    signal(SIGTERM, Quit);
    signal(SIGSTOP, Quit);
    if( timeout ) signal(SIGTSTP, Quit);
#ifdef SIGXCPU
    signal(SIGXCPU, Quit);
#endif
}


/* Parse_Arguments()                                                    */
/* Analyze command line arguments and set variables accordingly.        */

void
Parse_Arguments(argc, argv, geometry)
int     argc;
char    **argv;
char    **geometry;
{
    /* These variables are used to analyze the command line parameters. */
    int         option;
    extern char *optarg;
    double	argf;
    double	argi;
    double	buffer_scale = 1;

    /* Check the command line. */
    while( (option = getopt(argc,argv,"hrt:d:g:D:b:c:C:vRMa:m:B:l:T:") )
            != EOF)
    {
        switch( option )
        {
            case 'r':
                root = TRUE;
                break;

            case 't':
                timeout = atoi(optarg);
                break;

            case 'g':
                *geometry = optarg;
                break;

            case 'd':
                display.dname = optarg;
                break;

            case 'b':
                stars = atoi(optarg);
		if( max_stars < stars + MAX_COLLAPSAR )
		    max_stars = stars + MAX_COLLAPSAR;
                if( stars < 2 )
		{
		    Usage(*argv);
		    HandleError("You need at least two stars.", FATAL);
		}
		if( max_stars > MAX_NUM_STARS )
		{
		    Usage(*argv);
		    HandleError("Too many stars specified.", FATAL );
		}
                break;

            case 'D':
                delay = atoi(optarg) * 1000; /* convert to microseconds */
                break;

            case 'c':
                star_color = optarg;
                break;

            case 'C':
                bg_color = optarg;
                break;

	    case 'v':
		verbose++;
		break;

	    case 'R':
		rotate_colors++;
		break;

	    case 'M':
		multi_colors++;
		break;

	    case 'a':
		argf = atof(optarg);
		if( argf <= 0 )
		{
		    Usage(*argv);
		    HandleError("-a value must be greater than zero.", FATAL );
		}
		fv *= argf;
		fv_inv = 1/fv;
		fm /= argf*argf;
		collapsar = fm*COLLAPSAR;
		break;

	    case 'm':
		if( strcmp( optarg, "taylor3" ) == 0 )
		{
		    move_fna = move_stars_taylor3;
		    restart_fna = restart_taylor3;
		    init_fna = init_taylor3;
		    move_cost = 1;
		}
		else if( strcmp( optarg, "rk4" ) == 0 )
		{
		    move_fna = move_stars_rk4;
		    restart_fna = restart_rk4;
		    init_fna = init_rk4;
		    move_cost = 4;
		}
		else if( strcmp( optarg, "gpemce8" ) == 0 )
		{
		    move_fna = move_stars_gpemce8;
		    restart_fna = restart_gpemce8;
		    init_fna = init_gpemce8;
		    move_cost = 16;
		}
		else if( strcmp( optarg, "ab7" ) == 0 )
		{
		    move_fna = move_stars_ab7;
		    restart_fna = restart_ab7;
		    init_fna = init_ab7;
		    move_cost = 1.2;
		}
		else if( strcmp( optarg, "am7" ) == 0 )
		{
		    move_fna = move_stars_am7;
		    restart_fna = restart_am7;
		    init_fna = init_am7;
		    move_cost = 2.4;
		}
		else
		{
		    Usage(*argv);
		    HandleError("-m value must be greater than zero.", FATAL );
		}
		break;

	    case 'B':
		argf = atof(optarg);
		if( argf <= 0 )
		{
		    Usage(*argv);
		    HandleError("-B value must be greater than zero.", FATAL );
		}
		buffer_scale = argf;
		break;

	    case 'l':
		argf = atof(optarg);
		if( argf <= 0 )
		{
		    Usage(*argv);
		    HandleError("-l value must be greater than zero.", FATAL );
		}
		collide_dist = argf * argf;
		break;

	    case 'T':
		argi = atoi(optarg);
		if( argi < 0 )
		{
		    Usage(*argv);
		    HandleError("-T value must be greater than or equal to zero.", FATAL );
		}
		if( argi >= MAX_DISP_PT )
		{
		    Usage(*argv);
		    HandleError("-T value must is too large.", FATAL );
		}
		num_disp_pt = argi;
		break;

		
            case 'h':
                Usage(*argv);
		exit( 0 );
		break;

	    case '?':
                Usage(*argv);
                HandleError("The command line parameters were incorrect.",
                    FATAL);
                break;
        }
    }


    /*
     * calculate the X buffering value
     *
     * this is very approximate...  It really depends on the number
     * of FLOPS that the CPU cad do and such...  The assumption is that
     * the more stars that the user selects to have, the faster the
     * computer must be, so the more buffering we can do...
     */
    if( num_disp_pt )
	max_disp_skip = (stars*stars * buffer_scale) / (2.2*move_cost*fv);
    else
	max_disp_skip = (stars*stars * buffer_scale) / (1.9*move_cost*fv);

    if( max_disp_skip > 32*32 )
	max_disp_skip = 32*32;

    if( max_disp_skip < 2*2 )
	max_disp_skip = 2*2;



    /* make sure that the # of disp pts is large enough so that
     * if they add a *bunch* of stars, we still don't have to worry
     */
    if( num_disp_pt < (MAX_NUM_STARS + 1) * 2 )
	num_disp_pt = 0;
    

    /* The screen saver is incompatible with the */
    /* root window option. */
    if( timeout > 0 )
    {
        if( root )
        {
            printf("The root window option and the screen saver option are\n");
            printf("incompatible.  The screen saver takes precedence, so\n");
            printf("that's what you're getting.\n");
        }
        root = FALSE;
    }

    /* Open the display. */
    if( !(display.dpy = XOpenDisplay(display.dname)) )
    {
        HandleError("Cannot open display.\n", FATAL);
        exit(-1);
    }

    /* Record the screen number and root window. */
    display.screen = DefaultScreen(display.dpy);
    display.root = RootWindow(display.dpy, display.screen);

    /* Set the colors. */
    display.cmap = XDefaultColormap(display.dpy, display.screen);
    if( !display.cmap ) HandleError("There was no default colormap!", FATAL);

    if( star_color == NULL )
        display.star = WhitePixel(display.dpy, display.screen);
    else
        display.star = GetColor(&display, star_color, NULL);

    if( bg_color == NULL )
        display.bg = BlackPixel(display.dpy, display.screen);
    else
        display.bg = GetColor(&display, bg_color, &(display.bg_xcolor));

    /* Set up window size. */
    if( !root && !timeout )
    {
        winW = WINWIDTH;
        winH = WINHEIGHT;
        winX = (DisplayWidth(display.dpy, display.screen)
                - winW) / 2;
        winY = (DisplayHeight(display.dpy, display.screen)
                - winH) / 2;

        if( *geometry != NULL )
	{
	    int result;

            result = XParseGeometry(*geometry, &winX, &winY, &winW, &winH);
	    if( (result & XValue) && (result & XNegative) )
		winX += DisplayWidth(display.dpy, display.screen) - winW;
	    if( (result & YValue) && (result & YNegative) )
		winY += DisplayHeight(display.dpy, display.screen) - winH;

	    if( winW < WINMIN )
		winW = WINMIN;
	    if( winH < WINMIN )
		winH = WINMIN;
	}

	root_xoff = 0;
	root_yoff = 0;
    }
    else
    {
        winX = 0;
        winY = 0;
        if( *geometry != NULL )
            XParseGeometry(*geometry, &winX, &winY, &winW, &winH);
	root_xoff = winX;
	root_yoff = winY;
	
        winW = DisplayWidth(display.dpy, display.screen);
        winH = DisplayHeight(display.dpy, display.screen);
        winX = 0;
        winY = 0;
    }
    far_dist = (winW + winH)*1.5 + 750;
    half_winW = winW/2;
    half_winH = winH/2;
    center_x = half_winW + root_xoff;
    center_y = half_winH + root_yoff;
    min_x = 0 - center_x - 20;
    max_x = winW - center_x + 20;
    min_y = 0 - center_y - 20;
    max_y = winH - center_y + 20;

    /* If screen saving is on, then watch for events everywhere. */
    /* Traverse the window tree. */
    if( timeout ) Traverse_Tree(display.dpy, display.root);
}


/* Change_Screen_Saver()                                                */
/* Turn the server's screen saver on or off.                            */
/* This routine should be called with on=False before it is called      */
/* with on=True.                                                        */

void
Change_Screen_Saver(on)
int     on;     /* True or False */
{
    static int  timeout, interval, blanking, exposures;
    static int  set_yet = FALSE;

    if( on )
    {
        if( set_yet )
        {
            /* Restore the old settings. */
            XSetScreenSaver(display.dpy, timeout, interval, blanking,exposures);
            XActivateScreenSaver(display.dpy);
            XResetScreenSaver(display.dpy);
	    XFlush( display.dpy );
        }
    }
    else
    {
        /* Save the old settings and turn off the server's screen saver. */
        XGetScreenSaver(display.dpy, &timeout, &interval, &blanking,&exposures);
        XSetScreenSaver(display.dpy, 0, 0, DefaultBlanking, DefaultExposures);
        XResetScreenSaver(display.dpy);
        set_yet = TRUE;
    }
}


/* Traverse_Tree()                                                      */
/* Select some events from every single window that is a descendent      */
/* of "current".                                                        */

void
Traverse_Tree(display, current)
Display *display;
Window  current;
{
    Window              my_root, parent, *children;
    unsigned int        num_children;
    int         i;


    /* Watch for signs of life from the user in this window. */
    XSelectInput(display, current, ALIVE_MASK);

    /* Who are my kids? */
    XQueryTree(display, current, &my_root, &parent, &children,
        &num_children);

    /* Examine all of the children too. */
    for( i = 0; i < num_children; i++ )
        Traverse_Tree( display, children[i] );

    /* Let's not waste any memory. */
    if( num_children ) XFree( (char *) children );
}


/* Wait_For_Idleness()                                                  */
/* Wait for "timeout" seconds of user inactivity.                       */

void
Wait_For_Idleness()
{
    int         watching = TRUE;
    int         found;
    int         timer = 0;
    XEvent      event;

    while( watching )
    {
        if( timer >= timeout ) return;

        sleep(1);
        timer++;

        found = XCheckIfEvent(display.dpy, &event, Dummy_Predicate, NULL);
        if( found ) timer = 0;

        /* Flush events. */
        while( found )
        {
            switch( event.type )
            {
                /* Watch for events in new windows too. */
                case CreateNotify:
                    {
                        XCreateWindowEvent *ev = (XCreateWindowEvent *) &event;

                        /* What an insidious bug! */
                        /* I had always assumed that when I received a */
                        /* creation notification that the window in */
                        /* question would still exist.  However, the */
                        /* big screen saver window can sometimes be */
                        /* destroyed by the time we get to this point. */
                        /* We certainly don't want to do anything here */
                        /* with a window that no longer exists. */
                        if( display.win != ev->window )
                        {
                            XSelectInput(display.dpy, ev->window, ALIVE_MASK);
                            Traverse_Tree(display.dpy, ev->window);
                        }
                    }
                    break;
                default:
                    break;
            }
            /* Check for the existence of more events. */
            found = XCheckIfEvent(display.dpy, &event, Dummy_Predicate, NULL);
        }
    }
}


Bool
Dummy_Predicate(display, event, arg)
Display *display;
XEvent  *event;
char    *arg;
{
    return(True);
}


void
Create_Big_Window()
{
    /* Window Attributes */
    unsigned long               valuemask;
    XSetWindowAttributes        xswa;

    /* First time. */
    static int          first = TRUE;

    /* Cursor Stuff */
    static Cursor       cursor;

    /* Turn the cursor the color of the background. (Invisible) */
    if( first )
    {
        /* I don't care which cursor I get. */
        cursor = XCreateFontCursor(display.dpy, 0);

        XRecolorCursor(display.dpy, cursor, &(display.bg_xcolor),
                        &(display.bg_xcolor));
        first = FALSE;
    }

    /* Create a screen sized window. */
    xswa.cursor = cursor;
    xswa.background_pixel = display.bg;
    xswa.override_redirect = True;
    xswa.do_not_propagate_mask = KeyPressMask | KeyReleaseMask |
        ButtonPressMask | ButtonReleaseMask;
    valuemask = CWCursor | CWBackPixel | CWOverrideRedirect | CWDontPropagate;
    display.win = XCreateWindow(display.dpy, display.root,
        0, 0,
        winW, winH, 0,
        DefaultDepth(display.dpy, display.screen),
        InputOutput, DefaultVisual(display.dpy, display.screen),
        valuemask, &xswa);

    XMapWindow(display.dpy, display.win);

    /* Event Mask */
    XSelectInput(display.dpy, display.win, KeyPressMask);

    /* Set up the stars' graphics context. */
    display.star_gc = XCreateGC(display.dpy, display.win, 0, NULL);
    XSetForeground(display.dpy, display.star_gc, display.star);
    XSetBackground(display.dpy, display.star_gc, display.bg);

    /* Set up an erasing graphics context. */
    display.erase_gc = XCreateGC(display.dpy, display.win, 0, NULL);
    XCopyGC(display.dpy, display.star_gc, 0xffffffff, display.erase_gc);
    XSetForeground(display.dpy, display.erase_gc, display.bg);

    /* Clear the background. */
    XFillRectangle(display.dpy, display.win, display.erase_gc,
            0,0, winW, winH);
}



/* Create_Window()                                                      */
/* Create the window, making sure to set it up correctly.               */

void
Create_Window(geometry)
char    *geometry;
{
    XSetWindowAttributes xswa;
    XSizeHints          sizehint;
    XWMHints            wmhints;
    char	        wname[256];     /* Window Name */

    if( !root )
    {
        xswa.event_mask = 0;
        xswa.background_pixel = display.bg;
        display.win = XCreateWindow(display.dpy, display.root,
            winX, winY,
            winW, winH, 0,
            DefaultDepth(display.dpy, display.screen),
            InputOutput, DefaultVisual(display.dpy, display.screen),
            CWEventMask | CWBackPixel , &xswa);

        sizehint.x = winX;
        sizehint.y = winY;
        sizehint.width = winW;
        sizehint.height = winH;
        sizehint.min_width = WINMIN;
        sizehint.min_height = WINMIN;
        if( geometry != NULL )
            sizehint.flags = USPosition | USSize | PMinSize;
        else
            sizehint.flags = PPosition | PSize | PMinSize;
        XSetNormalHints(display.dpy, display.win, &sizehint);

        display.protocol_atom = XInternAtom(display.dpy, "WM_PROTOCOLS", False);
        display.kill_atom = XInternAtom(display.dpy, "WM_DELETE_WINDOW", False);
#ifdef X11R4
        XSetWMProtocols(display.dpy, display.win, &display.kill_atom, 1);
#endif

        /* Title */
        if( PATCHLEVEL )
          sprintf( wname, "XStar %s.%d", VERSION, PATCHLEVEL);
        else
          sprintf( wname, "XStar %s", VERSION);
        XChangeProperty(display.dpy, display.win,
            XA_WM_NAME, XA_STRING, 8, PropModeReplace, (uchar_t	*) wname,
            strlen(wname));

        /* Window Manager Hints (This is supposed to make input work.) */
        wmhints.flags = InputHint;
        wmhints.input = True;
        XSetWMHints(display.dpy, display.win, &wmhints);

        XMapWindow(display.dpy, display.win);
    }
    else
    {
        display.win = display.root;
	xswa.backing_store = Always;
	XChangeWindowAttributes(display.dpy, display.win,
				CWBackingStore, &xswa);
    }

    /* Event Mask */
    if( root )
        XSelectInput(display.dpy, display.win,
                        KeyPressMask | StructureNotifyMask | ExposureMask);
    else
        XSelectInput(display.dpy, display.win,
		     KeyPressMask | ButtonPressMask | StructureNotifyMask
		     | ExposureMask );

    /* Set up the stars' graphics context. */
    display.star_gc = XCreateGC(display.dpy, display.win, 0, NULL);
    XSetForeground(display.dpy, display.star_gc, display.star);
    XSetBackground(display.dpy, display.star_gc, display.bg);

    /* Set up an erasing graphics context. */
    display.erase_gc = XCreateGC(display.dpy, display.win, 0, NULL);
    XCopyGC(display.dpy, display.star_gc, 0xffffffff, display.erase_gc);
    XSetForeground(display.dpy, display.erase_gc, display.bg);

    /* Clear the background. */
    XSetWindowBackground(display.dpy, display.win, display.bg);
    XFillRectangle(display.dpy, display.win, display.erase_gc,
      0,0, winW, winH);
}



/*
** HandleEvent()
**
** process X events
*/

void
HandleEvent(event, sd)
XEvent  *event;
star_disp_type *sd;
{
    /* If the screen saver is on, then watch for signs of activity. */
    if( ((event->type == KeyPress) || (event->type == MotionNotify))
       && (timeout)
       )
	stop = TRUE;

    switch( event->type )
    {
        case ClientMessage: /* sent by f.delete from twm */
            {
                XClientMessageEvent     *ev = (XClientMessageEvent *) event;

                printf("Client message received.\n");
                if( ev->message_type == display.protocol_atom
		   && ev->data.l[0] == display.kill_atom )
                    stop = TRUE;
            }
            break;

        case ConfigureNotify:
            {
                XConfigureEvent *ev = (XConfigureEvent *) event;

                /* The screen saver should ignore resizes. */
                if( !timeout )
                {
                    winW = ev->width;
                    winH = ev->height;
		    far_dist = (winW + winH)*1.5 + 750;
		    half_winW = winW/2;
		    half_winH = winH/2;
		    center_x = half_winW + root_xoff;
		    center_y = half_winH + root_yoff;
		    min_x = 0 - center_x - 20;
		    max_x = winW - center_x + 20;
		    min_y = 0 - center_y - 20;
		    max_y = winH - center_y + 20;
                    winX = ev->x;
                    winY = ev->y;
                }
            };
            break;

        case KeyPress:
            {
                XKeyEvent *key_event = (XKeyEvent *) event;
                char buf[128] = "";
                KeySym ks;
                XComposeStatus status;

                XLookupString(key_event,buf,128,&ks,&status);
                if( buf[0]=='q' || buf[0]=='Q' )
                    stop = TRUE;
                if( buf[0]=='d' || buf[0]=='D' )
                    add = TRUE;
                if( buf[0]=='+' )
                    add_star = TRUE;
                if( buf[0]=='-' )
                    del_star = TRUE;
                if( buf[0]=='e' || buf[0]=='E' )
                    erase = TRUE;
                if( buf[0]=='n' || buf[0]=='N' )
                    new_start = TRUE;
                if( buf[0]=='p' || buf[0]=='P' )
		    pause_updt = TRUE;
		if( buf[0]=='m' || buf[0]=='M' )
		{
		    multi_colors = !multi_colors;
		    if( rotate_colors && multi_colors )
		    {
			rotate_colors = 0;
		    }
		    if( init_colors_done == 0 )
			init_colors( colors, color_gcs, sd, NUM_COLORS );

		    do_redraw = TRUE;
		}
		if( buf[0]=='r' || buf[0]=='R' )
		{
		    rotate_colors = !rotate_colors;
		    if( rotate_colors && multi_colors )
		    {
			multi_colors = 0;
		    }
		    if( init_colors_done == 0 )
			init_colors( colors, color_gcs, sd, NUM_COLORS );

		    do_redraw = TRUE;
		}
		if( buf[0]=='\f' )
		    do_redraw = TRUE;
            }
            break;

	case Expose:
            {
		XExposeEvent	*ev = (XExposeEvent *) event;

		sd->redraw_min_x = ev->x;
		sd->redraw_max_x = ev->x + ev->width;
		sd->redraw_min_y = ev->y;
		sd->redraw_max_y = ev->y + ev->height;

		screen_exposed = TRUE;
	    }
	    break;
	    
        default:
            break;
    }
}


#ifdef USE_USLEEP

static void alarmhandler()
{
}

void sleepms(msec)
int msec;
{
    struct itimerval value,ovalue;
    struct sigvec vec;
    long savemask, sigblock(), sigpause();

    vec.sv_handler = alarmhandler;
    vec.sv_mask = 0x0;
    vec.sv_flags = 0;
    sigvector(SIGALRM, &vec, &vec); /* Set up alarmhandler for SIGALRM */
    savemask = sigblock((long)(1L << (SIGALRM - 1)));

    value.it_interval.tv_sec = 0;
    value.it_interval.tv_usec = 0;
    value.it_value.tv_sec = msec/1000;
    value.it_value.tv_usec = (msec%1000)*1000;
    setitimer(ITIMER_REAL,&value,&ovalue);

    (void)sigpause(0L);
    (void)sigsetmask(savemask);

    sigvector(SIGALRM, &vec, NULL); /* Restore previous signal handler */
}

void usleep(us)
unsigned us;
{
    sleepms(us / 1000);
}

#else

/*  Put the process to sleep for a while.  */
void nap(sec,usec)
long sec, usec;
{
    static struct timeval tv;

    tv.tv_sec = sec;
    tv.tv_usec = usec;
    select(0, 0, 0, 0, &tv);
}

#endif


void
Usage(program)
char *program;
{
    printf("%s [options]  where options are listed below\n", program);
    printf("-h                display this message\n");
    printf("-b stars          number of stars  (2 to %d)\n",
	   MAX_NUM_STARS - MAX_COLLAPSAR );
    printf("-r                use root window\n");
    printf("-g geometry       window geometry\n");
    printf("-d host:display   X server to connect to\n");
    printf("-t timeout        screen saved after 'timeout' seconds\n");
    printf("-D delay          delay between updates (milliseconds)\n");
    printf("-c star_clr       star color\n");
    printf("-C bg_color       background color\n");
    printf("-R                Rotate the star colors  (Rainbow)\n");
    printf("-M                multiple colors, one per star\n");
    printf("-a float          accuracy of position calculations.   Larger values\n" );
    printf("                  increase accuracy.  Must be greater than zero.\n" );
    printf("-m [taylor3|rk4|gpemce8|ab7|am7]     method used to calculate points\n" );
    printf("-l float          distance between stars before they collide\n");
    printf("-T num_pts        number of points to display as star trails\n");
    printf("-B buf_factor     amount of buffering to be done with the X server.\n" );
    printf("-v                enable display of debug output\n");
    printf("\nDefaults are equivalent to the following command:\n");
    printf("%s -b %d -g =%dx%d -D %d -c White -C Black -a 1.0 -m ab7 -l %.1f -T %d -B 1.0\n",
	   program, STARS, WINWIDTH, WINHEIGHT, DELAY,
	   DEFAULT_COLLIDE, DEFAULT_DISP_PT );
    printf("\nCommands available while running are:\n" );
    printf("  d  create a collapsar (gravity well).\n");
    printf("  e  erase trails.\n");
    printf("  +  add a star to the system.\n");
    printf("  -  delete a star from the system.\n");
    printf("  n  get a new set of stars.\n");
    printf("  m  toggle multi-color mode.\n");
    printf("  r  toggle rainbow mode.\n");
    printf("  p  pause the updating.  Press p again to start.\n" );
    printf("  ^L redraw the screen.\n" );
    printf("  q  quit.\n\n");
}


void
HandleError(description, degree)
char    *description;
int     degree;
{
    fprintf(stderr, "An error has occurred.  The description is below...\n");
    fprintf(stderr, "%s\n", description);

    if( degree == FATAL )
    {
        fprintf(stderr, "Program aborting...\n");
        exit(-1);
    }
}

long
GetColor(display, color, final_color)
disp            *display;
char            *color;
XColor          *final_color;
{
    XColor      cdef;
    char        error_str[STD_STR];

    if( !XParseColor(display->dpy, display->cmap, color, &cdef) ||
        !XAllocColor(display->dpy, display->cmap, &cdef) )
    {
        sprintf(error_str, "Color \"%s\" wasn't found.", color);
        HandleError(error_str, FATAL);
    }

    /* Copy the final color. */
    if( final_color != NULL ) *final_color = cdef;

    return(cdef.pixel);
}


/* Check the return code of malloc(). */
void *
Safe_Malloc(bytes)
int bytes;
{
    void *pointer;

    pointer = (void *) malloc(bytes);
    if( NULL == pointer )
    {
        fprintf(stderr, "Error allocating %d bytes.\n", bytes);
        exit(-1);
    }

    return(pointer);
}

void *
Safe_Realloc(ptr, bytes)
void *ptr;
int bytes;
{
    void *pointer;

    pointer = (void *) realloc(ptr, bytes);
    if( NULL == pointer )
    {
        fprintf(stderr, "Error allocating %d bytes.\n", bytes);
        exit(-1);
    }

    return(pointer);
}
END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file init_stars.c ----------

filename="init_stars.c"

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



#include <time.h>

/*
 * Main driver routine for initializing the star system.
 */

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

    point_2d	*v_0 = vk[0];

    int		rnd_colors[ NUM_COLORS ];
    int		tmp_color;


    /* clean up from last time */
    sd->tstamp = time(NULL);
    if( verbose > 1 || (verbose == 1 && last_init == 0) )
    {
	fprintf( stderr, "                 started with                |             ended with\n" );
	fprintf( stderr, "            gravity                     no   | life\n" );
	fprintf( stderr, "stars/bin/++  well  size  cir lin bnce speed | span  collide  edge  ++    steps\n" );
    }

    if( verbose > 0 && last_init != 0 )
    {
	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 );

	if( verbose > 3 )
	{
	    prt_sys_const( max_stars, cur, m, v_0 );
	    prt_sys_const_delta( orig_sc, max_stars, cur, m, v_0 );
	}
    }

    if( last_init != 0 )
    {
	update_screen( sp, sd, &display );
	XFlush( display.dpy );
    }
    clear_disp_pt( sd );

    last_init = sd->tstamp;
    num_collide = 0;
    num_edge = 0;


    /* make sure the array is clean at the end... */
    for( i = 0; i < max_stars; i++ )
    {
	prev[i].x = cur[i].x = 0;
	prev[i].y = cur[i].y = 0;
	m[i] = 0;
	sd->last_disp[i].x = far_dist;
	sd->last_disp[i].y = far_dist;

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


#if 1
    /* setup up the system parameters */
    set_sys_param( sp, max_stars );

    /* set the mass, locations, velocities and accelerations */
    set_xmva( cur, prev, m, vk, ak, max_stars, sp );
#else
    /* use a predefined star system */
    init_sys_4( cur, prev, m, vk, ak, max_stars, sp );
/*    init_sys_8( cur, prev, m, vk, ak, max_stars, sp ); */
/*    init_sys_8b( cur, prev, m, vk, ak, max_stars, sp ); */
#endif
    
    /* remember the original system parameters */
    orig_sp = *sp;


    /* set up the star_disp variables */
    set_star_disp( sd, sp );
    set_buffer_factor( sd, sp );


    /* do any movement initialization that is required */
    calc_sys_const( &orig_sc, max_stars, cur, m, v_0 );
    if( (orig_sc.k + orig_sc.u)*fv*fv/fm > -1E-5 )
    {
	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] );
    }
    init_fna( max_stars, m, cur, v_0 );


    /*
     * set up a random star color system for multi-color mode
     */
    for( i = 0; i < NUM_COLORS; i++ )
	rnd_colors[ i ] = i;
    
    for( i = 0; i < NUM_COLORS-1; i++ )
    {
	j = i + lrand48() % (NUM_COLORS - i);
	tmp_color = rnd_colors[ i ];
	rnd_colors[ i ] = rnd_colors[ j ];
	rnd_colors[ j ] = tmp_color;
    }

    for( i = 0; i < max_stars; i++ )
	sd->star_color[ i ] = rnd_colors[ i % NUM_COLORS ];
    


    /*
     * print out the star configuration
     */

    if( verbose > 2 )
    {
	double	u, k;
	double	dx, dy;
	double	tdist;
	
	fprintf( stderr, "star   location         velocity     mass   tdist  kinetic   binding   -k/u\n" );
    
	for( i = 0; i < sp->live_stars; i++ )
	{
	    dx = cur[i].x;
	    dy = cur[i].y;
	    tdist = sqrt( dx*dx + dy*dy );
	    
	    k = kinetic_energy( i, max_stars, cur, m, v_0 );
	    u = binding_energy( i, max_stars, cur, m, v_0 );
	    
	    fprintf( stderr,
		    "%2d  %6.1f,%6.1f  %7.4f,%7.4f %6.3f  %6.2f %9.6f %9.6f %9.6f\n",
		    i, cur[i].x, cur[i].y, fv*v_0[i].x, fv*v_0[i].y, m[i]/fm,
		    tdist, k*fv*fv/fm, u*fv*fv/fm, -k/u );
	}
    }
#if 0
    /* dump the star system configuration for later reuse */
    dump_sys( cur, prev, m, vk, ak, max_stars, sp );
#endif
    

    /* clean up the display */
    sleep(1);
    XFillRectangle(display.dpy, display.win, display.erase_gc,
		   0,0, winW, winH);
    plot_collapsars( cur, m, sp, sd, &display );
    XFlush( display.dpy );
}



/*
 * add a new star to the system or reinitialize it
 */

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

void new_star( point_2d *cur, point_2d *prev, double *m, 
	      point_2d *vk[], point_2d *ak[],
	      int max_stars, sys_param_type *sp, star_disp_type *sd,
	      int add_collapsar )
{
    int		i, j;
    int		found;
    double	dx, dy;
    double	tdist;

    double	dist, cur_dist, max_dist;
    double	theta;
    


    static int	check = 0;


    /*
     * find a good slot to use.  (rotate through so we get different colors)
     */
    for( i = 0, found = -1; i < max_stars; i++ )
    {
	check++;
	if( check >= max_stars )
	    check = check % max_stars;
	
	if( m[check] <= 0 )
	{
	    found = check;
	    break;
	}
    }
    
    if( (sp->num_add)-- <= 0 || found == -1 )
    {
	init_stars( prev, cur, m, vk, ak, max_stars, sp, sd );
	return;
    }
    

    /*
     * create the star
     */

    max_dist = far_dist * far_dist;
    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;
	}

	if( add_collapsar && !sp->do_bounce )
	{
	    prev[found].x = cur[found].x = dist * .5 * cos( theta );
	    prev[found].y = cur[found].y = dist * .5 * sin( theta );
	}
	else
	{
	    prev[found].x = cur[found].x = dist * cos( theta );
	    prev[found].y = cur[found].y = dist * sin( theta );
	}
	
	
	/*
	 * don't put stars too close together
	 */
	cur_dist = far_dist * far_dist;
	for( j = 0; j < max_stars; j++ )
	{
	    if( found == j )
		continue;

	    if( m[j] <= 0 )
		continue;
	    
	    dx = prev[found].x - prev[j].x;
	    dy = prev[found].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 )
	  );
    
    if( add_collapsar )
    {
	m[found] = collapsar * 1.0001;
	vk[0][found].x = vk[0][found].y = ak[0][found].x = ak[0][found].y = 0;
    }
    else
    {
	if( sp->num_collapsar )
	    m[found] = fm*.1;
	else
	    m[found] = fm*.5;

	set_va( cur, prev, m, vk, ak, sp, found );
    }

    sd->last_disp[found].x = far_dist;
    sd->last_disp[found].y = far_dist;


    (sp->live_stars)++;
    (sd->live_erase)++;
    sd->tstamp = time(NULL);
    set_buffer_factor( sd, sp );

    if( verbose > 2 )
    {
	dx = cur[found].x;
	dy = cur[found].y;
	tdist = sqrt( dx*dx + dy*dy );
	
	fprintf( stderr,
		"%s:%d  star %2d @ (%6.1f,%6.1f)=|%7.4f,%7.4f|*%6.3f  tdist=%6.2f\n",
		__FILE__, __LINE__, found, cur[found].x, cur[found].y, fv*vk[0][found].x,
		fv*vk[0][found].y, m[found]/fm, tdist );
    }


    /* do any movement initialization that is required */
    init_fna( max_stars, m, cur, vk[0] );

    /* make sure collapsars get plotted... */
    if( add_collapsar )
	plot_collapsars( cur, m, sp, sd, &display );
}

END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file init_colors.c ----------

filename="init_colors.c"

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



/* straight out of Foley and Van Dam... */
static void hsv_to_rgb( double h, double s, double v,
		double *r, double *g, double *b )
{
    double f;
    double p;
    double q;
    double t;
    int sextant;
    
    if( h == 360.0 )
	h = 0.0;

    h /= 60.0;
    sextant = h;
    f = h - sextant;
    p = v * (1.0 - s);
    q = v * (1.0 - s * f);
    t = v * (1.0 - s * (1.0 - f));
    
    switch( sextant )
    {
    case 0 :
	*r = v;
	*g = t;
	*b = p;
	break;
	
    case 1 :
	*r = q;
	*g = v;
	*b = p;
	break;
	
    case 2 :
	*r = p;
	*g = v;
	*b = t;
	break;
	
    case 3 :
	*r = p;
	*g = q;
	*b = v;
	break;
	
    case 4 :
	*r = t;
	*g = p;
	*b = v;
	break;
	
    case 5 :
	*r = v;
	*g = p;
	*b = q;
	break;
    }
}


/*
 * initialize the color array to form a rainbow
 */

void init_colors( XColor *colors, GC *color_gcs, star_disp_type *sd, int num_elem )
{
    int			c;

    double		h, s = 1, v = 1;
    double		r, g, b;

    double		h_inc = 360. / NUM_COLORS;
    
    char        error_str[STD_STR];


    for( c = 0; c < NUM_COLORS; c++ )
    {
	h = c * h_inc;
	hsv_to_rgb( h, s, v, &r, &g, &b );


	/* allocate another color */
	colors[c].red = r * MAX_CVAL;
	colors[c].green = g * MAX_CVAL;
	colors[c].blue = b * MAX_CVAL;
	if( !XAllocColor(display.dpy, display.cmap, &colors[ c ]) )
	{
	    sprintf(error_str, "Color %d=(%d,%d,%d) couldn't be allocated.",
		    c, colors[c].red, colors[c].green,
		    colors[c].blue );
	    HandleError(error_str, FATAL);
	}
#if 0
	fprintf( stderr, "%s:%d  color[%d]=(%d,%d,%d)=%08X\n", __FILE__, __LINE__, c, colors[c].red, colors[c].green, colors[c].blue, colors[c].pixel );
#endif
	
	color_gcs[c] = XCreateGC(display.dpy, display.win, 0, NULL);
	XSetForeground( display.dpy, color_gcs[c], colors[c].pixel);
    }

    init_colors_done = 1;
}

END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file global_vars.c ----------

filename="global_vars.c"

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



/* X related */
int		winX, winY;
int		root_xoff, root_yoff;
double		center_x, center_y;
uint_t		winW, winH;
int		half_winW, half_winH;
int		max_x, min_x, max_y, min_y;
uint_t		far_dist;
disp		display;

/* animation related */
int     stars = STARS;          /* number of stars */
int     max_stars = STARS+MAX_COLLAPSAR;    /* maximum number of stars */
int     max_disp_skip;			/* max num of displays to skip */
int     delay = DELAY * 1000;   /* delay between updates, in microseconds */
int     timeout = 0;            /* time in seconds before screen saving */
char    *star_color = NULL;
char    *bg_color = NULL;
XColor	colors[NUM_COLORS];
GC	color_gcs[NUM_COLORS];
int	verbose = 0;
int	rotate_colors = 0;
int	init_colors_done = 0;
int	multi_colors = 0;
int	num_disp_pt = DEFAULT_DISP_PT;

int     stop = FALSE;
int     add  = FALSE;           /* flag to add a collapsar */
int     add_star  = FALSE;      /* flag to add a star */
int     del_star  = FALSE;      /* flag to delete a star */
int     erase = FALSE;          /* flag to erase trails */
int     new_start = FALSE;      /* flag to start a new system trails */
int	pause_updt = FALSE;	/* flag to pause the updates	*/
int	screen_exposed;		/* flag to update the screen	*/
int	do_redraw;		/* flag to redraw the screen	*/

int     root = FALSE;           /* display in root window */
time_t	last_init = 0;

double	collapsar = (FM*COLLAPSAR);
double	fv = FV;
double	fv_inv = 1./FV;
double	fm = FM;

sys_param_type	orig_sp;
sys_const orig_sc;
star_disp_type  orig_sd;

int num_collide;
int num_edge;

double	collide_dist = DEFAULT_COLLIDE*DEFAULT_COLLIDE;

void (*move_fna)( 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 )
    = move_stars_ab7;
    
void (*init_fna)( int n, double *m, point_2d *x, point_2d *v )
    = init_ab7;

void (*restart_fna)( int n, double *m, point_2d *x, point_2d *v )
    = restart_ab7;

double move_cost = 1.2;
END-OF-FILE

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

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

  chmod 660 $filename
fi

# ---------- file check_bounce.c ----------

filename="check_bounce.c"

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

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


/*
 * Check_bounce() is _really_ expensive.  Many banks charge $25 for the first
 * time, and many states require you to take 'classes' if you bounce checks.
 */

/*
** 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 code checks to see if any stars should bounce off of each other.  It
 * mostly works, but for multi-step routines like taylor3 or the Adam-Bashford
 * routines, it is fairly slow because it requires a movement restart.  There
 * are other inefficiencies in the code as well.  Fortunately, it only
 * calculates these expensive formulas when an actual collision is happening,
 * and that doesn't happen very often.
 *
 * This routine does not work in all cases though.  For some reason, if you
 * have more than two stars colliding at the same spot, the results are not
 * correct.  You can get this when you take a star circle, enable the bounce
 * code and start with no speed.  Then, all the stars drop toward the center
 * of the circle and when they get there, they should all bounce back to the
 * edge along the same path.  This only happens for n=2,4, for other values
 * n, the stars go off at random angles.
 *
 * [ I have checked into this a little bit more... It appears that the bounce
 *   code is working "correctly", if you take the "duration of the impulse" as
 *   being zero (or close too it) and you count stars with lower index numbers
 *   to be "slightly" ahead of the stars with higher index numbers.  Besides
 *   the fix described in the next paragraph, this problem might also be
 *   "fixed" by accumulating all the bounce effects and then somehow applying
 *   them after all stars have been checked.  I am not sure that this would
 *   be simpler than the following, better, fix... ]
 *
 * I have thought about how maybe we should bounce stars via the same method
 * that real things bounce.  That is, when objects get close enough together,
 * the atoms no longer look like neutral balls, and instead the electrical
 * repulsion between the electrons stars to take over.  If we created a very
 * short distance, very large magnitude negative force field around each
 * star, I think we could get rid of the above problem.  This would create
 * the non-zero duration of the impulse that real world items have.
 */

#define bounce_dist	(4)		/* large enough to be visible on scr */

/* #define	bounce_debug */


void check_bounce( 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 )
{
    register double	m_i, m_j;
    
    register int        i, j;		/* star index */
    
    point_2d	*v_1 = vk[1];
    
    
    
    /* bounce stars off of each other */
    for( i = 0; i < max_stars; i++ )
    {
	m_i = m[i];
	if( m_i <= 0.0 )
	    continue;
	
	
	for( j = i+1; j < max_stars; j++ )
	{
	    double dx, dy, n2;
	    point_2d tvi, tvj;
	    
	    m_j = m[j];
	    if( m_j <= 0.0 )
		continue;
	    
	    dx = prev[j].x - prev[i].x;
	    dy = prev[j].y - prev[i].y;
	    n2 = dx*dx + dy*dy;
	    
	    
	    /* this is probably too large of a value... */
	    if( (m[j] < collapsar && n2 < bounce_dist*bounce_dist )
	       || (m[j] >= collapsar
		   && n2 < bounce_dist*bounce_dist*bounce_dist
		   )
	       )
	    {
		/* collision:  perfectly elastic -> bounce */
		
		double n = sqrt(n2);
		
		
		/* projected velocities along alternate coordinate sys */
		point_2d	v1xp, v1yp, v2xp, v2yp;
		point_2d	xp, yp;
		
		/*
		 * figure out new coordinate system
		 *
		 * The x-prime axis is along the line connecting the
		 * two bodies.  The y-prime axis is "flat wall" that the
		 * bodies appear to bounce off of.  The momentum and
		 * velocity along the y-prime axis is not changed by
		 * the collision.
		 */
		xp.x = dx / n;
		xp.y = dy / n;
		
		yp.x = xp.y;
		yp.y = -xp.x;

		
		/* figure out projections of velocities along new axes */
		v1yp.x = (yp.x * v_1[i].x + yp.y * v_1[i].y) * yp.x;
		v1yp.y = (yp.x * v_1[i].x + yp.y * v_1[i].y) * yp.y;
		
		v2yp.x = (yp.x * v_1[j].x + yp.y * v_1[j].y) * yp.x;
		v2yp.y = (yp.x * v_1[j].x + yp.y * v_1[j].y) * yp.y;
		
		v1xp.x = (xp.x * v_1[i].x + xp.y * v_1[i].y) * xp.x;
		v1xp.y = (xp.x * v_1[i].x + xp.y * v_1[i].y) * xp.y;
		
		v2xp.x = (xp.x * v_1[j].x + xp.y * v_1[j].y) * xp.x;
		v2xp.y = (xp.x * v_1[j].x + xp.y * v_1[j].y) * xp.y;
		
		/*
		 * basically, we just need to treat the collision as a
		 * one dimensional collision along the x-prime axis.
		 */
		if( m_i >= collapsar )
		{
		    tvi = v_1[i];
		
		    tvj.x = v2yp.x - v2xp.x;
		    tvj.y = v2yp.y - v2xp.y;
		}
		else if( m_j >= collapsar )
		{
		    tvi.x = v1yp.x - v1xp.x;
		    tvi.y = v1yp.y - v1xp.y;

		    tvj = v_1[j];
		}
		else
		{
		    tvi.x = v1yp.x + (m_i * v1xp.x - m_j * v1xp.x + 2 * m_j * v2xp.x) / (m_i + m_j);
		    tvi.y = v1yp.y + (m_i * v1xp.y - m_j * v1xp.y + 2 * m_j * v2xp.y) / (m_i + m_j);
		
		    tvj.x = v2yp.x + (m_j * v2xp.x - m_i * v2xp.x + 2 * m_i * v1xp.x) / (m_i + m_j);
		    tvj.y = v2yp.y + (m_j * v2xp.y - m_i * v2xp.y + 2 * m_i * v1xp.y) / (m_i + m_j);
		}
		
		
		/*
		 * There has got to be a better way to see if the
		 * stars are moving away from each other than
		 * checking the distances like this.  This is way
		 * too expensive to calculate.
		 *
		 * the problem that I am fighting is that once inside the
		 * bounce distance, this code always says that the stars
		 * are colliding, even though after the first bounce, they
		 * are moving away from each-other.
		 *
		 * we know that the tvi and tvj vectors are collinear,
		 * and that should make things easier to calculate...
		 */
		dx = (prev[j].x + tvj.x) - (prev[i].x + tvi.x);
		dy = (prev[j].y + tvj.y) - (prev[i].y + tvi.y);
		
		if( n2 < (dx*dx + dy*dy) )
		{
#ifdef bounce_debug
		    static double	last_mag = 0, cur_mag = 0;
		    
		    fprintf( stderr, "\n%s:%d  i=(%g,%g)  j=(%g,%g)  delta=(%g,%g)  n=%g\n", __FILE__, __LINE__, prev[i].x, prev[i].y, prev[j].x, prev[j].y, dx, dy, n );
#define rtd(x)	((x)*360/(M_PI*2))
		    fprintf( stderr, "%s:%d  xp=(%g,%g) |%g| <%g\n",
			    __FILE__, __LINE__,
			    xp.x, xp.y,
			    sqrt((xp.x*xp.x) + (xp.y*xp.y)),
			    rtd( atan2( xp.y, xp.x ) ) );

		    fprintf( stderr, "%s:%d  yp=(%g,%g) |%g| <%g\n",
			    __FILE__, __LINE__,
			    yp.x, yp.y,
			    sqrt((yp.x*yp.x) + (yp.y*yp.y)),
			    rtd( atan2( yp.y, yp.x ) ) );

		    fprintf( stderr, "%s:%d  stars: %d,%d  v1=(%g,%g) |%g| <%g  v1yp=(%g,%g) v1xp=(%g,%g)  (%g,%g)\n",
			    __FILE__, __LINE__, i, j,
			    v_1[i].x, v_1[i].y,
			    sqrt((v_1[i].x*v_1[i].x) + (v_1[i].y*v_1[i].y)),
			    rtd( atan2( -v_1[i].y, -v_1[i].x ) ),
			    v1yp.x, v1yp.y, v1xp.x, v1xp.y,
			    v1xp.x + v1yp.x, v1yp.y + v1xp.y );
		    fprintf( stderr, "%s:%d  v1=(%g,%g) |%g| <%g\n",
			    __FILE__, __LINE__,
			    tvi.x, tvi.y,
			    sqrt((tvi.x*tvi.x) + (tvi.y*tvi.y)),
			    rtd( atan2( tvi.y, tvi.x ) ) );
		    fprintf( stderr, "%s:%d  stars: %d,%d  v2=(%g,%g) |%g| <%g  v2yp=(%g,%g) v2xp=(%g,%g)  (%g,%g)\n",
			    __FILE__, __LINE__, i, j,
			    v_1[j].x, v_1[j].y,
			    sqrt((v_1[j].x*v_1[j].x) + (v_1[j].y*v_1[j].y)),
			    rtd( atan2( -v_1[j].y, -v_1[j].x ) ),
			    v2yp.x, v2yp.y, v2xp.x, v2xp.y,
			    v2xp.x + v2yp.x, v2yp.y + v2xp.y );
		    fprintf( stderr, "%s:%d  v2=(%g,%g) |%g| <%g\n",
			    __FILE__, __LINE__,
			    tvj.x, tvj.y,
			    sqrt((tvj.x*tvj.x) + (tvj.y*tvj.y)),
			    rtd( atan2( tvj.y, tvj.x ) ) );
#endif
		    
		    
		    v_1[i] = tvi;
		    v_1[j] = tvj;

		    restart_fna( max_stars, m, prev, v_1 );

#ifdef bounce_debug
#if 0
		    cur_mag = sqrt((v_1[i].x*v_1[i].x) + (v_1[i].y*v_1[i].y));
		    
		    fprintf( stderr, "%s:%d  v1=(%g,%g) %g  diff=%.3g%%\n",
			    __FILE__, __LINE__,
			    v_1[i].x, v_1[i].y, cur_mag,
			    100.*(last_mag/cur_mag-1.) );
		    last_mag = cur_mag;
#endif
#endif
		}
#ifdef bounce_debug
		else
		    fprintf( stderr, "%s:%d  ignoring the bounce...\n", __FILE__, __LINE__ );
#endif
	    }
	}
	
    }
}
END-OF-FILE

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

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


