/* Next available MSG number is 114 */

/*    

   MAGNETS.CPP

   Copyright (C) 1996 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.
    
   .

      DESCRIPTION:


                  S Q U I R M I N G   M A G N E T S   !!!

        Magnetically perturbed pendulum chaos demonstrator for AutoCAD

        A sample ADS application.

        Designed and implemented by John Walker in October of 1989.

        The  physics  of  the  simulation  is  abstract,  not precise.
        First, I use arbitrary units for the mass and strength of  the
        magnets,  chosen so that values in the range from -10 to 10 do
        reasonable things.  Second,  I  treat  the  magnets  as  point
        sources  of  magnetic  field  lines--as  if they were magnetic
        monopoles (or  the  whole  simulation  were  of  electrostatic
        charges  rather than magnets).  Since real magnets always have
        two poles, and since magnets are of substantial size  compared
        to  the  scale of the simulation, there are substantial dipole
        and quadupole moments to the force exerted on the pendulum  in
        the actual apparatus.  None of these considerations are really
        relevant to any of the phenomena  of  attractors,  orbits,  or
        chaotic  processes,  only to the fidelity of the simulation to
        an actual pendulum.

        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 MAGNET and MAGEDIT
                    commands just like models entered manually.

        MAGNET      Creates  a  new  magnet.   You're invited to enter
                    its   position   (by   any  means  of  co-ordinate
                    specification),  and  its  strength  in  arbitrary
                    units.  If the strength is positive, it repels the
                    pendulum; if positive, it attracts it.

        MAGEDIT     Lets  you  modify  the  strength  of  an  existing
                    magnet.   Pick  a  single  magnet  by pointing.  A
                    dialogue is displayed with the properties  of  the
                    magnet.   Change them as you wish, then pick OK to
                    update  the  properties  of  the  magnet  in   the
                    database.   If  you pick Cancel, the magnet 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
                    seconds.  If you enter  a  positive  number,  it's
                    taken  as  a  step count.  Negative numbers (which
                    may include decimal fractions)  specify  simulated
                    time  in  seconds.  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.

        SETMAG      The  magnets  simulator has several variables that
                    control  its operation.  These variables are saved
                    with the drawing and may be inspected and modified
                    with   the  SETMAG  command.   SETMAG  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
                          magnets.  The  time  step  is  specified  in
                          terms  of seconds.  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   magnetic
                          attraction 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.

                       Dual pendulum offset?
                          To  observe  the  sensitive  dependence   on
                          initial conditions that characterises chaos,
                          you can simulate two pendulums  which  start
                          with  a  small  initial  displacement.   The
                          pendulums do not interact with one  another,
                          but  respond  to  the  magnets placed on the
                          baseplate.  If the offset is zero, a  single
                          pendulum  is  simulated, resulting in faster
                          execution.
*/

#include   <stdio.h>
#include   <stdlib.h>
#include   <string.h>
#include   <ctype.h>
#include   <math.h>
#include   <assert.h>

#include   "rxdefs.h"
#include   "adslib.h"

/*  Standard drawing object names  */

/* Utility frozen layer for information */
#define FrozenLayer /*MSG1*/"FROZEN-SOLID"
#define OrbitLayer  /*MSG2*/"PENDPATH"     /* Layer for pendulum 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 PendLength  20.0              /* Pendulum length */

/* 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 */

/*  Forward functions  */

Boolean        funcload    (void);
Boolean        funcunload  (void);
char          *alloc       (unsigned);
void           defmagblk   (void);
struct resbuf *entitem     (struct resbuf *, int);
void           entinfo     (ads_name, char *, ads_point, ads_real *, int *);
void           partext     (void);
void           sumvec      (ads_real *, ads_real *, ads_real, ads_real *);
void           varblockdef (void);
void           savemodes   (void);
Boolean        boolvar     (char *, char *);
void           varset      (void);
Boolean        initacad    (Boolean);
void           addmag      (char *, ads_point, ads_real);
void           magnet      (void);
void           demo        (void);
void           magedit     (void);
void           reset       (void);
void           run         (void);
void           setmag      (void);

