RunTimePerformance: cobra_nbody (initial).cobra

# The Computer Language Benchmarks Game
# http://shootout.alioth.debian.org/
#
# originally by Kevin Carson
# modified by Tupteq, Fredrik Johansson, and Daniel Nanz
# modified by Maciej Fijalkowski
# 2to3
#@number float
class Triple
    pro x from var as decimal 
    pro y from var as decimal 
    pro z from var as decimal
    get squared as decimal
        return .x*.x + .y*.y + .z*.z
    
    cue init(x as decimal, y as decimal, z as decimal)
        base.init
        .x = x
        .y = y
        .z = z
        
    def toString as String is override
        return "\[[.x], [.y], [.z]\]"

class Body
    pro coord from var as Triple
    pro velocity from var as Triple
    pro mass from var as decimal
    
    cue init(coord as Triple, velocity as Triple, mass as decimal)
        base.init
        .coord = coord
        .velocity = velocity
        .mass = mass  
        
class BodyPair
    pro body1 from var as Body
    pro body2 from var as Body  
    
    cue init(body1 as Body, body2 as Body)
        base.init
        .body1 = body1
        .body2 = body2

class Program
    def combinations(l) as List<of BodyPair>
        result = List<of BodyPair>()
        for x in l.count - 1
            ls = l[x+1:]
            for y in ls
                result.add(BodyPair(l[x],y))
        return result

    var bODIES
    cue init
        base.init
    
        .bODIES = {
            'sun': Body(Triple(0.0, 0.0, 0.0), Triple(0.0, 0.0, 0.0), .sOLAR_MASS),

            'jupiter': Body(Triple(4.84143144246472090,
                -1.16032004402742839,
                -0.103622044471123109),
                Triple(0.00166007664274403694 * .dAYS_PER_YEAR,
                0.00769901118419740425 * .dAYS_PER_YEAR,
                -0.0000690460016972063023 * .dAYS_PER_YEAR),
                0.000954791938424326609 * .sOLAR_MASS),

            'saturn': Body(Triple(8.34336671824457987,
                4.12479856412430479,
                -0.403523417114321381),
                Triple(-0.00276742510726862411 * .dAYS_PER_YEAR,
                0.00499852801234917238 * .dAYS_PER_YEAR,
                0.0000230417297573763929 * .dAYS_PER_YEAR),
                0.000285885980666130812 * .sOLAR_MASS),

            'uranus': Body(Triple(12.8943695621391310,
                -15.1111514016986312,
                -0.223307578892655734),
                Triple(0.00296460137564761618* .dAYS_PER_YEAR,
                0.00237847173959480950* .dAYS_PER_YEAR,
                -0.0000296589568540237556* .dAYS_PER_YEAR),
                0.0000436624404335156298* .sOLAR_MASS),

            'neptune': Body(Triple(15.3796971148509165,
                -25.9193146099879641,
                0.179258772950371181),
                Triple(0.00268067772490389322* .dAYS_PER_YEAR,
                0.00162824170038242295* .dAYS_PER_YEAR,
                -0.0000951592254519715870* .dAYS_PER_YEAR),
                0.0000515138902046611451* .sOLAR_MASS) 
        }
        .sYSTEM = [.bODIES['sun'], 
                    .bODIES['jupiter'],
                    .bODIES['saturn'],
                    .bODIES['uranus'],
                    .bODIES['neptune']
        ] 
        .pAIRS = .combinations(.sYSTEM)
    var sYSTEM 
    var _pI = 3.14159265358979323
    var sOLAR_MASS = 4 * 3.14159265358979323 *  3.14159265358979323
    var dAYS_PER_YEAR = 365.24

    var pAIRS


    def advance(dt as decimal, n as int)
        bodies=.sYSTEM
        pairs=.pAIRS

        for i in n
            for p in pairs
                x1, y1, z1 = p.body1.coord.x, p.body1.coord.y, p.body1.coord.z
                x2, y2, z2 = p.body2.coord.x, p.body2.coord.y, p.body2.coord.z
                m1, m2 = p.body1.mass, p.body2.mass
                v1, v2 =p.body1.velocity, p.body2.velocity
                dx = x1 - x2
                dy = y1 - y2
                dz = z1 - z2
                mag = dt * (Math.pow(Convert.toDouble(dx * dx + dy * dy + dz * dz), -1.5f) to decimal)
                b1m = m1 * mag
                b2m = m2 * mag
                v1.x = v1.x - dx * b2m
                v1.y = v1.y - dy * b2m
                v1.z = v1.z - dz * b2m
                v2.x = v2.x + dx * b1m
                v2.y = v2.y + dy * b1m
                v2.z = v2.z + dz * b1m
            for b in bodies
                r = b.coord
                vx,vy,vz = b.velocity.x, b.velocity.y, b.velocity.z
                r.x = r.x + dt * vx
                r.y = r.y + dt * vy
                r.z = r.z + dt * vz


    def report_energy
        bodies=.sYSTEM
        pairs=.pAIRS
        #print "Bodies Count:[bodies.count]\nPairs Count:[pairs.count]"
        e=0.0
        for p in pairs
            m = p.body1.mass * p.body2.mass
            #print m
            dx = p.body1.coord.x - p.body2.coord.x
            dy = p.body1.coord.y - p.body2.coord.y
            dz = p.body1.coord.z - p.body2.coord.z
            sqrs = Convert.toDouble(dx * dx + dy * dy + dz * dz) 
            e = e - (m / Math.pow(sqrs , 0.5f) to decimal)
        for b in bodies
            e = e + b.mass * b.velocity.squared / 2.0
        print e

    def offset_momentum(ref1 )
        bodies=.sYSTEM
        px=0.0 
        py=0.0
        pz=0.0

        for b in bodies
            vx ,vy ,vz  = b.velocity.x, b.velocity.y,b.velocity.z
            m = b.mass
            px = px - vx * m
            py =py - vy * m
            pz = pz - vz * m
        v, m = ref1.velocity, ref1.mass
        v.x = px / m
        v.y = py / m
        v.z = pz / m
        ref1.velocity = v
        print v

    def main
        sw = System.Diagnostics.Stopwatch()
        sw.start
        ref1 = 'sun'
        .report_energy
        .offset_momentum(  .bODIES[ref1])
        .report_energy
        .advance(0.01, 50000)
        .report_energy
        sw.stop
        print 'duration: [sw.elapsedMilliseconds] ms '