/* Next available MSG number is 141 */
/*

   GRAVITY.C

   Copyright (C) 1989, 1990, 1991, 1992, 1993, 1994 by Autodesk, Inc.
 
   Permission to use, copy, modify, and distribute this software in 
   object code form for any purpose and without fee is hereby granted, 
   provided that the above copyright notice appears in all copies and 
   that both that copyright notice and the limited warranty and 
   restricted rights notice below appear in all supporting 
   documentation.
 
   AUTODESK PROVIDES THIS PROGRAM "AS IS" AND WITH ALL FAULTS.  
   AUTODESK SPECIFICALLY DISCLAIMS ANY IMPLIED WARRANTY OF 
   MERCHANTABILITY OR FITNESS FOR A PARTICULAR USE.  AUTODESK, INC.
   DOES NOT WARRANT THAT THE OPERATION OF THE PROGRAM WILL BE 
   UNINTERRUPTED OR ERROR FREE.
 
   Use, duplication, or disclosure by the U.S. Government is subject to 
   restrictions set forth in FAR 52.227-19 (Commercial Computer 
   Software - Restricted Rights) and DFAR 252.227-7013(c)(1)(ii) 
   (Rights in Technical Data and Computer Software), as applicable.
    
   .

        N-Body Gravitational Interaction Simulator for AutoCAD

        A sample ADS application.

        Designed and implemented by John Walker in August of 1989.

           "Ubi materia, ibi geometria."
                      -- Johannes Kepler

           "The orbit of  any  one  planet  depends  on  the
           combined  motion  of  all  the  planets,  not  to
           mention the action of all of these on each other.
           But  to  consider simultaneously all these causes
           of motion and to define these  motions  by  exact
           laws   allowing   of   conventional   calculation
           exceeds, unless I am mistaken, the force  of  the
           entire human intellect."
                      -- Sir Isaac Newton, Principia


        The  physics  underlying  this simulation are explained in the
        the  chapter  "A  Cosmic  Ballet",  pp.   229-238  of  A.   K.
        Dewdney,  "The  Armchair Universe", W.  H.  Freeman: New York,
        1988.

        The following commands are implemented in this file:

        DEMO        Creates  one  of  a series of standard demo models
                    when executed in  a  new,  empty  drawing.   These
                    models  may be modified with the MASS and MASSEDIT
                    commands just like models entered manually.

        FRAME       Asks you to pick a mass entity.   Motion  is  then
                    displayed in that mass's reference frame (in other
                    words, that mass  appears  stationary  and  others
                    move  around  it).  The default reference frame is
                    "Inertial", an unaccelerated frame  at  rest  with
                    respect  to  the  distant  stars.   Specifying  no
                    entity to the  FRAME  command  re-establishes  the
                    inertial frame.

        MASS        Creates  a  new  mass.   You're invited to enter a
                    name for the mass, its position (specified by  any
                    means  of  co-ordinate  specification), a velocity
                    vector in units of  astronomical  units  per  year
                    (which  can  be either entered explicitly or drawn
                    by typing "V", then dragging the endpoint  of  the
                    velocity  vector to the correct position), and its
                    mass in units of Solar masses.

        MASSEDIT    Lets you modify the name, velocity, and mass of an
                    existing mass.  Pick a single mass by pointing.  A
                    dialogue is displayed with the properties  of  the
                    mass.   Change  them  as you wish, then pick OK to
                    update the properties of the mass in the database.
                    If you pick Cancel, the mass is not modified.

        RESET       When  a  simulation  is  run,  it  halts after the
                    specified number of steps or  simulated  time  has
                    elapsed.    A   subsequent  RUN  command  normally
                    resumes the simulation from the point at which the
                    last  stopped.  RESET erases all motion paths from
                    the screen and sets the  simulation  back  to  the
                    start.    RUN   after   a  RESET  will  begin  the
                    simulation from the initial state.

        RUN         Starts the simulation.  You're  asked  to  specify
                    the length of the simulation as either a number of
                    motion steps or by the number of simulated  years.
                    If  you  enter  a positive number, it's taken as a
                    step count.  Negative numbers (which  may  include
                    decimal   fractions)  specify  simulated  time  in
                    years.  The length in time of each simulation step
                    varies  based  upon the velocity and separation of
                    masses, so if objects approach  one  another  very
                    closely  a  very  large  number  of  steps will be
                    required to simulate a given time  interval.   The
                    calculation time per step is essentially constant,
                    so if you're investigating a system  with  unknown
                    behaviour  you'll probably want to RUN for a given
                    number of cycles, but when running a stable system
                    such as the Solar System, you can simulate a given
                    time span.  In  any  case,  you  can  terminate  a
                    simulation with the Control C key.  RUN  can  also
                    be  invoked as a function, (C:RUN <length>), where
                    <length> specifies  the  simulation  length  by  a
                    positive  or  negative  number as described above.
                    When called  as  a  function,  C:RUN  returns  the
                    simulation end time as its result.

        SETGRAV     The  Gravity  simulator has several variables that
                    control its operation.  These variables are  saved
                    with the drawing and may be inspected and modified
                    with the  SETGRAV  command.   SETGRAV  displays  a
                    dialogue  in  which  you  may  change  any  of the
                    following variables:

                       Output to display only?
                          This is set to Yes or  No  (True/False,  and
                          1/0  are  also  accepted).   If  set  to the
                          default value of Yes, motion paths are drawn
                          using  temporary  vectors within the current
                          viewport.  If the picture is  REDRAWn,  they
                          disappear.   If  set to No, motion paths are
                          added to the drawing as LINE entities.  This
                          causes paths to appear in all viewports, and
                          you can  ZOOM  on  sections  of  a  path  to
                          examine   it   in   greater  detail.   Paths
                          represented as  LINEs  are  saved  with  the
                          drawing.  Generating LINE entities for paths
                          is much slower than just drawing them on the
                          screen,  so  choose  this mode only when you
                          need it.

                       Display step number?
                          If this mode is  "Yes"  (the  default),  the
                          simulation  step  number is displayed in the
                          coordinate status  line  as  the  simulation
                          progresses.

                       Display time?
                          If "Yes" (the default) the simulated time in
                          years is displayed in the coordinate  status
                          line.

                       Step size?
                          This real number specifies the  factor  used
                          to  determine  the  size  of the integration
                          steps used in calculating the motion of  the
                          masses.   The  time  step  is  calculated by
                          dividing  the  distance  between   the   two
                          closest   masses  by  the  highest  relative
                          velocity  of  any  pair  of  masses.    This
                          quantity,  measured  in years, is multiplied
                          by the factor  given  by  this  variable  to
                          obtain  the length of the step.  The default
                          value  of  0.1  works  well  for  reasonably
                          well-behaved  simulations.  If you find that
                          a simulation is taking  too  long,  you  can
                          speed it up by increasing the step size, but
                          be aware that the results you see may not be
                          physically  accurate.   Increasing  the step
                          size magnifies the inaccuracies of  modeling
                          a  continuous process such as gravitation by
                          discrete  steps.   In   general,   you   can
                          increase the step size as long as you obtain
                          the same results.  When the  outcome  begins
                          to vary, you've set the step size too large.

                       Minimum step?
                          After  the  step  size  is   calculated   as
                          described  above  it  is  compared  with the
                          minimum step size and, if less, the  minimum
                          is  used.   The  minimum step size is set to
                          0.00001 years (about  5  minutes);  this  is
                          sufficient  to unstick many simulations that
                          involve a close encounter.  Amazingly,  even
                          a  minimum  this  short  can  yield  grossly
                          inaccurate results when solar masses execute
                          hairpin turns about one another

        UPDATE      A  simulation  normally  proceeds from the initial
                    positions, velocities, and masses of  the  objects
                    specified  by  the  MASS  command.  The simulation
                    keeps track of these values as it  progresses  but
                    does  not automatically adjust the mass objects in
                    the database.  If you want to  move  the  database
                    objects  to  the  positions at the end of the most
                    recent simulation (and adjust their velocities  to
                    the  corresponding  instantaneous velocities), use
                    the UPDATE command.  If you don't  do  an  UPDATE,
                    the masses will remain in their original positions
                    even if you  save  the  drawing  after  running  a
                    simulation.

*/

