<s> from visual import * <EOL> from visual . graph import * <EOL> from random import random <EOL> win = <NUM_LIT> <EOL> Natoms = <NUM_LIT:100> <EOL> L = <NUM_LIT:1.> <EOL> gray = ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) <EOL> Matom = <NUM_LIT> / <NUM_LIT> <EOL> Ratom = <NUM_LIT> <EOL> k = <NUM_LIT> <EOL> T = <NUM_LIT> <EOL> dt = <NUM_LIT> <EOL> scene = display ( title = "<STR_LIT>" , width = win , height = win , x = <NUM_LIT:0> , y = <NUM_LIT:0> , <EOL> center = ( L / <NUM_LIT> , L / <NUM_LIT> , L / <NUM_LIT> ) ) <EOL> deltav = <NUM_LIT> <EOL> vdist = gdisplay ( x = <NUM_LIT:0> , y = win , ymax = Natoms * deltav / <NUM_LIT> , <EOL> width = win , height = <NUM_LIT> * win , xtitle = '<STR_LIT:v>' , ytitle = '<STR_LIT>' ) <EOL> theory = gcurve ( color = color . cyan ) <EOL> dv = <NUM_LIT> <EOL> for v in arange ( <NUM_LIT:0.> , <NUM_LIT> + dv , dv ) : <EOL> theory . plot ( pos = ( v , <EOL> ( deltav / dv ) * Natoms * <NUM_LIT> * pi * ( ( Matom / ( <NUM_LIT> * pi * k * T ) ) ** <NUM_LIT> ) <EOL> * exp ( ( - <NUM_LIT:0.5> * Matom * v ** <NUM_LIT:2> ) / ( k * T ) ) * v ** <NUM_LIT:2> * dv ) ) <EOL> observation = ghistogram ( bins = arange ( <NUM_LIT:0.> , <NUM_LIT> , deltav ) , <EOL> accumulate = <NUM_LIT:1> , average = <NUM_LIT:1> , color = color . red ) <EOL> xaxis = curve ( pos = [ ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) , ( L , <NUM_LIT:0> , <NUM_LIT:0> ) ] , color = gray ) <EOL> yaxis = curve ( pos = [ ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) , ( <NUM_LIT:0> , L , <NUM_LIT:0> ) ] , color = gray ) <EOL> zaxis = curve ( pos = [ ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) , ( <NUM_LIT:0> , <NUM_LIT:0> , L ) ] , color = gray ) <EOL> xaxis2 = curve ( pos = [ ( L , L , L ) , ( <NUM_LIT:0> , L , L ) , ( <NUM_LIT:0> , <NUM_LIT:0> , L ) , ( L , <NUM_LIT:0> , L ) ] , color = gray ) <EOL> yaxis2 = curve ( pos = [ ( L , L , L ) , ( L , <NUM_LIT:0> , L ) , ( L , <NUM_LIT:0> , <NUM_LIT:0> ) , ( L , L , <NUM_LIT:0> ) ] , color = gray ) <EOL> zaxis2 = curve ( pos = [ ( L , L , L ) , ( L , L , <NUM_LIT:0> ) , ( <NUM_LIT:0> , L , <NUM_LIT:0> ) , ( <NUM_LIT:0> , L , L ) ] , color = gray ) <EOL> Atoms = [ ] <EOL> colors = [ color . red , color . green , color . blue , <EOL> color . yellow , color . cyan , color . magenta ] <EOL> poslist = [ ] <EOL> plist = [ ] <EOL> mlist = [ ] <EOL> rlist = [ ] <EOL> for i in range ( Natoms ) : <EOL> Lmin = <NUM_LIT> * Ratom <EOL> Lmax = L - Lmin <EOL> x = Lmin + ( Lmax - Lmin ) * random ( ) <EOL> y = Lmin + ( Lmax - Lmin ) * random ( ) <EOL> z = Lmin + ( Lmax - Lmin ) * random ( ) <EOL> r = Ratom <EOL> Atoms = Atoms + [ sphere ( pos = ( x , y , z ) , radius = r , color = colors [ i % <NUM_LIT:6> ] ) ] <EOL> mass = Matom * r ** <NUM_LIT:3> / Ratom ** <NUM_LIT:3> <EOL> pavg = sqrt ( <NUM_LIT> * mass * <NUM_LIT> * k * T ) <EOL> theta = pi * random ( ) <EOL> phi = <NUM_LIT:2> * pi * random ( ) <EOL> px = pavg * sin ( theta ) * cos ( phi ) <EOL> py = pavg * sin ( theta ) * sin ( phi ) <EOL> pz = pavg * cos ( theta ) <EOL> poslist . append ( ( x , y , z ) ) <EOL> plist . append ( ( px , py , pz ) ) <EOL> mlist . append ( mass ) <EOL> rlist . append ( r ) <EOL> pos = array ( poslist ) <EOL> p = array ( plist ) <EOL> m = array ( mlist ) <EOL> m . shape = ( Natoms , <NUM_LIT:1> ) <EOL> radius = array ( rlist ) <EOL> pos = pos + ( p / m ) * ( dt / <NUM_LIT> ) <EOL> while True : <EOL> rate ( <NUM_LIT:100> ) <EOL> observation . plot ( data = mag ( p / m ) ) <EOL> pos = pos + ( p / m ) * dt <EOL> r = pos - pos [ : , newaxis ] <EOL> rmag = sqrt ( sum ( square ( r ) , - <NUM_LIT:1> ) ) <EOL> hit = less_equal ( rmag , radius + radius [ : , None ] ) - identity ( Natoms ) <EOL> hitlist = sort ( nonzero ( hit . flat ) [ <NUM_LIT:0> ] ) . tolist ( ) <EOL> for ij in hitlist : <EOL> i , j = divmod ( ij , Natoms ) <EOL> hitlist . remove ( j * Natoms + i ) <EOL> ptot = p [ i ] + p [ j ] <EOL> mi = m [ i , <NUM_LIT:0> ] <EOL> mj = m [ j , <NUM_LIT:0> ] <EOL> vi = p [ i ] / mi <EOL> vj = p [ j ] / mj <EOL> ri = Atoms [ i ] . radius <EOL> rj = Atoms [ j ] . radius <EOL> a = mag ( vj - vi ) ** <NUM_LIT:2> <EOL> if a == <NUM_LIT:0> : continue <EOL> b = <NUM_LIT:2> * dot ( pos [ i ] - pos [ j ] , vj - vi ) <EOL> c = mag ( pos [ i ] - pos [ j ] ) ** <NUM_LIT:2> - ( ri + rj ) ** <NUM_LIT:2> <EOL> d = b ** <NUM_LIT:2> - <NUM_LIT> * a * c <EOL> if d < <NUM_LIT:0> : continue <EOL> deltat = ( - b + sqrt ( d ) ) / ( <NUM_LIT> * a ) <EOL> pos [ i ] = pos [ i ] - ( p [ i ] / mi ) * deltat <EOL> pos [ j ] = pos [ j ] - ( p [ j ] / mj ) * deltat <EOL> mtot = mi + mj <EOL> pcmi = p [ i ] - ptot * mi / mtot <EOL> pcmj = p [ j ] - ptot * mj / mtot <EOL> rrel = norm ( pos [ j ] - pos [ i ] ) <EOL> pcmi = pcmi - <NUM_LIT:2> * dot ( pcmi , rrel ) * rrel <EOL> pcmj = pcmj - <NUM_LIT:2> * dot ( pcmj , rrel ) * rrel <EOL> p [ i ] = pcmi + ptot * mi / mtot <EOL> p [ j ] = pcmj + ptot * mj / mtot <EOL> pos [ i ] = pos [ i ] + ( p [ i ] / mi ) * deltat <EOL> pos [ j ] = pos [ j ] + ( p [ j ] / mj ) * deltat <EOL> outside = less_equal ( pos , Ratom ) <EOL> p1 = p * outside <EOL> p = p - p1 + abs ( p1 ) <EOL> outside = greater_equal ( pos , L - Ratom ) <EOL> p1 = p * outside <EOL> p = p - p1 - abs ( p1 ) <EOL> for i in range ( Natoms ) : <EOL> Atoms [ i ] . pos = pos [ i ] </s>
