home *** CD-ROM | disk | FTP | other *** search
/ CDUTIL 13 / CDUTIL #13 Julio 1995.iso / windows / acadcom / ads / sample / gravity.c < prev    next >
Encoding:
Text File  |  1995-02-08  |  65.5 KB  |  1,824 lines
  1. /* Next available MSG number is 141 */
  2. /*
  3.  
  4.    GRAVITY.C
  5.  
  6.    Copyright (C) 1989, 1990, 1991, 1992, 1993, 1994 by Autodesk, Inc.
  7.  
  8.    Permission to use, copy, modify, and distribute this software in 
  9.    object code form for any purpose and without fee is hereby granted, 
  10.    provided that the above copyright notice appears in all copies and 
  11.    that both that copyright notice and the limited warranty and 
  12.    restricted rights notice below appear in all supporting 
  13.    documentation.
  14.  
  15.    AUTODESK PROVIDES THIS PROGRAM "AS IS" AND WITH ALL FAULTS.  
  16.    AUTODESK SPECIFICALLY DISCLAIMS ANY IMPLIED WARRANTY OF 
  17.    MERCHANTABILITY OR FITNESS FOR A PARTICULAR USE.  AUTODESK, INC.
  18.    DOES NOT WARRANT THAT THE OPERATION OF THE PROGRAM WILL BE 
  19.    UNINTERRUPTED OR ERROR FREE.
  20.  
  21.    Use, duplication, or disclosure by the U.S. Government is subject to 
  22.    restrictions set forth in FAR 52.227-19 (Commercial Computer 
  23.    Software - Restricted Rights) and DFAR 252.227-7013(c)(1)(ii) 
  24.    (Rights in Technical Data and Computer Software), as applicable.
  25.     
  26.    .
  27.  
  28.         N-Body Gravitational Interaction Simulator for AutoCAD
  29.  
  30.         A sample ADS application.
  31.  
  32.         Designed and implemented by John Walker in August of 1989.
  33.  
  34.            "Ubi materia, ibi geometria."
  35.                       -- Johannes Kepler
  36.  
  37.            "The orbit of  any  one  planet  depends  on  the
  38.            combined  motion  of  all  the  planets,  not  to
  39.            mention the action of all of these on each other.
  40.            But  to  consider simultaneously all these causes
  41.            of motion and to define these  motions  by  exact
  42.            laws   allowing   of   conventional   calculation
  43.            exceeds, unless I am mistaken, the force  of  the
  44.            entire human intellect."
  45.                       -- Sir Isaac Newton, Principia
  46.  
  47.  
  48.         The  physics  underlying  this simulation are explained in the
  49.         the  chapter  "A  Cosmic  Ballet",  pp.   229-238  of  A.   K.
  50.         Dewdney,  "The  Armchair Universe", W.  H.  Freeman: New York,
  51.         1988.
  52.  
  53.         The following commands are implemented in this file:
  54.  
  55.         DEMO        Creates  one  of  a series of standard demo models
  56.                     when executed in  a  new,  empty  drawing.   These
  57.                     models  may be modified with the MASS and MASSEDIT
  58.                     commands just like models entered manually.
  59.  
  60.         FRAME       Asks you to pick a mass entity.   Motion  is  then
  61.                     displayed in that mass's reference frame (in other
  62.                     words, that mass  appears  stationary  and  others
  63.                     move  around  it).  The default reference frame is
  64.                     "Inertial", an unaccelerated frame  at  rest  with
  65.                     respect  to  the  distant  stars.   Specifying  no
  66.                     entity to the  FRAME  command  re-establishes  the
  67.                     inertial frame.
  68.  
  69.         MASS        Creates  a  new  mass.   You're invited to enter a
  70.                     name for the mass, its position (specified by  any
  71.                     means  of  co-ordinate  specification), a velocity
  72.                     vector in units of  astronomical  units  per  year
  73.                     (which  can  be either entered explicitly or drawn
  74.                     by typing "V", then dragging the endpoint  of  the
  75.                     velocity  vector to the correct position), and its
  76.                     mass in units of Solar masses.
  77.  
  78.         MASSEDIT    Lets you modify the name, velocity, and mass of an
  79.                     existing mass.  Pick a single mass by pointing.  A
  80.                     dialogue is displayed with the properties  of  the
  81.                     mass.   Change  them  as you wish, then pick OK to
  82.                     update the properties of the mass in the database.
  83.                     If you pick Cancel, the mass is not modified.
  84.  
  85.         RESET       When  a  simulation  is  run,  it  halts after the
  86.                     specified number of steps or  simulated  time  has
  87.                     elapsed.    A   subsequent  RUN  command  normally
  88.                     resumes the simulation from the point at which the
  89.                     last  stopped.  RESET erases all motion paths from
  90.                     the screen and sets the  simulation  back  to  the
  91.                     start.    RUN   after   a  RESET  will  begin  the
  92.                     simulation from the initial state.
  93.  
  94.         RUN         Starts the simulation.  You're  asked  to  specify
  95.                     the length of the simulation as either a number of
  96.                     motion steps or by the number of simulated  years.
  97.                     If  you  enter  a positive number, it's taken as a
  98.                     step count.  Negative numbers (which  may  include
  99.                     decimal   fractions)  specify  simulated  time  in
  100.                     years.  The length in time of each simulation step
  101.                     varies  based  upon the velocity and separation of
  102.                     masses, so if objects approach  one  another  very
  103.                     closely  a  very  large  number  of  steps will be
  104.                     required to simulate a given time  interval.   The
  105.                     calculation time per step is essentially constant,
  106.                     so if you're investigating a system  with  unknown
  107.                     behaviour  you'll probably want to RUN for a given
  108.                     number of cycles, but when running a stable system
  109.                     such as the Solar System, you can simulate a given
  110.                     time span.  In  any  case,  you  can  terminate  a
  111.                     simulation with the Control C key.  RUN  can  also
  112.                     be  invoked as a function, (C:RUN <length>), where
  113.                     <length> specifies  the  simulation  length  by  a
  114.                     positive  or  negative  number as described above.
  115.                     When called  as  a  function,  C:RUN  returns  the
  116.                     simulation end time as its result.
  117.  
  118.         SETGRAV     The  Gravity  simulator has several variables that
  119.                     control its operation.  These variables are  saved
  120.                     with the drawing and may be inspected and modified
  121.                     with the  SETGRAV  command.   SETGRAV  displays  a
  122.                     dialogue  in  which  you  may  change  any  of the
  123.                     following variables:
  124.  
  125.                        Output to display only?
  126.                           This is set to Yes or  No  (True/False,  and
  127.                           1/0  are  also  accepted).   If  set  to the
  128.                           default value of Yes, motion paths are drawn
  129.                           using  temporary  vectors within the current
  130.                           viewport.  If the picture is  REDRAWn,  they
  131.                           disappear.   If  set to No, motion paths are
  132.                           added to the drawing as LINE entities.  This
  133.                           causes paths to appear in all viewports, and
  134.                           you can  ZOOM  on  sections  of  a  path  to
  135.                           examine   it   in   greater  detail.   Paths
  136.                           represented as  LINEs  are  saved  with  the
  137.                           drawing.  Generating LINE entities for paths
  138.                           is much slower than just drawing them on the
  139.                           screen,  so  choose  this mode only when you
  140.                           need it.
  141.  
  142.                        Display step number?
  143.                           If this mode is  "Yes"  (the  default),  the
  144.                           simulation  step  number is displayed in the
  145.                           coordinate status  line  as  the  simulation
  146.                           progresses.
  147.  
  148.                        Display time?
  149.                           If "Yes" (the default) the simulated time in
  150.                           years is displayed in the coordinate  status
  151.                           line.
  152.  
  153.                        Step size?
  154.                           This real number specifies the  factor  used
  155.                           to  determine  the  size  of the integration
  156.                           steps used in calculating the motion of  the
  157.                           masses.   The  time  step  is  calculated by
  158.                           dividing  the  distance  between   the   two
  159.                           closest   masses  by  the  highest  relative
  160.                           velocity  of  any  pair  of  masses.    This
  161.                           quantity,  measured  in years, is multiplied
  162.                           by the factor  given  by  this  variable  to
  163.                           obtain  the length of the step.  The default
  164.                           value  of  0.1  works  well  for  reasonably
  165.                           well-behaved  simulations.  If you find that
  166.                           a simulation is taking  too  long,  you  can
  167.                           speed it up by increasing the step size, but
  168.                           be aware that the results you see may not be
  169.                           physically  accurate.   Increasing  the step
  170.                           size magnifies the inaccuracies of  modeling
  171.                           a  continuous process such as gravitation by
  172.                           discrete  steps.   In   general,   you   can
  173.                           increase the step size as long as you obtain
  174.                           the same results.  When the  outcome  begins
  175.                           to vary, you've set the step size too large.
  176.  
  177.                        Minimum step?
  178.                           After  the  step  size  is   calculated   as
  179.                           described  above  it  is  compared  with the
  180.                           minimum step size and, if less, the  minimum
  181.                           is  used.   The  minimum step size is set to
  182.                           0.00001 years (about  5  minutes);  this  is
  183.                           sufficient  to unstick many simulations that
  184.                           involve a close encounter.  Amazingly,  even
  185.                           a  minimum  this  short  can  yield  grossly
  186.                           inaccurate results when solar masses execute
  187.                           hairpin turns about one another
  188.  
  189.         UPDATE      A  simulation  normally  proceeds from the initial
  190.                     positions, velocities, and masses of  the  objects
  191.                     specified  by  the  MASS  command.  The simulation
  192.                     keeps track of these values as it  progresses  but
  193.                     does  not automatically adjust the mass objects in
  194.                     the database.  If you want to  move  the  database
  195.                     objects  to  the  positions at the end of the most
  196.                     recent simulation (and adjust their velocities  to
  197.                     the  corresponding  instantaneous velocities), use
  198.                     the UPDATE command.  If you don't  do  an  UPDATE,
  199.                     the masses will remain in their original positions
  200.                     even if you  save  the  drawing  after  running  a
  201.                     simulation.
  202.  
  203. */
  204.  
  205. #include   <stdio.h>
  206. #include   <string.h>
  207. #ifndef Macintosh
  208. #include   <ctype.h>
  209. #endif
  210. #include   <math.h>
  211. #include    <assert.h>
  212. #include   "adslib.h"
  213.  
  214. /*  Standard drawing object names  */
  215.  
  216. /* Utility frozen layer for information */
  217. #define FrozenLayer /*MSG1*/"FROZEN-SOLID"
  218. #define OrbitLayer  /*MSG2*/"ORBITS"  /* Layer for orbital path traces */
  219.  
  220. /*  Data types  */
  221.  
  222. typedef enum {False = 0, True = 1} Boolean;
  223.  
  224. #define HANDLEN  18                   /* String long enough to hold a handle */
  225.  
  226. /* Definitions  to  wrap  around  submission  of  AutoCAD commands to
  227.    prevent their being echoed.  */
  228.  
  229. #define Cmdecho  False                /* Make True for debug command output */
  230.  
  231. #define CommandB()  { struct resbuf rBc, rBb, rBu, rBh; \
  232.         ads_getvar(/*MSG0*/"CMDECHO", &rBc); \
  233.         ads_getvar(/*MSG0*/"BLIPMODE", &rBb); \
  234.         ads_getvar(/*MSG0*/"HIGHLIGHT", &rBh); \
  235.         rBu.restype = RTSHORT; \
  236.         rBu.resval.rint = (int) Cmdecho; \
  237.         ads_setvar(/*MSG0*/"CMDECHO", &rBu); \
  238.         rBu.resval.rint = (int) False; \
  239.         ads_setvar(/*MSG0*/"BLIPMODE", &rBu); \
  240.         ads_setvar(/*MSG0*/"HIGHLIGHT", &rBu)
  241.  
  242. #define CommandE()  ads_setvar(/*MSG0*/"CMDECHO", &rBc); \
  243.                     ads_setvar(/*MSG0*/"BLIPMODE", &rBb); \
  244.                     ads_setvar(/*MSG0*/"HIGHLIGHT", &rBh); }
  245.  
  246. /*  Definitions  that permit you to push and pop system variables with
  247.     minimal complexity.  These don't work  (and  will  cause  horrible
  248.     crashes  if  used)  with  string  variables,  but since all string
  249.     variables are read-only, they cannot be saved and restored in  any
  250.     case.  */
  251.  
  252. #define PushVar(var, cell, newval, newtype) { struct resbuf cell, cNeW; \
  253.         ads_getvar(var, &cell); cNeW.restype = cell.restype;             \
  254.         cNeW.resval.newtype = newval; ads_setvar(var, &cNeW)
  255.  
  256. #define PopVar(var, cell) ads_setvar(var, &cell); }
  257.  
  258. /* Set point variable from three co-ordinates */
  259.  
  260. #define Spoint(pt, x, y, z)  pt[X] = (x);  pt[Y] = (y);  pt[Z] = (z)
  261.  
  262. /* Copy point  */
  263.  
  264. #define Cpoint(d, s)   d[X] = s[X];  d[Y] = s[Y];  d[Z] = s[Z]
  265.  
  266. /* Generation parameters for objects */
  267.  
  268. #define SphereLong  12                /* Longitudinal tabulations on sphere */
  269. #define SphereLat   12                /* Latitudinal tabulations on sphere */
  270.  
  271. /* Utility definition to get an  array's  element  count  (at  compile
  272.    time).   For  example:
  273.  
  274.        int  arr[] = {1,2,3,4,5};
  275.        ...
  276.        printf("%d", ELEMENTS(arr));
  277.  
  278.    would print a five.  ELEMENTS("abc") can also be used to  tell  how
  279.    many  bytes are in a string constant INCLUDING THE TRAILING NULL. */
  280.  
  281. #define ELEMENTS(array) (sizeof(array)/sizeof((array)[0]))
  282.  
  283. /* Utility definitions */
  284.  
  285. #ifndef abs
  286. #define abs(x)      ((x)<0 ? -(x) : (x))
  287. #endif  /* abs */
  288. #ifdef min
  289. #undef  min
  290. #endif
  291. #define min(a,b)    ((a)<(b) ? (a) : (b))
  292. #ifdef max
  293. #undef  max
  294. #endif
  295. #define max(a,b)    ((a)>(b) ? (a) : (b))
  296.  
  297. /*  Many C implementations may lack a cube root function.  Rather than
  298.     count  on  a system cbrt() function, we define our own in terms of
  299.     functions more likely to be available.  */
  300.  
  301. #define cuberoot(x) (exp(log(x) / 3.0))
  302.  
  303. /* All Function Prototypes for gravity.c */
  304.  
  305. void            main    _((int argc, char *argv[]));
  306. static Boolean  funcload _((void));
  307. static char *   alloc   _((unsigned nbytes));
  308. static void     defmassblk _((void));
  309. static struct resbuf *
  310.                 entitem _((struct resbuf *rchain, int gcode));
  311. static void     entinfo _((ads_name en,char *h,ads_point p,ads_real *r,int *c));
  312. static void     partext _((void));
  313. static void     addvec  _((ads_real *ap, ads_real *bp, ads_real *cp));
  314. static void     subvec  _((ads_real *ap, ads_real *bp, ads_real *cp));
  315. static ads_real sqabsv  _((ads_real *ap));
  316. static void     sumvec  _((ads_real *ap,ads_real *bp,ads_real t, ads_real *cp));
  317. static void     varblockdef _((void));
  318. static void     savemodes _((void));
  319. static Boolean  boolvar _((char *varname, char *s));
  320. static void     varset  _((void));
  321. static void     setframe _((void));
  322. static Boolean  initacad _((Boolean reset));
  323. static void     addmass _((char *mn,ads_point pos,ads_point vel,ads_real mass));
  324. static void     mass    _((void));
  325. static void     demo    _((void));
  326. static void     massedit _((void));
  327. static void     reset   _((void));
  328. static void     frame   _((void));
  329. static void     run     _((void));
  330. static void     setgrav _((void));
  331. static void     update  _((void));
  332. #ifdef  HIGHC
  333. static void     abort   _((void));
  334. #endif
  335.  
  336.  
  337. /*  This  program  works in a somewhat unconventional system of units.
  338.     Length is measured in astronomical units (the mean  distance  from
  339.     the  Earth  to the Sun), mass in units of the mass of the Sun, and
  340.     time in years.  The following definitions derive the value of  the
  341.     gravitational  constant  in this system of units from its handbook
  342.     definition in SI units.  */
  343.  
  344. #define G_SI        6.6732e-11        /* (Newton Metre^2) / Kilogram^2 */
  345. #define AU          149504094917.0    /* Metres / Astronomical unit */
  346. #define M_SUN       1.989e30          /* Kilograms / Mass of Sun */
  347. #define YEAR        (365.0 * 24 * 60 * 60) /* Seconds / Year */
  348.  
  349. /*  From Newton's second law, F = ma,
  350.  
  351.            Newton = kg m / sec^2
  352.  
  353.     the fundamental units of the gravitational constant are:
  354.  
  355.            G = N m^2 / kg^2
  356.              = (kg m / sec^2) m^2 / kg^2
  357.              = kg m^3 / sec^2 kg^2
  358.              = m^3 / sec^2 kg
  359.  
  360.     The conversion factor, therefore,  between  the  SI  gravitational
  361.     constant and its equivalent in our units is:
  362.  
  363.            K = AU^3 / YEAR^2 M_SUN
  364.  
  365. */
  366.  
  367. #define GRAV_CONV   ((AU * AU * AU) / ((YEAR * YEAR) * M_SUN))
  368.  
  369. /*  And finally the  gravitational  constant  itself  is  obtained  by
  370.     dividing the SI definition by this conversion factor.  */
  371.  
  372. #define GRAVCON     (G_SI / GRAV_CONV)
  373.  
  374. /*  We also want to come up with approximate sizes for the  objects we
  375.     create.   The  actual  sizes  based  on the density of the objects
  376.     result in everything looking like geometrical points  so,  in  the
  377.     rich  tradition  of  celestial  maps, we enormously exaggerate the
  378.     sizes of objects to render them visible.  Our magic number is  one
  379.     tenth of the cube root of the mass of the object.  */
  380.  
  381. #define DENSCON     0.1
  382.  
  383. static Boolean functional;            /* C:command is returning result */
  384. static ads_real numsteps = 50;        /* Default number of steps to run */
  385. static double simtime = 0.0;          /* Simulated time */
  386. static long stepno = 0;               /* Step number */
  387. static char fhandle[HANDLEN];         /* Handle of reference frame entity */
  388. static int framei = -1;               /* Reference frame object index */
  389.  
  390. /*  Command definition and dispatch table.  */
  391.  
  392. struct {
  393.         char *cmdname;
  394.         void (*cmdfunc)();
  395. } cmdtab[] = {
  396. /*        Name         Function  */
  397. {/*MSG3*/"DEMO",       demo},         /* Create demo case */
  398. {/*MSG4*/"FRAME",      frame},        /* Set local reference frame */
  399. {/*MSG5*/"MASS",       mass},         /* Create new mass */
  400. {/*MSG6*/"MASSEDIT",   massedit},     /* Edit existing mass */
  401. {/*MSG7*/"RESET",      reset},        /* Erase orbital paths */
  402. {/*MSG8*/"RUN",        run},          /* Run simulation */
  403. {/*MSG9*/"SETGRAV",    setgrav},      /* Set mode variables */
  404. {/*MSG10*/"UPDATE",     update}       /* Update masses to last state */
  405. };
  406.  
  407. /*  Particle structure.  */
  408.  
  409. typedef struct {
  410.         char partname[32];            /* Particle name */
  411.         ads_point position;           /* Location in space */
  412.         ads_point velocity;           /* Velocity vector */
  413.         ads_point acceleration;       /* Acceleration vector */
  414.         ads_point lastpos;            /* Last plotted position */
  415.         ads_real mass;                /* Mass */
  416.         ads_real radius;              /* Radius */
  417.         int colour;                   /* Entity colour */
  418.         char partblock[HANDLEN];      /* Block defining particle */
  419. } particle;
  420.  
  421. static particle pt;                   /* Static particle structure */
  422. static particle *ptab = NULL;         /* Particle table */
  423. static int nptab = 0;                 /* Number of in-memory particles */
  424.  
  425. /*  Particle definition block attributes.  */
  426.  
  427. #define ParticleBlock  /*MSG11*/"PARTICLE"  /* Particle block name */
  428. struct {                              /* Attribute tag table */
  429.         char *attname;                /* Attribute tag name */
  430.         ads_real *attvar;             /* Variable address */
  431.         char *attprompt;              /* Prompt for attribute */
  432. } partatt[] = {
  433.         {/*MSG12*/"VELOCITY_X", &pt.velocity[X], /*MSG13*/"Velocity (X)"},
  434.         {/*MSG14*/"VELOCITY_Y", &pt.velocity[Y], /*MSG15*/"Velocity (Y)"},
  435.         {/*MSG16*/"VELOCITY_Z", &pt.velocity[Z], /*MSG17*/"Velocity (Z)"},
  436.         {/*MSG18*/"MASS",       &pt.mass,        /*MSG19*/"Mass"}
  437.        };
  438.  
  439. /*  Modal  variables.   Default   initial   values   below   are   for
  440.     documentation  only.   The actual defaults are reset in initacad()
  441.     upon entry to the drawing editor, so that they will be placed with
  442.     the   default  values  in  the  modal  variable  block  if  it  is
  443.     subsequently created for a new drawing.  */
  444.  
  445. static Boolean drawonly = True,       /* DRAWONLY: Draw, but don't add entities
  446.                                                    if 1.  Make LINEs if 0. */
  447.                showstep = True,       /* SHOWSTEP: Show step in status line. */
  448.                showtime = True;       /* SHOWTIME: Show time in status line. */
  449. static ads_real stepsize = 0.1,       /* STEPSIZE: Motion step size factor. */
  450.                 stepmin = 0.00001;    /* STEPMIN:  Smallest step to use. */
  451.  
  452. /*  Modal attribute definition.  */
  453.  
  454. #define ModalBlock  /*MSG20*/"GRAVITY_MODES" /* Mode variable block name */
  455. struct {
  456.         char *attname;                /* Attribute tag name */
  457.         int attvari;                  /* Variable index */
  458.         char *attprompt;              /* Prompt for variable */
  459. } varatt[] = {
  460.         {/*MSG21*/"DRAWONLY", 1, /*MSG22*/"Output to display only? "},
  461.         {/*MSG23*/"SHOWSTEP", 3, /*MSG24*/"Display step number?    "},
  462.         {/*MSG25*/"SHOWTIME", 2, /*MSG26*/"Display time?           "},
  463.         {/*MSG27*/"STEPSIZE", 4, /*MSG28*/"Step size?              "},
  464.         {/*MSG29*/"STEPMIN",  5, /*MSG30*/"Minimum step (years)?   "}
  465.        };
  466.  
  467. /*-----------------------------------------------------------------------*/
  468. /* MAIN -- the main routine */
  469.  
  470. void main(argc, argv)
  471.   int argc;
  472.   char *argv[];
  473. {
  474.     int stat, cindex, scode = RSRSLT;
  475.     char errmsg[80];
  476.  
  477.     ads_init(argc, argv);             /* Initialise the application */
  478.  
  479.     /* Main dispatch loop. */
  480.  
  481.     while (True) {
  482.  
  483.         if ((stat = ads_link(scode)) < 0) {
  484.             sprintf(errmsg,
  485.                     /*MSG31*/"GRAVITY: bad status from ads_link() = %d\n",
  486.                     stat);
  487. #ifdef Macintosh
  488.             macalert(errmsg);
  489. #else
  490.             puts(errmsg);
  491.             fflush(stdout);
  492. #endif /* Macintosh */
  493.             exit(1);
  494.         }
  495.  
  496.         scode = RSRSLT;               /* Default return code */
  497.  
  498.         switch (stat) {
  499.  
  500.         case RQXLOAD:                 /* Load functions.  Called at the start
  501.                                          of the drawing editor.  Re-initialise
  502.                                          the application here. */
  503.             scode = -(funcload() ? RSRSLT : RSERR);
  504.             break;
  505.  
  506.         case RQXUNLD:                 /* Application unload request. */
  507.             break;
  508.  
  509.         case RQSUBR:                  /* Evaluate external lisp function */
  510.             cindex = ads_getfuncode();
  511.             functional = False;
  512.             if (!initacad(False)) {
  513.                 ads_printf(/*MSG32*/"\nUnable to initialise application.\n");
  514.             } else {
  515.  
  516.                 /* Execute the command from the command table with
  517.                    the index associated with this function. */
  518.  
  519.                 if (cindex > 0) {
  520.                     cindex--;
  521.                     assert(cindex < ELEMENTS(cmdtab));
  522.                     (*cmdtab[cindex].cmdfunc)();
  523.                 }
  524.             }
  525.             if (!functional)
  526.                 ads_retvoid();        /* Void result */
  527.             break;
  528.  
  529.         default:
  530.             break;
  531.         }
  532.     }
  533. }
  534.  
  535. /* FUNCLOAD  --  Load external functions into AutoLISP */
  536.  
  537. static Boolean funcload()
  538. {
  539.     char ccbuf[40];
  540.     int i;
  541.  
  542.     strcpy(ccbuf, /*MSG0*/"C:");
  543.     for (i = 0; i < ELEMENTS(cmdtab); i++) {
  544.         strcpy(ccbuf + 2, cmdtab[i].cmdname);
  545.         ads_defun(ccbuf, i + 1);
  546.     }
  547.  
  548.     return initacad(True);            /* Reset AutoCAD initialisation */
  549. }
  550.  
  551. /*  ALLOC  --  Allocate storage and fail if it can't be obtained.  */
  552.  
  553. static char *alloc(nbytes)
  554.   unsigned nbytes;
  555. {
  556.     char *cp;
  557.  
  558.     if ((cp = malloc(nbytes)) == NULL) {
  559.         ads_abort(/*MSG33*/"Out of memory");
  560.     }
  561.     return cp;
  562. }
  563.  
  564. /*  DEFMASSBLK  --  Create the mass definition block. */
  565.  
  566. static void defmassblk()
  567. {
  568.     ads_point centre, ucentre;
  569.     ads_real radius = 1;
  570.     ads_point ax, ax1;
  571.     struct resbuf rb1, rb2, ocolour;
  572.     ads_name e1, e2;
  573.     int i;
  574.     ads_name matbss;
  575.     ads_name ename;
  576.     ads_point atx;
  577.  
  578.     Spoint(ucentre, 4, 4, 0);
  579.     rb1.restype = rb2.restype = RTSHORT;
  580.     rb1.resval.rint = 1;              /* From UCS */
  581.     rb2.resval.rint = 0;              /* To world */
  582.     ads_trans(ucentre, &rb1, &rb2, False, centre);
  583.     CommandB();
  584.     ads_command(RTSTR, /*MSG0*/"_.UCS",
  585.                 RTSTR, /*MSG0*/"_X", RTREAL, 90.0, RTNONE);
  586.  
  587.     ads_getvar(/*MSG0*/"CECOLOR", &ocolour);
  588.     /* AutoCAD  currently returns  "human readable" colour strings
  589.        like "1 (red)" for the standard colours.  Trim  the  string
  590.        at  the  first space to guarantee we have a valid string to
  591.        restore the colour later.  */
  592.     if (strchr(ocolour.resval.rstring, ' ') != NULL)
  593.         *strchr(ocolour.resval.rstring, ' ') = EOS;
  594.  
  595.     ads_command(RTSTR, /*MSG0*/"_.COLOUR", RTSTR, /*MSG96*/"_BYBLOCK", RTNONE);
  596.     rb1.resval.rint = 0;              /* From world */
  597.     rb2.resval.rint = 1;              /* To new UCS */
  598.     ads_trans(centre, &rb1, &rb2, False, centre);
  599.     ax[X] = ax1[X] = centre[X];
  600.     ax[Y] = centre[Y] + radius;
  601.     ax[Z] = ax1[Z] = centre[Z];
  602.     ax1[Y] = centre[Y] - radius;
  603.     ads_command(RTSTR, /*MSG0*/"_.LINE", RT3DPOINT, ax, RT3DPOINT, ax1,
  604.                 RTSTR, "", RTNONE);
  605.     ads_entlast(e1);
  606.     ads_command(RTSTR, /*MSG0*/"_.Arc", RT3DPOINT, ax, RTSTR, /*MSG0*/"_E",
  607.                 RT3DPOINT, ax1, RTSTR, /*MSG0*/"_A", RTSTR, "<<180.0",
  608.                 RTNONE);
  609.     ads_entlast(e2);
  610.     PushVar(/*MSG0*/"SURFTAB1", stab1, SphereLong, rint);
  611.     PushVar(/*MSG0*/"SURFTAB2", stab2, SphereLat, rint);
  612.     ads_command(RTSTR, /*MSG0*/"_.REVSURF", RTLB, RTENAME, e2, RT3DPOINT, ax,
  613.                 RTLE, RTLB, RTENAME, e1, RT3DPOINT, centre, RTLE, RTSTR, "",
  614.                 RTSTR, "", RTNONE);
  615.     PopVar(/*MSG0*/"SURFTAB2", stab2);
  616.     PopVar(/*MSG0*/"SURFTAB1", stab1);
  617.     ads_command(RTSTR, /*MSG0*/"_.COLOUR",
  618.                 RTSTR, ocolour.resval.rstring, RTNONE);
  619.     free(ocolour.resval.rstring);
  620.     ads_entdel(e1);
  621.     ads_entdel(e2);
  622.     ads_entlast(e2);
  623.  
  624.     /* Create attributes */
  625.  
  626.     Spoint(atx, 4, 2.75, 0);
  627.     ads_ssadd(e2, NULL, matbss);
  628.  
  629.     PushVar(/*MSG0*/"AFLAGS", saflags, 1, rint);  /* Invisible */
  630.     ads_command(RTSTR, /*MSG0*/"_.ATTDEF", RTSTR, "",
  631.                 RTSTR, /*MSG34*/"PARTNAME",
  632.                 RTSTR, /*MSG35*/"Particle name", RTSTR, "", RT3DPOINT, atx,
  633.                 RTREAL, 0.2, RTREAL, 0.0, RTNONE);
  634.     ads_entlast(ename);
  635.     ads_ssadd(ename, matbss, matbss);
  636.  
  637.     for (i = 0; i < ELEMENTS(partatt); i++) {
  638.         atx[Y] -= 0.25;
  639.         ads_command(RTSTR, /*MSG0*/"_.ATTDEF",
  640.                     RTSTR, "", RTSTR, partatt[i].attname,
  641.                     RTSTR, partatt[i].attprompt, RTSTR, "0", RT3DPOINT, atx,
  642.                     RTREAL, 0.2, RTREAL, 0.0, RTNONE);
  643.         ads_entlast(ename);
  644.         ads_ssadd(ename, matbss, matbss);
  645.     }
  646.     PopVar(/*MSG0*/"AFLAGS", saflags);
  647.  
  648.     /* Collect sphere and attributes into a block. */
  649.  
  650.     ads_command(RTSTR, /*MSG0*/"_.BLOCK",
  651.                 RTSTR, ParticleBlock, RT3DPOINT, centre,
  652.                 RTPICKS, matbss, RTSTR, "", RTNONE);
  653.     ads_command(RTSTR, /*MSG0*/"_.UCS", RTSTR, /*MSG0*/"_Prev", RTNONE);
  654.     CommandE();
  655.     ads_ssfree(matbss);
  656. }
  657.  
  658. /*  ENTITEM  --  Search an entity buffer chain and return an item
  659.                  with the specified group code.  */
  660.  
  661. static struct resbuf *entitem(rchain, gcode)
  662.   struct resbuf *rchain;
  663.   int gcode;
  664. {
  665.     while ((rchain != NULL) && (rchain->restype != gcode))
  666.         rchain = rchain->rbnext;
  667.  
  668.     return rchain;
  669. }
  670.  
  671. /*  ENTINFO  --  Obtain information about a particle entity:
  672.  
  673.                        Handle
  674.                        Position
  675.                        Size (from scale factor of unit block)
  676.                        Colour
  677. */
  678.  
  679. static void entinfo(ename, h, p, r, c)
  680.   ads_name ename;
  681.   char *h;
  682.   ads_point p;
  683.   ads_real *r;
  684.   int *c;
  685. {
  686.     struct resbuf *rent, *rh;
  687.  
  688.     rent = ads_entget(ename);
  689.     if ((rh = entitem(rent, 5)) != NULL)
  690.         strcpy(h, rh->resval.rstring);
  691.     else
  692.         h[0] = EOS;
  693.     rh = entitem(rent, 10);
  694.     assert(rh != NULL);
  695.     Cpoint(p, rh->resval.rpoint);
  696.     rh = entitem(rent, 41);
  697.     assert(rh != NULL);
  698.     *r = rh->resval.rreal;
  699.     if ((rh = entitem(rent, 62)) != NULL) {
  700.         *c = rh->resval.rint;
  701.         if (*c == 0)                  /* Naked colour by block? */
  702.             *c = 7;                   /* Forbidden: make it white.  Q.C.D. */
  703.     } else {
  704.         /* This entity derives its colour from the layer.  Get
  705.            the colour from the layer table. */
  706.         *c = 7;
  707.         if ((rh = entitem(rent, 8)) != NULL) {
  708.             if ((rh = ads_tblsearch(/*MSG0*/"LAYER", rh->resval.rstring,
  709.                                     False)) != NULL) {
  710.                 struct resbuf *lh = entitem(rh, 62);
  711.                 if (lh != NULL) {
  712.                     int lc = abs(lh->resval.rint);
  713.                     if (lc > 0 && lc < 256) {
  714.                         *c = lc;
  715.                     }
  716.                 }
  717.                 ads_relrb(rh);
  718.             }
  719.         }
  720.     }
  721.     ads_relrb(rent);
  722. }
  723.  
  724. /*  PARTEXT  --  Extract particles present in drawing and build the
  725.                  in-memory particle table.  */
  726.  
  727. static void partext()
  728. {
  729.     long i, l;
  730.     struct resbuf rbet, rbb;
  731.     struct resbuf *rb, *rent, *vb;
  732.     ads_name ename, vname;
  733.  
  734.     if (ptab != NULL) {
  735.         free(ptab);
  736.     }
  737.     ptab = NULL;
  738.     nptab = 0;
  739.  
  740.     /* Build the SSGET entity buffer chain to filter for block
  741.        insertions of the named block on the named layer. */
  742.  
  743.     rbet.restype = 0;                 /* Entity type */
  744.     rbet.resval.rstring = /*MSG0*/"INSERT";
  745.     rbet.rbnext = &rbb;
  746.     rbb.restype = 2;                  /* Block name */
  747.     rbb.resval.rstring = ParticleBlock;
  748.     rbb.rbnext = NULL;
  749.  
  750.     if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbet, vname) != RTNORM) {
  751.         return;                       /* No definitions in database */
  752.     }
  753.  
  754.     /* We found one or more definitions.  Process  the  attributes
  755.        that  follow  each, plugging the values into material items
  756.        which get attached to the in-memory definition chain.  */
  757.  
  758.     if (ads_sslength(vname, &l) < 0)
  759.         l = 0;
  760.  
  761.     nptab = l;                        /* Save particle count */
  762.     ptab = (particle *) alloc(nptab * sizeof(particle));
  763.     for (i = 0; i < l; i++) {
  764.         ads_name pname;
  765.  
  766.         ads_ssname(vname, i, ename);
  767.         memcpy(pname, ename, sizeof ename);
  768.  
  769.         memset(&pt, 0, sizeof pt);
  770.         while (True) {
  771.             ads_entnext(ename, ename);
  772.             rent = ads_entget(ename);
  773.             if (rent == NULL) {
  774.                 ads_printf(/*MSG36*/"PARTEXT: Can't read attribute.\n");
  775.                 break;
  776.             }
  777.             rb = entitem(rent, 0);    /* Entity type */
  778.             if (rb != NULL) {
  779.                 if (strcmp(rb->resval.rstring, /*MSG0*/"ATTRIB") != 0)
  780.                     break;
  781.                 rb = entitem(rent, 2);  /* Attribute tag */
  782.                 vb = entitem(rent, 1);  /* Attribute value */
  783.                 if (rb != NULL && vb != NULL) {
  784.                     if (strcmp(rb->resval.rstring, /*MSG37*/"PARTNAME") == 0) {
  785.                         strcpy(pt.partname, vb->resval.rstring);
  786.                     } else {
  787.                         int j;
  788.  
  789.                         for (j = 0; j < ELEMENTS(partatt); j++) {
  790.                             if (strcmp(rb->resval.rstring,
  791.                                        partatt[j].attname) == 0) {
  792.                                 *partatt[j].attvar = atof(vb->resval.rstring);
  793.                                 break;
  794.                             }
  795.                         }
  796.                     }
  797.                 }
  798.             }
  799.             ads_relrb(rent);
  800.         }
  801.         entinfo(pname, pt.partblock, pt.position, &pt.radius, &pt.colour);
  802.         memcpy(&(ptab[(int) i]), &pt, sizeof(particle));
  803.     }
  804.     ads_ssfree(vname);                /* Release selection set */
  805. }
  806.  
  807. /*  ADDVEC  --  Add two vectors, a = b + c  */
  808.  
  809. static void addvec(ap, bp, cp)
  810.   ads_real *ap, *bp, *cp;
  811. {
  812.     ap[X] = bp[X] + cp[X];
  813.     ap[Y] = bp[Y] + cp[Y];
  814.     ap[Z] = bp[Z] + cp[Z];
  815. }
  816.  
  817. /*  SUBVEC  --  Subtract two vectors, a = b - c  */
  818.  
  819. static void subvec(ap, bp, cp)
  820.   ads_real *ap, *bp, *cp;
  821. {
  822.     ap[X] = bp[X] - cp[X];
  823.     ap[Y] = bp[Y] - cp[Y];
  824.     ap[Z] = bp[Z] - cp[Z];
  825. }
  826.  
  827. /*  SQABSV  --  Square of absolute value of a vector.  */
  828.  
  829. static ads_real sqabsv(ap)
  830.   ads_real *ap;
  831. {
  832.     return (ap[X] * ap[X] + ap[Y] * ap[Y] + ap[Z] * ap[Z]);
  833. }
  834.  
  835. /*  SUMVEC  --  Add a linear multiple to another vector, a = b + t * c.  */
  836.  
  837. static void sumvec(ap, bp, t, cp)
  838.   ads_real *ap, *bp, *cp;
  839.   ads_real t;
  840. {
  841.     ap[X] = bp[X] + t * cp[X];
  842.     ap[Y] = bp[Y] + t * cp[Y];
  843.     ap[Z] = bp[Z] + t * cp[Z];
  844. }
  845.  
  846. /*  VARBLOCKDEF  --  Create the block that carries the attributes that
  847.                      define the modal variables. */
  848.  
  849. static void varblockdef()
  850. {
  851.     int i;
  852.     ads_name varbss;
  853.     ads_name ename;
  854.     ads_point atx;
  855.  
  856.     atx[X] = 4.0;
  857.     atx[Y] = 4.0;
  858.     atx[Z] = 0.0;
  859.     ads_ssadd(NULL, NULL, varbss);
  860.  
  861.     CommandB();
  862.     PushVar(/*MSG0*/"AFLAGS", saflags, 1, rint);  /* Invisible */
  863.  
  864.     for (i = 0; i < ELEMENTS(varatt); i++) {
  865.         atx[Y] -= 0.25;
  866.         ads_command(RTSTR, /*MSG0*/"_.ATTDEF",
  867.                     RTSTR, "", RTSTR, varatt[i].attname,
  868.                     RTSTR, varatt[i].attprompt, RTSTR, "0", RT3DPOINT, atx,
  869.                     RTREAL, 0.2, RTREAL, 0.0, RTNONE);
  870.         ads_entlast(ename);
  871.         ads_ssadd(ename, varbss, varbss);
  872.     }
  873.     PopVar(/*MSG0*/"AFLAGS", saflags);
  874.  
  875.     ads_command(RTSTR, /*MSG0*/"_.BLOCK", RTSTR, ModalBlock, RTSTR, "4,4",
  876.                 RTPICKS, varbss, RTSTR, "", RTNONE);
  877.     CommandE();
  878.     ads_ssfree(varbss);
  879. }
  880.  
  881. /*  SAVEMODES  --  Save application modes as attributes of the mode
  882.                    block.  */
  883.  
  884. static void savemodes()
  885. {
  886.     int i;
  887.     struct resbuf rbb;
  888.     ads_name ename, vname;
  889.     long l;
  890.  
  891.     /* Build the SSGET entity buffer chain to filter for block
  892.        insertions of the named block on the named layer. */
  893.  
  894.     rbb.restype = 2;                  /* Block name */
  895.     rbb.resval.rstring = ModalBlock;
  896.     rbb.rbnext = NULL;
  897.  
  898.     if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {
  899.  
  900.         /* Delete all definitions found. */
  901.  
  902.         if (ads_sslength(vname, &l) < 0)
  903.             l = 0;
  904.  
  905.         if (l > 0) {
  906.             CommandB();
  907.             ads_command(RTSTR, /*MSG0*/"_.ERASE",
  908.                         RTPICKS, vname, RTSTR, "", RTNONE);
  909.             CommandE();
  910.         }
  911.         ads_ssfree(vname);            /* Release selection set */
  912.     }
  913.  
  914.     /* Now insert the modal variable block, attaching the
  915.        mode variables to its attributes. */
  916.  
  917.     CommandB();
  918.     ads_command(RTSTR, /*MSG0*/"_.INSERT", RTSTR, ModalBlock,
  919.                 RTSTR, "0,0", RTSTR, "1", RTSTR, "1", RTSTR, "0",
  920.                 RTNONE);
  921.     for (i = 0; i < ELEMENTS(varatt); i++) {
  922.         char attval[20];
  923.  
  924. #define YorN(x) ((x) ? /*MSG38*/"Yes" : /*MSG39*/"No")
  925.         switch (varatt[i].attvari) {
  926.         case 1:
  927.             strcpy(attval, YorN(drawonly));
  928.             break;
  929.         case 2:
  930.             strcpy(attval, YorN(showtime));
  931.             break;
  932.         case 3:
  933.             strcpy(attval, YorN(showstep));
  934.             break;
  935.         case 4:
  936.             sprintf(attval, "%.12g", stepsize);
  937.             break;
  938.         case 5:
  939.             sprintf(attval, "%.12g", stepmin);
  940.             break;
  941.         }
  942. #undef YorN
  943.         ads_command(RTSTR, attval, RTNONE);
  944.     }
  945.     ads_entlast(ename);
  946.     ads_command(RTSTR, /*MSG0*/"_.CHPROP", RTENAME, ename, RTSTR, "",
  947.                 RTSTR, /*MSG0*/"_layer", RTSTR, FrozenLayer, RTSTR, "",
  948.                 RTNONE);
  949.     CommandE();
  950. }
  951.  
  952. /*  BOOLVAR  --  Determine Boolean value from attribute string.  We
  953.                  recognise numbers (nonzero means True) and the
  954.                  strings "True", "False", "Yes", and "No",
  955.                  in either upper or lower case.  */
  956.  
  957. static Boolean boolvar(varname, s)
  958.   char *varname, *s;
  959. {
  960.     char ch = *s;
  961.  
  962.     if (islower(ch))
  963.         ch = toupper(ch);
  964.     if (isdigit(ch)) {
  965.         return (atoi(s) != 0) ? True : False;
  966.     }
  967.  
  968.     switch (ch) {
  969.     case /*MSG40*/'Y':
  970.     case /*MSG41*/'T':
  971.         return True;
  972.     case /*MSG42*/'N':
  973.     case /*MSG43*/'F':
  974.         return False;
  975.     default:
  976.         ads_printf(/*MSG44*/"Invalid Boolean value for %s: %s.\n", varname, s);
  977.         break;
  978.     }
  979.     return False;
  980. }
  981.  
  982. /*  VARSET  --  Set modal variables from the block attributes.  */
  983.  
  984. static void varset()
  985. {
  986.     struct resbuf rbb;
  987.     ads_name vname;
  988.     long l;
  989.  
  990.     /* Build the SSGET entity buffer chain to filter for block
  991.        insertions of the named block on the named layer. */
  992.  
  993.     rbb.restype = 2;                  /* Block name */
  994.     rbb.resval.rstring = ModalBlock;
  995.     rbb.rbnext = NULL;
  996.  
  997.     if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {
  998.  
  999.         if (ads_sslength(vname, &l) < 0)
  1000.             l = 0;
  1001.  
  1002.         if (l > 0) {
  1003.             ads_name ename;
  1004.  
  1005.             if (ads_ssname(vname, 0L, ename) == RTNORM) {
  1006.                 while (ads_entnext(ename, ename) == RTNORM) {
  1007.                     int i;
  1008.                     struct resbuf *rb = ads_entget(ename),
  1009.                                   *et, *at, *av;
  1010.  
  1011.                     if (rb == NULL)
  1012.                         break;
  1013.                     et = entitem(rb, 0);
  1014.                     assert(et != NULL);
  1015.                     if (strcmp(et->resval.rstring, /*MSG0*/"ATTRIB") == 0) {
  1016.                         at = entitem(rb, 2);  /* Attribute tag */
  1017.                         av = entitem(rb, 1);  /* Attribute value */
  1018.                         if (at == NULL || av == NULL)
  1019.                             break;
  1020.                         for (i = 0; i < ELEMENTS(varatt); i++) {
  1021.                             if (strcmp(at->resval.rstring,
  1022.                                        varatt[i].attname) == 0) {
  1023.                                 switch (varatt[i].attvari) {
  1024.                                 case 1:
  1025.                                     drawonly = boolvar(at->resval.rstring,
  1026.                                                        av->resval.rstring);
  1027.                                     break;
  1028.                                 case 2:
  1029.                                     showtime = boolvar(at->resval.rstring,
  1030.                                                        av->resval.rstring);
  1031.                                     break;
  1032.                                 case 3:
  1033.                                     showstep = boolvar(at->resval.rstring,
  1034.                                                        av->resval.rstring);
  1035.                                     break;
  1036.                                 case 4:
  1037.                                     stepsize = 0.1;    /* Default if error */
  1038.                                     sscanf(av->resval.rstring, "%lf",
  1039.                                            &stepsize);
  1040.                                     break;
  1041.                                 case 5:
  1042.                                     stepmin = 0.00001; /* Default if error */
  1043.                                     sscanf(av->resval.rstring, "%lf",
  1044.                                            &stepmin);
  1045.                                     break;
  1046.                                 }
  1047.                             }
  1048.                         }
  1049.                     } else if (strcmp(et->resval.rstring,
  1050.                                       /*MSG0*/"SEQEND") == 0) {
  1051.                         ads_relrb(rb);
  1052.                         break;
  1053.                     }
  1054.                     ads_relrb(rb);
  1055.                 }
  1056.             }
  1057.         }
  1058.         ads_ssfree(vname);            /* Release selection set */
  1059.     }
  1060. }
  1061.  
  1062. /*  SETFRAME  --  Set reference frame index from handle of
  1063.                   reference frame entity in effect.  */
  1064.  
  1065. static void setframe()
  1066. {
  1067.     int i;
  1068.  
  1069.     framei = -1;
  1070.     for (i = 0; i < nptab; i++) {
  1071.         if (strcmp(ptab[i].partblock, fhandle) == 0) {
  1072.             framei = i;
  1073.             break;
  1074.         }
  1075.     }
  1076. }
  1077.  
  1078. /*  INITACAD  --  Initialise the required modes in the AutoCAD
  1079.                   drawing.  */
  1080.  
  1081. static Boolean initacad(reset)
  1082.   Boolean reset;
  1083. {
  1084.     static Boolean initdone, initok;
  1085.  
  1086.     if (reset) {
  1087.         initdone = False;
  1088.         initok = True;
  1089.     } else {
  1090.         if (!initdone) {
  1091.             struct resbuf rb;
  1092.             struct resbuf *ep;
  1093.  
  1094.             initok = False;
  1095.  
  1096.             /* Reset the program modes to standard values upon
  1097.                entry to the drawing editor. */
  1098.  
  1099.             stepno = 0;               /* Step 0 */
  1100.             numsteps = 50;            /* Default number of steps */
  1101.             simtime = 0.0;            /* No elapsed time */
  1102.             stepsize = 0.1;           /* Default step size */
  1103.             stepmin = 0.00001;        /* Default minimum time: none */
  1104.             drawonly = True;          /* Draw using ads_grdraw() */
  1105.             showstep = True;          /* Show step number */
  1106.             showtime = True;          /* Show elapsed time */
  1107.             framei = -1;              /* No reference frame object */
  1108.             fhandle[0] = EOS;         /* Clear reference frame handle */
  1109.  
  1110.             /* Enable  handles  if  they aren't  already on.  We use
  1111.                handles to link  from  our  in-memory  table  to  the
  1112.                masses  in  the  AutoCAD  database,  and we also keep
  1113.                track of the reference frame entity  by  its  handle. */
  1114.  
  1115.             ads_getvar(/*MSG0*/"HANDLES", &rb);
  1116.             if (!rb.resval.rint) {
  1117.                 CommandB();
  1118.                 ads_command(RTSTR, /*MSG0*/"_.handles",
  1119.                             RTSTR, /*MSG0*/"_on", RTNONE);
  1120.                 CommandE();
  1121.                 ads_getvar(/*MSG0*/"HANDLES", &rb);
  1122.                 if (!rb.resval.rint) {
  1123.                     ads_printf(/*MSG45*/"Cannot enable handles.\n");
  1124.                     initdone = True;
  1125.                     return (initok = False);
  1126.                 }
  1127.             }
  1128.  
  1129.             /* Create required drawing objects. */
  1130.  
  1131.             /* If there isn't already a "frozen solid" layer, create one. */
  1132.  
  1133.             if ((ep = ads_tblsearch(/*MSG0*/"LAYER",
  1134.                                     FrozenLayer, False)) == NULL) {
  1135.                 CommandB();
  1136.                 ads_command(RTSTR, /*MSG0*/"_.LAYER",
  1137.                             RTSTR, /*MSG0*/"_NEW", RTSTR, FrozenLayer,
  1138.                             RTSTR, /*MSG0*/"_FREEZE",
  1139.                             RTSTR, FrozenLayer, RTSTR, "", RTNONE);
  1140.                 ads_command(RTSTR, /*MSG0*/"_.LAYER",
  1141.                             RTSTR, /*MSG0*/"_NEW", RTSTR, OrbitLayer,
  1142.                             RTSTR, "", RTNONE);
  1143.                 CommandE();
  1144.                 defmassblk();         /* Define the particle block */
  1145.             } else {
  1146.                 ads_relrb(ep);
  1147.             }
  1148.  
  1149.             /* Create the block definition for our modal variables. */
  1150.  
  1151.             if ((ep = ads_tblsearch(/*MSG0*/"BLOCK",
  1152.                                     ModalBlock, False)) == NULL) {
  1153.                 varblockdef();
  1154.                 savemodes();
  1155.             } else {
  1156.                 ads_relrb(ep);
  1157.                 varset();             /* Load modals from block */
  1158.             }
  1159.  
  1160.             initdone = initok = True;
  1161.         }
  1162.     }
  1163.     return initok;
  1164. }
  1165.  
  1166. /*  ADDMASS  --  Insert mass block with attributes.  */
  1167.  
  1168. static void addmass(mname, pos, vel, mass)
  1169.   char *mname;
  1170.   ads_point pos, vel;
  1171.   ads_real mass;
  1172. {
  1173.     ads_real size;
  1174.     char velx[32], vely[32], velz[32], smas[32];
  1175.  
  1176.     size = cuberoot(mass) * DENSCON;  /* Set size by cube root of mass */
  1177.     CommandB();
  1178.     sprintf(velx, "%.12g", vel[X]);
  1179.     sprintf(vely, "%.12g", vel[Y]);
  1180.     sprintf(velz, "%.12g", vel[Z]);
  1181.     sprintf(smas, "%.12g", mass);
  1182.     ads_command(RTSTR, /*MSG0*/"_.INSERT", RTSTR, ParticleBlock, RTPOINT, pos,
  1183.                 RTREAL, size, RTREAL, size, RTREAL, 0.0,
  1184.                 RTSTR, mname,
  1185.                 RTSTR, velx, RTSTR, vely, RTSTR, velz,
  1186.                 RTSTR, smas, RTNONE);
  1187.     CommandE();
  1188. }
  1189.  
  1190. /*  MASS  --  Add mass to database.  */
  1191.  
  1192. static void mass()
  1193. {
  1194.     int k;
  1195.     ads_point pos, vel;
  1196.     ads_real mass;
  1197.     char pname[134];
  1198.  
  1199.     if (ads_getstring(True, /*MSG46*/"\nParticle name (CR at end): ",
  1200.                       pname) != RTNORM)
  1201.         return;
  1202.     ads_initget(1 + 8 + 16, NULL);
  1203.     if (ads_getpoint(NULL, /*MSG47*/"\nPosition: ", pos) != RTNORM)
  1204.         return;
  1205.     ads_initget(1 + 8 + 16, /*MSG48*/"Vector");
  1206.     switch (ads_getpoint(NULL, /*MSG49*/"\nVelocity (or Vector) in AU/Year: ",
  1207.                          vel)) {
  1208.     case RTNORM:
  1209.         break;
  1210.     case RTKWORD:                     /* Vector keyword */
  1211.         if (ads_getpoint(pos, /*MSG50*/"\nVelocity vector: ", vel) != RTNORM)
  1212.             return;
  1213.         for (k = X; k <= Z; k++)
  1214.             vel[k] -= pos[k];
  1215.         break;
  1216.     default:
  1217.         return;
  1218.     }
  1219.     if (ads_getreal(/*MSG51*/"\nMass in suns: ", &mass) != RTNORM)
  1220.         return;
  1221.  
  1222.     addmass(pname, pos, vel, mass);
  1223. }
  1224.  
  1225. /*  DEMO  --  Load a canned demo world.  */
  1226.  
  1227. static void demo()
  1228. {
  1229.     int i, j;
  1230.  
  1231.     typedef struct {                  /* Demo mass table item */
  1232.        char *dmname;                  /* Mass name */
  1233.        ads_point dmpos;               /* Position */
  1234.        ads_point dmvel;               /* Velocity */
  1235.        ads_real dmmass;               /* Mass */
  1236.        int dmcolour;                  /* Colour */
  1237.     } dmass;
  1238.  
  1239.     typedef struct {
  1240.        char *dname;                   /* Name of demo case */
  1241.        ads_point dwind1,              /* Display window corners */
  1242.                  dwind2;
  1243.        ads_real dnstep;               /* Default number of steps */
  1244.        ads_real dssize;               /* Step size */
  1245.        ads_real dssmin;               /* Minimum step size */
  1246.        char *dframe;                  /* Reference frame */
  1247.        dmass *dmtab;                  /* Table of masses */
  1248.        int dmtabl;                    /* Number of masses in table */
  1249.     } demodesc;
  1250.  
  1251.     static dmass interlom[] = {
  1252.        {/*MSG52*/"Interloper", {-167.5, 168, 0}, {2,-1.5,0}, 2, 2},
  1253.        {/*MSG53*/"Star 3", {100, 50, 0}, {-1, 0, 0}, 20, 3},
  1254.        {/*MSG54*/"Star 2", {0, 100, 0}, {2, 0, 0}, 10, 1},
  1255.        {/*MSG55*/"Star 1", {0, 0, 0}, {-3, 0, 0}, 10, 5}
  1256.       };
  1257.     static dmass solarsm[] = {
  1258.        {/*MSG56*/"Sun", {0, 0, 0}, {0, 0, 0}, 1, 2},
  1259.        {/*MSG57*/"Mercury", {0.387, 0, 0}, {0, 10.0965, 0}, 1.666667e-7, 7},
  1260.        {/*MSG58*/"Venus", {0.723, 0, 0}, {0, 7.38628, 0}, 2.44786e-6, 7},
  1261.        {/*MSG59*/"Earth", {1, 0, 0}, {0, 6.21318531, 0}, 3.04044e-6, 5},
  1262.        {/*MSG60*/"Moon", {1,-0.00256952,0}, {0.215831,6.21318531,0},
  1263.         3.6944e-8, 7},
  1264.        {/*MSG61*/"Mars", {1.524, 0, 0}, {0, 5.0894, 0}, 3.22716e-7, 1},
  1265.        {/*MSG62*/"Jupiter", {5.203, 0, 0}, {0, 2.75456, 0}, 0.000954782, 7},
  1266.        {/*MSG63*/"Saturn", {9.539, 0, 0}, {0, 2.03534, 0}, 0.00285837, 2},
  1267.        {/*MSG64*/"Uranus", {19.182, 0, 0}, {0, 1.43423, 0}, 4.38596e-5, 4},
  1268.        {/*MSG65*/"Neptune", {30.058, 0, 0}, {0, 1.14527, 0}, 5.18135e-5, 4}
  1269.       };
  1270.     static dmass nemesis = {
  1271.        /*MSG66*/"Nemesis", {28.8879, 1.67301, 0}, {-65, 0, 0}, 8, 2
  1272.       };
  1273.     static demodesc dmt[] = {
  1274.        {/*MSG67*/"Interloper", {-176.3, -74, 0}, {176, 180, 0},
  1275.         3000, 0.1, 0.00001, NULL,
  1276.         interlom, ELEMENTS(interlom)},
  1277.        {/*MSG68*/"Solar system", {-3.35, -11.82, 0}, {31.9, 13.6, 0},
  1278.         -1, 100, 0, NULL,
  1279.         solarsm, ELEMENTS(solarsm)},
  1280.        {/*MSG69*/"Nemesis", {-1.59, -1.3, 0}, {2.67, 1.79, 0},
  1281.         -1, 100, 0, /*MSG70*/"Sun",
  1282.         solarsm, ELEMENTS(solarsm)}
  1283.       };
  1284.  
  1285.     if (nptab > 0) {
  1286.         ads_printf(/*MSG71*/"\n\
  1287. Demo can be created only in an empty drawing.\n");
  1288.         return;
  1289.     }
  1290.  
  1291.     partext();
  1292.     ads_textscr();                    /* Flip to text screen */
  1293.     ads_printf(/*MSG72*/"Available demonstration cases:\n");
  1294.     for (i = 0; i < ELEMENTS(dmt); i++) {
  1295.         ads_printf("\n%d.  %s", i + 1, dmt[i].dname);
  1296.     }
  1297.     ads_printf("\n");
  1298.     ads_initget(2 + 4, NULL);
  1299.     switch (ads_getint(/*MSG73*/"\nWhich demonstration? ", &i)) {
  1300.     case RTNORM:
  1301.         i--;
  1302.         if (i < ELEMENTS(dmt))
  1303.             break;
  1304.     default:
  1305.         return;
  1306.     }
  1307.  
  1308.     CommandB();
  1309.     ads_command(RTSTR, /*MSG0*/"_.ZOOM", RTSTR, /*MSG0*/"_Window",
  1310.                 RTPOINT, dmt[i].dwind1, RTPOINT, dmt[i].dwind2, RTNONE);
  1311.     CommandE();
  1312.     for (j = 0; j < dmt[i].dmtabl; j++) {
  1313.         CommandB();
  1314.         ads_command(RTSTR, /*MSG0*/"_.COLOUR",
  1315.                     RTSHORT, dmt[i].dmtab[j].dmcolour, RTNONE);
  1316.         CommandE();
  1317.         addmass(dmt[i].dmtab[j].dmname, dmt[i].dmtab[j].dmpos,
  1318.                 dmt[i].dmtab[j].dmvel, dmt[i].dmtab[j].dmmass);
  1319.     }
  1320.  
  1321.     /* Special case for Nemesis demo.  Load the standard Solar
  1322.        System definitions above, then jam the Nemesis mass
  1323.        on top of it.  Here comes trouble.... */
  1324.     if (i == 2) {
  1325.         CommandB();
  1326.         ads_command(RTSTR, /*MSG0*/"_.COLOUR",
  1327.                     RTSHORT, nemesis.dmcolour, RTNONE);
  1328.         CommandE();
  1329.         addmass(nemesis.dmname, nemesis.dmpos, nemesis.dmvel,
  1330.                 nemesis.dmmass);
  1331.     }
  1332.     /* We now return to our elegant program, already in progress. */
  1333.  
  1334.     numsteps = dmt[i].dnstep;
  1335.     stepsize = dmt[i].dssize;
  1336.     stepmin = dmt[i].dssmin;
  1337.     savemodes();
  1338.     CommandB();
  1339.     ads_command(RTSTR, /*MSG0*/"_.COLOUR", RTSTR, /*MSG127*/"_BYLAYER", RTNONE);
  1340.     CommandE();
  1341.     partext();
  1342.  
  1343.     /* If this demo requests a default reference frame, place it
  1344.        in effect. */
  1345.  
  1346.     if (dmt[i].dframe != NULL) {
  1347.         for (j = 0; j < nptab; j++) {
  1348.             if (strcmp(ptab[j].partname, dmt[i].dframe) == 0) {
  1349.                 strcpy(fhandle, ptab[j].partblock);
  1350.                 framei = j;
  1351.                 break;
  1352.             }
  1353.         }
  1354.         assert(j < nptab);
  1355.     }
  1356. }
  1357.  
  1358. /*  MASSEDIT  --  Edit properties of an existing mass.  */
  1359.  
  1360. static void massedit()
  1361. {
  1362.     int i = nptab + 1;
  1363.     ads_name ename;
  1364.     ads_point ppt;
  1365.  
  1366.     partext();                        /* Update in-memory mass table */
  1367.     if (ads_entsel(/*MSG74*/"Select mass to edit: ", ename, ppt) == RTNORM) {
  1368.         struct resbuf *eb = ads_entget(ename), *ha;
  1369.  
  1370.         if ((ha = entitem(eb, 5)) != NULL) {
  1371.             for (i = 0; i < nptab; i++) {
  1372.                 if (strcmp(ptab[i].partblock, ha->resval.rstring) == 0) {
  1373.                     break;
  1374.                 }
  1375.             }
  1376.         }
  1377.     }
  1378.     if (i < nptab) {
  1379.         if (stepno > 0) {
  1380.             reset();
  1381.         }
  1382.         CommandB();
  1383.         ads_command(RTSTR, /*MSG0*/"_.DDATTE", RTENAME, ename, RTNONE);
  1384.         CommandE();
  1385.         partext();
  1386.     } else {
  1387.         ads_printf(/*MSG75*/"\nNo mass selected.\n");
  1388.     }
  1389. }
  1390.  
  1391. /*  RESET  --  Reset simulation, erasing all orbital paths.  */
  1392.  
  1393. static void reset()
  1394. {
  1395.     struct resbuf rbb;
  1396.     ads_name vname;
  1397.     long l;
  1398.  
  1399.     /* Build the SSGET entity buffer chain to select all
  1400.        entities on the orbital path layer. */
  1401.  
  1402.     rbb.restype = 8;                  /* Layer name */
  1403.     rbb.resval.rstring = OrbitLayer;
  1404.     rbb.rbnext = NULL;
  1405.  
  1406.     if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {
  1407.  
  1408.         /* We found one or more definitions.  Delete them. */
  1409.  
  1410.         if (ads_sslength(vname, &l) < 0)
  1411.             l = 0;
  1412.  
  1413.         if (l > 0) {
  1414.             CommandB();
  1415.             ads_command(RTSTR, /*MSG0*/"_.ERASE",
  1416.                         RTPICKS, vname, RTSTR, "", RTNONE);
  1417.             CommandE();
  1418.         }
  1419.         ads_ssfree(vname);            /* Release selection set */
  1420.     }
  1421.     if (drawonly) {
  1422.         ads_grtext(-3, NULL, 0);      /* Clear time display */
  1423.         ads_redraw(NULL, 0);
  1424.     }
  1425.     stepno = 0;                       /* Reset step number */
  1426. }
  1427.  
  1428. /*  FRAME  --  Set reference frame for motion display.  */
  1429.  
  1430. static void frame()
  1431. {
  1432.     ads_name ename;
  1433.     ads_point ppt;
  1434.  
  1435.     fhandle[0] = EOS;                 /* Assume inertial frame */
  1436.     if (ads_entsel(/*MSG76*/"\
  1437. Select reference frame mass or RETURN for inertial: ",
  1438.                    ename, ppt) == RTNORM) {
  1439.         struct resbuf *eb = ads_entget(ename), *ha;
  1440.  
  1441.         if ((ha = entitem(eb, 5)) != NULL) {
  1442.             strcpy(fhandle, ha->resval.rstring);
  1443.         }
  1444.     }
  1445.     setframe();
  1446.     ads_printf(/*MSG77*/"\nReference frame: %s\n",
  1447.                framei < 0 ? /*MSG78*/"Inertial" :
  1448.                 ptab[framei].partname);
  1449. }
  1450.  
  1451. /*  RUN  --  Run simulation.  This can be called as an interactive
  1452.              command, RUN, or functionally as (C:RUN <steps/-time>).  */
  1453.  
  1454. static void run()
  1455. {
  1456.     int i, j, k, n = 0, lastcolour = -1;
  1457.     Boolean resume = (stepno > 0) ? True : False, timelimit;
  1458.     ads_real r, deltat, endtime = 0;
  1459.     char rpt[80];
  1460.     struct resbuf clayer, ocolour;
  1461.     struct resbuf *rb = ads_getargs();
  1462.  
  1463.     /* If an argument was supplied, set the simulation length
  1464.        from it rather than asking the user interactively. */
  1465.  
  1466.     if (rb != NULL) {
  1467.         Boolean argok = True;
  1468.  
  1469.         switch (rb->restype) {
  1470.         case RTREAL:
  1471.             numsteps = rb->resval.rreal;
  1472.             break;
  1473.         case RTSHORT:
  1474.             numsteps = rb->resval.rint;
  1475.             break;
  1476.         default:
  1477.             ads_fail(/*MSG79*/"RUN: Improper argument type");
  1478.             argok = False;
  1479.             break;
  1480.         }
  1481.         if (argok) {
  1482.             if (numsteps == 0) {
  1483.                 ads_fail(/*MSG80*/"RUN: Argument must be nonzero");
  1484.                 argok = False;
  1485.             } else {
  1486.                 if (rb->rbnext != NULL) {
  1487.                     ads_fail(/*MSG81*/"RUN: Too many arguments");
  1488.                     argok = False;
  1489.                 }
  1490.             }
  1491.         }
  1492.         if (!argok)
  1493.             return;
  1494.         functional = True;            /* Mark called as a function */
  1495.     } else {
  1496.  
  1497.         /* Ask user for length of simulation in years or number of
  1498.            integration steps to perform. */
  1499.  
  1500.         sprintf(rpt, /*MSG82*/"Steps or (-simulated years) <%.12g>: ",
  1501.                 numsteps);
  1502.         ads_initget(2, NULL);
  1503.         switch (ads_getreal(rpt, &r)) {
  1504.         case RTCAN:                   /* Control C */
  1505.             return;
  1506.         case RTNONE:                  /* Null input */
  1507.             break;
  1508.         case RTNORM:                  /* Number entered */
  1509.             numsteps = r;
  1510.             break;
  1511.         }
  1512.     }
  1513.  
  1514.     /* If we're starting a new simulation rather than resuming one
  1515.        already underway, reset  the  simulated  time  counter  and
  1516.        extract the current particle information from the database. */
  1517.  
  1518.     if (!resume) {
  1519.         simtime = 0.0;
  1520.         partext();
  1521.     }
  1522.  
  1523.     setframe();                       /* Activate reference frame */
  1524.     if (numsteps > 0) {
  1525.         timelimit = False;
  1526.         n = numsteps;
  1527.     } else {
  1528.         timelimit = True;
  1529.         endtime = simtime - numsteps;
  1530.     }
  1531.  
  1532.     /* Save original layer and colour */
  1533.  
  1534.     ads_getvar(/*MSG0*/"CLAYER", &clayer);
  1535.     ads_getvar(/*MSG0*/"CECOLOR", &ocolour);
  1536.     /* AutoCAD  currently returns  "human readable" colour strings
  1537.        like "1 (red)" for the standard colours.  Trim  the  string
  1538.        at  the  first space to guarantee we have a valid string to
  1539.        restore the colour later.  */
  1540.     if (strchr(ocolour.resval.rstring, ' ') != NULL)
  1541.         *strchr(ocolour.resval.rstring, ' ') = EOS;
  1542.  
  1543.     CommandB();
  1544.     ads_command(RTSTR, /*MSG0*/"_.UCS", RTSTR, /*MSG0*/"_World", RTNONE);
  1545.     if (!drawonly) {
  1546.         ads_command(RTSTR, /*MSG0*/"_.LAYER",
  1547.                     RTSTR, /*MSG0*/"_SET", RTSTR, OrbitLayer,
  1548.                     RTSTR, "", RTNONE);
  1549.     }
  1550.  
  1551.     /* Display frame of reference in mode line.  Do this  here  so
  1552.        it's  (a) outside the inner loop, and (b) after changing to
  1553.        the OrbitLayer (if !drawonly) wipes the mode line.  */
  1554.  
  1555.     sprintf(rpt, /*MSG84*/"Reference frame: %s",
  1556.             framei < 0 ? /*MSG85*/"Inertial" :
  1557.              ptab[framei].partname);
  1558.     ads_grtext(-1, rpt, False);
  1559.  
  1560.     /* Initialise last plotted position for each particle. */
  1561.  
  1562.     for (i = 0; i < nptab; i++) {
  1563.         Cpoint(ptab[i].lastpos, ptab[i].position);
  1564.     }
  1565.  
  1566.     while ((nptab > 1) &&
  1567.            (timelimit ? (simtime < endtime) : (n-- > 0)) &&
  1568.            !ads_usrbrk()) {
  1569.         ads_real mindist = 1E100,
  1570.                  maxvel = -1;
  1571.  
  1572.         /* Update acceleration for all bodies. */
  1573.  
  1574.         for (i = 0; i < nptab; i++) {
  1575.             Spoint(ptab[i].acceleration, 0, 0, 0);
  1576.             for (j = 0; j < nptab; j++) {
  1577.                 if (i != j) {
  1578.                     ads_real vdist = ads_distance(ptab[i].position,
  1579.                                                   ptab[j].position),
  1580.                              /* F = G (m1 m2) / r^2 */
  1581.                              force = -GRAVCON *
  1582.                               ((ptab[i].mass * ptab[j].mass) /
  1583.                                (vdist * vdist)),
  1584.                              /* F = ma, hence a = F/m */
  1585.                              accel = force / ptab[i].mass;
  1586.  
  1587.                     mindist = min(mindist, vdist);
  1588.                     if (vdist == 0.0) {
  1589.                         ads_printf(/*MSG86*/"Collision between %s and %s\n",
  1590.                                    ptab[i].partname, ptab[j].partname);
  1591.                     } else {
  1592.                         /* Update vector components of acceleration. */
  1593.                         for (k = X; k <= Z; k++) {
  1594.                             ptab[i].acceleration[k] +=
  1595.                                accel * (ptab[i].position[k] -
  1596.                                         ptab[j].position[k]) / vdist;
  1597.                         }
  1598.                     }
  1599.                 }
  1600.             }
  1601.         }
  1602.  
  1603.         /* Update velocity for all bodies. */
  1604.  
  1605.         for (i = 0; i < nptab; i++) {
  1606.             ads_point pacc;
  1607.             ads_real pvel;
  1608.  
  1609.             addvec(pacc, ptab[i].velocity, ptab[i].acceleration);
  1610.             pvel = sqabsv(pacc);
  1611.             maxvel = max(maxvel, pvel);
  1612.         }
  1613.         maxvel = sqrt(maxvel);
  1614.  
  1615.         deltat = stepsize * (mindist / maxvel);
  1616.         deltat = max(stepmin, deltat);
  1617.         for (i = 0; i < nptab; i++) {
  1618.             sumvec(ptab[i].velocity, ptab[i].velocity,
  1619.                    deltat, ptab[i].acceleration);
  1620.         }
  1621.  
  1622.         /* Update position for all bodies. */
  1623.  
  1624.         for (i = 0; i < nptab; i++) {
  1625.             sumvec(ptab[i].position, ptab[i].position,
  1626.                    deltat, ptab[i].velocity);
  1627.         }
  1628.  
  1629.         /* If we're viewing from one particle's frame of reference,
  1630.            translate the view to adjust  for  our  home  particle's
  1631.            most recent motion.  */
  1632.  
  1633.         if (framei >= 0) {
  1634.             ads_point delta;
  1635.  
  1636.             subvec(delta, ptab[framei].lastpos, ptab[framei].position);
  1637.             for (i = 0; i < nptab; i++) {
  1638.                 addvec(ptab[i].position, ptab[i].position, delta);
  1639.             }
  1640.         }
  1641.  
  1642.         /* Display motion since last update. */
  1643.  
  1644.         for (i = 0; i < nptab; i++) {
  1645.             if (drawonly) {
  1646.                 ads_grdraw(ptab[i].lastpos, ptab[i].position,
  1647.                            ptab[i].colour, False);
  1648.             } else {
  1649.                 if (lastcolour != ptab[i].colour) {
  1650.                     ads_command(RTSTR, /*MSG0*/"_.COLOUR",
  1651.                                 RTSHORT, ptab[i].colour,
  1652.                                 RTNONE);
  1653.                     lastcolour = ptab[i].colour;
  1654.                 }
  1655.                 ads_command(RTSTR, /*MSG0*/"_.LINE", RTPOINT, ptab[i].lastpos,
  1656.                             RTPOINT, ptab[i].position, RTSTR, "", RTNONE);
  1657.             }
  1658.             Cpoint(ptab[i].lastpos, ptab[i].position);
  1659.         }
  1660.  
  1661.         /* Update the step number, simulated time, and display
  1662.            the values on the status line if selected. */
  1663.  
  1664.         stepno++;
  1665.         simtime += deltat;
  1666.         rpt[0] = EOS;
  1667.         if (showtime) {
  1668.             sprintf(rpt, /*MSG87*/"Year %.2f", simtime);
  1669.         }
  1670.         if (showstep) {
  1671.             if (rpt[0] != EOS)
  1672.                 strcat(rpt, "  ");
  1673.             sprintf(rpt + strlen(rpt), /*MSG88*/"Step %ld", stepno);
  1674.         }
  1675.         if (rpt[0] != EOS)
  1676.             ads_grtext(-2, rpt, False);
  1677.     }
  1678.  
  1679.     /* Restore original layer, colour, and UCS. */
  1680.  
  1681.     if (!drawonly) {
  1682.         ads_command(RTSTR, /*MSG0*/"_.LAYER", RTSTR, /*MSG0*/"_SET",
  1683.                     RTSTR, clayer.resval.rstring, RTSTR, "", RTNONE);
  1684.         ads_command(RTSTR, /*MSG0*/"_.COLOUR",
  1685.                     RTSTR, ocolour.resval.rstring, RTNONE);
  1686.     }
  1687.     free(clayer.resval.rstring);
  1688.     free(ocolour.resval.rstring);
  1689.     ads_command(RTSTR, /*MSG0*/"_.UCS", RTSTR, /*MSG0*/"_Prev", RTNONE);
  1690.     CommandE();
  1691.  
  1692.     /* If called as a function, return the simulation end time.
  1693.        If called as a command, print the end time and step count. */
  1694.  
  1695.     if (functional)
  1696.         ads_retreal(simtime);
  1697.     else
  1698.         ads_printf(/*MSG89*/"\n\
  1699. End at step %ld, %.2f years.\n", stepno, simtime);
  1700. }
  1701.  
  1702. /*  SETGRAV  --  Set modal variables.  */
  1703.  
  1704. static void setgrav()
  1705. {
  1706.     struct resbuf rbb;
  1707.     ads_name vname;
  1708.     long l;
  1709.  
  1710.     /* Build the SSGET entity buffer chain to filter for block
  1711.        insertions of the named block on the named layer. */
  1712.  
  1713.     rbb.restype = 2;                  /* Block name */
  1714.     rbb.resval.rstring = ModalBlock;
  1715.     rbb.rbnext = NULL;
  1716.  
  1717.     if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {
  1718.  
  1719.         /* Delete all definitions found. */
  1720.  
  1721.         if (ads_sslength(vname, &l) < 0)
  1722.             l = 0;
  1723.  
  1724.         if (l > 0) {
  1725.             ads_name ename;
  1726.  
  1727.             if (ads_ssname(vname, 0L, ename) == RTNORM) {
  1728.                 CommandB();
  1729.                 ads_command(RTSTR, /*MSG0*/"_.DDATTE", RTENAME, ename, RTNONE);
  1730.                 CommandE();
  1731.                 varset();
  1732.             }
  1733.         }
  1734.         ads_ssfree(vname);            /* Release selection set */
  1735.     }
  1736. }
  1737.  
  1738. /*  UPDATE  --  Update masses to position them at the end of the current
  1739.                 point in the simulation.  */
  1740.  
  1741. static void update()
  1742. {
  1743.     int i;
  1744.  
  1745.     reset();                          /* Reset the simulation */
  1746.     for (i = 0; i < nptab; i++) {
  1747.         ads_name ename;
  1748.  
  1749.         if (ads_handent(ptab[i].partblock, ename) == RTNORM) {
  1750.             struct resbuf *rent, *rb, *vb;
  1751.  
  1752.             memcpy(&pt, &(ptab[i]), sizeof(particle));
  1753.             rent = ads_entget(ename);
  1754.             if (rent != NULL && ((vb = entitem(rent, 10)) != NULL)) {
  1755.                 Cpoint(vb->resval.rpoint, ptab[i].position);
  1756.                 if (ads_entmod(rent) != RTNORM) {
  1757.                     ads_printf(/*MSG90*/"\
  1758. UPDATE: Unable to update position.\n");
  1759.                 }
  1760.             }
  1761.             ads_relrb(rent);
  1762.  
  1763.             /* Update attributes */
  1764.  
  1765.             while (True) {
  1766.                 ads_entnext(ename, ename);
  1767.                 rent = ads_entget(ename);
  1768.                 if (rent == NULL) {
  1769.                     ads_printf(/*MSG91*/"UPDATE: Can't read attribute.\n");
  1770.                     break;
  1771.                 }
  1772.                 rb = entitem(rent, 0);  /* Entity type */
  1773.                 if (rb != NULL) {
  1774.                     if (strcmp(rb->resval.rstring, /*MSG0*/"ATTRIB") != 0)
  1775.                         break;
  1776.                     rb = entitem(rent, 2);  /* Attribute tag */
  1777.                     vb = entitem(rent, 1);  /* Attribute value */
  1778.                     if (rb != NULL && vb != NULL) {
  1779.                         int j;
  1780.  
  1781.                         for (j = 0; j < ELEMENTS(partatt); j++) {
  1782.                             if (strcmp(rb->resval.rstring,
  1783.                                        partatt[j].attname) == 0) {
  1784.                                 /* All right, we've found an attribute in
  1785.                                    the  table.   Edit  the new value into
  1786.                                    it, update the entity,  and  continue. */
  1787.                                 char *sp = vb->resval.rstring;
  1788.                                 char ebuf[32];
  1789.  
  1790.                                 sprintf(ebuf, "%.12g", *partatt[j].attvar);
  1791.                                 vb->resval.rstring = ebuf;
  1792.                                 if (ads_entmod(rent) != RTNORM) {
  1793.                                     ads_printf(
  1794.                                        /*MSG92*/"\
  1795. UPDATE:  Cannot update %s attribute.\n",
  1796.                                        partatt[j].attname);
  1797.                                 }
  1798.                                 /* Restore string buffer */
  1799.                                 vb->resval.rstring = sp;
  1800.                                 break;
  1801.                             }
  1802.                         }
  1803.                     }
  1804.                 }
  1805.                 ads_relrb(rent);
  1806.             }
  1807.         } else {
  1808.             ads_printf(/*MSG93*/"\nCannot retrieve %s.\n", ptab[i].partname);
  1809.         }
  1810.     }
  1811. }
  1812.  
  1813. #ifdef  HIGHC
  1814.  
  1815. /*  Earlier versions of High C put abort() in the same module with exit();
  1816.     ADS defines its own exit(), so we have to define our own abort(): */
  1817.  
  1818. static void abort(void)
  1819. {
  1820.     ads_abort("");
  1821. }
  1822.  
  1823. #endif  /* HIGHC */
  1824.