Here is the version of the galaxy program used to generate the movies of the last few days.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// PROGRAM galaxy
// model proposed by Schulman and Seiden
//
// ZFS - changed so that graphical output is handled in postscript
// included a minor fix to avoid stomping on memory and setting the time
// to 15, and increased the number of comments
//
// Usage
//
// Compile with:
// gcc galaxy.c
//
// Run with:
// ./a.out
//
// You will generate a collection of postscript files which show the
// evolution of your model galaxy
int main();
void parameter();
void initial();
void evolve();
void create(int r, int a, int *newactive, int newactive_r[], int newactive_a[]);
void plot_spiral();
float rand1();
// if printout is set to a positive integer this value is used to name
// the output postscript file
int printout=0;
int active_a[10001];
int active_r[10001];
int cell[51][301];
int dt;
int nactive;
int nring;
double p;
double s;
int t;
double two_pi;
double v;
double pi;
int main()
{
long int i;
parameter();
initial();
for(i=1;i<=1000;i++){
evolve();
printout=i;
plot_spiral();
t += dt;
}
return 0;
}
void parameter()
{
pi=3.141593;
nring = 50; // number of rings
nactive = 200; // number of initial active cells
v = 1; // circular velocity
p = 0.18; // star formation probability
dt = 10; // time step
s = 2*pi/6; // cell width
two_pi = 2*pi;
}
void initial()
{
int a, i, r, x, y, Itmp1_, Itmp2_;
float xwin, ywin;
float xtmp, ytmp;
for(Itmp1_ = 0; Itmp1_ <= 50; ++Itmp1_)
for(Itmp2_ = 0; Itmp2_ <= 300; ++Itmp2_)
cell[Itmp1_][Itmp2_] = 0;
// randomly activate cells
// this is done by choosing a ring by choosing points evenly on a grid
// and determining which ring the points lie in (to avoid overly sampling
// the small rings. The angle (a) about the ring can be chosen evenly for
// each ring.
i = 0;
do
{
do
{
x = (int)(nring*rand1()) + 1;
y = (int)(nring*rand1()) + 1;
xtmp=(float)x;
ytmp=(float)y;
r = (int)sqrt(xtmp*xtmp+ ytmp*ytmp) + 1;
} while(!(r <= nring));
a = (int)(6*r*rand1()) + 1; // array index corresponds to an angle
if(cell[r][a] == 0)
{
++i;
active_a[i] = a; // location of active region
active_r[i] = r;
// activate region, stars live for 15 time steps
cell[r][a] = 15;
}
} while(!(i == nactive));
t = 0; // initial time
}
void evolve()
{
int a, am, ap, i, newactive, newactive_a[10001], newactive_r[10001], r,
Itmp1_;
double angle, wt;
// number of active star clusters for next time step
newactive = 0;
for(i = 1; i <= nactive; ++i)
{
r = active_r[i];
a = active_a[i];
create(r, a+1, &newactive, newactive_r, newactive_a);
create(r, a-1, &newactive, newactive_r, newactive_a);
// angle is the angle of the current region, the angles of the next larger
// ring adjacent clusters are ap and next smaller ring clusters are am
//
// these angles are computed based on the ring velocity (which depends on the
// ring in question), the simulation time, and the overall velocity setting
angle = fmod((double)(a*s + v*t)/r, two_pi);
if(r < nring) // activate cells in next larger ring
{
wt = fmod((double)v*t/(r+1), two_pi);
ap = (int)((angle - wt)*(r+1)/s);
ap = ((ap) % (6*(r+1)));
create(r+1, ap, &newactive, newactive_r, newactive_a);
create(r+1, ap+1, &newactive, newactive_r, newactive_a);
}
if(r > 1) // activate cells in next smaller ring
{
wt = fmod((double)v*t/(r-1), two_pi);
am = (int)((angle - wt)*(r-1)/s);
am = ((am) % (6*(r-1)));
create(r-1, am, &newactive, newactive_r, newactive_a);
create(r-1, am+1, &newactive, newactive_r, newactive_a);
}
}
// copy the new actives into the actives lists - ready for the
// next time step
nactive = newactive;
for(Itmp1_ = 0; Itmp1_ <= 10000; ++Itmp1_) {
active_r[Itmp1_] = newactive_r[Itmp1_];
active_a[Itmp1_] = newactive_a[Itmp1_];
}
}
void create(int r, int a, int *newactive, int newactive_r[], int newactive_a[])
// create star clusters
{
if(a < 1)
a += 6*r;
// The original program had a bug in that a could be 301 on occasion (for 50th
// ring) and this then set the simulation time to '15' by reaching beyond the
// end of the cell matrix...this problem is avoided by the following test which
// imposes a circular boundary on the angle in the positive a direction (as
// was originally done in the negative a direction
if(a > 300)
a -= 6*r;
// Now see if a supplied region should become active because of its
// currenly proximal active region - this happens with a probability 'p'
if(rand1() < p && cell[r][a] != 15)
{
++(*newactive);
newactive_a[(*newactive)] = a;
newactive_r[(*newactive)] = r;
cell[r][a] = 15; // activate cell
}
}
void plot_spiral()
{
int a, r;
double plotsize, theta, x, y;
double theta2;
double scale=4.0;
double xoff=320.0;
double yoff=475.0;
char filename[20];
FILE * fd;
// Here we output postscript disks with radii which correspond to the
// age of a given star cluster. The location of a given cluster is
// determined based on the current simulation time and the velocity of a
// given ring.
sprintf(filename, "%07d.ps",printout);
if(printout){
fd=fopen(filename, "w");
fprintf(fd,"%%!PS\n");
fprintf(fd,"%%%%BoundingBox: 0 0 640 892\n");
fprintf(fd,"newpath\n");
fprintf(fd,"0.1 setlinewidth\n");
fprintf(fd,"0.0 0.0 0.0 setrgbcolor\n");
}
for(r = 1; r <= nring; ++r)
{
for(a = 1; a <= 6*r; ++a)
{
if(cell[r][a] > 0)
{
theta = (double)(a*s + (double)v*(double)t)/(double)r;
theta2 = fmod(theta, two_pi);
theta = theta2;
x = r*cos(theta);
y = r*sin(theta);
plotsize = cell[r][a]/15.0;
if(printout){
fprintf(fd, "newpath\n");
fprintf(fd, "%f %f\n", x*scale+xoff,y*scale+yoff);
fprintf(fd, "%f 0 360 arc closepath\n", 1.0*plotsize*scale);
fprintf(fd, "0.0 setgray fill\n");
fprintf(fd, "stroke\n");
}
// reduce star clusters lifetime
--(cell[r][a]);
}
}
}
if(printout){
fprintf(fd,"showpage\n");
fclose(fd);
}
}
// the code assumed that rand() returned a real number between 0 and 1,
// hence the following function which does this.
float rand1()
{
return((float)rand()/(float)RAND_MAX);
}
I found the original source code from here ... which comes from this compendium of code here from the book by Harvey Gould and Jan Tobochnik, 'Introduction to Computer Simulation Methods'. I have the first edition of this book (it has two volumes) and it is a really fantastic resource for those of us that did not pay enough attention in school but are interested in learning how physics operates.