Maxwell's Demon



Below are gifs generated with the included C code as various frame speeds and across multiple starting conditions:
20 particles 200 particles 2000 particles

Code Below:

#include
#include
#include
#include

// gcc max.c -o max
// ./max
// ffmpeg -framerate 20 -i %04d.pgm max.gif

#define MIN(X, Y) (((X) < (Y)) ? (X) : (Y))
#define MAX(X, Y) (((X) > (Y)) ? (X) : (Y))
#define PI 3.14159265359

typedef enum {FALSE, TRUE} bool;

typedef struct particle {
//location
double x;
double y;

//direction
double vel_x;
double vel_y;

//a fast or slow particle
bool fast;

} particle_t;

#define FRAMESCALE 2

#define NUMPARTICLES 200//num of particles
#define NUMFRAMES 2000 //number of frames, 20 frames a second
#define WIDTH 200
#define HEIGHT 100

#define MAXSPEED 4.5
#define MINSPEED 0.5
#define FASTCUTOFF 3.0 //pixel/sec

#define GATEX WIDTH/2
#define GATEHOLE HEIGHT/4
#define GATEWIDTH 3

#define FASTPIXELWIDTH 2
#define FASTGREY 200
#define SLOWPIXELWIDTH 5
#define SLOWGREY 100

#define TIMESTEP 1.2//how many seconds


unsigned char pixel[WIDTH][HEIGHT];
particle_t particles[NUMPARTICLES];
bool gate_closed = FALSE;

char filename[100];
FILE *f;


void init_particles(){
//initialize locations and motions of particles
for (int i = 0; i < NUMPARTICLES; i++){
double randx = (double)rand() * WIDTH / RAND_MAX;
double randy = (double)rand() * HEIGHT / RAND_MAX;

double randspeed = (double)rand() * (MAXSPEED - MINSPEED) / RAND_MAX + MINSPEED;
double randangle = (double)rand() * 2.0*PI / RAND_MAX;

bool fast = (randspeed >= FASTCUTOFF);

particles[i].x = randx;
particles[i].y = randy;

particles[i].vel_x = randspeed * cos(randangle);
particles[i].vel_y = randspeed * sin(randangle);

particles[i].fast = fast;

}

}


void print_particles(){
for (int i = 0; i < 3; i++){
printf("particle %d: %f,%f moving %f, %f fast: %d\n",
i, particles[i].x, particles[i].y, particles[i].vel_x, particles[i].vel_y, particles[i].fast);
}
}


void update_particles(){

for (int i = 0; i < NUMPARTICLES; i++){

//predict location
double x = particles[i].x;
double y = particles[i].y;
double upx = x + particles[i].vel_x*TIMESTEP;
double upy = y + particles[i].vel_y*TIMESTEP;

//check for outer box
if (upx <= 0) particles[i].vel_x *= -1;
if (upx >= WIDTH) particles[i].vel_x *= -1;

if (upy <= 0) particles[i].vel_y *= -1;
if (upy >= HEIGHT) particles[i].vel_y *= -1;

//check for gate
if ((x < (GATEX - GATEWIDTH) && upx >= (GATEX - GATEWIDTH)) || \
(x > (GATEX + GATEWIDTH) && upx <= (GATEX + GATEWIDTH))) {

if (gate_closed || upy < (HEIGHT/2 - GATEHOLE/2) || upy > (HEIGHT/2 + GATEHOLE/2)){
particles[i].vel_x *= -1;
}
}

//update
particles[i].x += particles[i].vel_x*TIMESTEP;
particles[i].y += particles[i].vel_y*TIMESTEP;
}

}


void draw_gatewall(){
for (int i=0; i if (gate_closed || i < (HEIGHT/2 - GATEHOLE/2) || i > (HEIGHT/2 + GATEHOLE/2)){
pixel[GATEX][i] = 0;
for (int j=1; j pixel[GATEX+j][i] = 0;
pixel[GATEX-j][i] = 0;
}
}
}
}


void demon_choice(){
//predict those that will cross and close/open based on best decision
int to_close = 0;
int or_not_to_close = 0;

for (int i = 0; i < NUMPARTICLES; i++){


//predict location
double x = particles[i].x;
double y = particles[i].y;
double upx = x + particles[i].vel_x*TIMESTEP;
double upy = y + particles[i].vel_y*TIMESTEP;

//check for crossing gate one way
if (x < (GATEX - GATEWIDTH) && upx >= (GATEX - GATEWIDTH)) {
if (upy > (HEIGHT/2 - GATEHOLE/2) && upy < (HEIGHT/2 + GATEHOLE/2)){
if (particles[i].fast){
to_close += 1;
} else {
or_not_to_close += 1;
}
}
}
//check the other way
else if (x > (GATEX + GATEWIDTH) && upx <= (GATEX + GATEWIDTH)) {
if (upy > (HEIGHT/2 - GATEHOLE/2) && upy < (HEIGHT/2 + GATEHOLE/2)){
if (particles[i].fast){
or_not_to_close += 1;
} else {
to_close += 1;
}
}
}
}

//open or close the gate!
if (to_close > or_not_to_close) gate_closed = TRUE;
if (to_close < or_not_to_close) gate_closed = FALSE;

}


void reset_pixels(){
//reset pixels in image to white
for (int i = 0; i < WIDTH; i++){
for (int j = 0; j < HEIGHT; j++){
pixel[i][j] = 255;
}
}
draw_gatewall();
}


void draw_circle(int x_center, int y_center, int radius, int pixelintensity){
for (int y = y_center - radius; y <= y_center + radius; y++){

double x_bound = sqrt( (radius*radius) - ((y-y_center)*(y-y_center)) );

for (int x = x_center - radius; x <= x_center + radius; x++){
if ( abs(x-x_center) < x_bound && 0 <= x && x < WIDTH && 0 <= y && y < HEIGHT ){
if ((int)pixel[x][y] < pixelintensity) pixel[x][y] = 0;
else pixel[x][y] -= pixelintensity;
}
}
}
}



void draw_particles(){
reset_pixels();

for (int i = 0; i < NUMPARTICLES; i++){

if (particles[i].fast){
draw_circle(particles[i].x, particles[i].y, FASTPIXELWIDTH, FASTGREY);
} else {
draw_circle(particles[i].x, particles[i].y, SLOWPIXELWIDTH, SLOWGREY);
}

}
}

void write_pixels(){
//write pixels to file
static int t = 1;

sprintf(filename, "%04d.pgm", t);
fprintf(stderr, "Writing %s...", filename);

f = fopen(filename, "w");
fprintf(f, "P5\n%d %d\n255\n", HEIGHT, WIDTH);

fwrite(pixel, 1, WIDTH*HEIGHT, f);
fclose(f);
fprintf(stderr, "DONE\n");

t++;
}

int main(){
srand(time(NULL)); //seed random num generator

init_particles();
print_particles();

draw_particles();
write_pixels();

for (int frames = 0; frames < NUMFRAMES*FRAMESCALE; frames++){
demon_choice();
update_particles();
if (frames%FRAMESCALE == 0){
draw_particles();
write_pixels();
}
}
}