home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume2 / molecule < prev    next >
Encoding:
Internet Message Format  |  1991-08-07  |  8.8 KB
  1. From: gfs@decwrl.dec.com@x.UUCP (Greg F. Shay)
  2. Newsgroups: comp.sources.misc
  3. Subject: v02i092: A simple 2-D particle system (animated atoms)
  4. Message-ID: <230@abvax.UUCP>
  5. Date: 5 Apr 88 13:34:17 GMT
  6. Approved: allbery@ncoast.UUCP
  7.  
  8. comp.sources.misc: Volume 2, Issue 92
  9. Submitted-By: "Greg F. Shay" <gfs@decwrl.dec.com@x.UUCP>
  10. Archive-Name: molecule
  11.  
  12. [Some people don't read Administrivia, it seems:  it's not shar'ed.  ++bsa]
  13.  
  14. For a request on comp.graphics:
  15.  
  16. With regard to the posting about a practical partical system, 
  17. here is a simple 2-d system I wrote for experimentation purposes a while
  18. ago that illustrates a few ideas in speeding up the computations and
  19. beating the O(n*n) computation complexity.
  20.  
  21. The system is a collection of simple body atoms that exhibit a force function
  22. on each other characteristic of a actual atoms:
  23. I use a piece of a cubic function to model this.  Note beyond a distance D
  24. there is no interaction.  At closer than D, atoms first attract (like
  25. surface tension phenomena), and then when too close, there is repulsion 
  26. because the atoms cannot occupy the same space.
  27.  
  28.    | +
  29.    | +  (repulsion)
  30.    | +    Force ^
  31.    | +
  32.    | +
  33.    |  +
  34.    |   +                    D
  35.    -----+-------------------++++----- Distance between atoms ->
  36.           +               +
  37.             +            +
  38.                +      +  (attraction)
  39.                    +
  40.  
  41. The program is rough, but it does work.  The first program computes the
  42. positions of the atoms for every time frame and outputs them into a file.
  43. The second reads the file and animates the data on a graphics screen.
  44. If the number of atoms is small, you can pipe the two together and compute
  45. in real time (I have even run the programs on different machines and ethernetted
  46. the intermediate data live.  If the computation gets too slow for good animation
  47. precompute and save the intermediate binary file.
  48.  
  49. The display program is presently written for a Masscomp system.  The basic
  50. calls are: mgibox which is used draw a solid box of 2x2 pixels.  Mgifb
  51. is flipping frame buffers at the next vertical retrace, and mgiclear is
  52. used to erase the back (non visible) frame buffer.  Adaptations can be
  53. made to other graphics output hardware.  If anyone writes a display routine
  54. for a SUN, let me have a copy please, I have not the time to do it myself.
  55.  
  56. The program computes forces, divides by mass to get acceleration, integrates 
  57. to get velocity, again to get distance.  The input file specifies the distance
  58. threshold D, mass, dt (delta time), number of iterations, and number of atoms,
  59. in that order. See the the program for more info.  If the time is too long
  60. or mass too small (forces are fixed), integration errors will get too large
  61. and the system will diverge (REALLY diverge, try it, set dt = 0.5, BOOM!).
  62.  
  63. A very interesting input atom set is a square lattice of atoms spaced about
  64. 0.25 or so apart (5 x 5 or so), the lattice will oscillate symmetrically before
  65. errors creep in and make the mass amorphous.  You can adjust the display program
  66. to magnify just the area where you put your lattice to show the details better.
  67. A VERY interesting modification I made once was to put in a 'cooling' function,
  68. where the reacting force was set to say 0.9700 of the sum of forces.  This
  69. causes the kinetic energy to gradually bleed off.  A cluster of atoms will
  70. eventually cool into a regular lattice, a crystal! Face centered cubic no less!
  71. Just like we learned in chemistry 101!
  72.  
  73. Now, the comments on computational speedup.  All the atoms are kept sorted
  74. by there y coordinate.  An index array is used to make swapping of floats
  75. unneccessary.  Making use of the distance D, for a given atom, all atoms
  76. ahead of it on the list by more than a distance D and all the atoms after
  77. it on the list by a distance D or more do not have to be considered.
  78. Only the atoms in the slab from (y-D,y+D) are considered.  Next, the
  79. x distance is tested to eliminate more.  Finally, dx^2+dy^2, the true
  80. distance is calculated only if inside the bounding box.
  81. After the forces are summed for all interacting atoms, the velocity and
  82. positions are updated.  A bubble sort is used because it has the property
  83. that if only a few items moved relative to each other, it is fairly quick.
  84. I know a better sort would be appropriate, but I find bubble's easy to write.
  85.  
  86. This speedup by trivial elimination is quite visible when running a realtime
  87. (computation on the fly) animation.  When there are many atoms clustered 
  88. together in a multi-body interaction, the program slows down.  When they 
  89. are separated, the program speeds up.  I would be curious to know how much
  90. of the time is spent in the sorting and how much in the floating computations.
  91. I use a 20mhz 68020+68881 system and it can perform a 6x6 (36 atom) lattice
  92. in real time over the ethernet quite well at 15-20fps or so.
  93.  
  94. Well, I hope someone has fun with this.  Drop me a line if anyone gets it
  95. running.
  96.         Greg Shay
  97.  
  98.     mandrill !
  99.     pyramid  !... abvax!gfs
  100.     decvax   !
  101.  
  102.  
  103.  
  104. /*--------Program molecule.c-------------------------*/
  105. #include <math.h>
  106.  
  107. main()
  108. {
  109.     float dx,dy,d;
  110.     float forcefunction();
  111.     float fx,fy,F;
  112.     int   aindex[100];
  113.     int   a,b,c;
  114.     int  start,end;
  115.     int  maxindex,itcount,itmax;
  116.     short i,j,k,l;
  117.     short idat[2],maxin,itinmax;
  118.     float mass,dt,thresh;
  119.     float atompos[100][2];
  120.     float atomv[100][2];
  121.     float atomf[100][2];
  122.     
  123. /* Read in initial positions of atoms*/
  124.  
  125. scanf("%f",&thresh);
  126. scanf("%f",&mass);
  127. scanf("%f",&dt);
  128. scanf("%d\n",&itmax);
  129. scanf("%d\n",&maxindex);
  130. maxindex = maxindex-1;
  131. maxin = (short)maxindex;
  132. itinmax = (short)itmax;
  133.  
  134. /*printf("%d \n %d \n",maxindex,itmax);*/
  135. write(1,&maxin,2);
  136. write(1,&itinmax,2);
  137.  
  138. for(a=0;a <= maxindex;a++)
  139.     {
  140.     for(i=0;i<2;i++)
  141.         {
  142.         scanf ("%f ",&atompos[a][i]);
  143.         }
  144.     for(i=0;i<2;i++)
  145.         {
  146.         scanf ("%f ",&atomv[a][i]);
  147.         }
  148.     scanf ("\n");
  149.     }
  150.  
  151.  
  152. /* Initialize index array */
  153. for (a=0;a <= maxindex;a++)
  154.     aindex[a]=a;
  155.  
  156.  
  157. /* Now start iteration loop */
  158. for (itcount =0;itcount<itmax;itcount++)
  159.     {
  160.  
  161.     bsort(maxindex,aindex,atompos);
  162.     for(b=0;b<=maxindex;b++) 
  163.         {
  164.         idat[0]=(short)(10*atompos[aindex[b]][0]);
  165.         idat[1]=(short)(10*atompos[aindex[b]][1]);
  166.         write(1,idat,4);
  167. /*printf("%d %d\n",(int)(10*atompos[aindex[b]][0]),(int)(10*atompos[aindex[b]][1])); */
  168.         }
  169.     
  170.     start = 0;
  171.     end   = 0;
  172.  
  173.     for(a=0;a<=maxindex;a++)
  174.         {
  175.         /* Do every atom */
  176.         /* Adjust start range pointer */
  177.         while (
  178.              (atompos[aindex[a]][1]-atompos[aindex[start]][1]) > thresh
  179.               )
  180.             start++;
  181.  
  182.         /* Adjust end range pointer */
  183.         while (
  184.            ((atompos[aindex[end]][1]-atompos[aindex[a]][1]) <= thresh)
  185.            &&
  186.            (end  < maxindex)
  187.               )
  188.             end++;
  189.  
  190.         /* Evaluate forces on this atom 'a' from every atom between*/
  191.         /*  start and end on the index list */
  192.      atomf[aindex[a]][0]=0.0;
  193.      atomf[aindex[a]][1]=0.0;
  194.      for(b=start;b<=end;b++)
  195.         if (b != a)
  196.         {
  197.         dx=atompos[aindex[b]][0]-atompos[aindex[a]][0];
  198.         if (dx <= thresh)
  199.             {
  200.             dy=atompos[aindex[b]][1]-atompos[aindex[a]][1];
  201.             d=sqrt(dx*dx+dy*dy);
  202.             F=forcefunction(d);
  203.             fx = F*(dx/d);
  204.             fy = F*(dy/d);
  205.  
  206.             atomf[aindex[a]][0]= fx+atomf[aindex[a]][0];
  207.             atomf[aindex[a]][1]= fy+atomf[aindex[a]][1];
  208.  
  209.             }
  210.         }
  211.  
  212.         }  /* Do all atoms */
  213. /* Now resolve forces */
  214. for(a=0;a<=maxindex;a++)
  215.     for(c=0;c<2;c++)
  216.     {
  217.     atomv[a][c]= dt*atomf[a][c]/mass+atomv[a][c];
  218.     atompos[a][c]=(dt*atomv[a][c])+atompos[a][c];
  219.     }
  220.  
  221. /* Now loop back for more iterations */
  222.         }  /* End of iteration loop */
  223. }  /* End of Main */
  224.  
  225.  
  226. bsort(limit,index,data)
  227.     int limit;
  228.     int index[];
  229.     float data[][2];
  230. {
  231.     short flag,j;
  232.     int   temp;
  233.  
  234.  
  235. do     {
  236.     flag = 0;
  237.     for (j=0;j<=limit-1;j++)
  238.         if (data[index[j]][1] > data[index[j+1]][1])
  239.             {
  240.             flag = 1;
  241.             temp = index[j+1];
  242.             index[j+1]=index[j];
  243.             index[j]=temp;
  244.             }
  245.     }
  246. while (flag == 1);
  247.  
  248. } /* end of bubble sort routine*/
  249.  
  250. float forcefunction(d)
  251.     float d;
  252. {
  253. /* return force depending on distance d */
  254.  
  255.     if (d > 1.0) return (0.);
  256.     d= (1.0-d);
  257.     return(60.0*d-160.0*d*d*d);
  258. }
  259. /*--------End Program molecule.c-------------------------*/
  260.  
  261.  
  262. /*-------- Program Display.c-------------------------------*/
  263. main()
  264. {
  265.     int i,j,k,x,y,sx,sy;
  266.     int fb,maxnum,iterno;
  267.     short idat[2],maxin,iterin;
  268.  
  269. system("cleargp");
  270. mginit (0,0);
  271. mgimodfunc (3,0,3);
  272. mgicm(2,4096*15);
  273.     
  274. mgihue(2);
  275.  
  276. fb = 1;
  277.  
  278. read(0,&maxin,2);
  279. read(0,&iterin,2);
  280. maxnum = (int)maxin;
  281. iterno=(int)iterin;
  282.  
  283. /*scanf("%d %d\n",&maxnum,&iterno);*/
  284.  
  285. for(k=0;k<iterno;k++)
  286.     {
  287.     mgifb(fb,3-fb);
  288.     fb = 3-fb;
  289.     mgiclearpln(2,-1,0);
  290.     for(i=0;i<=maxnum;i++)
  291.         {
  292.         /*scanf("%d %d\n",&x,&y);*/
  293.         read(0,idat,4);
  294.         mgibox(idat[0],idat[1],idat[0]+2,idat[1]+2);
  295.         }
  296.     }
  297. mgideagp (); 
  298. }
  299. /*----------_End program display.c --------------------------*/
  300. /* Sample input files                           */
  301. 1.0
  302. 1.0
  303. 0.07
  304. 1000
  305. 10
  306.  
  307. 16 16     0 0
  308. 10 10    0 0
  309. 14 14    1 1
  310. 24 24    0 0
  311. 23.5 23.5    0 0
  312. 18 18    -.25 .2
  313. 20 20    -.1 -.1
  314. 32 32    -.4 -.4
  315. 28 28    .4  -.3
  316. 40 40    -2 -2
  317.  
  318. /* Sample input files                           */
  319.  
  320. 1.0
  321. .2
  322. 0.01
  323. 800
  324. 10
  325.  
  326. 31 31.2    0 0
  327. 31 31    0 0
  328. 19 19    9 9.
  329. 19.5 19.5 9 9.
  330. 19.0 19.5 9 9.
  331. 31.2 31    0 0
  332. 31.2 31.5   0 0
  333. 31.4 32 0 0
  334. 30.5 31 0 0
  335. 30.4 31.5 0 0
  336. /*------------------------------------------*/
  337.