#include   <stdio.h>
#include   <string.h>
#ifndef Macintosh
#include   <ctype.h>
#endif
#include   <math.h>
#include    <assert.h>
#include   "adslib.h"

/*  Standard drawing object names  */

/* Utility frozen layer for information */
#define FrozenLayer /*MSG1*/"FROZEN-SOLID"
#define OrbitLayer  /*MSG2*/"ORBITS"  /* Layer for orbital path traces */

/*  Data types  */

typedef enum {False = 0, True = 1} Boolean;

#define HANDLEN  18                   /* String long enough to hold a handle */

/* Definitions  to  wrap  around  submission  of  AutoCAD commands to
   prevent their being echoed.  */

#define Cmdecho  False                /* Make True for debug command output */

#define CommandB()  { struct resbuf rBc, rBb, rBu, rBh; \
        ads_getvar(/*MSG0*/"CMDECHO", &rBc); \
        ads_getvar(/*MSG0*/"BLIPMODE", &rBb); \
        ads_getvar(/*MSG0*/"HIGHLIGHT", &rBh); \
        rBu.restype = RTSHORT; \
        rBu.resval.rint = (int) Cmdecho; \
        ads_setvar(/*MSG0*/"CMDECHO", &rBu); \
        rBu.resval.rint = (int) False; \
        ads_setvar(/*MSG0*/"BLIPMODE", &rBu); \
        ads_setvar(/*MSG0*/"HIGHLIGHT", &rBu)

#define CommandE()  ads_setvar(/*MSG0*/"CMDECHO", &rBc); \
                    ads_setvar(/*MSG0*/"BLIPMODE", &rBb); \
                    ads_setvar(/*MSG0*/"HIGHLIGHT", &rBh); }

/*  Definitions  that permit you to push and pop system variables with
    minimal complexity.  These don't work  (and  will  cause  horrible
    crashes  if  used)  with  string  variables,  but since all string
    variables are read-only, they cannot be saved and restored in  any
    case.  */

#define PushVar(var, cell, newval, newtype) { struct resbuf cell, cNeW; \
        ads_getvar(var, &cell); cNeW.restype = cell.restype;             \
        cNeW.resval.newtype = newval; ads_setvar(var, &cNeW)

#define PopVar(var, cell) ads_setvar(var, &cell); }

/* Set point variable from three co-ordinates */

#define Spoint(pt, x, y, z)  pt[X] = (x);  pt[Y] = (y);  pt[Z] = (z)

/* Copy point  */

#define Cpoint(d, s)   d[X] = s[X];  d[Y] = s[Y];  d[Z] = s[Z]

/* Generation parameters for objects */

#define SphereLong  12                /* Longitudinal tabulations on sphere */
#define SphereLat   12                /* Latitudinal tabulations on sphere */

/* Utility definition to get an  array's  element  count  (at  compile
   time).   For  example:

       int  arr[] = {1,2,3,4,5};
       ...
       printf("%d", ELEMENTS(arr));

   would print a five.  ELEMENTS("abc") can also be used to  tell  how
   many  bytes are in a string constant INCLUDING THE TRAILING NULL. */

#define ELEMENTS(array) (sizeof(array)/sizeof((array)[0]))

/* Utility definitions */

#ifndef abs
#define abs(x)      ((x)<0 ? -(x) : (x))
#endif  /* abs */
#ifdef min
#undef  min
#endif
#define min(a,b)    ((a)<(b) ? (a) : (b))
#ifdef max
#undef  max
#endif
#define max(a,b)    ((a)>(b) ? (a) : (b))

/*  Many C implementations may lack a cube root function.  Rather than
    count  on  a system cbrt() function, we define our own in terms of
    functions more likely to be available.  */

#define cuberoot(x) (exp(log(x) / 3.0))

/* All Function Prototypes for gravity.c */

void            main    _((int argc, char *argv[]));
static Boolean  funcload _((void));
static char *   alloc   _((unsigned nbytes));
static void     defmassblk _((void));
static struct resbuf *
                entitem _((struct resbuf *rchain, int gcode));
static void     entinfo _((ads_name en,char *h,ads_point p,ads_real *r,int *c));
static void     partext _((void));
static void     addvec  _((ads_real *ap, ads_real *bp, ads_real *cp));
static void     subvec  _((ads_real *ap, ads_real *bp, ads_real *cp));
static ads_real sqabsv  _((ads_real *ap));
static void     sumvec  _((ads_real *ap,ads_real *bp,ads_real t, ads_real *cp));
static void     varblockdef _((void));
static void     savemodes _((void));
static Boolean  boolvar _((char *varname, char *s));
static void     varset  _((void));
static void     setframe _((void));
static Boolean  initacad _((Boolean reset));
static void     addmass _((char *mn,ads_point pos,ads_point vel,ads_real mass));
static void     mass    _((void));
static void     demo    _((void));
static void     massedit _((void));
static void     reset   _((void));
static void     frame   _((void));
static void     run     _((void));
static void     setgrav _((void));
static void     update  _((void));
#ifdef  HIGHC
static void     abort   _((void));
#endif


/*  This  program  works in a somewhat unconventional system of units.
    Length is measured in astronomical units (the mean  distance  from
    the  Earth  to the Sun), mass in units of the mass of the Sun, and
    time in years.  The following definitions derive the value of  the
    gravitational  constant  in this system of units from its handbook
    definition in SI units.  */

#define G_SI        6.6732e-11        /* (Newton Metre^2) / Kilogram^2 */
#define AU          149504094917.0    /* Metres / Astronomical unit */
#define M_SUN       1.989e30          /* Kilograms / Mass of Sun */
#define YEAR        (365.0 * 24 * 60 * 60) /* Seconds / Year */

/*  From Newton's second law, F = ma,

           Newton = kg m / sec^2

    the fundamental units of the gravitational constant are:

           G = N m^2 / kg^2
             = (kg m / sec^2) m^2 / kg^2
             = kg m^3 / sec^2 kg^2
             = m^3 / sec^2 kg

    The conversion factor, therefore,  between  the  SI  gravitational
    constant and its equivalent in our units is:

           K = AU^3 / YEAR^2 M_SUN

*/

#define GRAV_CONV   ((AU * AU * AU) / ((YEAR * YEAR) * M_SUN))

/*  And finally the  gravitational  constant  itself  is  obtained  by
    dividing the SI definition by this conversion factor.  */

#define GRAVCON     (G_SI / GRAV_CONV)

/*  We also want to come up with approximate sizes for the  objects we
    create.   The  actual  sizes  based  on the density of the objects
    result in everything looking like geometrical points  so,  in  the
    rich  tradition  of  celestial  maps, we enormously exaggerate the
    sizes of objects to render them visible.  Our magic number is  one
    tenth of the cube root of the mass of the object.  */

#define DENSCON     0.1

static Boolean functional;            /* C:command is returning result */
static ads_real numsteps = 50;        /* Default number of steps to run */
static double simtime = 0.0;          /* Simulated time */
static long stepno = 0;               /* Step number */
static char fhandle[HANDLEN];         /* Handle of reference frame entity */
static int framei = -1;               /* Reference frame object index */

/*  Command definition and dispatch table.  */

struct {
        char *cmdname;
        void (*cmdfunc)();
} cmdtab[] = {
/*        Name         Function  */
{/*MSG3*/"DEMO",       demo},         /* Create demo case */
{/*MSG4*/"FRAME",      frame},        /* Set local reference frame */
{/*MSG5*/"MASS",       mass},         /* Create new mass */
{/*MSG6*/"MASSEDIT",   massedit},     /* Edit existing mass */
{/*MSG7*/"RESET",      reset},        /* Erase orbital paths */
{/*MSG8*/"RUN",        run},          /* Run simulation */
{/*MSG9*/"SETGRAV",    setgrav},      /* Set mode variables */
{/*MSG10*/"UPDATE",     update}       /* Update masses to last state */
};

/*  Particle structure.  */

