head	1.15;
access;
symbols;
locks
	ids:1.15; strict;
comment	@ * @;


1.15
date	2000.10.15.16.58.56;	author ids;	state Exp;
branches;
next	1.14;

1.14
date	99.05.12.20.17.16;	author ids;	state Exp;
branches;
next	1.13;

1.13
date	93.07.06.13.05.40;	author ids;	state Exp;
branches;
next	1.12;

1.12
date	93.06.26.20.57.05;	author ids;	state Exp;
branches;
next	1.11;

1.11
date	93.06.24.18.08.28;	author ids;	state Exp;
branches;
next	1.10;

1.10
date	93.06.22.12.46.38;	author ids;	state Exp;
branches;
next	1.9;

1.9
date	93.06.22.08.35.58;	author ids;	state Exp;
branches;
next	1.8;

1.8
date	93.06.21.20.53.17;	author ids;	state Exp;
branches;
next	1.7;

1.7
date	93.06.18.22.36.20;	author ids;	state Exp;
branches;
next	1.6;

1.6
date	93.06.18.22.24.49;	author ids;	state Exp;
branches;
next	1.5;

1.5
date	93.06.18.21.24.09;	author ids;	state Exp;
branches;
next	1.4;

1.4
date	93.06.11.23.05.35;	author ids;	state Exp;
branches;
next	1.3;

1.3
date	93.06.11.22.24.31;	author ids;	state Exp;
branches;
next	1.2;

1.2
date	93.06.11.20.46.31;	author ids;	state Exp;
branches;
next	1.1;

1.1
date	93.06.10.22.02.44;	author ids;	state Exp;
branches;
next	;


desc
@@


1.15
log
@changed to read data from a file rather than standard input.
@
text
@#include "gravity_support.h"
#include "GraphicsForGravity.h"
#include <stdio.h>
#include <math.h>
#include <string.h>


#define bool	int
#define TRUE	1
#define FALSE	0

#define G		6.673E-11   /* gravitational constant in SI units */

#define InitialTimeStep		1.0
/*
note: this is meant to be "very small" for the problem in question.
(Just in case by bad luck we're at the *start* (far point!) of a "comet-like"
 orbit...in that case all that cleverness in TIME_STEP_UP_BY comes to naught!)
*/

#define TIME_STEP_DOWN_BY	1.2
/*
for turning down the time step if there was an unsatisfactory update.
There's no point making it any smaller as it's MAX_REL_DISP that ultimately
controls accuracy. Any bigger and we waste time being "over-accurate" at
greater computational cost than just choosing a smaller MAX_REL_DISP.
*/

#define TIME_STEP_UP_BY		(1.0 + MAX_REL_DISP * (TIME_STEP_DOWN_BY - 1.0))
/*
for turning up the time step if there was a satisfactory update.
This choice ensures the time step never gets the chance to go spuriously high
at or around the far point of "comet-like" orbits. With a naive choice of
TIME_STEP_UP_BY, say just fixed at 1.01, the judging of safe time steps
by (temporarily very slow) changes in *position* is lulled into a false
sense of security. This definition avoids that problem. It also means at least
about a radian-worth of movement takes place before the time step
judders down; any finer-scale juddering would just be a waste.
*/


static double MAX_REL_DISP;	/* criterion for satisfactory update */
static int    PLOT_HOW_OFTEN;	/* one in this many states are plotted */
static int    PRINT_HOW_OFTEN;	/* one in this many states are printed */

static int state_count = 0;	/* out here for convenience */


/* private forward references */
static void dot (STATE *s, STATE_DOT *s_dot);
static bool WasSatisfactoryUpdate (STATE *CurrentState,
                                   double TimeStep,
                                   STATE *NextState);
static void ReadState (FILE *filehandle, STATE *AState);
static void PutState (STATE *AState, double ATimeStep);
static void PutEnsembleStatistics (STATE *AState);





void main(int argc, char *argv[])
{
  STATE  states[2];			/* TheState and NewState point here */
  STATE  *TheState, *NewState, *temp;	/* temp is just for swapping */
  double TimeStep;
  char   *filename;
  FILE   *filehandle;

  /*
   *	argument checking
   */
  if (argc != 2) {
    fprintf(stderr, "usage: %s <data-file-name>\n", argv[0]);
    exit(1);
  }/*endif*/

  filename=argv[1];

  /*
   *	do we have a readable file?
   */
  filehandle=fopen(filename, "r");
  if (filehandle == NULL) {
    fprintf(stderr, "%s: cannot open file \"%s\" for reading!\n", argv[0], filename);
    exit(1);
  }/*endif*/

  /*
   *	OK, we have a readable file - let's get the data!
   */
  if (fscanf(filehandle,
             "%lf %d %d%*[^\n]\n",
             &MAX_REL_DISP, &PLOT_HOW_OFTEN, &PRINT_HOW_OFTEN) != 3) {
    fprintf(stderr, "%s: bad input data\n", argv[0]);
    exit(1);
  }/*endif*/

  TheState = &states[0];
  NewState = &states[1];

  ReadState(filehandle, TheState);
  *NewState = *TheState;		/* to save having to update names etc */
  TimeStep = InitialTimeStep;

  InitializeGraphics(TheState);


  while (TRUE) {
    if (state_count % PLOT_HOW_OFTEN == 0) {
      DrawState(TheState);
    }/*endif*/;

    if (state_count % PRINT_HOW_OFTEN == 0) {
      PutState(TheState, TimeStep);
    }/*endif*/;

    while (! WasSatisfactoryUpdate(TheState, TimeStep, NewState)) {
      TimeStep /= TIME_STEP_DOWN_BY;
    }/*endwhile*/;
    temp = TheState; TheState = NewState; NewState = temp;
    TimeStep *= TIME_STEP_UP_BY;
    state_count++;
  }/*endwhile*/
}
             