static Boolean functional;            /* C:command is returning result */
static ads_real numsteps = 1000;      /* Default number of steps to run */
static double simtime = 0.0;          /* Simulated time */
static long stepno = 0;               /* Step number */

/*  Command definition and dispatch table.  */

struct {
        char *cmdname;
        void (*cmdfunc)();
} cmdtab[] = {
/*        Name         Function  */
{/*MSG3*/"DEMO",       demo},         /* Create demo case */
{/*MSG4*/"MAGNET",     magnet},       /* Create new magnet */
{/*MSG5*/"MAGEDIT",    magedit},      /* Edit existing magnet */
{/*MSG6*/"RESET",      reset},        /* Erase motion paths */
{/*MSG7*/"RUN",        run},          /* Run simulation */
{/*MSG8*/"SETMAG",     setmag}        /* Set mode variables */
};

/*  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 strength;            /* Magnetic strength */
        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 */

static int ready = 0;                 /* don't run till ready */

/*  Particle definition block attributes.  */

#define ParticleBlock  /*MSG9*/"PARTICLE"  /* Particle block name */
struct {                              /* Attribute tag table */
        char *attname;                /* Attribute tag name */
        ads_real *attvar;             /* Variable address */
        char *attprompt;              /* Prompt for attribute */
} partatt[] = {
        {/*MSG10*/"MAGSTRENGTH", &pt.strength, /*MSG11*/"Magnetic strength"}
       };

/*  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. */
                dualpend = 0.0;       /* DUALPEND: Offset of dual pendulum */
static int ntestp = 1;                /* Number of test particles */

/*  Modal attribute definition.  */

#define ModalBlock  /*MSG12*/"GRAVITY_MODES" /* Mode variable block name */
struct {
        char *attname;                /* Attribute tag name */
        int attvari;                  /* Variable index */
        char *attprompt;              /* Prompt for variable */
} varatt[] = {
        {/*MSG13*/"DRAWONLY", 1, /*MSG14*/"Output to display only? "},
        {/*MSG15*/"SHOWSTEP", 3, /*MSG16*/"Display step number?    "},
        {/*MSG17*/"SHOWTIME", 2, /*MSG18*/"Display time?           "},
        {/*MSG19*/"STEPSIZE", 4, /*MSG20*/"Step size?              "},
        {/*MSG21*/"DUALPEND", 5, /*MSG22*/"Dual pendulum offset?   "}
       };

// Initial entry point for this program
//
extern "C" AcRx::AppRetCode 
acrxEntryPoint(AcRx::AppMsgCode msg, void* ptr)
{
    int cindex;

    switch (msg) {
    case AcRx::kLoadDwgMsg:
        funcload();
        break;
    case AcRx::kInvkSubrMsg:
        cindex = ads_getfuncode();
        if (!initacad(False)) {
            ads_printf(/*MSG24*/"\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)();
            }
        }
        ads_retvoid();        /* Void result */
        break;
    }
    return AcRx::kRetOK;
}

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

static Boolean
funcload()
{
    char ccbuf[40];
    short 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 */
}

/* FUNCUNLOAD  --  Undefine this application's external functions  */

static Boolean
funcunload()
{
    short i;

    /* Undefine each function we defined */

    for (i = 0; i < ELEMENTS(cmdtab); i++) {
        ads_undef(cmdtab[i].cmdname,i);
    }

    return True;
}

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

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

    if ((cp = (char*)acad_malloc(nbytes)) == NULL) {
        ads_printf(/*MSG25*/"Out of memory");
    }
    return cp;
}

/*  DEFMAGBLK  --  Create the magnet definition block. */

static void
defmagblk()
{
    ads_point centre, ucentre;
    ads_real radius = 0.25;
    ads_point ax, ax1;
    struct resbuf rb1, rb2, ocolour;
    ads_name 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_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*/"_.ELEV", RTREAL, -0.5, RTREAL, 0.5, RTNONE);
    Spoint(ax, 0.0, 0.0, -0.5);
    ads_command(RTSTR, /*MSG0*/"_.COLOUR", RTSTR, /*MSG0*/"_BLUE", RTNONE);