typedef struct {
        char partname[32];            /* Particle name */
        ads_point position;           /* Location in space */
        ads_point velocity;           /* Velocity vector */
        ads_point acceleration;       /* Acceleration vector */
        ads_point lastpos;            /* Last plotted position */
        ads_real mass;                /* Mass */
        ads_real radius;              /* Radius */
        int colour;                   /* Entity colour */
        char partblock[HANDLEN];      /* Block defining particle */
} particle;

static particle pt;                   /* Static particle structure */
static particle *ptab = NULL;         /* Particle table */
static int nptab = 0;                 /* Number of in-memory particles */

/*  Particle definition block attributes.  */

#define ParticleBlock  /*MSG11*/"PARTICLE"  /* Particle block name */
struct {                              /* Attribute tag table */
        char *attname;                /* Attribute tag name */
        ads_real *attvar;             /* Variable address */
        char *attprompt;              /* Prompt for attribute */
} partatt[] = {
        {/*MSG12*/"VELOCITY_X", &pt.velocity[X], /*MSG13*/"Velocity (X)"},
        {/*MSG14*/"VELOCITY_Y", &pt.velocity[Y], /*MSG15*/"Velocity (Y)"},
        {/*MSG16*/"VELOCITY_Z", &pt.velocity[Z], /*MSG17*/"Velocity (Z)"},
        {/*MSG18*/"MASS",       &pt.mass,        /*MSG19*/"Mass"}
       };

/*  Modal  variables.   Default   initial   values   below   are   for
    documentation  only.   The actual defaults are reset in initacad()
    upon entry to the drawing editor, so that they will be placed with
    the   default  values  in  the  modal  variable  block  if  it  is
    subsequently created for a new drawing.  */

static Boolean drawonly = True,       /* DRAWONLY: Draw, but don't add entities
                                                   if 1.  Make LINEs if 0. */
               showstep = True,       /* SHOWSTEP: Show step in status line. */
               showtime = True;       /* SHOWTIME: Show time in status line. */
static ads_real stepsize = 0.1,       /* STEPSIZE: Motion step size factor. */
                stepmin = 0.00001;    /* STEPMIN:  Smallest step to use. */

/*  Modal attribute definition.  */

#define ModalBlock  /*MSG20*/"GRAVITY_MODES" /* Mode variable block name */
struct {
        char *attname;                /* Attribute tag name */
        int attvari;                  /* Variable index */
        char *attprompt;              /* Prompt for variable */
} varatt[] = {
        {/*MSG21*/"DRAWONLY", 1, /*MSG22*/"Output to display only? "},
        {/*MSG23*/"SHOWSTEP", 3, /*MSG24*/"Display step number?    "},
        {/*MSG25*/"SHOWTIME", 2, /*MSG26*/"Display time?           "},
        {/*MSG27*/"STEPSIZE", 4, /*MSG28*/"Step size?              "},
        {/*MSG29*/"STEPMIN",  5, /*MSG30*/"Minimum step (years)?   "}
       };

/*-----------------------------------------------------------------------*/
/* MAIN -- the main routine */

void main(argc, argv)
  int argc;
  char *argv[];
{
    int stat, cindex, scode = RSRSLT;
    char errmsg[80];

    ads_init(argc, argv);             /* Initialise the application */

    /* Main dispatch loop. */

    while (True) {

        if ((stat = ads_link(scode)) < 0) {
            sprintf(errmsg,
                    /*MSG31*/"GRAVITY: bad status from ads_link() = %d\n",
                    stat);
#ifdef Macintosh
            macalert(errmsg);
#else
            puts(errmsg);
            fflush(stdout);
#endif /* Macintosh */
            exit(1);
        }

        scode = RSRSLT;               /* Default return code */

        switch (stat) {

        case RQXLOAD:                 /* Load functions.  Called at the start
                                         of the drawing editor.  Re-initialise
                                         the application here. */
            scode = -(funcload() ? RSRSLT : RSERR);
            break;

        case RQXUNLD:                 /* Application unload request. */
            break;

        case RQSUBR:                  /* Evaluate external lisp function */
            cindex = ads_getfuncode();
            functional = False;
            if (!initacad(False)) {
                ads_printf(/*MSG32*/"\nUnable to initialise application.\n");
            } else {

                /* Execute the command from the command table with
                   the index associated with this function. */

                if (cindex > 0) {
                    cindex--;
                    assert(cindex < ELEMENTS(cmdtab));
                    (*cmdtab[cindex].cmdfunc)();
                }
            }
            if (!functional)
                ads_retvoid();        /* Void result */
            break;

        default:
            break;
        }
    }
}

/* FUNCLOAD  --  Load external functions into AutoLISP */

static Boolean funcload()
{
    char ccbuf[40];
    int i;

    strcpy(ccbuf, /*MSG0*/"C:");
    for (i = 0; i < ELEMENTS(cmdtab); i++) {
        strcpy(ccbuf + 2, cmdtab[i].cmdname);
        ads_defun(ccbuf, i + 1);
    }

    return initacad(True);            /* Reset AutoCAD initialisation */
}

/*  ALLOC  --  Allocate storage and fail if it can't be obtained.  */

static char *alloc(nbytes)
  unsigned nbytes;
{
    char *cp;

    if ((cp = malloc(nbytes)) == NULL) {
        ads_abort(/*MSG33*/"Out of memory");
    }
    return cp;
}

/*  DEFMASSBLK  --  Create the mass definition block. */

static void defmassblk()
{
    ads_point centre, ucentre;
    ads_real radius = 1;
    ads_point ax, ax1;
    struct resbuf rb1, rb2, ocolour;
    ads_name e1, e2;
    int i;
    ads_name matbss;
    ads_name ename;
    ads_point atx;

    Spoint(ucentre, 4, 4, 0);
    rb1.restype = rb2.restype = RTSHORT;
    rb1.resval.rint = 1;              /* From UCS */
    rb2.resval.rint = 0;              /* To world */
    ads_trans(ucentre, &rb1, &rb2, False, centre);
    CommandB();
    ads_command(RTSTR, /*MSG0*/"_.UCS",
                RTSTR, /*MSG0*/"_X", RTREAL, 90.0, RTNONE);

    ads_getvar(/*MSG0*/"CECOLOR", &ocolour);
    /* AutoCAD  currently returns  "human readable" colour strings
       like "1 (red)" for the standard colours.  Trim  the  string
       at  the  first space to guarantee we have a valid string to
       restore the colour later.  */
    if (strchr(ocolour.resval.rstring, ' ') != NULL)
        *strchr(ocolour.resval.rstring, ' ') = EOS;

    ads_command(RTSTR, /*MSG0*/"_.COLOUR", RTSTR, /*MSG96*/"_BYBLOCK", RTNONE);
    rb1.resval.rint = 0;              /* From world */
    rb2.resval.rint = 1;              /* To new UCS */
    ads_trans(centre, &rb1, &rb2, False, centre);
    ax[X] = ax1[X] = centre[X];
    ax[Y] = centre[Y] + radius;
    ax[Z] = ax1[Z] = centre[Z];
    ax1[Y] = centre[Y] - radius;
    ads_command(RTSTR, /*MSG0*/"_.LINE", RT3DPOINT, ax, RT3DPOINT, ax1,
                RTSTR, "", RTNONE);
    ads_entlast(e1);
    ads_command(RTSTR, /*MSG0*/"_.Arc", RT3DPOINT, ax, RTSTR, /*MSG0*/"_E",
                RT3DPOINT, ax1, RTSTR, /*MSG0*/"_A", RTSTR, "<<180.0",
                RTNONE);
    ads_entlast(e2);
    PushVar(/*MSG0*/"SURFTAB1", stab1, SphereLong, rint);
    PushVar(/*MSG0*/"SURFTAB2", stab2, SphereLat, rint);
    ads_command(RTSTR, /*MSG0*/"_.REVSURF", RTLB, RTENAME, e2, RT3DPOINT, ax,
                RTLE, RTLB, RTENAME, e1, RT3DPOINT, centre, RTLE, RTSTR, "",
                RTSTR, "", RTNONE);
    PopVar(/*MSG0*/"SURFTAB2", stab2);
    PopVar(/*MSG0*/"SURFTAB1", stab1);
    ads_command(RTSTR, /*MSG0*/"_.COLOUR",
                RTSTR, ocolour.resval.rstring, RTNONE);
    free(ocolour.resval.rstring);
    ads_entdel(e1);
    ads_entdel(e2);
    ads_entlast(e2);

    /* Create attributes */

    Spoint(atx, 4, 2.75, 0);
    ads_ssadd(e2, NULL, matbss);

    PushVar(/*MSG0*/"AFLAGS", saflags, 1, rint);  /* Invisible */
    ads_command(RTSTR, /*MSG0*/"_.ATTDEF", RTSTR, "",
                RTSTR, /*MSG34*/"PARTNAME",
                RTSTR, /*MSG35*/"Particle name", RTSTR, "", RT3DPOINT, atx,
                RTREAL, 0.2, RTREAL, 0.0, RTNONE);
    ads_entlast(ename);
    ads_ssadd(ename, matbss, matbss);

    for (i = 0; i < ELEMENTS(partatt); i++) {
        atx[Y] -= 0.25;
        ads_command(RTSTR, /*MSG0*/"_.ATTDEF",
                    RTSTR, "", RTSTR, partatt[i].attname,
                    RTSTR, partatt[i].attprompt, RTSTR, "0", RT3DPOINT, atx,
                    RTREAL, 0.2, RTREAL, 0.0, RTNONE);
        ads_entlast(ename);
        ads_ssadd(ename, matbss, matbss);
    }
    PopVar(/*MSG0*/"AFLAGS", saflags);

    /* Collect sphere and attributes into a block. */

    ads_command(RTSTR, /*MSG0*/"_.BLOCK",
                RTSTR, ParticleBlock, RT3DPOINT, centre,
                RTPICKS, matbss, RTSTR, "", RTNONE);
    ads_command(RTSTR, /*MSG0*/"_.UCS", RTSTR, /*MSG0*/"_Prev", RTNONE);
    CommandE();
    ads_ssfree(matbss);
}