/* the procedures at the core of the simulation... */

static void dot (STATE *s, STATE_DOT *s_dot)
/* finds the time derivative of the state. Treats first arg as read-only. */
{
  double mass_i, mass_j, multiplier;
  int    i,j,d;
  VECTOR i_to_j;

  s_dot->n = s->n;

  for (i=0; i<s->n; i++) {
    for (d=0; d<NDIMS; d++) {
      s_dot->bodydata[i].position_dot[d] = s->bodydata[i].velocity[d];
    }/*endfor*/
  }/*endfor*/

  for (i=0; i<s->n; i++) {
    CreateZeroVector (s_dot->bodydata[i].velocity_dot);
    mass_i = s->bodydata[i].mass;
    for (j=0; j<i; j++) {
      mass_j = s->bodydata[j].mass;

      if ((mass_i != 0.0) || (mass_j != 0.0)) {
        for (d=0; d<NDIMS; d++) {
          i_to_j[d] = s->bodydata[j].position[d] - s->bodydata[i].position[d];
        }/*endfor*/;
        multiplier = G / pow(lengthsquared(i_to_j), NDIMS/2.0);

        /* see if body i pulls on body j */
        if (mass_i != 0.0) {
          for (d=0; d<NDIMS; d++) {
            s_dot->bodydata[j].velocity_dot[d] -=
                    multiplier * mass_i * i_to_j[d];
          }/*endfor*/
        }/*endif*/

        /* see if body i gets pulled by body j */
        if (mass_j != 0.0) {
          for (d=0; d<NDIMS; d++) {
            s_dot->bodydata[i].velocity_dot[d] +=
                    multiplier * mass_j * i_to_j[d];
          }/*endfor*/
        }/*endif*/

      }/*endif*/
    }/*endfor*/
  }/*endfor*/
}





static bool WasSatisfactoryUpdate (STATE *CurrentState,
                                   double TimeStep,
                                   STATE *NextState)
/*
attempts update, and reports whether satisfactory.
treats first arg as read-only; assumes "static" parts of *NextState
are already filled in.


We use the Runge-Kutta method, as follows.

We define these intermediate estimates of the mean "speed" (change against time)
that the state changes at, between now and now + TimeStep:

  speed0 = dot(CurrentState)

  speed1 = dot(CurrentState + FRAC00 * TimeStep * speed0)

  speed2 = dot(CurrentState + FRAC10 * TimeStep * speed0
                            + FRAC11 * TimeStep * speed1)

  speed3 = dot(CurrentState + FRAC20 * TimeStep * speed0
                            + FRAC21 * TimeStep * speed1
                            + FRAC22 * TimeStep * speed2)

And so on, say up to speed_n-1.
We then weight these into a final estimate of the mean speed:

  speed  =   FRAC_n-1_0   * speed0
           + FRAC_n-1_1   * speed1
           + FRAC_n-1_2   * speed2
             .
             .
             .
           + FRAC_n-1_n-1 * speed_n-1

And finally, we say NextState = CurrentState + TimeStep * speed.

This agrees with the Taylor series expansion for the state up to some
particular high order in TimeStep (this order is not necessarily the "n" of
speed_n-1, etc, above, which we call the complexity),
leaving errors of the next order higher still
(or *relative* errors [error/TimeStep] of the same particular order).
However, it takes far less computation than any explicit evaluation
of the Taylor series would.

The magic numbers FRACij are given in the static local data structure below.
*/