#define PcX   3
    ads_command(RTSTR, /*MSG0*/"_.Circle",
                RT3DPOINT, ax, RTREAL, PendLength / PcX,
                RTNONE);
    ads_command(RTSTR, /*MSG0*/"_.ELEV", RTREAL, 0.0, RTREAL, 0.0, RTNONE);
#define PlX   40.0
    Spoint(ax, -PendLength / PlX, -PendLength / PlX, 0.0);
    Spoint(ax1, PendLength / PlX, PendLength / PlX, 0.0);
    ads_command(RTSTR, /*MSG*/"_.COLOUR", RTSTR, /*MSG0*/"_WHITE", RTNONE);
    ads_command(RTSTR, /*MSG*/"_.LINE", RT3DPOINT, ax, RT3DPOINT, ax1,
                RTSTR, "", RTNONE);
    Spoint(ax, -PendLength / PlX, PendLength / PlX, 0.0);
    Spoint(ax1, PendLength / PlX, -PendLength / PlX, 0.0);
    ads_command(RTSTR, /*MSG0*/"_.LINE", RT3DPOINT, ax, RT3DPOINT, ax1,
                RTSTR, "", RTNONE);

    ads_command(RTSTR, /*MSG0*/"_.COLOUR", RTSTR, /*MSG0*/"_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*/"_.ELEV", RTREAL, 0.0, RTREAL, 0.1,
                RTNONE);
    ads_command(RTSTR, /*MSG0*/"_.Circle", RT3DPOINT, ax, RTREAL, radius,
                RTNONE);
    ads_command(RTSTR, /*MSG0*/"_.ELEV", RTREAL, 0.0, RTREAL, 0.0, RTNONE);
    ads_command(RTSTR, /*MSG0*/"_.COLOUR",
                RTSTR, ocolour.resval.rstring, RTNONE);
    acad_free(ocolour.resval.rstring);
    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, /*MSG28*/"PARTNAME",
                RTSTR, /*MSG29*/"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 magnet and attributes into a block. */

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

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

static struct resbuf* 
entitem(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(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) {
        acad_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 magnet count */
    ntestp = (dualpend == 0.0) ? 1 : 2;  /* Number of test particles */
    ptab = (particle *) alloc((nptab + ntestp) * 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(/*MSG30*/"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, /*MSG31*/"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));
    }

    /* Initialise the test particle(s) to their initial state. */

    strcpy(ptab[nptab].partname, /*MSG32*/"Pendulum");
    Spoint(ptab[nptab].position, 4.0, 4.0,
           (PendLength + 1.0) - sqrt(PendLength * PendLength - 4.0 * 4.0));
    Spoint(ptab[nptab].velocity, 0.0, 0.0, 0.0);
    Spoint(ptab[nptab].acceleration, 0.0, 0.0, 0.0);
    Cpoint(ptab[nptab].lastpos, ptab[nptab].position);
    ptab[nptab].strength = 1.0;
    ptab[nptab].radius = 0.25;
    ptab[nptab].colour = 1;           /* Red */

    if (ntestp > 1) {
        double mactemp = (4.0 * 4.0 + (4.0 + dualpend) * (4.0 + dualpend));

        strcpy(ptab[nptab + 1].partname, /*MSG33*/"Displaced pendulum");
        Spoint(ptab[nptab + 1].position, 4.0, 4.0 + dualpend,
               (PendLength + 1.0) - sqrt(PendLength * PendLength - mactemp));
        Spoint(ptab[nptab + 1].velocity, 0.0, 0.0, 0.0);
        Spoint(ptab[nptab + 1].acceleration, 0.0, 0.0, 0.0);
        Cpoint(ptab[nptab + 1].lastpos, ptab[nptab].position);
        ptab[nptab + 1].strength = 1.0;
        ptab[nptab + 1].radius = 0.25;
        ptab[nptab + 1].colour = 3;   /* Green */
    }

    ads_ssfree(vname);                /* Release selection set */
}

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

