Wed May 18 08:58:33 PDT 2016

Galaxy Simulation Program

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.

Comments are closed


If you would like to get in touch with me, please mail zfs at themolecularuniverse.com

recent comments

Posted by ZFS | Permanent link | File under: chemistry, general
[StumbleUpon] [Digg] [Reddit] [Facebook] [Google]