{
  #define RK_IS_BUTCHER_ORDER_6


  #if defined(RK_IS_CLASSICAL_ORDER_4)
    #define RungeKuttaComplexity 4
    static double FRAC[RungeKuttaComplexity][RungeKuttaComplexity] = {
      { 0.5 },
      { 0,     0.5 },
      { 0,     0,     1 },
      { 1.0/6, 2.0/6, 2.0/6, 1.0/6 }
    };
  #elif defined(RK_IS_BUTCHER_ORDER_5)
    #define RungeKuttaComplexity 6
    static double FRAC[RungeKuttaComplexity][RungeKuttaComplexity] = {
      { 0.25 },
      { 0.125,   0.125 },
      { 0,      -0.5,    1 },
      { 3.0/16,  0,      0,        9.0/16 },
      {-3.0/7,   2.0/7,  12.0/7,  -12.0/7,  8.0/7 },
      { 7.0/90,  0,      32.0/90,  12.0/90, 32.0/90, 7.0/90 }
    };
  #elif defined(RK_IS_BUTCHER_ORDER_6)
    #define RungeKuttaComplexity 7
    static double FRAC[RungeKuttaComplexity][RungeKuttaComplexity] = {
      { 1.0/3 },
      { 0,        2.0/3 },
      { 1.0/12,   1.0/3,  -1.0/12 },
      {-1.0/16,   1.125,  -3.0/16,  -0.375 },
      { 0,        1.125,  -0.375,   -0.75,     0.5 },
      { 9.0/44,  -9.0/11,  63.0/44,  18.0/11,  0,      -16.0/11 },
      { 11.0/120, 0,       27.0/40,  27.0/40, -4.0/15, -4.0/15, 11.0/120 }
    };
  #endif


  STATE_DOT speed[RungeKuttaComplexity];
  VECTOR    current_i_to_j, new_i_to_j, change_in_i_to_j;
  int       i,j,d;


  /* compute NextState using Runge-Kutta as explained above */

  for (i=0; i<=RungeKuttaComplexity; i++) {
    copy_state(CurrentState, NextState);
    for (j=0; j<i; j++) {
      incr_state(NextState, FRAC[i-1][j] * TimeStep, &speed[j]);
    }/*endfor*/

    if (i<RungeKuttaComplexity) {
      dot(NextState, &speed[i]);
    }/*endif*/
  }/*endfor*/


  /* OK, that's NextState computed...now see if it's satisfactory! */

  for (i=0; i<CurrentState->n; i++) {
    for (j=0; j<i; j++) {
      if (CurrentState->bodydata[i].mass != 0 ||
          CurrentState->bodydata[j].mass != 0) {

        /* check non-trivial pair of bodies aren't swinging by too fast */
        for (d=0; d<NDIMS; d++) {
          current_i_to_j[d] = CurrentState->bodydata[j].position[d] -
                              CurrentState->bodydata[i].position[d];
          new_i_to_j[d]     = NextState->bodydata[j].position[d] -
                              NextState->bodydata[i].position[d];
          change_in_i_to_j[d] = new_i_to_j[d] - current_i_to_j[d];
        }/*endfor*/;
        if (lengthsquared(change_in_i_to_j) > pow(MAX_REL_DISP, 2.0) *
            lengthsquared(current_i_to_j)) {
          return FALSE;
        }/*endif*/

      }/*endif*/
    }/*endfor*/
  }/*endfor*/;

  return TRUE;
}







/* I/O functions. */

static void ReadState (FILE *filehandle, STATE *AState)
/* specialized to 3-D for now */
{
  AState->Time = 0.0;
  AState->n    = 0;

  while (AState->n < NMAX
         &&
         fscanf(filehandle,
                "%20s %lf  %lf %lf %lf  %lf %lf %lf%*[^\n]\n",
                /* NOTE: 20 in %20s is really NAMELENGTH */
                AState->bodydata[AState->n].name,
                &AState->bodydata[AState->n].mass,

                &AState->bodydata[AState->n].position[0],
                &AState->bodydata[AState->n].position[1],
                &AState->bodydata[AState->n].position[2],

                &AState->bodydata[AState->n].velocity[0],
                &AState->bodydata[AState->n].velocity[1],
                &AState->bodydata[AState->n].velocity[2]) == 8) {
    (AState->n)++;
  }/*endwhile*/
}





static void PutState (STATE *AState, double ATimeStep)
/* treats first arg as read-only */
{
  int i, d;

  printf("STATE NO. %d OF THIS RUN\n", state_count);

  printf("Time=%18g\n", AState->Time);
  printf("Time step to be tried next = %18g\n", ATimeStep);
  printf("Max relative displacement tolerated = %18g\n", MAX_REL_DISP);
  printf("[note: plotting one state in %d, printing one state in %d]\n",
           PLOT_HOW_OFTEN, PRINT_HOW_OFTEN);
  printf("-------------------- ------------------ ------------------ ------------------\n");
  printf("        BODY                MASS             POSITION           VELOCITY     \n");
  printf("-------------------- ------------------ ------------------ ------------------\n");

  for (i=0; i<AState->n; i++) {
    printf("%-20s %18g ", AState->bodydata[i].name, AState->bodydata[i].mass);
    for (d=0; d<NDIMS; d++) {
      printf("%18g %18g\n                                        ",
             AState->bodydata[i].position[d], AState->bodydata[i].velocity[d]);
    }/*endfor*/;
    printf("\n");
  }/*endfor*/;
  PutEnsembleStatistics(AState);
  printf("\n\n\n\n");
}





