sys/src/games/mole.c

   1 #include <u.h>
   2 #include <libc.h>
   3 #include <draw.h>
   4 #include <event.h>
   5
   6 int N = 20, refresh = 0;
   7
   8 double dt = 0.01,
   9         xmin = -40,
  10         xmax = 40,
  11         ymin = -40,
  12         ymax = 40,
  13         v0 = 0.1;
  14
  15 #define mini(a,b) (((a)<(b))?(a):(b))
  16
  17 typedef struct Particle Particle;
  18 struct Particle {
  19         double x, y;
  20         double vx, vy;
  21         double ax, ay;
  22         double prevx, prevy;
  23         Image* col;
  24 };
  25
  26 int colors[] = {
  27               DBlack,
  28               DRed,
  29               DGreen,
  30               DBlue,
  31               DCyan,
  32               DMagenta,
  33               DDarkyellow,
  34               DDarkgreen,
  35               DPalegreen,
  36               DMedgreen,
  37               DDarkblue,
  38               DPalebluegreen,
  39               DPaleblue,
  40               DBluegreen,
  41               DGreygreen,
  42               DPalegreygreen,
  43               DYellowgreen,
  44               DMedblue,
  45               DGreyblue,
  46               DPalegreyblue,
  47               DPurpleblue
  48 };
  49
  50 Particle *A, *B;
  51 Particle *prev, *cur;
  52
  53 void
  54 reset(void)
  55 {
  56         int j, grid = sqrt(N)+0.5;
  57         Particle *p;
  58         draw(screen, screen->r, display->white, 0, ZP);
  59         for(j=0;j<N;j++) {
  60                 p = prev+j;
  61                 p->x = 2*(j%grid)+frand()/2;
  62                 p->y = 2*(j/grid)+frand()/2;
  63                 p->vx = 1.*v0*frand();
  64                 p->vy = 1.*v0*frand();
  65                 p->prevx = p->x - p->vx * dt;
  66                 p->prevy = p->y - p->vy * dt;
  67                 p->col = allocimage(display, Rect(0,0,1,1), screen->chan, 1, colors[rand()%(sizeof(colors)/sizeof(int))]);
  68                 if(!p->col) sysfatal("allocimage");
  69         }
  70 }
  71
  72 void
  73 usage(void)
  74 {
  75         print("USAGE: mole options\n");
  76         print(" -N number of particles [20]\n");
  77         print(" -x left boundary [-40]\n");
  78         print(" -X right boundary [40]\n");
  79         print(" -y top boundary [-40]\n");
  80         print(" -Y bottom boundary [40]\n");
  81         print(" -t time step [0.01]\n");
  82         print(" -v maximum start velocity [0.1]\n");
  83         print(" -F clear every <n> frames (0:off) [0]\n");
  84         exits("usage");
  85 }
  86
  87 void
  88 main(int argc, char** argv)
  89 {
  90         int i, j;
  91         Particle *p, *q;
  92         double dx, dx1, dx2, dy, dy1, dy2, R, F;
  93         char* f;
  94
  95         #define FARG(c, v) case c: if(!(f=ARGF())) usage(); v = atof(f); break;
  96         ARGBEGIN {
  97         case 'N': if(!(f=ARGF())) usage(); N = atoi(f); break;
  98         case 'F': if(!(f=ARGF())) usage(); refresh = atoi(f); break;
  99         FARG('v', v0);
 100         FARG('x', xmin);
 101         FARG('X', xmax);
 102         FARG('y', ymin);
 103         FARG('Y', ymax);
 104         FARG('t', dt);
 105         default: usage();
 106         } ARGEND;
 107
 108         A = calloc(sizeof(Particle), N);
 109         B = calloc(sizeof(Particle), N);
 110         prev = A;
 111         cur = B;
 112         srand(time(0));
 113         initdraw(0, 0, "Molecular Dynamics");
 114         einit(Emouse | Ekeyboard);
 115         reset();
 116
 117         i=0;
 118         while(1) {
 119                 if(refresh && ((++i)%refresh)==0) draw(screen, screen->r, display->white, 0, ZP);
 120                 memset(cur, 0, sizeof(Particle) * N);
 121                 for(p=prev;p<prev+N;p++) {
 122                         for(q=prev;q<p;q++) {
 123                                 dx1 = fabs(p->x - q->x);
 124                                 dx2 = xmax - xmin - dx1;
 125                                 dx = mini(dx1, dx2);
 126                                 dy1 = fabs(p->y - q->y);
 127                                 dy2 = ymax - ymin - dy1;
 128                                 dy = mini(dy1, dy2);
 129                                 R = dx*dx + dy*dy;
 130                                 if(R >= 9) continue;
 131                                 R = 1/sqrt(R);
 132                                 double R2, R4, R6, R12;
 133                                 R2 = R * R;
 134                                 R4 = R2 * R2;
 135                                 R6 = R4 * R2;
 136                                 R12 = R6 * R6;
 137                                 F = 24*(2*R12 - R6);
 138                                 if(p->x < q->x) dx = -dx;
 139                                 if(p->y < q->y) dy = -dy;
 140                                 if(dx1 > dx2) dx = -dx;
 141                                 if(dy1 > dy2) dy = -dy;
 142                                 dx *= F;
 143                                 dy *= F;
 144                                 (p-prev+cur)->ax += dx;
 145                                 (p-prev+cur)->ay += dy;
 146                                 (q-prev+cur)->ax -= dx;
 147                                 (q-prev+cur)->ay -= dy;
 148                         }
 149                 }
 150                 for(j=0;j<N;j++) {
 151                         int x, y;
 152                         p = prev+j;
 153                         q = cur+j;
 154                         q->x = 2*p->x - p->prevx + q->ax * dt*dt;
 155                         q->y = 2*p->y - p->prevy + q->ay * dt*dt;
 156                         q->vx = (q->x - p->prevx) / (2*dt);
 157                         q->vy = (q->y - p->prevy) / (2*dt);
 158                         q->prevx = p->x;
 159                         q->prevy = p->y;
 160                         if(q->x > xmax) {q->x -= xmax - xmin; q->prevx -= xmax - xmin;}
 161                         if(q->x < xmin) {q->x += xmax - xmin; q->prevx += xmax - xmin;}
 162                         if(q->y > ymax) {q->y -= ymax - ymin; q->prevy -= ymax - ymin;}
 163                         if(q->y < ymin) {q->y += ymax - ymin; q->prevy += ymax - ymin;}
 164                         q->col = p->col;
 165                         x = (screen->r.max.x - screen->r.min.x) * (q->x - xmin) / (xmax - xmin) + screen->r.min.x;
 166                         y = (screen->r.max.y - screen->r.min.y) * (q->y - ymin) / (ymax - ymin) + screen->r.min.y;
 167                         draw(screen, Rect(x, y, x+1, y+1), p->col, 0, ZP);
 168                 }
 169
 170                 Particle* tmp = prev;
 171                 prev = cur;
 172                 cur = tmp;
 173                 flushimage(display, 1);
 174
 175
 176                 if(ecankbd()) {
 177                         switch(ekbd()) {
 178                                 case 'q': exits(0); break;
 179                                 case 'r': reset(); break;
 180                                 case 'f': draw(screen, screen->r, display->white, 0, ZP); break;
 181                         }
 182                 }
 183         }
 184 }
 185
 186 void
 187 eresized(int new)
 188 {
 189         if(new) getwindow(display, Refnone);
 190 }