/*  ENTITEM  --  Search an entity buffer chain and return an item
                 with the specified group code.  */

static struct resbuf *entitem(rchain, gcode)
  struct resbuf *rchain;
  int gcode;
{
    while ((rchain != NULL) && (rchain->restype != gcode))
        rchain = rchain->rbnext;

    return rchain;
}

/*  ENTINFO  --  Obtain information about a particle entity:

                       Handle
                       Position
                       Size (from scale factor of unit block)
                       Colour
*/

static void entinfo(ename, h, p, r, c)
  ads_name ename;
  char *h;
  ads_point p;
  ads_real *r;
  int *c;
{
    struct resbuf *rent, *rh;

    rent = ads_entget(ename);
    if ((rh = entitem(rent, 5)) != NULL)
        strcpy(h, rh->resval.rstring);
    else
        h[0] = EOS;
    rh = entitem(rent, 10);
    assert(rh != NULL);
    Cpoint(p, rh->resval.rpoint);
    rh = entitem(rent, 41);
    assert(rh != NULL);
    *r = rh->resval.rreal;
    if ((rh = entitem(rent, 62)) != NULL) {
        *c = rh->resval.rint;
        if (*c == 0)                  /* Naked colour by block? */
            *c = 7;                   /* Forbidden: make it white.  Q.C.D. */
    } else {
        /* This entity derives its colour from the layer.  Get
           the colour from the layer table. */
        *c = 7;
        if ((rh = entitem(rent, 8)) != NULL) {
            if ((rh = ads_tblsearch(/*MSG0*/"LAYER", rh->resval.rstring,
                                    False)) != NULL) {
                struct resbuf *lh = entitem(rh, 62);
                if (lh != NULL) {
                    int lc = abs(lh->resval.rint);
                    if (lc > 0 && lc < 256) {
                        *c = lc;
                    }
                }
                ads_relrb(rh);
            }
        }
    }
    ads_relrb(rent);
}

/*  PARTEXT  --  Extract particles present in drawing and build the
                 in-memory particle table.  */

static void partext()
{
    long i, l;
    struct resbuf rbet, rbb;
    struct resbuf *rb, *rent, *vb;
    ads_name ename, vname;

    if (ptab != NULL) {
        free(ptab);
    }
    ptab = NULL;
    nptab = 0;

    /* Build the SSGET entity buffer chain to filter for block
       insertions of the named block on the named layer. */

    rbet.restype = 0;                 /* Entity type */
    rbet.resval.rstring = /*MSG0*/"INSERT";
    rbet.rbnext = &rbb;
    rbb.restype = 2;                  /* Block name */
    rbb.resval.rstring = ParticleBlock;
    rbb.rbnext = NULL;

    if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbet, vname) != RTNORM) {
        return;                       /* No definitions in database */
    }

    /* We found one or more definitions.  Process  the  attributes
       that  follow  each, plugging the values into material items
       which get attached to the in-memory definition chain.  */

    if (ads_sslength(vname, &l) < 0)
        l = 0;

    nptab = l;                        /* Save particle count */
    ptab = (particle *) alloc(nptab * sizeof(particle));
    for (i = 0; i < l; i++) {
        ads_name pname;

        ads_ssname(vname, i, ename);
        memcpy(pname, ename, sizeof ename);

        memset(&pt, 0, sizeof pt);
        while (True) {
            ads_entnext(ename, ename);
            rent = ads_entget(ename);
            if (rent == NULL) {
                ads_printf(/*MSG36*/"PARTEXT: Can't read attribute.\n");
                break;
            }
            rb = entitem(rent, 0);    /* Entity type */
            if (rb != NULL) {
                if (strcmp(rb->resval.rstring, /*MSG0*/"ATTRIB") != 0)
                    break;
                rb = entitem(rent, 2);  /* Attribute tag */
                vb = entitem(rent, 1);  /* Attribute value */
                if (rb != NULL && vb != NULL) {
                    if (strcmp(rb->resval.rstring, /*MSG37*/"PARTNAME") == 0) {
                        strcpy(pt.partname, vb->resval.rstring);
                    } else {
                        int j;

                        for (j = 0; j < ELEMENTS(partatt); j++) {
                            if (strcmp(rb->resval.rstring,
                                       partatt[j].attname) == 0) {
                                *partatt[j].attvar = atof(vb->resval.rstring);
                                break;
                            }
                        }
                    }
                }
            }
            ads_relrb(rent);
        }
        entinfo(pname, pt.partblock, pt.position, &pt.radius, &pt.colour);
        memcpy(&(ptab[(int) i]), &pt, sizeof(particle));
    }
    ads_ssfree(vname);                /* Release selection set */
}

/*  ADDVEC  --  Add two vectors, a = b + c  */

static void addvec(ap, bp, cp)
  ads_real *ap, *bp, *cp;
{
    ap[X] = bp[X] + cp[X];
    ap[Y] = bp[Y] + cp[Y];
    ap[Z] = bp[Z] + cp[Z];
}

/*  SUBVEC  --  Subtract two vectors, a = b - c  */

static void subvec(ap, bp, cp)
  ads_real *ap, *bp, *cp;
{
    ap[X] = bp[X] - cp[X];
    ap[Y] = bp[Y] - cp[Y];
    ap[Z] = bp[Z] - cp[Z];
}

/*  SQABSV  --  Square of absolute value of a vector.  */

static ads_real sqabsv(ap)
  ads_real *ap;
{
    return (ap[X] * ap[X] + ap[Y] * ap[Y] + ap[Z] * ap[Z]);
}

/*  SUMVEC  --  Add a linear multiple to another vector, a = b + t * c.  */

static void sumvec(ap, bp, t, cp)
  ads_real *ap, *bp, *cp;
  ads_real t;
{
    ap[X] = bp[X] + t * cp[X];
    ap[Y] = bp[Y] + t * cp[Y];
    ap[Z] = bp[Z] + t * cp[Z];
}

/*  VARBLOCKDEF  --  Create the block that carries the attributes that
                     define the modal variables. */