static void PutEnsembleStatistics (STATE *AState)
/* treats arg as read-only. Potential energy is specialized to 3-D for now. */
{
  double total_mass;
  VECTOR total_moment;
  VECTOR total_momentum;
  double total_potential_energy;
  double total_kinetic_energy;
  VECTOR centre_of_mass_position;
  VECTOR centre_of_mass_velocity;
  double centre_of_mass_virtual_kinetic_energy;
  int    i, j, d;
  VECTOR i_to_j;

  total_mass = 0.0;
  CreateZeroVector(total_moment);
  CreateZeroVector(total_momentum);
  total_potential_energy = 0.0;
  total_kinetic_energy = 0.0;
  for (i=0; i<AState->n; i++) {
    if (AState->bodydata[i].mass != 0.0) {
      total_mass += AState->bodydata[i].mass;
      for (d=0; d<NDIMS; d++) {
        total_moment[d]   += AState->bodydata[i].mass * AState->bodydata[i].position[d];
        total_momentum[d] += AState->bodydata[i].mass * AState->bodydata[i].velocity[d];
      }/*endfor*/;
      for (j=i+1; j<AState->n; j++) {
        for (d=0; d<NDIMS; d++) {
          i_to_j[d] = AState->bodydata[j].position[d] - AState->bodydata[i].position[d];
        }/*endfor*/;
        total_potential_energy +=
         - G * AState->bodydata[i].mass * AState->bodydata[j].mass / sqrt(lengthsquared(i_to_j));
      }/*endfor*/;
      total_kinetic_energy +=
       0.5 * AState->bodydata[i].mass * lengthsquared(AState->bodydata[i].velocity);
    }/*endif*/
  }/*endfor*/

  for (d=0; d<NDIMS; d++) {
    centre_of_mass_position[d] = total_moment[d] / total_mass;
    centre_of_mass_velocity[d] = total_momentum[d] / total_mass;
  }/*endfor*/;
  centre_of_mass_virtual_kinetic_energy =
   lengthsquared(total_momentum) / (2.0 * total_mass);

  printf("ENSEMBLE\n");
  printf(" (AS CENTRE OF MASS) ");
  printf("%18g ", total_mass);
  for (d=0; d<NDIMS; d++) {
    printf("%18g %18g\n                                        ",
            centre_of_mass_position[d], centre_of_mass_velocity[d]);
  }/*endfor*/;
  printf("\n");
  
  printf("FURTHER ENSEMBLE STATISTICS:\n");
  printf(" centre-of-mass virtual kinetic energy = %18g\n\n",
           centre_of_mass_virtual_kinetic_energy);
  printf(" IN STANDARD COORDINATE FRAME:\n");
  printf("  potential energy = %18g\n", total_potential_energy);
  printf("  kinetic energy   = %18g\n", total_kinetic_energy);
  printf("  total energy     = %18g\n\n", total_potential_energy + total_kinetic_energy);
  printf(" IN ENSEMBLE'S REST FRAME:\n");
  printf("  potential energy = %18g\n", total_potential_energy);
  printf("  kinetic energy   = %18g\n", total_kinetic_energy -
                 centre_of_mass_virtual_kinetic_energy);
  printf("  total energy     = %18g\n\n",
            total_potential_energy + total_kinetic_energy -
                 centre_of_mass_virtual_kinetic_energy);
}
@


1.14
log
@at last! we seem to be able to dispense with the rigmarole of
 "plus give proper ANSI headers to..." as a hanger-on of "#include <math.h>".
 (It was things like sqrt, pow, sin, cos, atan2 that used to cause all the
  trouble - presumably they once had non-ANSI headers by default or something.)

