Wiki

RunTimePerformance: cobra_nbody (initial).cobra

File cobra_nbody (initial).cobra, 5.7 KB (added by torial, 9 months ago)

Initial Cobra version (unoptimized)

Line 
1# The Computer Language Benchmarks Game
2# http://shootout.alioth.debian.org/
3#
4# originally by Kevin Carson
5# modified by Tupteq, Fredrik Johansson, and Daniel Nanz
6# modified by Maciej Fijalkowski
7# 2to3
8#@number float
9class Triple
10    pro x from var as decimal 
11    pro y from var as decimal 
12    pro z from var as decimal
13    get squared as decimal
14        return .x*.x + .y*.y + .z*.z
15   
16    cue init(x as decimal, y as decimal, z as decimal)
17        base.init
18        .x = x
19        .y = y
20        .z = z
21       
22    def toString as String is override
23        return "\[[.x], [.y], [.z]\]"
24
25class Body
26    pro coord from var as Triple
27    pro velocity from var as Triple
28    pro mass from var as decimal
29   
30    cue init(coord as Triple, velocity as Triple, mass as decimal)
31        base.init
32        .coord = coord
33        .velocity = velocity
34        .mass = mass 
35       
36class BodyPair
37    pro body1 from var as Body
38    pro body2 from var as Body 
39   
40    cue init(body1 as Body, body2 as Body)
41        base.init
42        .body1 = body1
43        .body2 = body2
44
45class Program
46    def combinations(l) as List<of BodyPair>
47        result = List<of BodyPair>()
48        for x in l.count - 1
49            ls = l[x+1:]
50            for y in ls
51                result.add(BodyPair(l[x],y))
52        return result
53
54    var bODIES
55    cue init
56        base.init
57   
58        .bODIES = {
59            'sun': Body(Triple(0.0, 0.0, 0.0), Triple(0.0, 0.0, 0.0), .sOLAR_MASS),
60
61            'jupiter': Body(Triple(4.84143144246472090,
62                -1.16032004402742839,
63                -0.103622044471123109),
64                Triple(0.00166007664274403694 * .dAYS_PER_YEAR,
65                0.00769901118419740425 * .dAYS_PER_YEAR,
66                -0.0000690460016972063023 * .dAYS_PER_YEAR),
67                0.000954791938424326609 * .sOLAR_MASS),
68
69            'saturn': Body(Triple(8.34336671824457987,
70                4.12479856412430479,
71                -0.403523417114321381),
72                Triple(-0.00276742510726862411 * .dAYS_PER_YEAR,
73                0.00499852801234917238 * .dAYS_PER_YEAR,
74                0.0000230417297573763929 * .dAYS_PER_YEAR),
75                0.000285885980666130812 * .sOLAR_MASS),
76
77            'uranus': Body(Triple(12.8943695621391310,
78                -15.1111514016986312,
79                -0.223307578892655734),
80                Triple(0.00296460137564761618* .dAYS_PER_YEAR,
81                0.00237847173959480950* .dAYS_PER_YEAR,
82                -0.0000296589568540237556* .dAYS_PER_YEAR),
83                0.0000436624404335156298* .sOLAR_MASS),
84
85            'neptune': Body(Triple(15.3796971148509165,
86                -25.9193146099879641,
87                0.179258772950371181),
88                Triple(0.00268067772490389322* .dAYS_PER_YEAR,
89                0.00162824170038242295* .dAYS_PER_YEAR,
90                -0.0000951592254519715870* .dAYS_PER_YEAR),
91                0.0000515138902046611451* .sOLAR_MASS)
92        }
93        .sYSTEM = [.bODIES['sun'], 
94                    .bODIES['jupiter'],
95                    .bODIES['saturn'],
96                    .bODIES['uranus'],
97                    .bODIES['neptune']
98        ] 
99        .pAIRS = .combinations(.sYSTEM)
100    var sYSTEM
101    var _pI = 3.14159265358979323
102    var sOLAR_MASS = 4 * 3.14159265358979323 *  3.14159265358979323
103    var dAYS_PER_YEAR = 365.24
104
105    var pAIRS
106
107
108    def advance(dt as decimal, n as int)
109        bodies=.sYSTEM
110        pairs=.pAIRS
111
112        for i in n
113            for p in pairs
114                x1, y1, z1 = p.body1.coord.x, p.body1.coord.y, p.body1.coord.z
115                x2, y2, z2 = p.body2.coord.x, p.body2.coord.y, p.body2.coord.z
116                m1, m2 = p.body1.mass, p.body2.mass
117                v1, v2 =p.body1.velocity, p.body2.velocity
118                dx = x1 - x2
119                dy = y1 - y2
120                dz = z1 - z2
121                mag = dt * (Math.pow(Convert.toDouble(dx * dx + dy * dy + dz * dz), -1.5f) to decimal)
122                b1m = m1 * mag
123                b2m = m2 * mag
124                v1.x = v1.x - dx * b2m
125                v1.y = v1.y - dy * b2m
126                v1.z = v1.z - dz * b2m
127                v2.x = v2.x + dx * b1m
128                v2.y = v2.y + dy * b1m
129                v2.z = v2.z + dz * b1m
130            for b in bodies
131                r = b.coord
132                vx,vy,vz = b.velocity.x, b.velocity.y, b.velocity.z
133                r.x = r.x + dt * vx
134                r.y = r.y + dt * vy
135                r.z = r.z + dt * vz
136
137
138    def report_energy
139        bodies=.sYSTEM
140        pairs=.pAIRS
141        #print "Bodies Count:[bodies.count]\nPairs Count:[pairs.count]"
142        e=0.0
143        for p in pairs
144            m = p.body1.mass * p.body2.mass
145            #print m
146            dx = p.body1.coord.x - p.body2.coord.x
147            dy = p.body1.coord.y - p.body2.coord.y
148            dz = p.body1.coord.z - p.body2.coord.z
149            sqrs = Convert.toDouble(dx * dx + dy * dy + dz * dz)
150            e = e - (m / Math.pow(sqrs , 0.5f) to decimal)
151        for b in bodies
152            e = e + b.mass * b.velocity.squared / 2.0
153        print e
154
155    def offset_momentum(ref1 )
156        bodies=.sYSTEM
157        px=0.0 
158        py=0.0
159        pz=0.0
160
161        for b in bodies
162            vx ,vy ,vz  = b.velocity.x, b.velocity.y,b.velocity.z
163            m = b.mass
164            px = px - vx * m
165            py =py - vy * m
166            pz = pz - vz * m
167        v, m = ref1.velocity, ref1.mass
168        v.x = px / m
169        v.y = py / m
170        v.z = pz / m
171        ref1.velocity = v
172        print v
173
174    def main
175        sw = System.Diagnostics.Stopwatch()
176        sw.start
177        ref1 = 'sun'
178        .report_energy
179        .offset_momentum(  .bODIES[ref1])
180        .report_energy
181        .advance(0.01, 50000)
182        .report_energy
183        sw.stop
184        print 'duration: [sw.elapsedMilliseconds] ms '
185