static void varblockdef()
{
    int i;
    ads_name varbss;
    ads_name ename;
    ads_point atx;

    atx[X] = 4.0;
    atx[Y] = 4.0;
    atx[Z] = 0.0;
    ads_ssadd(NULL, NULL, varbss);

    CommandB();
    PushVar(/*MSG0*/"AFLAGS", saflags, 1, rint);  /* Invisible */

    for (i = 0; i < ELEMENTS(varatt); i++) {
        atx[Y] -= 0.25;
        ads_command(RTSTR, /*MSG0*/"_.ATTDEF",
                    RTSTR, "", RTSTR, varatt[i].attname,
                    RTSTR, varatt[i].attprompt, RTSTR, "0", RT3DPOINT, atx,
                    RTREAL, 0.2, RTREAL, 0.0, RTNONE);
        ads_entlast(ename);
        ads_ssadd(ename, varbss, varbss);
    }
    PopVar(/*MSG0*/"AFLAGS", saflags);

    ads_command(RTSTR, /*MSG0*/"_.BLOCK", RTSTR, ModalBlock, RTSTR, "4,4",
                RTPICKS, varbss, RTSTR, "", RTNONE);
    CommandE();
    ads_ssfree(varbss);
}

/*  SAVEMODES  --  Save application modes as attributes of the mode
                   block.  */

static void savemodes()
{
    int i;
    struct resbuf rbb;
    ads_name ename, vname;
    long l;

    /* Build the SSGET entity buffer chain to filter for block
       insertions of the named block on the named layer. */

    rbb.restype = 2;                  /* Block name */
    rbb.resval.rstring = ModalBlock;
    rbb.rbnext = NULL;

    if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {

        /* Delete all definitions found. */

        if (ads_sslength(vname, &l) < 0)
            l = 0;

        if (l > 0) {
            CommandB();
            ads_command(RTSTR, /*MSG0*/"_.ERASE",
                        RTPICKS, vname, RTSTR, "", RTNONE);
            CommandE();
        }
        ads_ssfree(vname);            /* Release selection set */
    }

    /* Now insert the modal variable block, attaching the
       mode variables to its attributes. */

    CommandB();
    ads_command(RTSTR, /*MSG0*/"_.INSERT", RTSTR, ModalBlock,
                RTSTR, "0,0", RTSTR, "1", RTSTR, "1", RTSTR, "0",
                RTNONE);
    for (i = 0; i < ELEMENTS(varatt); i++) {
        char attval[20];

#define YorN(x) ((x) ? /*MSG38*/"Yes" : /*MSG39*/"No")
        switch (varatt[i].attvari) {
        case 1:
            strcpy(attval, YorN(drawonly));
            break;
        case 2:
            strcpy(attval, YorN(showtime));
            break;
        case 3:
            strcpy(attval, YorN(showstep));
            break;
        case 4:
            sprintf(attval, "%.12g", stepsize);
            break;
        case 5:
            sprintf(attval, "%.12g", stepmin);
            break;
        }
#undef YorN
        ads_command(RTSTR, attval, RTNONE);
    }
    ads_entlast(ename);
    ads_command(RTSTR, /*MSG0*/"_.CHPROP", RTENAME, ename, RTSTR, "",
                RTSTR, /*MSG0*/"_layer", RTSTR, FrozenLayer, RTSTR, "",
                RTNONE);
    CommandE();
}

/*  BOOLVAR  --  Determine Boolean value from attribute string.  We
                 recognise numbers (nonzero means True) and the
                 strings "True", "False", "Yes", and "No",
                 in either upper or lower case.  */

static Boolean boolvar(varname, s)
  char *varname, *s;
{
    char ch = *s;

    if (islower(ch))
        ch = toupper(ch);
    if (isdigit(ch)) {
        return (atoi(s) != 0) ? True : False;
    }

    switch (ch) {
    case /*MSG40*/'Y':
    case /*MSG41*/'T':
        return True;
    case /*MSG42*/'N':
    case /*MSG43*/'F':
        return False;
    default:
        ads_printf(/*MSG44*/"Invalid Boolean value for %s: %s.\n", varname, s);
        break;
    }
    return False;
}

/*  VARSET  --  Set modal variables from the block attributes.  */

static void varset()
{
    struct resbuf rbb;
    ads_name vname;
    long l;

    /* Build the SSGET entity buffer chain to filter for block
       insertions of the named block on the named layer. */

    rbb.restype = 2;                  /* Block name */
    rbb.resval.rstring = ModalBlock;
    rbb.rbnext = NULL;

    if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {

        if (ads_sslength(vname, &l) < 0)
            l = 0;

        if (l > 0) {
            ads_name ename;

            if (ads_ssname(vname, 0L, ename) == RTNORM) {
                while (ads_entnext(ename, ename) == RTNORM) {
                    int i;
                    struct resbuf *rb = ads_entget(ename),
                                  *et, *at, *av;

                    if (rb == NULL)
                        break;
                    et = entitem(rb, 0);
                    assert(et != NULL);
                    if (strcmp(et->resval.rstring, /*MSG0*/"ATTRIB") == 0) {
                        at = entitem(rb, 2);  /* Attribute tag */
                        av = entitem(rb, 1);  /* Attribute value */
                        if (at == NULL || av == NULL)
                            break;
                        for (i = 0; i < ELEMENTS(varatt); i++) {
                            if (strcmp(at->resval.rstring,
                                       varatt[i].attname) == 0) {
                                switch (varatt[i].attvari) {
                                case 1:
                                    drawonly = boolvar(at->resval.rstring,
                                                       av->resval.rstring);
                                    break;
                                case 2:
                                    showtime = boolvar(at->resval.rstring,
                                                       av->resval.rstring);
                                    break;
                                case 3:
                                    showstep = boolvar(at->resval.rstring,
                                                       av->resval.rstring);
                                    break;
                                case 4:
                                    stepsize = 0.1;    /* Default if error */
                                    sscanf(av->resval.rstring, "%lf",
                                           &stepsize);
                                    break;
                                case 5:
                                    stepmin = 0.00001; /* Default if error */
                                    sscanf(av->resval.rstring, "%lf",
                                           &stepmin);
                                    break;
                                }
                            }
                        }
                    } else if (strcmp(et->resval.rstring,
                                      /*MSG0*/"SEQEND") == 0) {
                        ads_relrb(rb);
                        break;
                    }
                    ads_relrb(rb);
                }
            }
        }
        ads_ssfree(vname);            /* Release selection set */
    }
}

/*  SETFRAME  --  Set reference frame index from handle of
                  reference frame entity in effect.  */

static void setframe()
{
    int i;

    framei = -1;
    for (i = 0; i < nptab; i++) {
        if (strcmp(ptab[i].partblock, fhandle) == 0) {
            framei = i;
            break;
        }
    }
}

/*  INITACAD  --  Initialise the required modes in the AutoCAD
                  drawing.  */

