1 #include <u.h> 2 #include <libc.h> 3 #include <draw.h> 4 5 #define NPJW 48 6 #define NSPLIT 6 7 #define NDOT 14 8 #define MAXTRACKS 128 9 #define MINV 6 10 #define PDUP 1 /* useful for debugging */ 11 12 #define nels(a) (sizeof(a) / sizeof(a[0])) 13 #define imin(a,b) (((a) <= (b)) ? (a) : (b)) 14 #define imax(a,b) (((a) >= (b)) ? (a) : (b)) 15 #define sqr(x) ((x) * (x)) 16 17 typedef struct vector vector; 18 struct vector 19 { 20 double x; 21 double y; 22 }; 23 24 typedef struct dot Dot; 25 struct dot 26 { 27 Point pos; 28 Point ivel; 29 vector vel; 30 double mass; 31 int charge; 32 int height; /* precalculated for speed */ 33 int width; /* precalculated for speed */ 34 int facei; 35 int spin; 36 Image *face; 37 Image *faces[4]; 38 Image *mask; 39 Image *masks[4]; 40 Image *save; 41 int ntracks; 42 Point track[MAXTRACKS]; 43 int next_clear; 44 int next_write; 45 }; 46 47 Dot dot[NDOT]; 48 Dot *edot = &dot[NDOT]; 49 int total_spin = 3; 50 vector no_gravity; 51 vector gravity = 52 { 53 0.0, 54 1.0, 55 }; 56 57 /* static Image *track; */ 58 static Image *im; 59 static int track_length; 60 static int track_width; 61 static int iflag; 62 static double k_floor = 0.9; 63 Image *screen; 64 65 #include "andrew.bits" 66 #include "bart.bits" 67 #include "bwk.bits" 68 #include "dmr.bits" 69 #include "doug.bits" 70 #include "gerard.bits" 71 #include "howard.bits" 72 #include "ken.bits" 73 #include "philw.bits" 74 #include "pjw.bits" 75 #include "presotto.bits" 76 #include "rob.bits" 77 #include "sean.bits" 78 #include "td.bits" 79 80 Image* 81 eallocimage(Display *d, Rectangle r, ulong chan, int repl, int col) 82 { 83 Image *i; 84 i = allocimage(d, r, chan, repl, col); 85 if(i == nil) { 86 fprint(2, "cannot allocate image\n"); 87 exits("allocimage"); 88 } 89 return i; 90 } 91 92 /* 93 * p1 dot p2 94 * inner product 95 */ 96 static 97 double 98 dotprod(vector *p1, vector *p2) 99 { 100 return p1->x * p2->x + p1->y * p2->y; 101 } 102 103 /* 104 * r = p1 + p2 105 */ 106 static 107 void 108 vecaddpt(vector *r, vector *p1, vector *p2) 109 { 110 r->x = p1->x + p2->x; 111 r->y = p1->y + p2->y; 112 } 113 114 /* 115 * r = p1 - p2 116 */ 117 static 118 void 119 vecsub(vector *r, vector *p1, vector *p2) 120 { 121 r->x = p1->x - p2->x; 122 r->y = p1->y - p2->y; 123 } 124 125 /* 126 * r = v * s 127 * multiplication of a vector by a scalar 128 */ 129 static 130 void 131 scalmult(vector *r, vector *v, double s) 132 { 133 r->x = v->x * s; 134 r->y = v->y * s; 135 } 136 137 static 138 double 139 separation(Dot *p, Dot *q) 140 { 141 return sqrt((double)(sqr(q->pos.x - p->pos.x) + sqr(q->pos.y - p->pos.y))); 142 } 143 144 static 145 int 146 close_enough(Dot *p, Dot *q, double *s) 147 { 148 int sepx; 149 int width; 150 int sepy; 151 152 sepx = p->pos.x - q->pos.x; 153 width = p->width; 154 155 if (sepx < -width || sepx > width) 156 return 0; 157 158 sepy = p->pos.y - q->pos.y; 159 160 if (sepy < -width || sepy > width) 161 return 0; 162 163 if ((*s = separation(p, q)) > (double)width) 164 return 0; 165 166 return 1; 167 } 168 169 static 170 void 171 ptov(vector *v, Point *p) 172 { 173 v->x = (double)p->x; 174 v->y = (double)p->y; 175 } 176 177 static 178 void 179 vtop(Point *p, vector *v) 180 { 181 p->x = (int)(v->x + 0.5); 182 p->y = (int)(v->y + 0.5); 183 } 184 185 static 186 void 187 exchange_spin(Dot *p, Dot *q) 188 { 189 int tspin; 190 191 if (p->spin == 0 && q->spin == 0) 192 return; 193 194 tspin = p->spin + q->spin; 195 196 p->spin = imax(nrand(imin(3, tspin + 1)), tspin + 1 - 3); 197 198 if (p->spin == 0) 199 { 200 p->face = p->faces[0]; 201 p->mask = p->masks[0]; 202 } 203 204 q->spin = tspin - p->spin; 205 206 if (q->spin == 0) 207 { 208 q->face = q->faces[0]; 209 q->mask = q->masks[0]; 210 } 211 } 212 213 static 214 void 215 exchange_charge(Dot *p, Dot *q) 216 { 217 if (p->charge == 0 && q->charge == 0) 218 { 219 switch (nrand(16)) 220 { 221 case 0: 222 p->charge = 1; 223 q->charge = -1; 224 break; 225 226 case 1: 227 p->charge = -1; 228 q->charge = 1; 229 break; 230 } 231 } 232 } 233 234 /* 235 * The aim is to completely determine the final 236 * velocities of the two interacting particles as 237 * a function of their initial velocities, masses 238 * and point of collision. We start with the two 239 * quantities that we know are conserved in 240 * elastic collisions -- momentum and kinetic energy. 241 * 242 * Let the masses of the particles be m1 and m2; 243 * their initial velocities V1 and V2 (vector quantities 244 * represented by upper case variables, scalars lower 245 * case); their final velocities V1' and V2'. 246 * 247 * Conservation of momentum gives us: 248 * 249 * (1) m1 * V1 + m2 * V2 = m1 * V1' + m2 * V2' 250 * 251 * and conservation of kinetic energy gives: 252 * 253 * (2) 1/2 * m1 * |V1| * |V1| + 1/2 * m2 * |V2| * |V2| = 254 * 1/2 * m1 * |V1'| * |V1'| + 1/2 * m2 * |V2'| * |V2'| 255 * 256 * Then, decomposing (1) into its 2D scalar components: 257 * 258 * (1a) m1 * v1x + m2 * v2x = m1 * v1x' + m2 * v2x' 259 * (1b) m1 * v1y + m2 * v2y = m1 * v1y' + m2 * v2y' 260 * 261 * and simplifying and expanding (2): 262 * 263 * (2a) m1 * (v1x * v1x + v1y * v1y) + 264 * m2 * (v2x * v2x + v2y * v2y) 265 * = 266 * m1 * (v1x' * v1x' + v1y' * v1y') + 267 * m2 * (v2x' * v2x' + v2y' * v2y') 268 * 269 * We know m1, m2, V1 and V2 which leaves four unknowns: 270 * 271 * v1x', v1y', v2x' and v2y' 272 * 273 * but we have just three equations: 274 * 275 * (1a), (1b) and (2a) 276 * 277 * To remove this extra degree of freedom we add the assumption that 278 * the forces during the collision act instantaneously and exactly 279 * along the (displacement) vector joining the centres of the two objects. 280 * 281 * And unfortunately, I've forgotten the extra equation that this 282 * adds to the system(!), but it turns out that this extra constraint 283 * does allow us to solve the augmented system of equations. 284 * (I guess I could reverse engineer it from the code ...) 285 */ 286 static 287 void 288 rebound(Dot *p, Dot *q, double s) 289 { 290 vector pposn; 291 vector qposn; 292 vector pvel; 293 vector qvel; 294 vector newqp; /* vector difference of the positions */ 295 vector newqpu; /* unit vector parallel to newqp */ 296 vector newqv; /* " " " " velocities */ 297 double projp; /* projection of p's velocity in newqp direction */ 298 double projq; /* " " q's velocity in newqp direction */ 299 double pmass; /* mass of p */ 300 double qmass; /* mass of q */ 301 double tmass; /* sum of mass of p and q */ 302 double dmass; /* difference of mass of p and q */ 303 double newp; 304 double newq; 305 306 if (s == 0.0) 307 return; 308 309 ptov(&pposn, &p->pos); 310 ptov(&qposn, &q->pos); 311 pvel = p->vel; 312 qvel = q->vel; 313 314 /* 315 * Transform so that p is stationary at (0,0,0); 316 */ 317 vecsub(&newqp, &qposn, &pposn); 318 vecsub(&newqv, &qvel, &pvel); 319 320 scalmult(&newqpu, &newqp, -1.0 / s); 321 322 if (dotprod(&newqv, &newqpu) <= 0.0) 323 return; 324 325 projp = dotprod(&pvel, &newqpu); 326 projq = dotprod(&qvel, &newqpu); 327 pmass = p->mass; 328 qmass = q->mass; 329 tmass = pmass + qmass; 330 dmass = pmass - qmass; 331 newp = (2.0 * qmass * projq + dmass * projp) / tmass; 332 newq = (2.0 * pmass * projp - dmass * projq) / tmass; 333 scalmult(&newqp, &newqpu, projp - newp); 334 scalmult(&newqv, &newqpu, projq - newq); 335 vecsub(&pvel, &pvel, &newqp); 336 vecsub(&qvel, &qvel, &newqv); 337 338 p->vel = pvel; 339 vtop(&p->ivel, &p->vel); 340 q->vel = qvel; 341 vtop(&q->ivel, &q->vel); 342 343 exchange_spin(p, q); 344 345 if (iflag) 346 exchange_charge(p, q); 347 } 348 349 static 350 int 351 rrand(int lo, int hi) 352 { 353 return lo + nrand(hi + 1 - lo); 354 } 355 356 static 357 void 358 drawdot(Dot *d) 359 { 360 Rectangle r; 361 char buf[10]; 362 363 r = d->face->r; 364 r = rectsubpt(r, d->face->r.min); 365 r = rectaddpt(r, d->pos); 366 r = rectaddpt(r, screen->r.min); 367 368 draw(screen, r, d->face, d->mask, d->face->r.min); 369 370 if(PDUP > 1) { /* assume debugging */ 371 sprint(buf, "(%d,%d)", d->pos.x, d->pos.y); 372 string(screen, r.min, display->black, ZP, font, buf); 373 } 374 } 375 376 static 377 void 378 undraw(Dot *d) 379 { 380 Rectangle r; 381 r = d->face->r; 382 r = rectsubpt(r, d->face->r.min); 383 r = rectaddpt(r, d->pos); 384 r = rectaddpt(r, screen->r.min); 385 386 draw(screen, r, display->black, d->mask, d->face->r.min); 387 388 /* 389 if (track_width > 0) 390 { 391 if (d->ntracks >= track_length) 392 { 393 bitblt(&screen, d->track[d->next_clear], track, track->r, D ^ ~S); 394 d->next_clear = (d->next_clear + 1) % track_length; 395 d->ntracks--; 396 } 397 398 d->track[d->next_write] = Pt(d->pos.x + d->width / 2 - track_width / 2, d->pos.y + d->height / 2 - track_width / 2); 399 bitblt(&screen, d->track[d->next_write], track, track->r, D ^ ~S); 400 d->next_write = (d->next_write + 1) % track_length; 401 d->ntracks++; 402 } 403 */ 404 } 405 406 static void 407 copy(Image *a, Image *b) 408 { 409 if(PDUP==1 || eqrect(a->r, b->r)) 410 draw(a, a->r, b, nil, b->r.min); 411 else { 412 int x, y; 413 for(x = a->r.min.x; x < a->r.max.x; x++) 414 for(y = a->r.min.y; y < a->r.max.y; y++) 415 draw(a, Rect(x,y,x+1,y+1), b, nil, Pt(x/PDUP,y/PDUP)); 416 } 417 } 418 419 static void 420 transpose(Image *b) 421 { 422 int x, y; 423 424 for(x = b->r.min.x; x < b->r.max.x; x++) 425 for(y = b->r.min.y; y < b->r.max.y; y++) 426 draw(im, Rect(y,x,y+1,x+1), b, nil, Pt(x,y)); 427 428 draw(b, b->r, im, nil, ZP); 429 } 430 431 static void 432 reflect1_lr(Image *b) 433 { 434 int x; 435 for(x = 0; x < PDUP*NPJW; x++) 436 draw(im, Rect(x, 0, x, PDUP*NPJW), b, nil, Pt(b->r.max.x-x,0)); 437 draw(b, b->r, im, nil, ZP); 438 } 439 440 static void 441 reflect1_ud(Image *b) 442 { 443 int y; 444 for(y = 0; y < PDUP*NPJW; y++) 445 draw(im, Rect(0, y, PDUP*NPJW, y), b, nil, Pt(0,b->r.max.y-y)); 446 draw(b, b->r, im, nil, ZP); 447 } 448 449 static 450 void 451 rotate_clockwise(Image *b) 452 { 453 transpose(b); 454 reflect1_lr(b); 455 } 456 457 static 458 void 459 spin(Dot *d) 460 { 461 int i; 462 463 if (d->spin > 0) 464 { 465 i = (d->facei + d->spin) % nels(d->faces); 466 d->face = d->faces[i]; 467 d->mask = d->masks[i]; 468 d->facei = i; 469 } 470 sleep(0); 471 } 472 473 static 474 void 475 restart(Dot *d) 476 { 477 static int beam; 478 479 d->charge = 0; 480 481 d->pos.x = 0; 482 d->vel.x = 20.0 + (double)rrand(-2, 2); 483 484 if (beam ^= 1) 485 { 486 d->pos.y = screen->r.max.y - d->height; 487 d->vel.y = -20.0 + (double)rrand(-2, 2); 488 } 489 else 490 { 491 d->pos.y = 0; 492 d->vel.y = 20.0 + (double)rrand(-2, 2); 493 } 494 495 vtop(&d->ivel, &d->vel); 496 } 497 498 static 499 void 500 upd(Dot *d) 501 { 502 Dot *dd; 503 504 spin(d); 505 506 /* 507 * Bounce off others? 508 */ 509 for (dd = d + 1; dd != edot; dd++) 510 { 511 double s; 512 513 if (close_enough(d, dd, &s)) 514 rebound(d, dd, s); 515 } 516 517 if (iflag) 518 { 519 /* 520 * Going off-screen? 521 */ 522 if 523 ( 524 d->pos.x + d->width < 0 525 || 526 d->pos.x >= Dx(screen->r) 527 || 528 d->pos.y + d->height < 0 529 || 530 d->pos.y >= Dy(screen->r) 531 ) 532 restart(d); 533 534 /* 535 * Charge. 536 */ 537 if (d->charge != 0) 538 { 539 /* 540 * TODO: This is wrong -- should 541 * generate closed paths. Can't 542 * simulate effect of B field just 543 * by adding in an orthogonal 544 * velocity component... (sigh) 545 */ 546 vector f; 547 double s; 548 549 f.x = -d->vel.y; 550 f.y = d->vel.x; 551 552 if ((s = sqrt(sqr(f.x) + sqr(f.y))) != 0.0) 553 { 554 scalmult(&f, &f, -((double)d->charge) / s); 555 556 vecaddpt(&d->vel, &f, &d->vel); 557 vtop(&d->ivel, &d->vel); 558 } 559 } 560 } 561 else 562 { 563 /* 564 * Bounce off left or right border? 565 */ 566 if (d->pos.x < 0) 567 { 568 d->vel.x = fabs(d->vel.x); 569 vtop(&d->ivel, &d->vel); 570 } 571 else if (d->pos.x + d->width >= Dx(screen->r)) 572 { 573 d->vel.x = -fabs(d->vel.x); 574 vtop(&d->ivel, &d->vel); 575 } 576 577 /* 578 * Bounce off top or bottom border? 579 * (bottom is slightly inelastic) 580 */ 581 if (d->pos.y < 0) 582 { 583 d->vel.y = fabs(d->vel.y); 584 vtop(&d->ivel, &d->vel); 585 } 586 else if (d->pos.y + d->height >= Dy(screen->r)) 587 { 588 if (gravity.y == 0.0) 589 d->vel.y = -fabs(d->vel.y); 590 else if (d->ivel.y >= -MINV && d->ivel.y <= MINV) 591 { 592 /* 593 * y-velocity is too small -- 594 * give it a random kick. 595 */ 596 d->vel.y = (double)-rrand(30, 40); 597 d->vel.x = (double)rrand(-10, 10); 598 } 599 else 600 d->vel.y = -fabs(d->vel.y) * k_floor; 601 vtop(&d->ivel, &d->vel); 602 } 603 } 604 605 if (gravity.x != 0.0 || gravity.y != 0.0) 606 { 607 vecaddpt(&d->vel, &gravity, &d->vel); 608 vtop(&d->ivel, &d->vel); 609 } 610 611 d->pos = addpt(d->pos, d->ivel); 612 613 drawdot(d); 614 } 615 616 static 617 void 618 setup(Dot *d, char *who, uchar *face, int n_els) 619 { 620 int i, j, k, n; 621 int repl; 622 uchar mask; 623 int nbits, bits; 624 uchar tmpface[NPJW*NPJW]; 625 uchar tmpmask[NPJW*NPJW]; 626 static Image *im; /* not the global */ 627 static Image *imask; 628 629 if(im == nil) 630 { 631 im = eallocimage(display, Rect(0,0,NPJW,NPJW), CMAP8, 0, DNofill); 632 imask = eallocimage(display, Rect(0,0,NPJW,NPJW), CMAP8, 0, DNofill); 633 } 634 635 repl = (NPJW*NPJW)/n_els; 636 if(repl > 8) { 637 fprint(2, "can't happen --rsc\n"); 638 exits("repl"); 639 } 640 nbits = 8/repl; 641 mask = (1<<(nbits))-1; 642 643 if(0) print("converting %s... n_els=%d repl=%d mask=%x nbits=%d...\n", 644 who, n_els, repl, mask, nbits); 645 n = 0; 646 for (i = 0; i < n_els; i++) 647 { 648 for(j = repl; j--; ) { 649 bits = (face[i] >> (j*nbits)) & mask; 650 tmpface[n] = 0; 651 tmpmask[n] = 0; 652 for(k = 0; k < repl; k++) 653 { 654 tmpface[n] |= (mask-bits) << (k*nbits); 655 tmpmask[n] |= (bits==mask ? 0 : mask) << (k*nbits); 656 } 657 n++; 658 } 659 660 } 661 662 if(n != sizeof tmpface) { 663 fprint(2, "can't happen2 --rsc\n"); 664 exits("n!=tmpface"); 665 } 666 667 loadimage(im, im->r, tmpface, n); 668 loadimage(imask, imask->r, tmpmask, n); 669 670 for (i = 0; i < nels(d->faces); i++) 671 { 672 d->faces[i] = eallocimage(display, Rect(0,0,PDUP*NPJW, PDUP*NPJW), 673 screen->chan, 0, DNofill); 674 d->masks[i] = eallocimage(display, Rect(0,0,PDUP*NPJW, PDUP*NPJW), 675 GREY1, 0, DNofill); 676 677 switch (i) { 678 case 0: 679 copy(d->faces[i], im); 680 copy(d->masks[i], imask); 681 break; 682 683 case 1: 684 copy(d->faces[i], im); 685 copy(d->masks[i], imask); 686 rotate_clockwise(d->faces[i]); 687 rotate_clockwise(d->masks[i]); 688 break; 689 690 default: 691 copy(d->faces[i], d->faces[i-2]); 692 copy(d->masks[i], d->masks[i-2]); 693 reflect1_lr(d->faces[i]); 694 reflect1_ud(d->faces[i]); 695 reflect1_lr(d->masks[i]); 696 reflect1_ud(d->masks[i]); 697 break; 698 } 699 } 700 d->face = d->faces[0]; 701 d->mask = d->masks[0]; 702 703 d->height = Dy(im->r); 704 d->width = Dx(im->r); 705 706 d->mass = 1.0; 707 708 d->spin = nrand(imin(3, total_spin + 1)); 709 total_spin -= d->spin; 710 711 d->pos.x = nrand(screen->r.max.x - d->width); 712 d->pos.y = nrand(screen->r.max.y - d->height); 713 714 d->vel.x = (double)rrand(-20, 20); 715 d->vel.y = (double)rrand(-20, 20); 716 vtop(&d->ivel, &d->vel); 717 718 drawdot(d); 719 } 720 721 int 722 msec(void) 723 { 724 static int fd; 725 int n; 726 char buf[64]; 727 728 if(fd <= 0) 729 fd = open("/dev/msec", OREAD); 730 if(fd < 0) 731 return 0; 732 if(seek(fd, 0, 0) < 0) 733 return 0; 734 if((n=read(fd, buf, sizeof(buf)-1)) < 0) 735 return 0; 736 buf[n] = 0; 737 return atoi(buf); 738 } 739 740 /* 741 * debugging: make del pause so that we can 742 * inspect window. 743 */ 744 jmp_buf j; 745 static void 746 myhandler(void *v, char *s) 747 { 748 if(strcmp(s, "interrupt") == 0) 749 notejmp(v, j, -1); 750 noted(NDFLT); 751 } 752 753 void 754 main(int argc, char *argv[]) 755 { 756 int c; 757 long now, then; 758 759 ARGBEGIN 760 { 761 case 'i': 762 iflag = 1; 763 gravity = no_gravity; 764 track_length = 64; 765 track_width = 2; 766 break; 767 768 case 'k': 769 k_floor = atof(ARGF()); 770 break; 771 772 case 'n': 773 track_length = atoi(ARGF()); 774 break; 775 776 case 'w': 777 track_width = atoi(ARGF()); 778 break; 779 780 case 'x': 781 gravity.x = atof(ARGF()); 782 break; 783 784 case 'y': 785 gravity.y = atof(ARGF()); 786 break; 787 788 default: 789 fprint(2, "Usage: %s [-i] [-k k_floor] [-n track_length] [-w track_width] [-x gravityx] [-y gravityy]\n", argv0); 790 exits("usage"); 791 } ARGEND 792 793 if (track_length > MAXTRACKS) 794 track_length = MAXTRACKS; 795 796 if (track_length == 0) 797 track_width = 0; 798 799 srand(time(0)); 800 801 initdraw(0,0,0); 802 im = eallocimage(display, Rect(0, 0, PDUP*NPJW, PDUP*NPJW), CMAP8, 0, DNofill); 803 804 draw(screen, screen->r, display->black, nil, ZP); 805 806 /* track = balloc(Rect(0, 0, track_width, track_width), 0); */ 807 808 edot = &dot[0]; 809 810 setup(edot++, "andrew", andrewbits, nels(andrewbits)); 811 setup(edot++, "bart", bartbits, nels(bartbits)); 812 setup(edot++, "bwk", bwkbits, nels(bwkbits)); 813 setup(edot++, "dmr", dmrbits, nels(dmrbits)); 814 setup(edot++, "doug", dougbits, nels(dougbits)); 815 setup(edot++, "gerard", gerardbits, nels(gerardbits)); 816 setup(edot++, "howard", howardbits, nels(howardbits)); 817 setup(edot++, "ken", kenbits, nels(kenbits)); 818 setup(edot++, "philw", philwbits, nels(philwbits)); 819 setup(edot++, "pjw", pjwbits, nels(pjwbits)); 820 setup(edot++, "presotto", presottobits, nels(presottobits)); 821 setup(edot++, "rob", robbits, nels(robbits)); 822 setup(edot++, "sean", seanbits, nels(seanbits)); 823 setup(edot++, "td", tdbits, nels(tdbits)); 824 825 if(PDUP > 1) { /* assume debugging */ 826 setjmp(j); 827 read(0, &c, 1); 828 829 /* make DEL pause so that we can inspect window */ 830 notify(myhandler); 831 } 832 SET(c); 833 USED(c); 834 835 #define DELAY 100 836 for (then = msec();; then = msec()) 837 { 838 Dot *d; 839 840 for (d = dot; d != edot; d++) 841 undraw(d); 842 for (d = dot; d != edot; d++) 843 upd(d); 844 draw(screen, screen->r, screen, nil, screen->r.min); 845 flushimage(display, 1); 846 now = msec(); 847 if(now - then < DELAY) 848 sleep(DELAY - (now - then)); 849 } 850 } 851