<s> from visual import * <EOL> from time import clock <EOL> from random import random <EOL> win = <NUM_LIT> <EOL> Nstars = <NUM_LIT:20> <EOL> G = <NUM_LIT> <EOL> Msun = <NUM_LIT> <EOL> Rsun = <NUM_LIT> <EOL> L = <NUM_LIT> <EOL> vsun = <NUM_LIT> * sqrt ( G * Msun / Rsun ) <EOL> scene = display ( title = "<STR_LIT>" , width = win , height = win , <EOL> range = <NUM_LIT:2> * L , forward = ( - <NUM_LIT:1> , - <NUM_LIT:1> , - <NUM_LIT:1> ) ) <EOL> xaxis = curve ( pos = [ ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) , ( L , <NUM_LIT:0> , <NUM_LIT:0> ) ] , color = ( <NUM_LIT:0.5> , <NUM_LIT:0.5> , <NUM_LIT:0.5> ) ) <EOL> yaxis = curve ( pos = [ ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) , ( <NUM_LIT:0> , L , <NUM_LIT:0> ) ] , color = ( <NUM_LIT:0.5> , <NUM_LIT:0.5> , <NUM_LIT:0.5> ) ) <EOL> zaxis = curve ( pos = [ ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) , ( <NUM_LIT:0> , <NUM_LIT:0> , L ) ] , color = ( <NUM_LIT:0.5> , <NUM_LIT:0.5> , <NUM_LIT:0.5> ) ) <EOL> Stars = [ ] <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 ( Nstars ) : <EOL> x = - L + <NUM_LIT:2> * L * random ( ) <EOL> y = - L + <NUM_LIT:2> * L * random ( ) <EOL> z = - L + <NUM_LIT:2> * L * random ( ) <EOL> r = Rsun / <NUM_LIT:2> + Rsun * random ( ) <EOL> Stars = Stars + [ sphere ( pos = ( x , y , z ) , radius = r , color = colors [ i % <NUM_LIT:6> ] , <EOL> make_trail = True , interval = <NUM_LIT:10> ) ] <EOL> mass = Msun * r ** <NUM_LIT:3> / Rsun ** <NUM_LIT:3> <EOL> px = mass * ( - vsun + <NUM_LIT:2> * vsun * random ( ) ) <EOL> py = mass * ( - vsun + <NUM_LIT:2> * vsun * random ( ) ) <EOL> pz = mass * ( - vsun + <NUM_LIT:2> * vsun * random ( ) ) <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 = ( Nstars , <NUM_LIT:1> ) <EOL> radius = array ( rlist ) <EOL> vcm = sum ( p ) / sum ( m ) <EOL> p = p - m * vcm <EOL> dt = <NUM_LIT> <EOL> Nsteps = <NUM_LIT:0> <EOL> pos = pos + ( p / m ) * ( dt / <NUM_LIT> ) <EOL> time = clock ( ) <EOL> Nhits = <NUM_LIT:0> <EOL> while True : <EOL> rate ( <NUM_LIT:100> ) <EOL> r = pos - pos [ : , newaxis ] <EOL> for n in range ( Nstars ) : <EOL> r [ n , n ] = <NUM_LIT> <EOL> rmag = sqrt ( sum ( square ( r ) , - <NUM_LIT:1> ) ) <EOL> hit = less_equal ( rmag , radius + radius [ : , newaxis ] ) - identity ( Nstars ) <EOL> hitlist = sort ( nonzero ( hit . flat ) [ <NUM_LIT:0> ] ) . tolist ( ) <EOL> F = G * m * m [ : , newaxis ] * r / rmag [ : , : , newaxis ] ** <NUM_LIT:3> <EOL> for n in range ( Nstars ) : <EOL> F [ n , n ] = <NUM_LIT:0> <EOL> p = p + sum ( F , <NUM_LIT:1> ) * dt <EOL> pos = pos + ( p / m ) * dt <EOL> for i in range ( Nstars ) : <EOL> Stars [ i ] . pos = pos [ i ] <EOL> for ij in hitlist : <EOL> i , j = divmod ( ij , Nstars ) <EOL> if not Stars [ i ] . visible : continue <EOL> if not Stars [ j ] . visible : continue <EOL> newpos = ( pos [ i ] * m [ i , <NUM_LIT:0> ] + pos [ j ] * m [ j , <NUM_LIT:0> ] ) / ( m [ i , <NUM_LIT:0> ] + m [ j , <NUM_LIT:0> ] ) <EOL> newmass = m [ i , <NUM_LIT:0> ] + m [ j , <NUM_LIT:0> ] <EOL> newp = p [ i ] + p [ j ] <EOL> newradius = Rsun * ( ( newmass / Msun ) ** ( <NUM_LIT:1.> / <NUM_LIT> ) ) <EOL> iset , jset = i , j <EOL> if radius [ j ] > radius [ i ] : <EOL> iset , jset = j , i <EOL> Stars [ iset ] . radius = newradius <EOL> m [ iset , <NUM_LIT:0> ] = newmass <EOL> pos [ iset ] = newpos <EOL> p [ iset ] = newp <EOL> Stars [ jset ] . trail_object . visible = False <EOL> Stars [ jset ] . visible = <NUM_LIT:0> <EOL> p [ jset ] = vector ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) <EOL> m [ jset , <NUM_LIT:0> ] = Msun * <NUM_LIT> <EOL> Nhits = Nhits + <NUM_LIT:1> <EOL> pos [ jset ] = ( <NUM_LIT> * L * Nhits , <NUM_LIT:0> , <NUM_LIT:0> ) <EOL> Nsteps += <NUM_LIT:1> </s>