static Boolean initacad(reset)
  Boolean reset;
{
    static Boolean initdone, initok;

    if (reset) {
        initdone = False;
        initok = True;
    } else {
        if (!initdone) {
            struct resbuf rb;
            struct resbuf *ep;

            initok = False;

            /* Reset the program modes to standard values upon
               entry to the drawing editor. */

            stepno = 0;               /* Step 0 */
            numsteps = 50;            /* Default number of steps */
            simtime = 0.0;            /* No elapsed time */
            stepsize = 0.1;           /* Default step size */
            stepmin = 0.00001;        /* Default minimum time: none */
            drawonly = True;          /* Draw using ads_grdraw() */
            showstep = True;          /* Show step number */
            showtime = True;          /* Show elapsed time */
            framei = -1;              /* No reference frame object */
            fhandle[0] = EOS;         /* Clear reference frame handle */

            /* Enable  handles  if  they aren't  already on.  We use
               handles to link  from  our  in-memory  table  to  the
               masses  in  the  AutoCAD  database,  and we also keep
               track of the reference frame entity  by  its  handle. */

            ads_getvar(/*MSG0*/"HANDLES", &rb);
            if (!rb.resval.rint) {
                CommandB();
                ads_command(RTSTR, /*MSG0*/"_.handles",
                            RTSTR, /*MSG0*/"_on", RTNONE);
                CommandE();
                ads_getvar(/*MSG0*/"HANDLES", &rb);
                if (!rb.resval.rint) {
                    ads_printf(/*MSG45*/"Cannot enable handles.\n");
                    initdone = True;
                    return (initok = False);
                }
            }

            /* Create required drawing objects. */

            /* If there isn't already a "frozen solid" layer, create one. */

            if ((ep = ads_tblsearch(/*MSG0*/"LAYER",
                                    FrozenLayer, False)) == NULL) {
                CommandB();
                ads_command(RTSTR, /*MSG0*/"_.LAYER",
                            RTSTR, /*MSG0*/"_NEW", RTSTR, FrozenLayer,
                            RTSTR, /*MSG0*/"_FREEZE",
                            RTSTR, FrozenLayer, RTSTR, "", RTNONE);
                ads_command(RTSTR, /*MSG0*/"_.LAYER",
                            RTSTR, /*MSG0*/"_NEW", RTSTR, OrbitLayer,
                            RTSTR, "", RTNONE);
                CommandE();
                defmassblk();         /* Define the particle block */
            } else {
                ads_relrb(ep);
            }

            /* Create the block definition for our modal variables. */

            if ((ep = ads_tblsearch(/*MSG0*/"BLOCK",
                                    ModalBlock, False)) == NULL) {
                varblockdef();
                savemodes();
            } else {
                ads_relrb(ep);
                varset();             /* Load modals from block */
            }

            initdone = initok = True;
        }
    }
    return initok;
}

/*  ADDMASS  --  Insert mass block with attributes.  */

static void addmass(mname, pos, vel, mass)
  char *mname;
  ads_point pos, vel;
  ads_real mass;
{
    ads_real size;
    char velx[32], vely[32], velz[32], smas[32];

    size = cuberoot(mass) * DENSCON;  /* Set size by cube root of mass */
    CommandB();
    sprintf(velx, "%.12g", vel[X]);
    sprintf(vely, "%.12g", vel[Y]);
    sprintf(velz, "%.12g", vel[Z]);
    sprintf(smas, "%.12g", mass);
    ads_command(RTSTR, /*MSG0*/"_.INSERT", RTSTR, ParticleBlock, RTPOINT, pos,
                RTREAL, size, RTREAL, size, RTREAL, 0.0,
                RTSTR, mname,
                RTSTR, velx, RTSTR, vely, RTSTR, velz,
                RTSTR, smas, RTNONE);
    CommandE();
}

/*  MASS  --  Add mass to database.  */

static void mass()
{
    int k;
    ads_point pos, vel;
    ads_real mass;
    char pname[134];

    if (ads_getstring(True, /*MSG46*/"\nParticle name (CR at end): ",
                      pname) != RTNORM)
        return;
    ads_initget(1 + 8 + 16, NULL);
    if (ads_getpoint(NULL, /*MSG47*/"\nPosition: ", pos) != RTNORM)
        return;
    ads_initget(1 + 8 + 16, /*MSG48*/"Vector");
    switch (ads_getpoint(NULL, /*MSG49*/"\nVelocity (or Vector) in AU/Year: ",
                         vel)) {
    case RTNORM:
        break;
    case RTKWORD:                     /* Vector keyword */
        if (ads_getpoint(pos, /*MSG50*/"\nVelocity vector: ", vel) != RTNORM)
            return;
        for (k = X; k <= Z; k++)
            vel[k] -= pos[k];
        break;
    default:
        return;
    }
    if (ads_getreal(/*MSG51*/"\nMass in suns: ", &mass) != RTNORM)
        return;

    addmass(pname, pos, vel, mass);
}

/*  DEMO  --  Load a canned demo world.  */

static void demo()
{
    int i, j;

    typedef struct {                  /* Demo mass table item */
       char *dmname;                  /* Mass name */
       ads_point dmpos;               /* Position */
       ads_point dmvel;               /* Velocity */
       ads_real dmmass;               /* Mass */
       int dmcolour;                  /* Colour */
    } dmass;

    typedef struct {
       char *dname;                   /* Name of demo case */
       ads_point dwind1,              /* Display window corners */
                 dwind2;
       ads_real dnstep;               /* Default number of steps */
       ads_real dssize;               /* Step size */
       ads_real dssmin;               /* Minimum step size */
       char *dframe;                  /* Reference frame */
       dmass *dmtab;                  /* Table of masses */
       int dmtabl;                    /* Number of masses in table */
    } demodesc;

    static dmass interlom[] = {
       {/*MSG52*/"Interloper", {-167.5, 168, 0}, {2,-1.5,0}, 2, 2},
       {/*MSG53*/"Star 3", {100, 50, 0}, {-1, 0, 0}, 20, 3},
       {/*MSG54*/"Star 2", {0, 100, 0}, {2, 0, 0}, 10, 1},
       {/*MSG55*/"Star 1", {0, 0, 0}, {-3, 0, 0}, 10, 5}
      };
    static dmass solarsm[] = {
       {/*MSG56*/"Sun", {0, 0, 0}, {0, 0, 0}, 1, 2},
       {/*MSG57*/"Mercury", {0.387, 0, 0}, {0, 10.0965, 0}, 1.666667e-7, 7},
       {/*MSG58*/"Venus", {0.723, 0, 0}, {0, 7.38628, 0}, 2.44786e-6, 7},
       {/*MSG59*/"Earth", {1, 0, 0}, {0, 6.21318531, 0}, 3.04044e-6, 5},
       {/*MSG60*/"Moon", {1,-0.00256952,0}, {0.215831,6.21318531,0},
        3.6944e-8, 7},
       {/*MSG61*/"Mars", {1.524, 0, 0}, {0, 5.0894, 0}, 3.22716e-7, 1},
       {/*MSG62*/"Jupiter", {5.203, 0, 0}, {0, 2.75456, 0}, 0.000954782, 7},
       {/*MSG63*/"Saturn", {9.539, 0, 0}, {0, 2.03534, 0}, 0.00285837, 2},
       {/*MSG64*/"Uranus", {19.182, 0, 0}, {0, 1.43423, 0}, 4.38596e-5, 4},
       {/*MSG65*/"Neptune", {30.058, 0, 0}, {0, 1.14527, 0}, 5.18135e-5, 4}
      };
    static dmass nemesis = {
       /*MSG66*/"Nemesis", {28.8879, 1.67301, 0}, {-65, 0, 0}, 8, 2
      };
    static demodesc dmt[] = {
       {/*MSG67*/"Interloper", {-176.3, -74, 0}, {176, 180, 0},
        3000, 0.1, 0.00001, NULL,
        interlom, ELEMENTS(interlom)},
       {/*MSG68*/"Solar system", {-3.35, -11.82, 0}, {31.9, 13.6, 0},
        -1, 100, 0, NULL,
        solarsm, ELEMENTS(solarsm)},
       {/*MSG69*/"Nemesis", {-1.59, -1.3, 0}, {2.67, 1.79, 0},
        -1, 100, 0, /*MSG70*/"Sun",
        solarsm, ELEMENTS(solarsm)}
      };

    if (nptab > 0) {
        ads_printf(/*MSG71*/"\n\
Demo can be created only in an empty drawing.\n");
        return;
    }

    partext();
    ads_textscr();                    /* Flip to text screen */
    ads_printf(/*MSG72*/"Available demonstration cases:\n");
    for (i = 0; i < ELEMENTS(dmt); i++) {
        ads_printf("\n%d.  %s", i + 1, dmt[i].dname);
    }
    ads_printf("\n");
    ads_initget(2 + 4, NULL);
    switch (ads_getint(/*MSG73*/"\nWhich demonstration? ", &i)) {
    case RTNORM:
        i--;
        if (i < ELEMENTS(dmt))
            break;
    default:
        return;
    }

    CommandB();
    ads_command(RTSTR, /*MSG0*/"_.ZOOM", RTSTR, /*MSG0*/"_Window",
                RTPOINT, dmt[i].dwind1, RTPOINT, dmt[i].dwind2, RTNONE);
    CommandE();
    for (j = 0; j < dmt[i].dmtabl; j++) {
        CommandB();
        ads_command(RTSTR, /*MSG0*/"_.COLOUR",
                    RTSHORT, dmt[i].dmtab[j].dmcolour, RTNONE);
        CommandE();
        addmass(dmt[i].dmtab[j].dmname, dmt[i].dmtab[j].dmpos,
                dmt[i].dmtab[j].dmvel, dmt[i].dmtab[j].dmmass);
    }

    /* Special case for Nemesis demo.  Load the standard Solar
       System definitions above, then jam the Nemesis mass
       on top of it.  Here comes trouble.... */
    if (i == 2) {
        CommandB();
        ads_command(RTSTR, /*MSG0*/"_.COLOUR",
                    RTSHORT, nemesis.dmcolour, RTNONE);
        CommandE();
        addmass(nemesis.dmname, nemesis.dmpos, nemesis.dmvel,
                nemesis.dmmass);
    }
    /* We now return to our elegant program, already in progress. */

    numsteps = dmt[i].dnstep;
    stepsize = dmt[i].dssize;
    stepmin = dmt[i].dssmin;
    savemodes();
    CommandB();
    ads_command(RTSTR, /*MSG0*/"_.COLOUR", RTSTR, /*MSG127*/"_BYLAYER", RTNONE);
    CommandE();
    partext();

    /* If this demo requests a default reference frame, place it
       in effect. */

    if (dmt[i].dframe != NULL) {
        for (j = 0; j < nptab; j++) {
            if (strcmp(ptab[j].partname, dmt[i].dframe) == 0) {
                strcpy(fhandle, ptab[j].partblock);
                framei = j;
                break;
            }
        }
        assert(j < nptab);
    }
}

