The magazine of the Melbourne PC User Group
Dancing Stars - Part 6
 Let's Do It In C++  
Ken Holmes

Continued from PC Update for May 2002
I hope the series has stirred some interest and provided a small insight into the workings of the mysterious gravity force in shaping the universe. This final Part will deal with the same situation as in Part 5 but using C language. Listing 6 was written in the DOS-based Version 3.0 of Borland Turbo C++ but is more like ordinary C as it does not use more modern features. Real programmers will not be impressed by having all code in the main module, but chopping it up into many small procedures requires them to be declared, defined and called which would consume more space. With a reasonable short program, I prefer to see it in one piece and feel it is easier to follow; others may feel the opposite.
 
An exception is the procedure, called "setscreen" here, to bring up the graphics screen. This is used for all graphics programs, as with Basic's much more expedient SCREEN command. It also needs the egavgaf.obj file to obtain the Borland Graphics Interface and the program needs to be handled as a Project. QBasic has a severe limitation of 64 kB of memory, including both code and data, but by using the Huge memory model and far variables in C we could handle many more stars. With 30 stars, the five 30 x 30 arrays of 4 byte floats need 18,000 bytes which QBasic can handle but this rises as the square of the number of stars so it would be in trouble by about 35 stars. 
Depending on the compiler and by specific allocation of memory, many more stars can be handled. Here, we will settle for 60 as this sufficiently clutters the screen. You will also note that the include files are at the top of the Listing instead of in a separate .h file. You may wish to re-arrange the program into separate files and procedures or, if using a Windows (R) -based programming language, follow your usual practices.

Note in the Listing that we use a 0.01 time interval and do 100 gravity calculations in the innermost loop, ie. one time unit, before doing a graphics plotting calculation. Where velocities reach 2 pixels per time unit we are doing a calculation at each 0.02 pixels along a stars path which should provide a reasonably accurate calculation of encounters as close as .2 pixels. With a star following a highly curved path and too high time intervals, the calculation lets the star "skid out" and gain apparent energy. This would accentuate any tendency for stars to "evaporate" from the cluster. It illustrates that, the closer you investigate things, the more you find to investigate.

As an additional feature, we plot, at the top and at the start of each 640 time units, a small circle showing the velocity of each star against its distance from the centre of the group. On average, these should hover about an inverse squared relationship; there is a suggestion of this in Figure 11 but masked by a considerable scatter in velocity as stars are at different phases of interaction with neighbours. Note some stars approaching maxd with high velocity, in the process of escaping. At the start of this 4th Run of 640 time units, there were 47 stars left within maxd.


Figure 11
With a much larger number of stars, it would be interesting to plot a histogram of the number of stars in different ranges of distance; one relevant characteristic of a real cluster must be its density distribution. Another feature is that we check at each 640 time units for stars that have moved away by more than maxd and assume that they won't come back; they are subsequently excluded to cut down on unnecessary calculation and speed up the program. When the cluster evaporates down to 30 stars, for example, it runs four times as fast.
In extended runs evaporation distils down to 3 or 4 stars in close orbits and it doesn't seem possible to get a sustained group with a reasonable number of stars.
 
However, we have come to realise that nothing is permanent and the universe continues to evolve; even real global clusters are probably constantly losing stars and may be no more in a few hundred billion years. Smaller bodies such as planets and brown dwarfs would likely be the first to go, much as do smaller molecules evaporate first from such mixtures as petrol.
 
'----------------------------------------------------------

'Listing 6        'Global cluster'



#include <graphics.h>

#include <stdlib.h>   //rand,random,randomize

#include <math.h> //sqrt

#include <stdio.h>    //sprintf

#include <conio.h>    //kbhit,getch



float c[60],x[60],y[60],z[60],dc[60];  //colr, posn & distance

float vx[60], vy[60], vz[60], vi[60];   //velocity comps & vel

float far dsq[60][60], a[60][60]; //dist. sqrd & accel. i to j

float far ax[60][60], ay[60][60], az[60][60];  // x/y/z accels

float atx[60], aty[60], atz[60]; // x/y/z comps of total accel

float xm = 480, ym = 240;                //x & y mean position

float eyexr=400, eyexl =240, eyey =0, eyez =-600; //eyes' posn

float xl, xr, yb, zf, df = 1000, vf = 300, maxd = 500;

float g = 1, ti=.01;                  //gravity, time interval

int n, m, r, i, j, ns = 60, nr, rm;

float tx, tvx, ty, tvy, tz, tvz, xd, yd, zd, d;

