<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
"http://www.w3.org/TR/html4/loose.dtd">
<html>
<head>
<meta http-equiv=Content-Type content="text/html; charset=utf8">
<title>/usr/web/sources/contrib/aiju/mole.c - Plan 9 from Bell Labs</title>
<!-- THIS FILE IS AUTOMATICALLY GENERATED. -->
<!-- EDIT sources.tr INSTEAD. -->
</meta>
</head>
<body>
<p style="margin-top: 0; margin-bottom: 0.17in"></p>
<p style="line-height: 1.2em; margin-left: 1.00in; text-indent: 0.00in; margin-right: 1.00in; margin-top: 0; margin-bottom: 0; text-align: center;">
<span style="font-size: 10pt"><a href="/plan9/">Plan 9 from Bell Labs</a>’s /usr/web/sources/contrib/aiju/mole.c</span></p>
<p style="margin-top: 0; margin-bottom: 0.17in"></p>
<p style="margin-top: 0; margin-bottom: 0.17in"></p>
<center><font size=-1>
Copyright © 2009 Alcatel-Lucent.<br />
Distributed under the
<a href="/plan9/license.html">Lucent Public License version 1.02</a>.
<br />
<a href="/plan9/download.html">Download the Plan 9 distribution.</a>
</font>
</center>
<p style="margin-top: 0; margin-bottom: 0.17in"></p>
<table width="100%" cellspacing=0 border=0><tr><td align="center">
<table cellspacing=0 cellpadding=5 bgcolor="#eeeeff"><tr><td align="left">
<pre>
<!-- END HEADER -->
#include <u.h>
#include <libc.h>
#include <draw.h>
#include <event.h>
int N = 20, refresh = 0;
double dt = 0.01,
xmin = -40,
xmax = 40,
ymin = -40,
ymax = 40,
v0 = 0.1;
#define mini(a,b) (((a)<(b))?(a):(b))
typedef struct Particle Particle;
struct Particle {
double x, y;
double vx, vy;
double ax, ay;
double prevx, prevy;
Image* col;
};
int colors[] = {
DBlack,
DRed,
DGreen,
DBlue,
DCyan,
DMagenta,
DDarkyellow,
DDarkgreen,
DPalegreen,
DMedgreen,
DDarkblue,
DPalebluegreen,
DPaleblue,
DBluegreen,
DGreygreen,
DPalegreygreen,
DYellowgreen,
DMedblue,
DGreyblue,
DPalegreyblue,
DPurpleblue
};
Particle *A, *B;
Particle *prev, *cur;
void
reset(void)
{
int j, grid = sqrt(N)+0.5;
Particle *p;
draw(screen, screen->r, display->white, 0, ZP);
for(j=0;j<N;j++) {
p = prev+j;
p->x = 2*(j%grid)+frand()/2;
p->y = 2*(j/grid)+frand()/2;
p->vx = 1.*v0*frand();
p->vy = 1.*v0*frand();
p->prevx = p->x - p->vx * dt;
p->prevy = p->y - p->vy * dt;
p->col = allocimage(display, Rect(0,0,1,1), screen->chan, 1, colors[rand()%(sizeof(colors)/sizeof(int))]);
if(!p->col) sysfatal("allocimage");
}
}
void
usage(void)
{
print("USAGE: mole options\n");
print(" -N number of particles [20]\n");
print(" -x left boundary [-40]\n");
print(" -X right boundary [40]\n");
print(" -y top boundary [-40]\n");
print(" -Y bottom boundary [40]\n");
print(" -t time step [0.01]\n");
print(" -v maximum start velocity [0.1]\n");
print(" -F clear every <n> frames (0:off) [0]\n");
exits("usage");
}
void
main(int argc, char** argv)
{
int i, j;
Particle *p, *q;
double dx, dx1, dx2, dy, dy1, dy2, R, F;
char* f;
#define FARG(c, v) case c: if(!(f=ARGF())) usage(); v = atof(f); break;
ARGBEGIN {
case 'N': if(!(f=ARGF())) usage(); N = atoi(f); break;
case 'F': if(!(f=ARGF())) usage(); refresh = atoi(f); break;
FARG('v', v0);
FARG('x', xmin);
FARG('X', xmax);
FARG('y', ymin);
FARG('Y', ymax);
FARG('t', dt);
default: usage();
} ARGEND;
A = calloc(sizeof(Particle), N);
B = calloc(sizeof(Particle), N);
prev = A;
cur = B;
srand(time(0));
initdraw(0, 0, "Molecular Dynamics");
einit(Emouse | Ekeyboard);
reset();
i=0;
while(1) {
if(refresh && ((++i)%refresh)==0) draw(screen, screen->r, display->white, 0, ZP);
memset(cur, 0, sizeof(Particle) * N);
for(p=prev;p<prev+N;p++) {
for(q=prev;q<p;q++) {
dx1 = fabs(p->x - q->x);
dx2 = xmax - xmin - dx1;
dx = mini(dx1, dx2);
dy1 = fabs(p->y - q->y);
dy2 = ymax - ymin - dy1;
dy = mini(dy1, dy2);
R = dx*dx + dy*dy;
if(R >= 9) continue;
R = 1/sqrt(R);
double R2, R4, R6, R12;
R2 = R * R;
R4 = R2 * R2;
R6 = R4 * R2;
R12 = R6 * R6;
F = 24*(2*R12 - R6);
if(p->x < q->x) dx = -dx;
if(p->y < q->y) dy = -dy;
if(dx1 > dx2) dx = -dx;
if(dy1 > dy2) dy = -dy;
dx *= F;
dy *= F;
(p-prev+cur)->ax += dx;
(p-prev+cur)->ay += dy;
(q-prev+cur)->ax -= dx;
(q-prev+cur)->ay -= dy;
}
}
for(j=0;j<N;j++) {
int x, y;
p = prev+j;
q = cur+j;
q->x = 2*p->x - p->prevx + q->ax * dt*dt;
q->y = 2*p->y - p->prevy + q->ay * dt*dt;
q->vx = (q->x - p->prevx) / (2*dt);
q->vy = (q->y - p->prevy) / (2*dt);
q->prevx = p->x;
q->prevy = p->y;
if(q->x > xmax) {q->x -= xmax - xmin; q->prevx -= xmax - xmin;}
if(q->x < xmin) {q->x += xmax - xmin; q->prevx += xmax - xmin;}
if(q->y > ymax) {q->y -= ymax - ymin; q->prevy -= ymax - ymin;}
if(q->y < ymin) {q->y += ymax - ymin; q->prevy += ymax - ymin;}
q->col = p->col;
x = (screen->r.max.x - screen->r.min.x) * (q->x - xmin) / (xmax - xmin) + screen->r.min.x;
y = (screen->r.max.y - screen->r.min.y) * (q->y - ymin) / (ymax - ymin) + screen->r.min.y;
draw(screen, Rect(x, y, x+1, y+1), p->col, 0, ZP);
}
Particle* tmp = prev;
prev = cur;
cur = tmp;
flushimage(display, 1);
if(ecankbd()) {
switch(ekbd()) {
case 'q': exits(0); break;
case 'r': reset(); break;
case 'f': draw(screen, screen->r, display->white, 0, ZP); break;
}
}
}
}
void
eresized(int new)
{
if(new) getwindow(display, Refnone);
}
<!-- BEGIN TAIL -->
</pre>
</td></tr></table>
</td></tr></table>
<p style="margin-top: 0; margin-bottom: 0.17in"></p>
<p style="line-height: 1.2em; margin-left: 1.00in; text-indent: 0.00in; margin-right: 1.00in; margin-top: 0; margin-bottom: 0; text-align: center;">
<span style="font-size: 10pt"></span></p>
<p style="margin-top: 0; margin-bottom: 0.50in"></p>
<p style="margin-top: 0; margin-bottom: 0.33in"></p>
<center><table border="0"><tr>
<td valign="middle"><a href="http://www.alcatel-lucent.com/"><img border="0" src="/plan9/img/logo_ft.gif" alt="Bell Labs" />
</a></td>
<td valign="middle"><a href="http://www.opensource.org"><img border="0" alt="OSI certified" src="/plan9/img/osi-certified-60x50.gif" />
</a></td>
<td><img style="padding-right: 45px;" alt="Powered by Plan 9" src="/plan9/img/power36.gif" />
</td>
</tr></table></center>
<p style="margin-top: 0; margin-bottom: 0.17in"></p>
<center>
<span style="font-size: 10pt">(<a href="/plan9/">Return to Plan 9 Home Page</a>)</span>
</center>
<p style="margin-top: 0; margin-bottom: 0.17in"></p>
<center><font size=-1>
<span style="font-size: 10pt"><a href="http://www.lucent.com/copyright.html">Copyright</a></span>
<span style="font-size: 10pt">© 2009 Alcatel-Lucent.</span>
<span style="font-size: 10pt">All Rights Reserved.</span>
<br />
<span style="font-size: 10pt">Comments to</span>
<span style="font-size: 10pt"><a href="mailto:webmaster@plan9.bell-labs.com">webmaster@plan9.bell-labs.com</a>.</span>
</font></center>
</body>
</html>