static void
sumvec(ads_real* ap, ads_real* bp, ads_real t, ads_real* cp)
{
    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) ? /*MSG34*/"Yes" : /*MSG35*/"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", dualpend);
            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(char* varname, char* s)
{
    char ch = *s;

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

    switch (ch) {
    case /*MSG36*/'Y':
    case /*MSG37*/'T':
        return True;
    case /*MSG38*/'N':
    case /*MSG39*/'F':
        return False;
    default:
        ads_printf(/*MSG40*/"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:
                                    dualpend = 0.0;    /* Default if error */
                                    sscanf(av->resval.rstring, "%lf",
                                           &dualpend);
                                    break;
                                }
                            }
                        }
                    } else if (strcmp(et->resval.rstring,
                                      /*MSG0*/"SEQEND") == 0) {
                        ads_relrb(rb);
                        break;
                    }
                    ads_relrb(rb);
                }
            }
        }
        ads_ssfree(vname);            /* Release selection set */
    }
}

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

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

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

            initok = False;

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

            stepno = 0;               /* Step 0 */
            numsteps = 1000;          /* Default number of steps */
            simtime = 0.0;            /* No elapsed time */
            stepsize = 0.1;           /* Default step size */
            dualpend = 0.0;           /* Default 2nd pendulum offset: none */
            drawonly = True;          /* Draw using ads_grdraw() */
            showstep = True;          /* Show step number */
            showtime = True;          /* Show elapsed time */

            /* Enable  handles  if  they aren't  already on.  We use
               handles to link  from  our  in-memory  table  to  the
               magnets 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(/*MSG41*/"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();
                defmagblk();          /* 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 */
            }

            /* Establish default initial view. */

            CommandB();
            Spoint(pzero, 0.0, 0.0, 0.0);
            ads_command(RTSTR, /*MSG0*/"_.ZOOM",
                        RTSTR, /*MSG0*/"_C", RT3DPOINT, pzero,
                        RTREAL, 2.0 * (PendLength / PcX), RTNONE);
            CommandE();

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

/*  ADDMAG  --  Insert magnet block with attributes.  */

static void
addmag(char* mname, ads_point pos, ads_real strength)
{
    ads_real size;
    char smas[32];

    size = 1.0;
    CommandB();
    sprintf(smas, "%.12g", strength);
    ads_command(RTSTR, /*MSG0*/"_.INSERT", RTSTR, ParticleBlock, RTPOINT, pos,
                RTREAL, size, RTREAL, size, RTREAL, 0.0,
                RTSTR, mname,
                RTSTR, smas, RTNONE);
    CommandE();
}

/*  MAGNET  --  Add magnet to database.  */

static void
magnet()
{
    static int magseq = 0;
    ads_point pos;
    ads_real strength;
    char pname[134];

    sprintf(pname, /*MSG42*/"Magnet %d", ++magseq);
    ads_initget(1 + 8 + 16, NULL);
    if (ads_getpoint(NULL, /*MSG43*/"\nPosition: ", pos) != RTNORM)
        return;
    ads_initget(2, NULL);             /* No zero-strength magnets */
    switch (ads_getreal(/*MSG44*/"\nStrength of magnet <1>: ", &strength)) {
    case RTNONE:
        strength = 1.0;
    case RTNORM:
        break;
    default:
        return;
    }
    addmag(pname, pos, strength);
	ready++;
}

/*  DEMO  --  Load canned demo case.  */

