sys/src/games/mole.c

   1 #include <u.h>
   2 #include <libc.h>
   3 #include <draw.h>
   4 #include <event.h>
   5
   6 enum {
   7         Kdel = 0x7f
   8 };
   9
  10 int N = 49,
  11         pathlen = 1000,
  12         nosnake;
  13
  14 double dt = 0.01,
  15         xmin = -40,
  16         xmax = 40,
  17         ymin = -40,
  18         ymax = 40,
  19         v0 = 0.1;
  20
  21 #define mini(a,b) (((a)<(b))?(a):(b))
  22
  23 typedef struct Particle Particle;
  24 struct Particle {
  25         double x, y;
  26         double vx, vy;
  27         double ax, ay;
  28         double prevx, prevy;
  29         Image* col;
  30 };
  31
  32 typedef struct Path Path;
  33 struct Path {
  34         int *x, *y;
  35 };
  36
  37 int colors[] = {
  38               DBlack,
  39               DRed,
  40               DGreen,
  41               DBlue,
  42               DCyan,
  43               DMagenta,
  44               DDarkyellow,
  45               DDarkgreen,
  46               DPalegreen,
  47               DMedgreen,
  48               DDarkblue,
  49               DPalebluegreen,
  50               DPaleblue,
  51               DBluegreen,
  52               DGreygreen,
  53               DPalegreygreen,
  54               DYellowgreen,
  55               DMedblue,
  56               DGreyblue,
  57               DPalegreyblue,
  58               DPurpleblue
  59 };
  60
  61 Particle        *A, *B;
  62 Particle        *prev, *cur;
  63 Path            *paths;
  64
  65 void
  66 reset(void)
  67 {
  68         int j, grid = sqrt(N)+0.5;
  69         Particle *p;
  70         draw(screen, screen->r, display->white, 0, ZP);
  71         for(j=0;j<N;j++) {
  72                 p = prev+j;
  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");
  81         }
  82 }
  83
  84 void
  85 reverse(void)
  86 {
  87         Particle *p, *q;
  88         Path *pa;
  89         int i;
  90
  91         draw(screen, screen->r, display->white, 0, ZP);
  92         for(i=0;i<N;i++){
  93                 pa=paths+i;
  94                 memset(pa->x, 0, sizeof(int) * pathlen);
  95                 memset(pa->y, 0, sizeof(int) * pathlen);
  96                 p=prev+i;
  97                 q=cur+i;
  98                 p->vx = -q->vx;
  99                 p->vy = -q->vy;
 100                 p->prevx = p->x;
 101                 p->prevy = p->y;
 102                 p->x = q->x;
 103                 p->y = q->y;
 104         }
 105 }
 106
 107 void
 108 drawpath(Path *p, Image *col, int i)
 109 {
 110         int j;
 111
 112         if((j = i+1) == pathlen)
 113                 j = 0;
 114         draw(screen, Rect(p->x[i], p->y[i], p->x[i]+1, p->y[i]+1), col, 0, ZP);
 115         if(nosnake)
 116                 return;
 117         draw(screen, Rect(p->x[j], p->y[j], p->x[j]+1, p->y[j]+1), display->white, 0, ZP);
 118 }
 119
 120 void
 121 usage(void)
 122 {
 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");
 132         exits("usage");
 133 }
 134
 135 void
 136 main(int argc, char** argv)
 137 {
 138         int i, j;
 139         Particle *p, *q;
 140         Path *pa;
 141         double dx, dx1, dx2, dy, dy1, dy2, R, F;
 142         char* f;
 143
 144         #define FARG(c, v) case c: if(!(f=ARGF())) usage(); v = atof(f); break;
 145         ARGBEGIN {
 146         case 'N': if(!(f=ARGF())) usage(); N = atoi(f); break;
 147         case 'P': if(!(f=ARGF())) usage(); pathlen = atoi(f); break;
 148         FARG('v', v0);
 149         FARG('x', xmin);
 150         FARG('X', xmax);
 151         FARG('y', ymin);
 152         FARG('Y', ymax);
 153         FARG('t', dt);
 154         default: usage();
 155         } ARGEND;
 156
 157         if(pathlen == 0) {
 158                 nosnake = 1;
 159                 pathlen = 1000;
 160         }
 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);
 167         }
 168         prev = A;
 169         cur = B;
 170         srand(time(0));
 171         initdraw(0, 0, "Molecular Dynamics");
 172         einit(Emouse | Ekeyboard);
 173         reset();
 174
 175         for(i=0;; i++) {
 176                 if(i == pathlen)
 177                         i = 0;
 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;
 183                                 dx = mini(dx1, dx2);
 184                                 dy1 = fabs(p->y - q->y);
 185                                 dy2 = ymax - ymin - dy1;
 186                                 dy = mini(dy1, dy2);
 187                                 R = dx*dx + dy*dy;
 188                                 if(R >= 9) continue;
 189                                 R = 1/sqrt(R);
 190                                 double R2, R4, R6, R12;
 191                                 R2 = R * R;
 192                                 R4 = R2 * R2;
 193                                 R6 = R4 * R2;
 194                                 R12 = R6 * R6;
 195                                 F = 24*(2*R12 - R6);
 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;
 200                                 dx *= F;
 201                                 dy *= F;
 202                                 (p-prev+cur)->ax += dx;
 203                                 (p-prev+cur)->ay += dy;
 204                                 (q-prev+cur)->ax -= dx;
 205                                 (q-prev+cur)->ay -= dy;
 206                         }
 207                 }
 208                 for(j=0;j<N;j++) {
 209                         pa = paths+j;
 210                         p = prev+j;
 211                         q = cur+j;
 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);
 216                         q->prevx = p->x;
 217                         q->prevy = p->y;
 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;}
 222                         q->col = p->col;
 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);
 226                 }
 227
 228                 Particle* tmp = prev;
 229                 prev = cur;
 230                 cur = tmp;
 231                 flushimage(display, 1);
 232
 233
 234                 if(ecankbd()) {
 235                         switch(ekbd()) {
 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;
 240                         }
 241                 }
 242         }
 243 }
 244
 245 void
 246 eresized(int new)
 247 {
 248         if(new) getwindow(display, Refnone);
 249 }