Note: we do still have to define pi by hand. There's always something!
@
text
@d54 1
a54 1
static void GetState (STATE *AState);
d67 2
d70 24
a93 1
  if (scanf("%lf %d %d%*[^\n]\n",
d102 1
a102 1
  GetState(TheState);
d324 1
a324 1
static void GetState (STATE *AState)
d332 13
a344 12
         scanf("%20s %lf  %lf %lf %lf  %lf %lf %lf%*[^\n]\n",
               /* NOTE: 20 in %20s is really NAMELENGTH */
               AState->bodydata[AState->n].name,
               &AState->bodydata[AState->n].mass,

               &AState->bodydata[AState->n].position[0],
               &AState->bodydata[AState->n].position[1],
               &AState->bodydata[AState->n].position[2],

               &AState->bodydata[AState->n].velocity[0],
               &AState->bodydata[AState->n].velocity[1],
               &AState->bodydata[AState->n].velocity[2]) == 8) {
@


1.13
log
@improved checking of whether an update was satisfactory: doesn't waste time
 checking pairs of test bodies. (So you can now start test bodies on top of
 each other, if you want.)
@
text
@a4 3
/* plus give proper ANSI headers to sqrt and pow */
extern double sqrt(double);
extern double pow (double, double);
@


1.12
log
@new cleaner C conventions used, bearing a much stronger resemblance to
 Modula's modularization.
@
text
@d271 6
a276 2
      for (d=0; d<NDIMS; d++) {
        current_i_to_j[d] = CurrentState->bodydata[j].position[d] -
d278 9
a286 6
        new_i_to_j[d] = NextState->bodydata[j].position[d] -
                          NextState->bodydata[i].position[d];
        change_in_i_to_j[d] = new_i_to_j[d] - current_i_to_j[d];
      }/*endfor*/;
      if (lengthsquared(change_in_i_to_j) > MAX_REL_DISP * MAX_REL_DISP * lengthsquared(current_i_to_j)) {
        return FALSE;
@


1.11
log
@MAJOR UPGRADE:
 * now uses Runge-Kutta order 6: a vast improvement in speed and/or accuracy.
 * much cleaner implementation of the Runge-Kutta algorithm, ready to work
   for any order, given appropriate magic numbers for that order.
 * popular magic numbers for orders 4, 5 and 6 present in the source-code,
   with conditional compilation to choose between them (though I think
   order 6 is the best for most uses, hence the choice of it as the
   "active" one).
moved some low-level stuff to gravity_support
minor C tidyups.
@
text
@a0 1
#include "gravity.h"
d2 1
d9 1
a9 1
#include "GraphicsForGravity.h"
d11 3
d49 1
a49 1

d52 8
a60 3
static int state_count = 0;	/* out here for convenience */


d111 2
a112 1
void dot (STATE *s, STATE_DOT *s_dot)
d163 3
a165 3
bool WasSatisfactoryUpdate (STATE *CurrentState,
                            double TimeStep,
                            STATE *NextState)
d167 5
d295 1
a295 1
void GetState (STATE *AState)
d323 2
a324 1
void PutState (STATE *AState, double ATimeStep)
d355 2
a356 1
void PutEnsembleStatistics (STATE *AState)
@


1.10
log
@minor C tidyups.
@
text
@a11 2
#define TRUE		1
#define FALSE		0
a13 1

d154 3
a156 3
BOOLEAN WasSatisfactoryUpdate (STATE *CurrentState,
                               double TimeStep,
                               STATE *NextState)
d158 1
a158 3
We use the Runge-Kutta order-4 method, as follows.
(Well, the imperative equivalent: there's no need to have all these variables
 actually existing as we can systematically re-use some of them.)
d163 3
a165 4
  speed1 = dot(CurrentState)
  speed2 = dot(CurrentState + FRAC1 * TimeStep * speed1)
  speed3 = dot(CurrentState + FRAC2 * TimeStep * speed2)
  speed4 = dot(CurrentState + FRAC3 * TimeStep * speed3)
d167 2
a168 1
We weight these into a final estimate of the mean speed:
d170 14
a183 1
  speed  = WEIGHT1*speed1 + WEIGHT2*speed2 + WEIGHT3*speed3 + WEIGHT4*speed4
d187 6
a192 3
This agrees with the Taylor series expansion for the state up to order 4 in
TimeStep, leaving errors of order 5 (or *relative* errors [error/TimeStep]
of order 4). However, it takes far less computation than any explicit evaluation
d195 1
a195 2
The magic numbers are FRAC1 = 0.5, FRAC2 = 0.5, FRAC3 = 1.0,
WEIGHT1 = 1/6, WEIGHT2 = 2/6, WEIGHT3 = 2/6, WEIGHT4 = 1/6.
d199 33
a231 7
  #define FRAC1   0.5
  #define FRAC2   0.5
  #define FRAC3   1.0
  #define WEIGHT1 (1.0/6.0)
  #define WEIGHT2 (2.0/6.0)
  #define WEIGHT3 (2.0/6.0)
  #define WEIGHT4 (1.0/6.0)
d233 2
a234 2
  STATE     IntermediateState;
  STATE_DOT speed;
d238 12
a249 8
  /*
  get IntermediateState adequately initialized by doing what incr_state
  doesn't bother with for efficiency reasons.
  NOTE: even here we don't bother with the names of the bodies.
  */
  IntermediateState.n = CurrentState->n;
  for (i=0; i<CurrentState->n; i++) {
    IntermediateState.bodydata[i].mass = CurrentState->bodydata[i].mass;
a251 14
  dot(CurrentState, &speed);		/* speed = speed1 above */
  incr_state(CurrentState, WEIGHT1 * TimeStep, &speed, NextState);

  incr_state(CurrentState, FRAC1 * TimeStep,   &speed, &IntermediateState);
  dot(&IntermediateState, &speed);	/* speed = speed2 above */
  incr_state(NextState,    WEIGHT2 * TimeStep, &speed, NextState);

  incr_state(CurrentState, FRAC2 * TimeStep,   &speed, &IntermediateState);
  dot(&IntermediateState, &speed);	/* speed = speed3 above */
  incr_state(NextState,    WEIGHT3 * TimeStep, &speed, NextState);

  incr_state(CurrentState, FRAC3 * TimeStep,   &speed, &IntermediateState);
  dot(&IntermediateState, &speed);      /* speed = speed4 above */
  incr_state(NextState,    WEIGHT4 * TimeStep, &speed, NextState);
d279 1
a279 1
/* I/O functions */
a406 27
}





/* low-level functions. */

void incr_state (STATE *state, double delta_t, STATE_DOT *speed,
                 STATE *newstate)
{
  int i,d;

  newstate->Time = state->Time + delta_t;

  for (i=0; i<state->n; i++) {
    for (d=0; d<NDIMS; d++) {
      newstate->bodydata[i].position[d] =
         state->bodydata[i].position[d] + delta_t *
         speed->bodydata[i].position_dot[d];
    }/*endfor*/
    for (d=0; d<NDIMS; d++) {
      newstate->bodydata[i].velocity[d] =
         state->bodydata[i].velocity[d] + delta_t *
         speed->bodydata[i].velocity_dot[d];
    }/*endfor*/
  }/*endfor*/
@


1.9
log
@better procedure names, following the "IC DoC conventions"
(Put___ for standard output, Draw___ for graphics).
@
text
@d45 3
a47 3
double MAX_REL_DISP;	/* criterion for satisfactory update */
int    PLOT_HOW_OFTEN;	/* one in this many states are plotted graphically */
int    PRINT_HOW_OFTEN;	/* one in this many states are printed */
d53 1
a53 1
int state_count = 0;	/* out here for convenience */
@


1.8
log
@added graphics plotting capability
moved some stuff into gravity_support.
@
text
@d84 1
a84 1
      PlotGraphically(TheState);
d88 1
a88 1
      Print(TheState, TimeStep);
d283 1
a283 1
void Print (STATE *AState, double ATimeStep)
d306 1
a306 1
  PrintEnsembleStatistics(AState);
d314 1
a314 1
void PrintEnsembleStatistics (STATE *AState)
@


1.7
log
@purely aesthetic change: order of procedures is now more "top-down",
 matching the order in gravity.h.
@
text
@d2 1
d9 1
d45 3
a47 2
double MAX_REL_DISP;		/* criterion for satisfactory update */
int    PRINT_HOW_OFTEN;		/* one in this many states are displayed */
d66 2
a67 1
  if (scanf("%lf %d%*[^\n]\n", &MAX_REL_DISP, &PRINT_HOW_OFTEN) != 2) {
d79 1
d81 1
d83 4
d291 3
a293 1
  printf("[note: only 1 computed state in %d is being displayed]\n", PRINT_HOW_OFTEN);
a407 29
}





void CreateZeroVector (VECTOR AVector)
{
  int d;

  for (d=0; d<NDIMS; d++) {
    AVector[d] = 0.0;
  }/*endfor*/
}





double lengthsquared (VECTOR v)
{
  double total;
  int    d;

  total = 0.0;
  for (d=0; d<NDIMS; d++) {
    total += v[d] * v[d];
  }/*endfor*/;
  return total;
@


1.6
log
@inner loop optimization in "dot":
compute the i_to_j vectors only for the "lower triangle" of combinations of
 i and j, since the other combinations are just (the negations of) those
 same vectors again.
This gives a near doubling of speed! (And yes, I checked!)
@
text
@d94 1
a94 126
void GetState (STATE *AState)
/* specialized to 3-D for now */
{
  AState->Time = 0.0;
  AState->n    = 0;

  while (AState->n < NMAX
         &&
         scanf("%20s %lf  %lf %lf %lf  %lf %lf %lf%*[^\n]\n",
               /* NOTE: 20 in %20s is really NAMELENGTH */
               AState->bodydata[AState->n].name,
               &AState->bodydata[AState->n].mass,

               &AState->bodydata[AState->n].position[0],
               &AState->bodydata[AState->n].position[1],
               &AState->bodydata[AState->n].position[2],

               &AState->bodydata[AState->n].velocity[0],
               &AState->bodydata[AState->n].velocity[1],
               &AState->bodydata[AState->n].velocity[2]) == 8) {
    (AState->n)++;
  }/*endwhile*/
}





void Print (STATE *AState, double ATimeStep)
{
  int i, d;

  printf("STATE NO. %d OF THIS RUN\n", state_count);

  printf("Time=%18g\n", AState->Time);
  printf("Time step to be tried next = %18g\n", ATimeStep);
  printf("[note: only 1 computed state in %d is being displayed]\n", PRINT_HOW_OFTEN);
  printf("-------------------- ------------------ ------------------ ------------------\n");
  printf("        BODY                MASS             POSITION           VELOCITY     \n");
  printf("-------------------- ------------------ ------------------ ------------------\n");

  for (i=0; i<AState->n; i++) {
    printf("%-20s %18g ", AState->bodydata[i].name, AState->bodydata[i].mass);
    for (d=0; d<NDIMS; d++) {
      printf("%18g %18g\n                                        ",
             AState->bodydata[i].position[d], AState->bodydata[i].velocity[d]);
    }/*endfor*/;
    printf("\n");
  }/*endfor*/;
  PrintEnsembleStatistics(AState);
  printf("\n\n\n\n");
}





void PrintEnsembleStatistics (STATE *AState)
{
  double total_mass;
  VECTOR total_moment;
  VECTOR total_momentum;
  double total_potential_energy;
  double total_kinetic_energy;
  VECTOR centre_of_mass_position;
  VECTOR centre_of_mass_velocity;
  double centre_of_mass_virtual_kinetic_energy;
  int    i, j, d;
  VECTOR i_to_j;

  total_mass = 0.0;
  CreateZeroVector(total_moment);
  CreateZeroVector(total_momentum);
  total_potential_energy = 0.0;
  total_kinetic_energy = 0.0;
  for (i=0; i<AState->n; i++) {
    if (AState->bodydata[i].mass != 0.0) {
      total_mass += AState->bodydata[i].mass;
      for (d=0; d<NDIMS; d++) {
        total_moment[d]   += AState->bodydata[i].mass * AState->bodydata[i].position[d];
        total_momentum[d] += AState->bodydata[i].mass * AState->bodydata[i].velocity[d];
      }/*endfor*/;
      for (j=i+1; j<AState->n; j++) {
        for (d=0; d<NDIMS; d++) {
          i_to_j[d] = AState->bodydata[j].position[d] - AState->bodydata[i].position[d];
        }/*endfor*/;
        total_potential_energy +=
         - G * AState->bodydata[i].mass * AState->bodydata[j].mass / sqrt(lengthsquared(i_to_j));
      }/*endfor*/;
      total_kinetic_energy +=
       0.5 * AState->bodydata[i].mass * lengthsquared(AState->bodydata[i].velocity);
    }/*endif*/
  }/*endfor*/

  for (d=0; d<NDIMS; d++) {
    centre_of_mass_position[d] = total_moment[d] / total_mass;
    centre_of_mass_velocity[d] = total_momentum[d] / total_mass;
  }/*endfor*/;
  centre_of_mass_virtual_kinetic_energy =
   lengthsquared(total_momentum) / (2.0 * total_mass);

  printf("ENSEMBLE\n");
  printf(" (AS CENTRE OF MASS) ");
  printf("%18g ", total_mass);
  for (d=0; d<NDIMS; d++) {
    printf("%18g %18g\n                                        ",
            centre_of_mass_position[d], centre_of_mass_velocity[d]);
  }/*endfor*/;
  printf("\n");
  
  printf("FURTHER ENSEMBLE STATISTICS:\n");
  printf(" centre-of-mass virtual kinetic energy = %18g\n\n",
           centre_of_mass_virtual_kinetic_energy);
  printf(" IN STANDARD COORDINATE FRAME:\n");
  printf("  potential energy = %18g\n", total_potential_energy);
  printf("  kinetic energy   = %18g\n", total_kinetic_energy);
  printf("  total energy     = %18g\n\n", total_potential_energy + total_kinetic_energy);
  printf(" IN ENSEMBLE'S REST FRAME:\n");
  printf("  potential energy = %18g\n", total_potential_energy);
  printf("  kinetic energy   = %18g\n", total_kinetic_energy -
                 centre_of_mass_virtual_kinetic_energy);
  printf("  total energy     = %18g\n\n",
            total_potential_energy + total_kinetic_energy -
                 centre_of_mass_virtual_kinetic_energy);
}

a95 3



d243 133
a375 1
/* low-level procedures */
@


1.5
log
@MAJOR REVISION:
 * Runge-Kutta order 4 algorithm installed. (huge improvement in accuracy!)
 * better choice of time-step control constants.
 * proper ANSI headers for sqrt, pow.
 * improved comments, documenting rationale for choice of algorithm,
   time-step control constants, etc.
 * minor tidyups.
@
text
@d226 1
a226 1
  double massj, multiplier;
d240 5
a244 3
    for (j=0; j<s->n; j++) {
      massj = s->bodydata[j].mass;
      if ((j != i) && (massj != 0.0)) {
d248 18
a265 4
        multiplier = G * massj / pow(lengthsquared(i_to_j), NDIMS/2.0);
        for (d=0; d<NDIMS; d++) {
          s_dot->bodydata[i].velocity_dot[d] += multiplier * i_to_j[d];
        }/*endfor*/
@


1.4
log
@made MAX_REL_DISP and PRINT_HOW_OFTEN be read from the standard input
 before reading the state. (previously they were constants)
@
text
@d4 3
d9 1
d12 29
a41 3
#define G		6.673E-11   /* gravitational constant in SI units */
#define TIME_STEP_DOWN_BY	2.0 /* for unsatisfactory updates */
#define TIME_STEP_UP_BY		1.1 /* for satisfactory updates */
d59 2
a60 2
  State  states[2];			/* TheState and NewState point here */
  State  *TheState, *NewState, *temp;	/* temp is just for swapping */
a62 1

d73 1
a73 1
  TimeStep = 1.0;
d94 1
a94 1
void GetState (State *AState)
d120 3
a122 1
void Print (State *AState, double ATimeStep)
d149 3
a151 1
void PrintEnsembleStatistics (State *AState)
d181 1
a181 1
         (- G * AState->bodydata[i].mass * AState->bodydata[j].mass / length(i_to_j));
d184 1
a184 1
       0.5 * AState->bodydata[i].mass * pow(length(AState->bodydata[i].velocity), 2.0);
d193 1
a193 1
   pow(length(total_momentum), 2.0) / (2.0 * total_mass);
d224 1
a224 3
BOOLEAN WasSatisfactoryUpdate (State *CurrentState,
                               double TimeStep,
                               State *NextState)
d226 3
a228 11
  int    i;
  VECTOR HalfWayPositions[NMAX];
  int    j;
  double multiplier;
  VECTOR acceleration;
  double massj;
  VECTOR current_i_to_j;
  VECTOR halfway_i_to_j;
  VECTOR new_i_to_j;
  VECTOR change_in_i_to_j;
  int    d;
d230 1
a230 1
  NextState->Time = CurrentState->Time + TimeStep;
d232 1
a232 1
  for (i=0; i<CurrentState->n; i++) {
d234 1
a234 3
      HalfWayPositions[i][d] =
       CurrentState->bodydata[i].position[d] +
        0.5 * TimeStep * CurrentState->bodydata[i].velocity[d];
d236 1
a236 1
  }/*endfor*/;
d238 4
a241 4
  for (i=0; i<CurrentState->n; i++) {
    CreateZeroVector (acceleration);
    for (j=0; j<CurrentState->n; j++) {
      massj = CurrentState->bodydata[j].mass;
d244 1
a244 1
          halfway_i_to_j[d] = HalfWayPositions[j][d] - HalfWayPositions[i][d];
d246 1
a246 1
        multiplier = G * massj / pow(length(halfway_i_to_j), (double) NDIMS);
d248 1
a248 1
          acceleration[d] += multiplier * halfway_i_to_j[d];
d251 77
a327 12
    }/*endfor*/; /* acceleration of i now computed */
    for (d=0; d<NDIMS; d++) {
      NextState->bodydata[i].velocity[d] = 
       CurrentState->bodydata[i].velocity[d] +
        TimeStep * acceleration[d];
    }/*endfor*/;
    for (d=0; d<NDIMS; d++) {
      NextState->bodydata[i].position[d] =
       CurrentState->bodydata[i].position[d] +
        TimeStep * (CurrentState->bodydata[i].velocity[d] +
                     NextState->bodydata[i].velocity[d]) / 2.0;
    }/*endfor*/;
d329 1
a329 1
    /* now check if latest computed position is satisfactory */
d331 1
d340 1
a340 1
      if (length(change_in_i_to_j) > MAX_REL_DISP * length(current_i_to_j)) {
d357 25
d395 1
a395 1
double length (VECTOR v)
d404 1
a404 1
  return sqrt(total);
@


1.3
log
@now reads state from standard input instead of embedded in source code.
@
text
@a9 1
#define MAX_REL_DISP		0.1 /* to test if satisfactory update */
a11 1
#define PRINT_HOW_OFTEN		10000/* one in this many states are displayed */
d13 2
d19 1
d27 1
a27 1
void main()
d33 6
d45 1
@


1.2
log
@efficiency improvement using pointers even for value-passing
(when the thing being passed is large, ie, a state).
minor I/O and low-level function tidyups.
@
text
@d35 1
a35 1
  Initialize(TheState);
d57 2
a58 1
void Initialize (State *AState)
d61 1
d63 16
a78 27
  /*.........DATA FROM .dat FILES STARTS HERE.........*/
  AState->n = 2;

  strncpy(AState->bodydata[0].name, "Sun", NAMELENGTH+1);
  AState->bodydata[0].mass = 2.E30;
  Create3DVector(AState->bodydata[0].position, 0.0, 0.0, 0.0);
  Create3DVector(AState->bodydata[0].velocity, 0.0, 0.0, 0.0);

  strncpy(AState->bodydata[1].name, "Earth", NAMELENGTH+1);
  AState->bodydata[1].mass = 6.E24;
  Create3DVector(AState->bodydata[1].position, 1.5E11, 0.0, 0.0);
  Create3DVector(AState->bodydata[1].velocity, 0.0, 25000.0, 0.0);
  /*..........DATA FROM .dat FILES ENDS HERE..........*/
}



void Create3DVector (VECTOR AVector, double x,double y,double z)
{
  if (NDIMS == 3) {
    AVector[0] = x;
    AVector[1] = y;
    AVector[2] = z;
  } else {
    fprintf(stderr, "ERROR!!! Cannot create 3D vector in %d-space.\n, NDIMS");
    exit(1);
  }/*endif*/
a79 1

@


1.1
log
@Initial revision
@
text
@d13 1
a13 1
#define PRINT_HOW_OFTEN		1000/* one in this many states are displayed */
d28 2
a29 1
  State  TheState, NewState;
d32 5
a36 1
  Initialize(&TheState);
d44 1
a44 1
    while (! WasSatisfactoryUpdate(TheState, TimeStep, &NewState)) {
d47 1
a47 1
    TheState = NewState;
d93 1
a93 1
void Print (State AState, double ATimeStep)
d99 2
a100 2
  printf("Time=%g\n", AState.Time);
  printf("Time step to be tried next = %g\n", ATimeStep);
d106 2
a107 2
  for (i=0; i<AState.n; i++) {
    printf("%20s %18g ", AState.bodydata[i].name, AState.bodydata[i].mass);
d110 1
a110 1
             AState.bodydata[i].position[d], AState.bodydata[i].velocity[d]);
d120 1
a120 1
void PrintEnsembleStatistics (State AState)
d138 3
a140 3
  for (i=0; i<AState.n; i++) {
    if (AState.bodydata[i].mass != 0.0) {
      total_mass += AState.bodydata[i].mass;
d142 2
a143 2
        total_moment[d]   += AState.bodydata[i].mass * AState.bodydata[i].position[d];
        total_momentum[d] += AState.bodydata[i].mass * AState.bodydata[i].velocity[d];
d145 1
a145 1
      for (j=i+1; j<AState.n; j++) {
d147 1
a147 1
          i_to_j[d] = AState.bodydata[j].position[d] - AState.bodydata[i].position[d];
d150 1
a150 1
         (- G * AState.bodydata[i].mass * AState.bodydata[j].mass / length(i_to_j));
d153 1
a153 1
       0.5 * AState.bodydata[i].mass * ToPower(length(AState.bodydata[i].velocity), 2);
d162 1
a162 1
   ToPower(length(total_momentum), 2) / (2.0 * total_mass);
d174 1
a174 1
  printf(" centre-of-mass virtual kinetic energy = %g\n\n",
d177 3
a179 3
  printf("  potential energy = %g\n", total_potential_energy);
  printf("  kinetic energy   = %g\n", total_kinetic_energy);
  printf("  total energy     = %g\n\n", total_potential_energy + total_kinetic_energy);
d181 2
a182 2
  printf("  potential energy = %g\n", total_potential_energy);
  printf("  kinetic energy   = %g\n", total_kinetic_energy -
d184 1
a184 1
  printf("  total energy     = %g\n\n",
d193 1
a193 1
BOOLEAN WasSatisfactoryUpdate (State  CurrentState,
a195 1
/* attempts update, and reports whether satisfactory */
d209 1
a209 2
  NextState->Time = CurrentState.Time + TimeStep;
  NextState->n = CurrentState.n;
d211 1
a211 3
  for (i=0; i<CurrentState.n; i++) {
    strncpy(NextState->bodydata[i].name, CurrentState.bodydata[i].name, NAMELENGTH+1);
    NextState->bodydata[i].mass = CurrentState.bodydata[i].mass;
d214 2
a215 2
       CurrentState.bodydata[i].position[d] +
        0.5 * TimeStep * CurrentState.bodydata[i].velocity[d];
d219 1
a219 1
  for (i=0; i<CurrentState.n; i++) {
d221 2
a222 2
    for (j=0; j<CurrentState.n; j++) {
      massj = CurrentState.bodydata[j].mass;
d227 1
a227 1
        multiplier = G * massj / ToPower(length(halfway_i_to_j), NDIMS);
d235 1
a235 1
       CurrentState.bodydata[i].velocity[d] +
d240 2
a241 2
       CurrentState.bodydata[i].position[d] +
        TimeStep * (CurrentState.bodydata[i].velocity[d] +
d249 2
a250 2
        current_i_to_j[d] = CurrentState.bodydata[j].position[d] -
                              CurrentState.bodydata[i].position[d];
a294 16
}





double ToPower (double r, int c)
{
  double product;
  int    count;

  product = 1.0;
  for (count=0; count<c; count++) {
    product *= r;
  }/*endfor*/;
  return product;
@