home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / X / mit / demos / xgas / dynamics.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-07-22  |  7.4 KB  |  292 lines
  1. /*
  2.  * dynamics.c
  3.  *   xgas: Copyright 1990 Larry Medwin: @(#)dynamics.c    1.5 2/9/90
  4.  *   Larry Medwin -- Dec. 15, 1989
  5.  *   fixed: lm 3-24-91
  6.  */
  7.  
  8. #include "xgas.h"
  9. static void inertia();
  10. static void collide();
  11. static double hit();
  12. static double hitHole();
  13.  
  14. /* DYNAMICS */
  15. void dynamics( mol, data )
  16.     Molecule *mol;
  17.     LabData *data;
  18. {
  19.     /* Move molecule */
  20.     while( mol->collisionTime <= data->timestepSize) {
  21.     collide( mol, data );
  22.     findNextCollision( mol, data);
  23.     }
  24.     inertia( mol, data );
  25. }
  26.  
  27. /* INERTIA() */
  28. static void inertia( mol, data)
  29.     Molecule *mol;
  30.     LabData *data;
  31. {
  32.     /* Solve equations of motion */
  33.     mol->pos.x = mol->xCoeff[0] += mol->xCoeff[1] * data->timestepSize;
  34.     mol->pos.y = mol->yCoeff[0] += mol->yCoeff[1] * data->timestepSize;
  35.     mol->collisionTime -= data->timestepSize;
  36. }
  37.  
  38. /*
  39.  * COLLIDE()
  40.  * bounce a molecule off a wall with given temperature
  41.  *   recompute its new velocity,
  42.  *   including next collision info.
  43.  */
  44. static void collide( mol, data)
  45.     Molecule *mol;
  46.     LabData *data;
  47. {
  48.     float vMagnitude, hypot;
  49.     float x, y;
  50.     float theta, thetaRand;
  51.     WallType *wallParams;
  52.     int corner;
  53.  
  54.  
  55.     /* Find out about the wall we're about to hit */
  56.     wallParams = &WallParam[ mol->collisionWall->type ];
  57.  
  58.     /*
  59.      * Are we colliding with a corner?
  60.      */
  61.     if ((corner = whichCorner( mol->collisionPos.x, mol->collisionPos.y,
  62.         mol->thisBox, data)) != 0)
  63.     wallParams = &WallParam[ corner ];
  64.  
  65.     /* Reflect the trajectory */
  66.     mol->xCoeff[1] *= wallParams->wallReflect.x;
  67.     mol->yCoeff[1] *= wallParams->wallReflect.y;
  68.  
  69.     /* Rotate the trajectory into the canonical orientation */
  70.     x = mol->xCoeff[1] * wallParams->toRotate[0][0] +
  71.     mol->yCoeff[1] * wallParams->toRotate[0][1];
  72.     y = mol->xCoeff[1] * wallParams->toRotate[1][0] +
  73.     mol->yCoeff[1] * wallParams->toRotate[1][1];
  74.  
  75.     /* Find reflection angle theta */
  76.     if ( fabs ( x ) > SMALL ) theta = atan ( y / x );
  77.     else if ( y > 0.0 ) theta = M_PI_2;
  78.     else theta = -M_PI_2;
  79.  
  80.     /* Randomize angle by the weight "randomBounce" */
  81.     if (corner) thetaRand = frand( 0.0, M_PI_2);
  82.     else thetaRand = frand( -M_PI_2, M_PI_2);
  83.  
  84.     theta = (1.0 - data->randomBounce) * theta
  85.         + data->randomBounce * thetaRand;
  86.  
  87.     /* Molecule approaches equilibrium according to weight "equilibrium" */
  88.     mol->temperature = (1.0 - data->equilibrium) * mol->temperature
  89.     + data->equilibrium * data->chamber[mol->thisBox].temperature;
  90.  
  91.     /* Find velocity vector magnitude */
  92.     vMagnitude = vEquilibrium( mol->temperature);
  93.  
  94.     /* ... and components */
  95.     x = vMagnitude * cos( theta );
  96.     y = vMagnitude * sin( theta );
  97.  
  98.     /* Rotate the trajectory back to its original orientation */
  99.     mol->xCoeff[1] = x * wallParams->fromRotate[0][0] +
  100.             y * wallParams->fromRotate[0][1];
  101.     mol->yCoeff[1] = x * wallParams->fromRotate[1][0] +
  102.             y * wallParams->fromRotate[1][1];
  103. }
  104.  
  105. /*
  106.  * FIND TIME AND POSITION OF NEXT COLLISION
  107.  *   given the equations of motion for the molecule
  108.  */
  109. void findNextCollision( mol, data)
  110.     Molecule *mol;
  111.     LabData *data;
  112. {
  113.     int i;
  114.     int box;
  115.     double deltaT;
  116.     Coord lastCollision;
  117.  
  118.     /* Save position of last collision */
  119.     lastCollision.x = mol->collisionPos.x;
  120.     lastCollision.y = mol->collisionPos.y;
  121.  
  122.     box = mol->thisBox;
  123.  
  124.     /* Do we go through the hole? (walls[0]) */
  125.     if( hitHole( mol, &data->chamber[box].walls[0], data) > 0.0) {
  126.  
  127.         /* Move into the other box */
  128.         box = 1 - box;
  129.     }
  130.  
  131.     /*
  132.      * Now check for collisions with the walls
  133.      * Update collisionPos in hit()
  134.      */
  135.     for( i=1; i<NWALLS ; i++ ) {
  136.     if((mol->collisionWall != &data->chamber[box].walls[i])
  137.       && (( deltaT = hit( mol, &data->chamber[box].walls[i]))
  138.          > 0.0)) {
  139.       break;
  140.     }
  141.     }
  142.  
  143.     if( deltaT == -1) {
  144.     error("In findNextCollision(): couldn't find a wall to hit.",
  145.         data->time);
  146.     }
  147.  
  148.     /* Correct the equations of motion for the particle */
  149.     mol->xCoeff[1]
  150.     = (mol->collisionPos.x - lastCollision.x) / deltaT;
  151.     mol->yCoeff[1]
  152.     = (mol->collisionPos.y - lastCollision.y) / deltaT;
  153.     mol->xCoeff[0] = lastCollision.x;
  154.     mol->yCoeff[0] = lastCollision.y;
  155.  
  156.     /* Update collision info */
  157.     mol->collisionTime = deltaT;
  158.     mol->collisionWall = &data->chamber[box].walls[i];
  159.     mol->thisBox = box;
  160. }
  161.  
  162. /*
  163.  * HIT
  164.  * Return collision time if the molecule will hit this wall
  165.  *   otherwise -1.0
  166.  */
  167. static double hit( mol, wall)
  168.     Molecule *mol;
  169.     Wall *wall;
  170. {
  171.     double deltaT;
  172.     int xPos, yPos;
  173.  
  174.     /* Find intersection of trajectory and wall */
  175.  
  176.     /* Horizontal wall */
  177.     if( wall->type == TOP || wall->type == BOTTOM) {
  178.  
  179.         if( fabs( mol->yCoeff[1]) >= SMALL) {
  180.         deltaT = (wall->end[0].y-mol->collisionPos.y) / mol->yCoeff[1];
  181.  
  182.         /* Are we headed in that direction? */
  183.         if( deltaT < 0.0)
  184.         return -1.0;
  185.  
  186.         xPos = mol->collisionPos.x + mol->xCoeff[1] * deltaT;
  187.  
  188.         /* Hits wall between endpoints */
  189.         if( xPos >= wall->end[0].x && xPos <= wall->end[1].x) {
  190.         mol->collisionPos.x = xPos;
  191.         mol->collisionPos.y = wall->end[0].y;
  192.         return deltaT;
  193.         }
  194.         else return -1.0;
  195.     }
  196.     else return -1.0;
  197.     }
  198.  
  199.     /* Vertical wall */
  200.     else if ( wall->type == LEFT || wall->type == RIGHT) {
  201.         if( fabs( mol->xCoeff[1]) >= SMALL) {
  202.         deltaT = (wall->end[0].x - mol->collisionPos.x) / mol->xCoeff[1];
  203.  
  204.         /* Are we headed in that direction? */
  205.         if( deltaT < 0.0)
  206.         return -1.0;
  207.  
  208.         yPos = mol->collisionPos.y + mol->yCoeff[1] * deltaT;
  209.  
  210.         /* Hits wall between endpoints */
  211.         if( yPos >= wall->end[0].y && yPos <= wall->end[1].y) {
  212.         mol->collisionPos.x = wall->end[0].x;
  213.         mol->collisionPos.y = yPos;
  214.         return deltaT;
  215.         }
  216.         else return -1.0;
  217.     }
  218.     else return -1.0;
  219.     }
  220.     else {
  221.     error(" In hit, illegal wall type.", 0);
  222.     return -1.0;
  223.     }
  224. }
  225.  
  226. /*
  227.  * HITHOLE
  228.  * Return collision time if the molecule will go through hole
  229.  *   otherwise -1.0
  230.  */
  231. static double hitHole( mol, wall, data)
  232.     Molecule *mol;
  233.     Wall *wall;
  234.     LabData *data;
  235. {
  236.     double deltaT;
  237.     int xPos, yPos;
  238.  
  239.     /* Find intersection of trajectory and hole */
  240.  
  241.     /* OK to divide? */
  242.     if( fabs( mol->xCoeff[1]) < SMALL)
  243.     return -1.0;
  244.  
  245.     deltaT = (wall->end[0].x - mol->collisionPos.x) / mol->xCoeff[1];
  246.  
  247.     /* Are we headed in that direction? */
  248.     /* ... or already on wall or corner? use <=, not <   lm 4/4/91 */
  249.     if( deltaT <= 0.0)
  250.     return -1.0;
  251.  
  252.     /* FInd intersection with barrier wall */
  253.     yPos = mol->collisionPos.y + mol->yCoeff[1] * deltaT;
  254.  
  255.     /* Does it go through hole between endpoints? */
  256.     if( !( yPos > wall->end[0].y && yPos < wall->end[1].y))
  257.     return -1.0;
  258.  
  259.     /*
  260.      * Check that x-intersection with horizontal wall will not
  261.      *   be truncated to corner position.   This happens once
  262.      *   every 10^7 collisions.
  263.      */
  264.  
  265.     /* This isn't a problem for molecules moving in -x direction */
  266.     if( mol->xCoeff[1] < 0.0)
  267.     return deltaT;
  268.  
  269.     /* OK to divide? */
  270.     if( fabs( mol->yCoeff[1]) < SMALL)
  271.     return -1.0;
  272.  
  273.     /* Moving in +y or -y direction? */
  274.     if( mol->yCoeff[1] > 0.0)
  275.     deltaT = (data->heightMM - mol->collisionPos.y) / mol->yCoeff[1];
  276.     else
  277.     deltaT = (0 - mol->collisionPos.y) / mol->yCoeff[1];
  278.  
  279.     /* Where do we hit the horizontal wall? */
  280.     xPos = mol->collisionPos.x + mol->xCoeff[1] * deltaT;
  281.  
  282.     /*
  283.      *  If xPos truncated so that it lies at the x-position of the corner,
  284.      *    then the molecule really doesn't go through the hole.
  285.      */
  286.     if( xPos == data->chamber[0].walls[0].end[0].x)
  287.     return -1.0;
  288.  
  289.     /* It went through! */
  290.     return deltaT;
  291. }
  292.