static void
demo()
{
    int i, j;

    typedef struct {                  /* Demo magnet table item */
       char *dmname;                  /* Magnet name */
       ads_point dmpos;               /* Position */
       ads_real dmstrength;           /* Magnetic strength */
       int dmcolour;                  /* Colour */
    } dmag;

    typedef struct {
       char *dname;                   /* Name of demo case */
       ads_real dnstep;               /* Default number of steps */
       ads_real dssize;               /* Step size */
       ads_real dualp;                /* Dual pendulum offset */
       dmag *dmtab;                   /* Table of magnets */
       int dmtabl;                    /* Number of magnets in table */
    } demodesc;

    static dmag interlom[] = {
       {/*MSG45*/"Demo magnet 1", {3, 3, 0}, 1, 2},
       {/*MSG46*/"Demo magnet 2", {0, 3, 0}, 2, 3},
       {/*MSG47*/"Demo magnet 3", {-1, -2, 0}, 1, 1},
       {/*MSG48*/"Demo magnet 3", {2, -3, 0}, 1, 5}
      };
    static demodesc dmt[] = {
       {/*MSG49*/"Capt. Magneto",
        3000, 0.1, 0.001,
        interlom, ELEMENTS(interlom)}
      };

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

#ifndef MULTIPLE_DEMOS
    i = 0;                            /* Always use first demo */
#else
    ads_textscr();                    /* Flip to text screen */
    ads_printf(/*MSG51*/"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(/*MSG52*/"\nWhich demonstration? ", &i)) {
    case RTNORM:
        i--;
        if (i < ELEMENTS(dmt))
            break;
    default:
        return;
    }
#endif

    for (j = 0; j < dmt[i].dmtabl; j++) {
        CommandB();
        ads_command(RTSTR, /*MSG0*/"_.COLOUR",
                    RTSHORT, dmt[i].dmtab[j].dmcolour, RTNONE);
        CommandE();
        addmag(dmt[i].dmtab[j].dmname, dmt[i].dmtab[j].dmpos,
               dmt[i].dmtab[j].dmstrength);
    }

    numsteps = dmt[i].dnstep;
    stepsize = dmt[i].dssize;
    dualpend = dmt[i].dualp;
    savemodes();
    CommandB();
    ads_command(RTSTR, /*MSG0*/"_.COLOUR", RTSTR, /*MSG0*/"_BYLAYER", RTNONE);
    CommandE();
    partext();
	ready++;
}

/*  MAGEDIT  --  Edit properties of an existing magnet.  */

static void
magedit()
{
    int i;
    ads_name ename;
    ads_point ppt;

    partext();                        /* Update in-memory magnet table */
    i = nptab; 
    if (ads_entsel(/*MSG53*/"Select magnet 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(/*MSG54*/"\nNo magnet selected.\n");
    }
}

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

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 */
}

