21 #define mini(a,b) (((a)<(b))?(a):(b))
23 typedef struct Particle Particle;
32 typedef struct Path Path;
68 int j, grid = sqrt(N)+0.5;
70 draw(screen, screen->r, display->white, 0, ZP);
73 p->x = 2*(j%grid)+frand()/2;
74 p->y = 2*(j/grid)+frand()/2;
75 p->vx = 1.*v0*frand();
76 p->vy = 1.*v0*frand();
77 p->prevx = p->x - p->vx * dt;
78 p->prevy = p->y - p->vy * dt;
79 p->col = allocimage(display, Rect(0,0,1,1), screen->chan, 1, colors[rand()%(sizeof(colors)/sizeof(int))]);
80 if(!p->col) sysfatal("allocimage");
91 draw(screen, screen->r, display->white, 0, ZP);
94 memset(pa->x, 0, sizeof(int) * pathlen);
95 memset(pa->y, 0, sizeof(int) * pathlen);
108 drawpath(Path *p, Image *col, int i)
112 if((j = i+1) == pathlen)
114 draw(screen, Rect(p->x[i], p->y[i], p->x[i]+1, p->y[i]+1), col, 0, ZP);
117 draw(screen, Rect(p->x[j], p->y[j], p->x[j]+1, p->y[j]+1), display->white, 0, ZP);
123 print("USAGE: mole options\n");
124 print(" -N number of particles [49]\n");
125 print(" -x left boundary [-40]\n");
126 print(" -X right boundary [40]\n");
127 print(" -y top boundary [-40]\n");
128 print(" -Y bottom boundary [40]\n");
129 print(" -t time step [0.01]\n");
130 print(" -v maximum start velocity [0.1]\n");
131 print(" -P path length [1000]\n");
136 main(int argc, char** argv)
141 double dx, dx1, dx2, dy, dy1, dy2, R, F;
144 #define FARG(c, v) case c: if(!(f=ARGF())) usage(); v = atof(f); break;
146 case 'N': if(!(f=ARGF())) usage(); N = atoi(f); break;
147 case 'P': if(!(f=ARGF())) usage(); pathlen = atoi(f); break;
161 A = calloc(sizeof(Particle), N);
162 B = calloc(sizeof(Particle), N);
163 paths = calloc(sizeof(Path), N);
164 for(pa=paths;pa<paths+N;pa++){
165 pa->x = calloc(sizeof(int), pathlen);
166 pa->y = calloc(sizeof(int), pathlen);
171 initdraw(0, 0, "Molecular Dynamics");
172 einit(Emouse | Ekeyboard);
178 memset(cur, 0, sizeof(Particle) * N);
179 for(p=prev;p<prev+N;p++) {
180 for(q=prev;q<p;q++) {
181 dx1 = fabs(p->x - q->x);
182 dx2 = xmax - xmin - dx1;
184 dy1 = fabs(p->y - q->y);
185 dy2 = ymax - ymin - dy1;
190 double R2, R4, R6, R12;
196 if(p->x < q->x) dx = -dx;
197 if(p->y < q->y) dy = -dy;
198 if(dx1 > dx2) dx = -dx;
199 if(dy1 > dy2) dy = -dy;
202 (p-prev+cur)->ax += dx;
203 (p-prev+cur)->ay += dy;
204 (q-prev+cur)->ax -= dx;
205 (q-prev+cur)->ay -= dy;
212 q->x = 2*p->x - p->prevx + q->ax * dt*dt;
213 q->y = 2*p->y - p->prevy + q->ay * dt*dt;
214 q->vx = (q->x - p->prevx) / (2*dt);
215 q->vy = (q->y - p->prevy) / (2*dt);
218 if(q->x > xmax) {q->x -= xmax - xmin; q->prevx -= xmax - xmin;}
219 if(q->x < xmin) {q->x += xmax - xmin; q->prevx += xmax - xmin;}
220 if(q->y > ymax) {q->y -= ymax - ymin; q->prevy -= ymax - ymin;}
221 if(q->y < ymin) {q->y += ymax - ymin; q->prevy += ymax - ymin;}
223 pa->x[i] = (screen->r.max.x - screen->r.min.x) * (q->x - xmin) / (xmax - xmin) + screen->r.min.x;
224 pa->y[i] = (screen->r.max.y - screen->r.min.y) * (q->y - ymin) / (ymax - ymin) + screen->r.min.y;
225 drawpath(pa, p->col, i);
228 Particle* tmp = prev;
231 flushimage(display, 1);
236 case 'q': case Kdel: exits(0); break;
237 case 'r': reset(); break;
238 case 'R': reverse(); break;
239 case 'f': draw(screen, screen->r, display->white, 0, ZP); break;
248 if(new) getwindow(display, Refnone);