Pastey nbody update
- Posted by mattlewis (admin) Oct 05, 2012
-- The Computer Language Shootout Benchmarks -- http://shootout.alioth.debian.org/ -- -- Converted to Euphoria by Jason Gade -- Optimized by Matt Lewis -- run: eui nbody.ex N without warning without type_check include get.e constant PI = 3.141592653589793, SOLAR_MASS = 4 * PI * PI, DAYS_PER_YEAR = 365.24 -- struct planet constant name = 1, x = 2, y = 3, z = 4, vx = 5, vy = 6, vz = 7, mass = 8 -- end struct sequence bodies bodies = { { "Sun", 0, 0, 0, 0, 0, 0, SOLAR_MASS}, { "Jupiter", 4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01, 1.66007664274403694e-03 * DAYS_PER_YEAR, 7.69901118419740425e-03 * DAYS_PER_YEAR, -6.90460016972063023e-05 * DAYS_PER_YEAR, 9.54791938424326609e-04 * SOLAR_MASS }, { "Saturn", 8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01, -2.76742510726862411e-03 * DAYS_PER_YEAR, 4.99852801234917238e-03 * DAYS_PER_YEAR, 2.30417297573763929e-05 * DAYS_PER_YEAR, 2.85885980666130812e-04 * SOLAR_MASS }, { "Uranus", 1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01, 2.96460137564761618e-03 * DAYS_PER_YEAR, 2.37847173959480950e-03 * DAYS_PER_YEAR, -2.96589568540237556e-05 * DAYS_PER_YEAR, 4.36624404335156298e-05 * SOLAR_MASS }, { "Neptune", 1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01, 2.68067772490389322e-03 * DAYS_PER_YEAR, 1.62824170038242295e-03 * DAYS_PER_YEAR, -9.51592254519715870e-05 * DAYS_PER_YEAR, 5.15138902046611451e-05 * SOLAR_MASS } } constant NBODIES = length(bodies) sequence NAME = {}, X = {}, Y = {}, Z = {}, VX = {}, VY = {}, VZ = {}, MASS = {} for i = 1 to NBODIES do for f = name to mass do switch f do case name then NAME = append( NAME, bodies[i][f] ) case x then X &= bodies[i][f] case y then Y &= bodies[i][f] case z then Z &= bodies[i][f] case vx then VX &= bodies[i][f] case vy then VY &= bodies[i][f] case vz then VZ &= bodies[i][f] case mass then MASS &= bodies[i][f] end switch end for end for procedure advance(atom dt) atom dx, dy, dz, distance, mag, mass, distance_2, mass_mag for i = 1 to NBODIES do for j = i + 1 to NBODIES do dx = X[i] - X[j] dy = Y[i] - Y[j] dz = Z[i] - Z[j] distance_2 = dx*dx distance_2 += dy*dy distance_2 += dz*dz distance = sqrt(distance_2) distance *= distance_2 mag = dt / distance mass = MASS[j] mass_mag = mass * mag VX[i] -= dx * mass_mag VY[i] -= dy * mass_mag VZ[i] -= dz * mass_mag mass = MASS[i] mass_mag = mass * mag VX[j] += dx * mass_mag VY[j] += dy * mass_mag VZ[j] += dz * mass_mag end for X[i] += dt * VX[i] Y[i] += dt * VY[i] Z[i] += dt * VZ[i] end for end procedure -- advance function energy() atom e, dx, dy, dz, distance, mass_i e = 0.0 for i = 1 to NBODIES do mass_i = MASS[i] e += 0.5 * mass_i * (VX[i]*VX[i] + VY[i]*VY[i] + VZ[i]*VZ[i]) for j = i + 1 to NBODIES do dx = X[i] - X[j] dy = Y[i] - Y[j] dz = Z[i] - Z[j] distance = sqrt(dx*dx + dy*dy + dz*dz) e -= (mass_i*MASS[j])/distance end for end for return e end function -- energy procedure offset_momentum() atom px, py, pz px = 0.0 py = 0.0 pz = 0.0 for i = 1 to NBODIES do px += VX[i] * MASS[i] py += VY[i] * MASS[i] pz += VZ[i] * MASS[i] end for VX[1] = - px / SOLAR_MASS VY[1] = - py / SOLAR_MASS VZ[1] = - pz / SOLAR_MASS end procedure -- offset_momentum procedure main(sequence argv) object n if length(argv) > 2 then n = value(argv[3]) n = n[2] else n = 1000 end if offset_momentum() printf(1, "%.9f\n", energy()) for i = 1 to n do advance(0.01) end for printf(1, "%.9f\n", energy()) end procedure -- main main(command_line())