/*  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, endtime = 0;
    char rpt[80];
    struct resbuf clayer, ocolour;
    struct resbuf *rb = ads_getargs();

	if (!ready) {
		ads_printf( "\nSet up some magnets first or run Demo.\n");
		return;
	}
		
    /* 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(/*MSG55*/"RUN: Improper argument type");
            argok = False;
            break;
        }
        if (argok) {
            if (numsteps == 0) {
                ads_fail(/*MSG56*/"RUN: Argument must be nonzero");
                argok = False;
            } else {
                if (rb->rbnext != NULL) {
                    ads_fail(/*MSG57*/"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, /*MSG58*/"Steps or (-simulated seconds) <%.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) {
        ads_point bpx;

        Spoint(bpx, 0.0, 0.0, 0.0);
        simtime = 0.0;
        partext();
        ads_initget(1 + 8, NULL);
        if (ads_getpoint(bpx, /*MSG59*/"Pendulum start position: ", 
				ptab[nptab].position) < 0)
		{
            Spoint(ptab[nptab].position, 4.0, 4.0, (PendLength + 1.0) -
                                                   sqrt(PendLength *
                                                        PendLength -
                                                        4.0 * 4.0));
        } else {
            ptab[nptab].position[Z] = (PendLength + 1.0) -
                                      sqrt(PendLength * PendLength -
                                           (ptab[nptab].position[X] *
                                            ptab[nptab].position[X] +
                                            ptab[nptab].position[Y] *
                                            ptab[nptab].position[Y]));
        }
        Cpoint(ptab[nptab].lastpos, ptab[nptab].position);
        if (ntestp > 1) {
            Cpoint(ptab[nptab + 1].position, ptab[nptab].position);
            ptab[nptab + 1].position[Y] += dualpend;
            ptab[nptab + 1].position[Z] = (PendLength + 1.0) -
                                          sqrt(PendLength * PendLength -
                                               (ptab[nptab + 1].position[X] *
                                                ptab[nptab].position[X] +
                                                ptab[nptab + 1].position[Y] *
                                                ptab[nptab].position[Y]));
            Cpoint(ptab[nptab + 1].lastpos, ptab[nptab + 1].position);
        }
    }

    if (numsteps > 0) {
        timelimit = False;
        n = (int)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);
    }

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

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

    while ((timelimit ? (simtime < endtime) : (n-- > 0)) &&
           !ads_usrbrk()) {

        /* Update acceleration for all test bodies. */

        for (i = nptab; i < (nptab + ntestp); i++) {
            ads_real dist, pacce;
            ads_point pzero;

            /* Initialise acceleration to that resulting from the
               restoring force exerted by the current position of
               the pendulum. */

            Spoint(pzero, 0.0, 0.0, 0.0);
            dist = ads_distance(ptab[i].position, pzero);
            pacce = - ptab[i].strength * (dist / PendLength);
            Spoint(ptab[i].acceleration,
                   pacce * (ptab[i].position[X] / dist),
                   pacce * (ptab[i].position[Y] / dist),
                   0.0);
            for (j = 0; j < nptab; j++) {
                if (i != j) {
                    ads_real vdist = ads_distance(ptab[i].position,
                                                  ptab[j].position),
                             /* F = K (m1 m2) / r^2 */
                             force = ((ptab[i].strength * ptab[j].strength) /
                                      (vdist * vdist)),
                             /* F = ma, hence a = F/m.  We take the mass to
                                be 1, so the acceleration is just equal to
                                the force. */
                             accel = force;

                    if (vdist == 0.0) {
                        ads_printf(/*MSG61*/"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 = nptab; i < (nptab + ntestp); i++) {
            sumvec(ptab[i].velocity, ptab[i].velocity,
                   stepsize, ptab[i].acceleration);
        }

        /* Update position for all bodies. */

        for (i = nptab; i < (nptab + ntestp); i++) {
            ads_real pdx;

            sumvec(ptab[i].position, ptab[i].position,
                   stepsize, ptab[i].velocity);
            pdx = ptab[i].position[X] * ptab[i].position[X] +
                  ptab[i].position[Y] * ptab[i].position[Y];
            if (pdx > (PendLength * PendLength)) {
                ads_real patx;

                ads_printf(/*MSG62*/"\
Boink!!  Hit the stop--magnets too strong!\n");
                Spoint(ptab[i].acceleration, 0.0, 0.0, 0.0);
                Spoint(ptab[i].velocity, 0.0, 0.0, 0.0);
                patx = atan2(ptab[i].position[Y], ptab[i].position[X]);
                ptab[i].position[X] = PendLength * cos(patx);
                ptab[i].position[Y] = PendLength * sin(patx);
                pdx = PendLength * PendLength;
            }
            ptab[i].position[Z] = (PendLength + 1.0) -
                                  sqrt(PendLength * PendLength - pdx);
        }

        /* Display motion since last update. */

        for (i = nptab; i < (nptab + ntestp); 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 += stepsize;
        rpt[0] = EOS;
        if (showtime) {
            sprintf(rpt, /*MSG63*/"Time %.2f", simtime);
        }
        if (showstep) {
            if (rpt[0] != EOS)
                strcat(rpt, "  ");
            sprintf(rpt + strlen(rpt), /*MSG64*/"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);
    }
    acad_free(clayer.resval.rstring);
    acad_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(/*MSG65*/"\nEnd at step %ld, %.2f seconds.\n",
                   stepno, simtime);
}

/*  SETMAG  --  Set modal variables.  */

static void
setmag()
{
    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 */
    }
}