/*  MASSEDIT  --  Edit properties of an existing mass.  */

static void massedit()
{
    int i = nptab + 1;
    ads_name ename;
    ads_point ppt;

    partext();                        /* Update in-memory mass table */
    if (ads_entsel(/*MSG74*/"Select mass to edit: ", ename, ppt) == RTNORM) {
        struct resbuf *eb = ads_entget(ename), *ha;

        if ((ha = entitem(eb, 5)) != NULL) {
            for (i = 0; i < nptab; i++) {
                if (strcmp(ptab[i].partblock, ha->resval.rstring) == 0) {
                    break;
                }
            }
        }
    }
    if (i < nptab) {
        if (stepno > 0) {
            reset();
        }
        CommandB();
        ads_command(RTSTR, /*MSG0*/"_.DDATTE", RTENAME, ename, RTNONE);
        CommandE();
        partext();
    } else {
        ads_printf(/*MSG75*/"\nNo mass selected.\n");
    }
}

/*  RESET  --  Reset simulation, erasing all orbital paths.  */

static void reset()
{
    struct resbuf rbb;
    ads_name vname;
    long l;

    /* Build the SSGET entity buffer chain to select all
       entities on the orbital path layer. */

    rbb.restype = 8;                  /* Layer name */
    rbb.resval.rstring = OrbitLayer;
    rbb.rbnext = NULL;

    if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {

        /* We found one or more definitions.  Delete them. */

        if (ads_sslength(vname, &l) < 0)
            l = 0;

        if (l > 0) {
            CommandB();
            ads_command(RTSTR, /*MSG0*/"_.ERASE",
                        RTPICKS, vname, RTSTR, "", RTNONE);
            CommandE();
        }
        ads_ssfree(vname);            /* Release selection set */
    }
    if (drawonly) {
        ads_grtext(-3, NULL, 0);      /* Clear time display */
        ads_redraw(NULL, 0);
    }
    stepno = 0;                       /* Reset step number */
}

/*  FRAME  --  Set reference frame for motion display.  */

static void frame()
{
    ads_name ename;
    ads_point ppt;

    fhandle[0] = EOS;                 /* Assume inertial frame */
    if (ads_entsel(/*MSG76*/"\
Select reference frame mass or RETURN for inertial: ",
                   ename, ppt) == RTNORM) {
        struct resbuf *eb = ads_entget(ename), *ha;

        if ((ha = entitem(eb, 5)) != NULL) {
            strcpy(fhandle, ha->resval.rstring);
        }
    }
    setframe();
    ads_printf(/*MSG77*/"\nReference frame: %s\n",
               framei < 0 ? /*MSG78*/"Inertial" :
                ptab[framei].partname);
}

/*  RUN  --  Run simulation.  This can be called as an interactive
             command, RUN, or functionally as (C:RUN <steps/-time>).  */