char msg[40];                 //for formatting screen printing

void setscreen();        //Declare VGA 640x480 setup procedure

void main() {                                    //main module

 setscreen(); randomize();    //set VGA and randomise on timer

 for (i = 0; i < ns; i++) {          //allocate posns and vels

  c[i] = i % 15 + 1;                  //individual star colour

  x[i] = rand() % 160 - 80;  y[i] = rand() % 160 - 80;

  z[i] = rand() % 160 - 80;//ie. in x/y/z range +-80/+-80/+-80

  dc[i]=sqrt(x[i]*x[i]+y[i]*y[i]+z[i]*z[i]);//dist from centre

  if(dc[i] < 20) dc[i] = 20;   //avoid excess veloc. at centre

  rm = int(df / sqrt(dc[i]));          //random number maximum

  vx[i] = (float(random(rm) -rm/2))/ vf;  //ie. velocity comps

  vy[i] = (float(random(rm) -rm/2))/ vf;  //...in range +-rm/2

  vz[i] = (float(random(rm) -rm/2))/ vf;

 }

 do { r++;                       //start of outer program loop

  cleardevice();      //each time, clear screen and zero totals

  nr = 0; tx = 0; tvx = 0; ty = 0; tvy = 0; tz = 0; tvz = 0;

  for (i = 0; i < ns; i++) {     //get 'creep' of central group

   if (dc[i] < maxd) {  //if star was still in range to maxd...

    dc[i] =sqrt(x[i]*x[i] +y[i]*y[i]+ z[i]*z[i]); //check again

    if (dc[i] < maxd) { nr++;     //count number STILL in range

     setcolor(c[i]); //set star colour & plot velocity versus..

     circle (dc[i], 130 -vi[i]*100, 3);//..distance with circle

     tx += x[i]; ty += y[i]; tz += z[i];  //total dist & vel...

     tvx += vx[i]; tvy += vy[i]; tvz += vz[i]; //...comps of nr

    } //end of  if (dc <...  ie. STILL in range

   } //end of  if (dc <...   ie. previously in range

  } //end of for (i = 0;...

  setcolor (15); line (320, 0, 320, 380);     //draw background

  line (0, 30, 640, 30); line (0, 130, 640, 130);

  line (0, 380, 640, 380); rectangle (0, 0, 639, 479);

  line (maxd, 40, maxd, 120);

  sprintf(msg, "Run %d,  %d stars left", r, nr);

  outtextxy(140, 10, msg);

  sprintf(msg, "distance  maxd"); outtextxy(maxd-70, 70, msg);

  sprintf(msg, "velocity = 1");

  outtextxy(530, 20, msg); outtextxy(530, 370, msg);

  sprintf(msg, "velocity = 0");

  outtextxy(530, 120, msg); outtextxy(530, 470, msg);

  sprintf(msg, "M"); outtextxy(325, 200, msg);

  sprintf(msg, "I"); outtextxy(325, 220, msg);

  sprintf(msg, "R"); outtextxy(325, 240, msg);

  outtextxy(325, 260, msg); outtextxy(325, 300, msg);

  sprintf(msg, "O"); outtextxy (325, 280, msg);

  tx /= nr; ty /= nr; tz /= nr;        //average position of nr

  tvx /= nr; tvy /= nr; tvz /= nr;     //average velocity of nr

  for (i = 0; i < ns; i++) {             //adjust only stars...

   if (dc[i] < maxd) {                   //...still in range...

    x[i] -= tx; y[i] -= ty; z[i] -= tz;    //...for position...

    vx[i] -= tvx; vy[i] -= tvy; vz[i] -= tvz; //...and velocity

   } //end of  if (dc <...  ie. STILL in range

  } //end of for (i = 0;....

  for (n = 0; n < 640; n++) {  //inner loop - 640 units of time

   if (kbhit()) n = 641;  //any key escapes inner & outer loops

   for (m = 0; m < 100; m++) {  //innermost loop, 100 x ti(.01)

    for (i = 0; i < ns - 1; i++) {          //for each star ...

     if (dc[i] < maxd) {                 //...still in range...

      for (j = i+1; j < ns; j++) {      //...with all higher...

       if (dc[j] < maxd) {   //...numbered stars still in range

    xd=x[j]-x[i]; yd=y[j]-y[i]; zd=z[j]-z[i];// x/y/z diffs

    dsq[i][j]= xd*xd + yd*yd + zd*zd;//separation dist sqrd

    if(dsq[i][j] < .1) dsq[i][j] = .1; //stop excess accels

    d = sqrt(dsq[i][j]);//separation distance between i & j

    a[i][j] = g / dsq[i][j]; //gravity accel of i towards j

    ax[i][j] = a[i][j] * xd / d; //x comp of i towards j...

    ax[j][i] = -ax[i][j];//...and j to i (equal & opposite)

    ay[i][j] = a[i][j] * yd / d; //y comp of i towards j...

    ay[j][i] = -ay[i][j];//...and j to i (equal & opposite)

    az[i][j] = a[i][j] * zd / d; //z comp of i towards j...

    az[j][i] = -az[i][j];//...and j to i (equal & opposite)

       } // end of if (dc[j] < maxd)

      } //end for (j = i+1;...

     } // end of if (dc[i] < maxd)

    } //end for (i = 0;...

    for (i = 0; i < ns; i++) {      //now process every star...

     if (dc[i] < maxd) {                  //..still in range...

      atx[i] = 0; aty[i] = 0; atz[i] = 0;//zero the total accel

      for (j = 0; j < ns; j++) {     //add accel towards all...

       if (dc[j] < maxd) {   //...OTHER stars still in range...

    if (i != j ) { atx[i] += ax[i][j];     //cumul. x accel

               aty[i] += ay[i][j];     //cumul. y accel

               atz[i] += az[i][j];     //cumul. z accel

    } //end of if (i != j )  ie. NOT towards self

       } // end of if (dc[j] < maxd)   ie. OTHER star

      } //end for (j = 0; j < ns;

      vx[i] += atx[i] * ti;      //adjust x velocity comp. of i

      x[i] +=  vx[i] * ti;       //adjust x position comp. of i

      vy[i] +=  aty[i] * ti;     //adjust y velocity comp. of i

      y[i] += vy[i] * ti;        //adjust y position comp. of i

      vz[i] += atz[i] * ti;      //adjust z velocity comp. of i

      z[i] +=  vz[i] * ti;       //adjust z position comp. of i

     } // end of if (dc[i] < maxd)  ie. every star

    }  //end for (i = 0; i < ns...

   } //end of for (m = 1;... innermost loop 100 x ti

   for (i = 0; i < ns; i++) {       //now stereo plot each star

    if (dc[i] < maxd) {                   //..still in range...

     zf = eyez / (eyez - z[i]);      //z factor for perspective

     xr = eyexr +(x[i]+xm -eyexr) * zf;//screen x for right eye

     xl = eyexl +(x[i]+xm -eyexl) * zf; //screen x for left eye

     yb = eyey + (y[i]+ym -eyey) * zf; //screen y for both eyes

     if (xr > 320) putpixel(xr, yb, c[i]);//plot right eye view

     if (xl > 320) putpixel(640-xl, yb, c[i]);//left - mirrored

     vi[i] =sqrt(vx[i]*vx[i] +vy[i]*vy[i] +vz[i]*vz[i]);//veloc

     putpixel (n, 480 - vi[i] * 100, c[i]); //veloc versus time

    } // end of if (dc[i] < maxd)

   } //end for (i = 0; i < ns   ie. stars in range plotted

  } //end of for (n = 0;...    inner loop (640 x unit time)

 } while(n < 641); //end do {...  outer loop

 closegraph();//shut down graphics when escaped from outer loop

}; //end main()                                  //program ends



void setscreen(void){ //define procedure to get graphics screen

 int gdriver = VGA, gmode=VGAHI, errorcode;   //set VGA 640x480

 if(registerfarbgidriver(EGAVGA_driver_far)<0) exit(1);//driver

 initgraph(&gdriver, &gmode, "..//bgi");  //initialize graphics

 errorcode = graphresult(); //read result of initialization

 if (errorcode != grOk)          //if an error occurred

 {  printf("Graphics error; %s\n", grapherrormsg(errorcode));

    printf("Press any key to halt;");    getch();

    exit(1);   }               //return with error code

};  //end setscreen()





Click to download DSTARS6.ZIP for source listing

and executable version of this program.

Conclusion

That concludes the series; I hope it has illustrated to non-programmers what (some) programming looks like and has inveigled learners to try their hand at something different to their previous experience. 

Note: This series consists of Parts 1 to 6 in PC Update issues from December 2001 to June 2002.  Click to download DSTARS.ZIP file (225 KB) which contains all text, source and executable copies of programs used in this series.

Reprinted from the June 2002 issue of PC Update, the magazine of Melbourne PC User Group, Australia

[About Melbourne PC User Group]