static void run()
{
    int i, j, k, n = 0, lastcolour = -1;
    Boolean resume = (stepno > 0) ? True : False, timelimit;
    ads_real r, deltat, endtime = 0;
    char rpt[80];
    struct resbuf clayer, ocolour;
    struct resbuf *rb = ads_getargs();

    /* If an argument was supplied, set the simulation length
       from it rather than asking the user interactively. */

    if (rb != NULL) {
        Boolean argok = True;

        switch (rb->restype) {
        case RTREAL:
            numsteps = rb->resval.rreal;
            break;
        case RTSHORT:
            numsteps = rb->resval.rint;
            break;
        default:
            ads_fail(/*MSG79*/"RUN: Improper argument type");
            argok = False;
            break;
        }
        if (argok) {
            if (numsteps == 0) {
                ads_fail(/*MSG80*/"RUN: Argument must be nonzero");
                argok = False;
            } else {
                if (rb->rbnext != NULL) {
                    ads_fail(/*MSG81*/"RUN: Too many arguments");
                    argok = False;
                }
            }
        }
        if (!argok)
            return;
        functional = True;            /* Mark called as a function */
    } else {

        /* Ask user for length of simulation in years or number of
           integration steps to perform. */

        sprintf(rpt, /*MSG82*/"Steps or (-simulated years) <%.12g>: ",
                numsteps);
        ads_initget(2, NULL);
        switch (ads_getreal(rpt, &r)) {
        case RTCAN:                   /* Control C */
            return;
        case RTNONE:                  /* Null input */
            break;
        case RTNORM:                  /* Number entered */
            numsteps = r;
            break;
        }
    }

    /* If we're starting a new simulation rather than resuming one
       already underway, reset  the  simulated  time  counter  and
       extract the current particle information from the database. */

    if (!resume) {
        simtime = 0.0;
        partext();
    }

    setframe();                       /* Activate reference frame */
    if (numsteps > 0) {
        timelimit = False;
        n = numsteps;
    } else {
        timelimit = True;
        endtime = simtime - numsteps;
    }

    /* Save original layer and colour */

    ads_getvar(/*MSG0*/"CLAYER", &clayer);
    ads_getvar(/*MSG0*/"CECOLOR", &ocolour);
    /* AutoCAD  currently returns  "human readable" colour strings
       like "1 (red)" for the standard colours.  Trim  the  string
       at  the  first space to guarantee we have a valid string to
       restore the colour later.  */
    if (strchr(ocolour.resval.rstring, ' ') != NULL)
        *strchr(ocolour.resval.rstring, ' ') = EOS;

    CommandB();
    ads_command(RTSTR, /*MSG0*/"_.UCS", RTSTR, /*MSG0*/"_World", RTNONE);
    if (!drawonly) {
        ads_command(RTSTR, /*MSG0*/"_.LAYER",
                    RTSTR, /*MSG0*/"_SET", RTSTR, OrbitLayer,
                    RTSTR, "", RTNONE);
    }

    /* Display frame of reference in mode line.  Do this  here  so
       it's  (a) outside the inner loop, and (b) after changing to
       the OrbitLayer (if !drawonly) wipes the mode line.  */

    sprintf(rpt, /*MSG84*/"Reference frame: %s",
            framei < 0 ? /*MSG85*/"Inertial" :
             ptab[framei].partname);
    ads_grtext(-1, rpt, False);

    /* Initialise last plotted position for each particle. */

    for (i = 0; i < nptab; i++) {
        Cpoint(ptab[i].lastpos, ptab[i].position);
    }

    while ((nptab > 1) &&
           (timelimit ? (simtime < endtime) : (n-- > 0)) &&
           !ads_usrbrk()) {
        ads_real mindist = 1E100,
                 maxvel = -1;

        /* Update acceleration for all bodies. */

        for (i = 0; i < nptab; i++) {
            Spoint(ptab[i].acceleration, 0, 0, 0);
            for (j = 0; j < nptab; j++) {
                if (i != j) {
                    ads_real vdist = ads_distance(ptab[i].position,
                                                  ptab[j].position),
                             /* F = G (m1 m2) / r^2 */
                             force = -GRAVCON *
                              ((ptab[i].mass * ptab[j].mass) /
                               (vdist * vdist)),
                             /* F = ma, hence a = F/m */
                             accel = force / ptab[i].mass;

                    mindist = min(mindist, vdist);
                    if (vdist == 0.0) {
                        ads_printf(/*MSG86*/"Collision between %s and %s\n",
                                   ptab[i].partname, ptab[j].partname);
                    } else {
                        /* Update vector components of acceleration. */
                        for (k = X; k <= Z; k++) {
                            ptab[i].acceleration[k] +=
                               accel * (ptab[i].position[k] -
                                        ptab[j].position[k]) / vdist;
                        }
                    }
                }
            }
        }

        /* Update velocity for all bodies. */

        for (i = 0; i < nptab; i++) {
            ads_point pacc;
            ads_real pvel;

            addvec(pacc, ptab[i].velocity, ptab[i].acceleration);
            pvel = sqabsv(pacc);
            maxvel = max(maxvel, pvel);
        }
        maxvel = sqrt(maxvel);

        deltat = stepsize * (mindist / maxvel);
        deltat = max(stepmin, deltat);
        for (i = 0; i < nptab; i++) {
            sumvec(ptab[i].velocity, ptab[i].velocity,
                   deltat, ptab[i].acceleration);
        }

        /* Update position for all bodies. */

        for (i = 0; i < nptab; i++) {
            sumvec(ptab[i].position, ptab[i].position,
                   deltat, ptab[i].velocity);
        }

        /* If we're viewing from one particle's frame of reference,
           translate the view to adjust  for  our  home  particle's
           most recent motion.  */

        if (framei >= 0) {
            ads_point delta;

            subvec(delta, ptab[framei].lastpos, ptab[framei].position);
            for (i = 0; i < nptab; i++) {
                addvec(ptab[i].position, ptab[i].position, delta);
            }
        }

        /* Display motion since last update. */

        for (i = 0; i < nptab; i++) {
            if (drawonly) {
                ads_grdraw(ptab[i].lastpos, ptab[i].position,
                           ptab[i].colour, False);
            } else {
                if (lastcolour != ptab[i].colour) {
                    ads_command(RTSTR, /*MSG0*/"_.COLOUR",
                                RTSHORT, ptab[i].colour,
                                RTNONE);
                    lastcolour = ptab[i].colour;
                }
                ads_command(RTSTR, /*MSG0*/"_.LINE", RTPOINT, ptab[i].lastpos,
                            RTPOINT, ptab[i].position, RTSTR, "", RTNONE);
            }
            Cpoint(ptab[i].lastpos, ptab[i].position);
        }

        /* Update the step number, simulated time, and display
           the values on the status line if selected. */

        stepno++;
        simtime += deltat;
        rpt[0] = EOS;
        if (showtime) {
            sprintf(rpt, /*MSG87*/"Year %.2f", simtime);
        }
        if (showstep) {
            if (rpt[0] != EOS)
                strcat(rpt, "  ");
            sprintf(rpt + strlen(rpt), /*MSG88*/"Step %ld", stepno);
        }
        if (rpt[0] != EOS)
            ads_grtext(-2, rpt, False);
    }

    /* Restore original layer, colour, and UCS. */

    if (!drawonly) {
        ads_command(RTSTR, /*MSG0*/"_.LAYER", RTSTR, /*MSG0*/"_SET",
                    RTSTR, clayer.resval.rstring, RTSTR, "", RTNONE);
        ads_command(RTSTR, /*MSG0*/"_.COLOUR",
                    RTSTR, ocolour.resval.rstring, RTNONE);
    }
    free(clayer.resval.rstring);
    free(ocolour.resval.rstring);
    ads_command(RTSTR, /*MSG0*/"_.UCS", RTSTR, /*MSG0*/"_Prev", RTNONE);
    CommandE();

    /* If called as a function, return the simulation end time.
       If called as a command, print the end time and step count. */

    if (functional)
        ads_retreal(simtime);
    else
        ads_printf(/*MSG89*/"\n\
End at step %ld, %.2f years.\n", stepno, simtime);
}

/*  SETGRAV  --  Set modal variables.  */

static void setgrav()
{
    struct resbuf rbb;
    ads_name vname;
    long l;

    /* Build the SSGET entity buffer chain to filter for block
       insertions of the named block on the named layer. */

    rbb.restype = 2;                  /* Block name */
    rbb.resval.rstring = ModalBlock;
    rbb.rbnext = NULL;

    if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {

        /* Delete all definitions found. */

        if (ads_sslength(vname, &l) < 0)
            l = 0;

        if (l > 0) {
            ads_name ename;

            if (ads_ssname(vname, 0L, ename) == RTNORM) {
                CommandB();
                ads_command(RTSTR, /*MSG0*/"_.DDATTE", RTENAME, ename, RTNONE);
                CommandE();
                varset();
            }
        }
        ads_ssfree(vname);            /* Release selection set */
    }
}

/*  UPDATE  --  Update masses to position them at the end of the current
                point in the simulation.  */

static void update()
{
    int i;

    reset();                          /* Reset the simulation */
    for (i = 0; i < nptab; i++) {
        ads_name ename;

        if (ads_handent(ptab[i].partblock, ename) == RTNORM) {
            struct resbuf *rent, *rb, *vb;

            memcpy(&pt, &(ptab[i]), sizeof(particle));
            rent = ads_entget(ename);
            if (rent != NULL && ((vb = entitem(rent, 10)) != NULL)) {
                Cpoint(vb->resval.rpoint, ptab[i].position);
                if (ads_entmod(rent) != RTNORM) {
                    ads_printf(/*MSG90*/"\
UPDATE: Unable to update position.\n");
                }
            }
            ads_relrb(rent);

            /* Update attributes */

            while (True) {
                ads_entnext(ename, ename);
                rent = ads_entget(ename);
                if (rent == NULL) {
                    ads_printf(/*MSG91*/"UPDATE: Can't read attribute.\n");
                    break;
                }
                rb = entitem(rent, 0);  /* Entity type */
                if (rb != NULL) {
                    if (strcmp(rb->resval.rstring, /*MSG0*/"ATTRIB") != 0)
                        break;
                    rb = entitem(rent, 2);  /* Attribute tag */
                    vb = entitem(rent, 1);  /* Attribute value */
                    if (rb != NULL && vb != NULL) {
                        int j;

                        for (j = 0; j < ELEMENTS(partatt); j++) {
                            if (strcmp(rb->resval.rstring,
                                       partatt[j].attname) == 0) {
                                /* All right, we've found an attribute in
                                   the  table.   Edit  the new value into
                                   it, update the entity,  and  continue. */
                                char *sp = vb->resval.rstring;
                                char ebuf[32];

                                sprintf(ebuf, "%.12g", *partatt[j].attvar);
                                vb->resval.rstring = ebuf;
                                if (ads_entmod(rent) != RTNORM) {
                                    ads_printf(
                                       /*MSG92*/"\
UPDATE:  Cannot update %s attribute.\n",
                                       partatt[j].attname);
                                }
                                /* Restore string buffer */
                                vb->resval.rstring = sp;
                                break;
                            }
                        }
                    }
                }
                ads_relrb(rent);
            }
        } else {
            ads_printf(/*MSG93*/"\nCannot retrieve %s.\n", ptab[i].partname);
        }
    }
}

#ifdef  HIGHC

/*  Earlier versions of High C put abort() in the same module with exit();
    ADS defines its own exit(), so we have to define our own abort(): */

static void abort(void)
{
    ads_abort("");
}

#endif  /* HIGHC */
