Commit | Line | Data |
---|---|---|
7f918cf1 CE |
1 | (* From the SML/NJ benchmark suite. *) |
2 | (* vector-sig.sml | |
3 | * | |
4 | * COPYRIGHT (c) 1993, AT&T Bell Laboratories. | |
5 | * | |
6 | * The abstract interface of vectors and matrices in some dimension. | |
7 | *) | |
8 | ||
9 | signature VECTOR = | |
10 | sig | |
11 | type 'a vec | |
12 | ||
13 | val dim : int (* dimension of the vectors *) | |
14 | ||
15 | val tabulate : (int -> 'a) -> 'a vec | |
16 | ||
17 | val equal : real vec * real vec -> bool | |
18 | val zerov : real vec | |
19 | val addv : (real vec * real vec) -> real vec | |
20 | val subv : (real vec * real vec) -> real vec | |
21 | val dotvp : (real vec * real vec) -> real | |
22 | val crossvp : (real vec * real vec) -> real vec | |
23 | val addvs : (real vec * real) -> real vec | |
24 | val mulvs : (real vec * real) -> real vec | |
25 | val divvs : (real vec * real) -> real vec | |
26 | ||
27 | val mapv : ('a -> 'b) -> 'a vec -> 'b vec | |
28 | val map3v : (('a * 'b * 'c) -> 'd) -> ('a vec * 'b vec * 'c vec) -> 'd vec | |
29 | val foldv : ('a * 'b -> 'b) -> 'a vec -> 'b -> 'b | |
30 | val format : {lp : string, sep : string, rp : string, cvt : 'a -> string} | |
31 | -> 'a vec -> string | |
32 | val explode : 'a vec -> 'a list | |
33 | val implode : 'a list -> 'a vec | |
34 | ||
35 | type matrix (* matrices are always real valued *) | |
36 | ||
37 | val zerom : matrix | |
38 | val addm : (matrix * matrix) -> matrix | |
39 | val outvp : (real vec * real vec) -> matrix | |
40 | ||
41 | end | |
42 | (* space.sml | |
43 | * | |
44 | * COPYRIGHT (c) 1993, AT&T Bell Laboratories. | |
45 | * | |
46 | * The quad/oct-tree representation of space. | |
47 | *) | |
48 | ||
49 | signature SPACE = | |
50 | sig | |
51 | ||
52 | structure V : VECTOR | |
53 | ||
54 | datatype body = Body of { | |
55 | mass : real, | |
56 | pos : real V.vec ref, | |
57 | vel : real V.vec ref, | |
58 | acc : real V.vec ref, | |
59 | phi : real ref | |
60 | } | |
61 | ||
62 | datatype cell | |
63 | = BodyCell of body | |
64 | | Cell of node Array.array | |
65 | ||
66 | and node | |
67 | = Empty | |
68 | | Node of { | |
69 | mass : real ref, | |
70 | pos : real V.vec ref, | |
71 | cell : cell | |
72 | } | |
73 | ||
74 | datatype space = Space of { | |
75 | rmin : real V.vec, | |
76 | rsize : real, | |
77 | root : node | |
78 | } | |
79 | ||
80 | val nsub : int (* number of sub cells / cell (2 ^ V.dim) *) | |
81 | ||
82 | val putCell : (cell * int * node) -> unit | |
83 | val getCell : (cell * int) -> node | |
84 | val mkCell : unit -> cell | |
85 | val mkBodyNode : body -> node | |
86 | val mkCellNode : cell -> node | |
87 | val eqBody : body * body -> bool | |
88 | ||
89 | (* debugging code *) | |
90 | val dumpTree : node -> unit | |
91 | val prBody : body -> string | |
92 | val prNode : node -> string | |
93 | ||
94 | end; (* SPACE *) | |
95 | ||
96 | functor Space (V : VECTOR) : SPACE = | |
97 | struct | |
98 | ||
99 | structure V = V | |
100 | ||
101 | datatype body = Body of { | |
102 | mass : real, | |
103 | pos : real V.vec ref, | |
104 | vel : real V.vec ref, | |
105 | acc : real V.vec ref, | |
106 | phi : real ref | |
107 | } | |
108 | ||
109 | datatype cell | |
110 | = BodyCell of body | |
111 | | Cell of node Array.array | |
112 | ||
113 | and node | |
114 | = Empty | |
115 | | Node of { | |
116 | mass : real ref, | |
117 | pos : real V.vec ref, | |
118 | cell : cell | |
119 | } | |
120 | ||
121 | datatype space = Space of { | |
122 | rmin : real V.vec, | |
123 | rsize : real, | |
124 | root : node | |
125 | } | |
126 | ||
127 | fun eqBody(Body{mass,pos,vel,acc,phi}, | |
128 | Body{mass=m1,pos=p1,vel=v1,acc=a1,phi=h1}) = | |
129 | (Real.==(mass, m1) andalso Real.==(!phi, !h1) | |
130 | andalso V.equal(!pos, !p1) andalso V.equal(!vel, !v1) | |
131 | andalso V.equal(!acc, !a1)) | |
132 | ||
133 | (* number of sub cells per cell (2 ^ V.dim) *) | |
134 | val nsub = Word.toInt(Word.<<(0w1, Word.fromInt V.dim)) | |
135 | ||
136 | fun putCell (Cell a, i, nd) = Array.update(a, i, nd) | |
137 | fun getCell (Cell a, i) = Array.sub(a, i) | |
138 | fun mkCell () = Cell(Array.array(nsub, Empty)) | |
139 | fun mkBodyNode (body as Body{pos, mass, ...}) = Node{ | |
140 | cell = BodyCell body, | |
141 | mass = ref mass, | |
142 | pos = ref (!pos) | |
143 | } | |
144 | fun mkCellNode cell = Node{cell = cell, mass = ref 0.0, pos = ref V.zerov} | |
145 | ||
146 | (* debugging code *) | |
147 | local | |
148 | val rfmt = Real.toString | |
149 | val vfmt = V.format{lp="[", rp="]", sep=",", cvt = rfmt} | |
150 | in | |
151 | fun prBody (Body{mass, pos, vel, acc, phi}) = String.concat [ | |
152 | "B{m=", rfmt mass, | |
153 | ", p=", vfmt(!pos), | |
154 | ", v=", vfmt(!vel), | |
155 | ", a=", vfmt(!acc), | |
156 | ", phi=", rfmt(!phi), "}" | |
157 | ] | |
158 | fun prNode Empty = "Empty" | |
159 | | prNode (Node{mass, pos, cell}) = let | |
160 | val cell = (case cell | |
161 | of (Cell _) => "Cell" | |
162 | | (BodyCell b) => (*prBody b*) "Body" | |
163 | (* end case *)) | |
164 | in | |
165 | String.concat [ | |
166 | "N{m=", rfmt(!mass), | |
167 | ", p=", vfmt(!pos), | |
168 | cell, "}" | |
169 | ] | |
170 | end | |
171 | end | |
172 | ||
173 | fun dumpTree tree = let | |
174 | fun printf items = TextIO.output(TextIO.stdOut, String.concat items) | |
175 | fun indent i = StringCvt.padLeft #" " (i+1) "" | |
176 | fun dump (node, l) = let | |
177 | fun dump' (Node{cell=Cell a, ...}) = let | |
178 | fun dump'' i = (dump(Array.sub(a, i), l+1); dump''(i+1)) | |
179 | in | |
180 | (dump'' 0) handle _ => () | |
181 | end | |
182 | | dump' _ = () | |
183 | in | |
184 | printf [ | |
185 | StringCvt.padLeft #" " 2 (Int.toString l), | |
186 | indent l, | |
187 | prNode node, "\n" | |
188 | ]; | |
189 | dump' node | |
190 | end | |
191 | in | |
192 | dump (tree, 0) | |
193 | end | |
194 | ||
195 | end; (* Space *) | |
196 | ||
197 | (* load.sml | |
198 | * | |
199 | * COPYRIGHT (c) 1993, AT&T Bell Laboratories. | |
200 | * | |
201 | * Code to build the tree from a list of bodies. | |
202 | *) | |
203 | ||
204 | signature LOAD = | |
205 | sig | |
206 | ||
207 | structure S : SPACE | |
208 | structure V : VECTOR | |
209 | sharing S.V = V | |
210 | ||
211 | val makeTree : (S.body list * real V.vec * real) -> S.space | |
212 | ||
213 | end; (* LOAD *) | |
214 | ||
215 | functor Load (S : SPACE) : LOAD = | |
216 | struct | |
217 | ||
218 | structure S = S | |
219 | structure V = S.V | |
220 | ||
221 | exception NotIntCoord | |
222 | ||
223 | fun rshift (n, k) = Word.toInt(Word.~>>(Word.fromInt n, Word.fromInt k)) | |
224 | ||
225 | val IMAX = 0x20000000 (* 2^29 *) | |
226 | val IMAXrs1 = rshift(IMAX, 1) | |
227 | val rIMAX = real IMAX | |
228 | ||
229 | (* compute integerized coordinates. Raises the NotIntCoord exception, | |
230 | * if rp is out of bounds. | |
231 | *) | |
232 | fun intcoord (rp, rmin, rsize) = let | |
233 | val xsc = V.divvs (V.subv(rp, rmin), rsize) | |
234 | fun cvt x = if ((0.0 <= x) andalso (x < 1.0)) | |
235 | then floor(rIMAX * x) | |
236 | else raise NotIntCoord | |
237 | in | |
238 | V.mapv cvt xsc | |
239 | end | |
240 | ||
241 | (* determine which subcell to select. *) | |
242 | fun subindex (iv, l) = let | |
243 | fun aux (v, (i, k)) = if (Word.andb(Word.fromInt v, Word.fromInt l) <> 0w0) | |
244 | then (i + rshift(S.nsub, k+1), k+1) | |
245 | else (i, k+1) | |
246 | in | |
247 | #1 (V.foldv aux iv (0, 0)) | |
248 | end | |
249 | ||
250 | (* enlarge cubical "box", salvaging existing tree structure. *) | |
251 | fun expandBox (nd as S.Body{pos, ...}, box as S.Space{rmin, rsize, root}) = ( | |
252 | (intcoord (!pos, rmin, rsize); box) | |
253 | handle NotIntCoord => let | |
254 | val rmid = V.addvs (rmin, 0.5 * rsize) | |
255 | val rmin' = V.map3v (fn (x,y,z) => | |
256 | if x < y then z - rsize else z) (!pos, rmid, rmin) | |
257 | val rsize' = 2.0 * rsize | |
258 | fun mksub (v, r) = let | |
259 | val x = intcoord (v, rmin', rsize') | |
260 | val k = subindex (x, IMAXrs1) | |
261 | val cell = S.mkCell () | |
262 | in | |
263 | S.putCell (cell, k, r); cell | |
264 | end | |
265 | val box = (case root | |
266 | of S.Empty => S.Space{rmin=rmin', rsize=rsize', root=root} | |
267 | | _ => S.Space{ | |
268 | rmin = rmin', | |
269 | rsize = rsize', | |
270 | root = S.mkCellNode (mksub (rmid, root)) | |
271 | } | |
272 | (* end case *)) | |
273 | in | |
274 | expandBox (nd, box) | |
275 | end) | |
276 | ||
277 | ||
278 | (* insert a single node into the tree *) | |
279 | fun loadTree (body as S.Body{pos=posp, ...}, S.Space{rmin, rsize, root}) = let | |
280 | val xp = intcoord (!posp, rmin, rsize) | |
281 | fun insert (S.Empty, _) = S.mkBodyNode body | |
282 | | insert (n as S.Node{cell=S.BodyCell _, pos=posq, ...}, l) = let | |
283 | val xq = intcoord (!posq, rmin, rsize) | |
284 | val k = subindex (xq, l) | |
285 | val a = S.mkCell() | |
286 | in | |
287 | S.putCell(a, k, n); | |
288 | insert (S.mkCellNode a, l) | |
289 | end | |
290 | | insert (n as S.Node{cell, ...}, l) = let | |
291 | val k = subindex (xp, l) | |
292 | val subtree = insert (S.getCell (cell, k), rshift(l, 1)) | |
293 | in | |
294 | S.putCell (cell, k, subtree); | |
295 | n | |
296 | end | |
297 | in | |
298 | S.Space{rmin = rmin, rsize = rsize, root = insert (root, IMAXrs1)} | |
299 | end | |
300 | ||
301 | (* descend tree finding center-of-mass coordinates. *) | |
302 | fun hackCofM S.Empty = () | |
303 | | hackCofM (S.Node{cell = S.BodyCell _, ...}) = () | |
304 | | hackCofM (S.Node{cell = S.Cell subcells, mass, pos}) = let | |
305 | fun sumMass (i, totMass, cofm) = if (i < S.nsub) | |
306 | then (case Array.sub(subcells, i) | |
307 | of S.Empty => sumMass (i+1, totMass, cofm) | |
308 | | (nd as S.Node{mass, pos, ...}) => let | |
309 | val _ = hackCofM nd | |
310 | val m = !mass | |
311 | in | |
312 | sumMass (i+1, totMass + m, V.addv(cofm, V.mulvs(!pos, m))) | |
313 | end | |
314 | (* end case *)) | |
315 | else ( | |
316 | mass := totMass; | |
317 | pos := V.divvs(cofm, totMass)) | |
318 | in | |
319 | sumMass (0, 0.0, V.zerov) | |
320 | end | |
321 | ||
322 | (* initialize tree structure for hack force calculation. *) | |
323 | fun makeTree (bodies, rmin, rsize) = let | |
324 | fun build ([], space) = space | |
325 | | build ((body as S.Body{mass, ...}) :: r, space) = | |
326 | if Real.==(mass, 0.0) then build (r, space) | |
327 | else let | |
328 | val box = expandBox (body, space) | |
329 | val box = loadTree(body, box) | |
330 | in build (r, box) | |
331 | end | |
332 | val (space as S.Space{root, ...}) = | |
333 | build (bodies, S.Space{rmin=rmin, rsize=rsize, root=S.Empty}) | |
334 | in | |
335 | hackCofM root; | |
336 | space | |
337 | end | |
338 | ||
339 | end; (* functor Load *) | |
340 | (* grav.sml | |
341 | * | |
342 | * COPYRIGHT (c) 1993, AT&T Bell Laboratories. | |
343 | * | |
344 | * Gravity module for hierarchical N-body code; routines to compute gravity. | |
345 | *) | |
346 | ||
347 | signature GRAV = | |
348 | sig | |
349 | ||
350 | structure S : SPACE | |
351 | structure V : VECTOR | |
352 | sharing S.V = V | |
353 | ||
354 | val hackGrav : {body:S.body, root:S.node, rsize:real, tol:real, eps : real} | |
355 | -> {n2bterm:int, nbcterm:int, skipSelf:bool} | |
356 | ||
357 | end; (* GRAV *) | |
358 | ||
359 | functor Grav (S : SPACE) : GRAV = | |
360 | struct | |
361 | ||
362 | structure S = S | |
363 | structure V = S.V | |
364 | ||
365 | fun walk {acc0, phi0, pos0, pskip, eps, rsize, tol, root} = let | |
366 | val skipSelf = ref false | |
367 | val nbcterm = ref 0 and n2bterm = ref 0 | |
368 | val tolsq = (tol * tol) | |
369 | (* compute a single body-body or body-cell interaction. *) | |
370 | fun gravsub (S.Empty, phi0, acc0, _) = (phi0, acc0) | |
371 | | gravsub (p as S.Node{mass, pos, cell, ...}, phi0, acc0, memo) = let | |
372 | val (dr, drsq) = (case memo | |
373 | of NONE => let | |
374 | val dr = V.subv(!pos, pos0) | |
375 | in | |
376 | (dr, V.dotvp(dr, dr) + (eps*eps)) | |
377 | end | |
378 | | SOME(dr', drsq') => (dr', drsq' + (eps*eps)) | |
379 | (* end case *)) | |
380 | val phii = !mass / (Math.sqrt drsq) | |
381 | in | |
382 | case cell | |
383 | of (S.Cell _) => nbcterm := !nbcterm + 1 | |
384 | | _ => n2bterm := !n2bterm + 1 | |
385 | (* end case *); | |
386 | (phi0 - phii, V.addv(acc0, V.mulvs(dr, phii / drsq))) | |
387 | end (* gravsub *) | |
388 | (* recursive routine to do hackwalk operation. This combines the | |
389 | * subdivp and walksub routines from the C version. | |
390 | *) | |
391 | fun walksub (p, dsq, phi0, acc0) = ( | |
392 | (*print(implode[" walksub: dsq = ", makestring dsq, ", ", S.prNode p, "\n"]);*) | |
393 | case p | |
394 | of S.Empty => (phi0, acc0) | |
395 | | (S.Node{cell = S.BodyCell body, ...}) => | |
396 | if S.eqBody(body, pskip) | |
397 | then (skipSelf := true; (phi0, acc0)) | |
398 | else gravsub (p, phi0, acc0, NONE) | |
399 | | (S.Node{cell = S.Cell a, pos, ...}) => let | |
400 | val dr = V.subv(!pos, pos0) | |
401 | val drsq = V.dotvp(dr, dr) | |
402 | in | |
403 | if ((tolsq * drsq) < dsq) | |
404 | then let (* open p up *) | |
405 | fun loop (i, phi0, acc0) = if (i < S.nsub) | |
406 | then let | |
407 | val (phi0', acc0') = walksub ( | |
408 | Array.sub(a, i), dsq/4.0, phi0, acc0) | |
409 | in | |
410 | loop (i+1, phi0', acc0') | |
411 | end | |
412 | else (phi0, acc0) | |
413 | in | |
414 | loop (0, phi0, acc0) | |
415 | end | |
416 | else gravsub (p, phi0, acc0, SOME(dr, drsq)) | |
417 | end | |
418 | (* end case *)) | |
419 | val (phi0, acc0) = walksub (root, rsize*rsize, phi0, acc0) | |
420 | in | |
421 | { phi0 = phi0, acc0 = acc0, | |
422 | nbcterm = !nbcterm, n2bterm = !n2bterm, skip = !skipSelf | |
423 | } | |
424 | end (* walk *) | |
425 | ||
426 | (* evaluate grav field at a given particle. *) | |
427 | fun hackGrav {body as S.Body{pos, phi, acc, ...}, root, rsize, eps, tol} = let | |
428 | val {phi0, acc0, nbcterm, n2bterm, skip} = walk { | |
429 | acc0 = V.zerov, phi0 = 0.0, pos0 = !pos, pskip = body, | |
430 | eps = eps, rsize = rsize, tol = tol, root = root | |
431 | } | |
432 | in | |
433 | phi := phi0; | |
434 | acc := acc0; | |
435 | (** | |
436 | app (fn (fmt, items) => print(Format.format fmt items)) [ | |
437 | ("pos = [%f %f %f]\n", map Format.REAL (V.explode(!pos))), | |
438 | ("pos = [%f %f %f]\n", map Format.REAL (V.explode acc0)), | |
439 | ("phi = %f\n", [Format.REAL phi0]) | |
440 | ]; | |
441 | raise Fail ""; | |
442 | **) | |
443 | {nbcterm=nbcterm, n2bterm=n2bterm, skipSelf=skip} | |
444 | end (* hackgrav *) | |
445 | ||
446 | end; (* Grav *) | |
447 | (* data-io.sml | |
448 | * | |
449 | * COPYRIGHT (c) 1993, AT&T Bell Laboratories. | |
450 | * | |
451 | * I/O routines for export version of hierarchical N-body code. | |
452 | *) | |
453 | ||
454 | signature DATA_IO = | |
455 | sig | |
456 | ||
457 | structure S : SPACE | |
458 | ||
459 | val inputData : string -> { | |
460 | nbody : int, | |
461 | bodies : S.body list, | |
462 | tnow : real, | |
463 | headline : string | |
464 | } | |
465 | ||
466 | (* output routines *) | |
467 | val initOutput : { | |
468 | outfile : string, headline : string, nbody : int, tnow : real, | |
469 | dtime : real, eps : real, tol : real, dtout : real, tstop : real | |
470 | } -> unit | |
471 | val output : { | |
472 | nbody : int, bodies : S.body list, n2bcalc : int, nbccalc : int, | |
473 | selfint : int, tnow : real | |
474 | } -> unit | |
475 | val stopOutput : unit -> unit | |
476 | ||
477 | end; | |
478 | ||
479 | functor DataIO (S : SPACE) : DATA_IO = | |
480 | struct | |
481 | ||
482 | structure SS = Substring | |
483 | structure S = S | |
484 | structure V = S.V | |
485 | ||
486 | val atoi = valOf o Int.scan StringCvt.DEC SS.getc | |
487 | ||
488 | (* NOTE: this really out to be implemented using the lazy IO streams, | |
489 | * but SML/NJ doesn't implement these correctly yet. | |
490 | *) | |
491 | fun inputData fname = let | |
492 | val strm = TextIO.openIn fname | |
493 | val buf = ref(SS.full "") | |
494 | fun getLn () = (case (TextIO.inputLine strm) | |
495 | of NONE => raise Fail "inputData: EOF" | |
496 | | SOME s => buf := SS.full s | |
497 | (* end case *)) | |
498 | fun skipWS () = let | |
499 | val buf' = SS.dropl Char.isSpace (!buf) | |
500 | in | |
501 | if (SS.isEmpty buf') | |
502 | then (getLn(); skipWS()) | |
503 | else buf' | |
504 | end | |
505 | fun readInt () = let | |
506 | val (n, ss) = atoi (skipWS ()) | |
507 | in | |
508 | buf := ss; n | |
509 | end | |
510 | fun readReal () = let | |
511 | val (r, ss) = valOf (Real.scan SS.getc (skipWS())) | |
512 | in | |
513 | buf := ss; r | |
514 | end | |
515 | val nbody = readInt() | |
516 | val _ = if (nbody < 1) | |
517 | then raise Fail "absurd nbody" | |
518 | else () | |
519 | val ndim = readInt() | |
520 | val _ = if (ndim <> V.dim) | |
521 | then raise Fail "absurd ndim" | |
522 | else () | |
523 | val tnow = readReal() | |
524 | fun iter f = let | |
525 | fun loop (0, l) = l | |
526 | | loop (n, l) = loop (n-1, f() :: l) | |
527 | in | |
528 | fn n => loop (n, []) | |
529 | end | |
530 | fun readVec () = V.implode (rev (iter readReal ndim)) | |
531 | val massList = iter readReal nbody | |
532 | val posList = iter readVec nbody | |
533 | val velList = iter readVec nbody | |
534 | fun mkBodies ([], [], [], l) = l | |
535 | | mkBodies (m::r1, p::r2, v::r3, l) = let | |
536 | val b = S.Body{ | |
537 | mass = m, | |
538 | pos = ref p, | |
539 | vel = ref v, | |
540 | acc = ref V.zerov, | |
541 | phi = ref 0.0 | |
542 | } | |
543 | in | |
544 | mkBodies(r1, r2, r3, b::l) | |
545 | end | |
546 | in | |
547 | TextIO.closeIn strm; | |
548 | { nbody = nbody, | |
549 | bodies = mkBodies (massList, posList, velList, []), | |
550 | tnow = tnow, | |
551 | headline = concat["Hack code: input file ", fname, "\n"] | |
552 | } | |
553 | end | |
554 | ||
555 | local | |
556 | val timer = ref (Timer.startCPUTimer ()) | |
557 | in | |
558 | fun initTimer () = timer := Timer.startCPUTimer() | |
559 | fun cputime () = let | |
560 | val {usr, sys, ...} = Timer.checkCPUTimer(!timer) | |
561 | val totTim = usr | |
562 | in | |
563 | (Time.toReal totTim) / 60.0 | |
564 | end | |
565 | end | |
566 | ||
567 | type out_state = { | |
568 | tout : real, | |
569 | dtout : real, | |
570 | dtime : real, | |
571 | strm : TextIO.outstream | |
572 | } | |
573 | val outState = ref (NONE : out_state option) | |
574 | ||
575 | fun fprintf (strm, items) = TextIO.output(strm, String.concat items) | |
576 | fun printf items = fprintf(TextIO.stdOut, items) | |
577 | fun pad n s = StringCvt.padLeft #" " n s | |
578 | fun fmtInt (wid, i) = pad wid (Int.toString i) | |
579 | fun fmtReal (wid, prec, r) = pad wid (Real.fmt (StringCvt.FIX(SOME prec)) r) | |
580 | fun fmtRealE (wid, prec, r) = pad wid (Real.fmt (StringCvt.SCI(SOME prec)) r) | |
581 | local | |
582 | fun itemFmt r = fmtReal (9, 4, r) | |
583 | val fmt = V.format{lp="", sep="", rp="", cvt=itemFmt} | |
584 | in | |
585 | fun printvec (init, vec) = printf [ | |
586 | "\t ", pad 9 init, fmt vec, "\n" | |
587 | ] | |
588 | end (* local *) | |
589 | ||
590 | fun stopOutput () = (case (! outState) | |
591 | of NONE => () | |
592 | | (SOME{strm, ...}) => (TextIO.closeOut strm; outState := NONE) | |
593 | (* end case *)) | |
594 | ||
595 | fun initOutput {outfile, headline, nbody, tnow, dtime, eps, tol, dtout, tstop} = ( | |
596 | initTimer(); | |
597 | printf ["\n\t\t", headline, "\n\n"]; | |
598 | printf (map (pad 12) ["nbody", "dtime", "eps", "tol", "dtout", "tstop"]); | |
599 | printf ["\n"]; | |
600 | printf [fmtInt(12, nbody), fmtReal(12, 5, dtime)]; | |
601 | printf [ | |
602 | fmtInt(12, nbody), fmtReal(12, 5, dtime), | |
603 | fmtReal(12, 4, eps), fmtReal(12, 2, tol), | |
604 | fmtReal(12, 3, dtout), fmtReal(12, 2, tstop), "\n\n" | |
605 | ]; | |
606 | case outfile | |
607 | of "" => stopOutput() | |
608 | | _ => outState := SOME{ | |
609 | dtime = dtime, | |
610 | tout = tnow, | |
611 | dtout = dtout, | |
612 | strm = TextIO.openOut outfile | |
613 | } | |
614 | (* end case *)) | |
615 | ||
616 | (* compute set of dynamical diagnostics. *) | |
617 | fun diagnostics bodies = let | |
618 | fun loop ([], arg) = { | |
619 | mtot = #totM arg, (* total mass *) | |
620 | totKE = #totKE arg, (* total kinetic energy *) | |
621 | totPE = #totPE arg, (* total potential energy *) | |
622 | cOfMPos = #cOfMPos arg, (* center of mass: position *) | |
623 | cOfMVel = #cOfMVel arg, (* center of mass: velocity *) | |
624 | amVec = #amVec arg (* angular momentum vector *) | |
625 | } | |
626 | | loop (S.Body{ | |
627 | mass, pos=ref pos, vel=ref vel, acc=ref acc, phi=ref phi | |
628 | } :: r, arg) = let | |
629 | val velsq = V.dotvp(vel, vel) | |
630 | val halfMass = 0.5 * mass | |
631 | val posXmass = V.mulvs(pos, mass) | |
632 | in | |
633 | loop ( r, { | |
634 | totM = (#totM arg) + mass, | |
635 | totKE = (#totKE arg) + halfMass * velsq, | |
636 | totPE = (#totPE arg) + halfMass * phi, | |
637 | keTen = V.addm(#keTen arg, V.outvp(V.mulvs(vel, halfMass), vel)), | |
638 | peTen = V.addm(#peTen arg, V.outvp(posXmass, acc)), | |
639 | cOfMPos = V.addv(#cOfMPos arg, posXmass), | |
640 | cOfMVel = V.addv(#cOfMVel arg, V.mulvs(vel, mass)), | |
641 | amVec = V.addv(#amVec arg, V.mulvs(V.crossvp(pos, vel), mass)) | |
642 | }) | |
643 | end | |
644 | in | |
645 | loop (bodies, { | |
646 | totM = 0.0, totKE = 0.0, totPE = 0.0, | |
647 | keTen = V.zerom, peTen = V.zerom, | |
648 | cOfMPos = V.zerov, cOfMVel = V.zerov, | |
649 | amVec = V.zerov | |
650 | }) | |
651 | end (* diagnostics *) | |
652 | ||
653 | fun outputData (strm, tnow, nbody, bodies) = let | |
654 | fun outInt i = fprintf(strm, [" ", Int.toString i, "\n"]) | |
655 | fun outReal r = fprintf(strm, [" ", fmtRealE(21, 14, r), "\n"]) | |
656 | fun prReal r = fprintf(strm, [" ", fmtRealE(21, 14, r)]) | |
657 | fun outVec v = let | |
658 | fun out [] = TextIO.output(strm, "\n") | |
659 | | out (x::r) = (prReal x; out r) | |
660 | in | |
661 | out(V.explode v) | |
662 | end | |
663 | in | |
664 | outInt nbody; | |
665 | outInt V.dim; | |
666 | outReal tnow; | |
667 | app (fn (S.Body{mass, ...}) => outReal mass) bodies; | |
668 | app (fn (S.Body{pos, ...}) => outVec(!pos)) bodies; | |
669 | app (fn (S.Body{vel, ...}) => outVec(!vel)) bodies; | |
670 | printf ["\n\tparticle data written\n"] | |
671 | end; | |
672 | ||
673 | fun output {nbody, bodies, n2bcalc, nbccalc, selfint, tnow} = let | |
674 | val nttot = n2bcalc + nbccalc | |
675 | val nbavg = floor(real n2bcalc / real nbody) | |
676 | val ncavg = floor(real nbccalc / real nbody) | |
677 | val data = diagnostics bodies | |
678 | in | |
679 | printf ["\n"]; | |
680 | printf (map (pad 9) [ | |
681 | "tnow", "T+U", "T/U", "nttot", "nbavg", "ncavg", "selfint", | |
682 | "cputime" | |
683 | ]); | |
684 | printf ["\n"]; | |
685 | printf [ | |
686 | fmtReal(9, 3, tnow), fmtReal(9, 4, #totKE data + #totPE data), | |
687 | fmtReal(9, 4, #totKE data / #totPE data), fmtInt(9, nttot), | |
688 | fmtInt(9, nbavg), fmtInt(9, ncavg), fmtInt(9, selfint), | |
689 | fmtReal(9, 2, cputime()), "\n\n" | |
690 | ]; | |
691 | printvec ("cm pos", #cOfMPos data); | |
692 | printvec ("cm vel", #cOfMVel data); | |
693 | printvec ("am pos", #amVec data); | |
694 | case !outState | |
695 | of NONE => () | |
696 | | (SOME{tout, dtout, dtime, strm}) => | |
697 | if ((tout - 0.01 * dtime) <= tnow) | |
698 | then ( | |
699 | outputData (strm, tnow, nbody, bodies); | |
700 | outState := SOME{ | |
701 | tout=tout+dtout, dtout=dtout, dtime=dtime, strm=strm | |
702 | }) | |
703 | else () | |
704 | (* end case *) | |
705 | end | |
706 | ||
707 | end; (* DataIO *) | |
708 | ||
709 | (* getparam.sml | |
710 | * | |
711 | * COPYRIGHT (c) 1993, AT&T Bell Laboratories. | |
712 | *) | |
713 | ||
714 | structure GetParam : sig | |
715 | ||
716 | exception EOF | |
717 | ||
718 | val initParam : (string list * string list) -> unit | |
719 | val getParam : string -> string | |
720 | val getIParam : string -> int | |
721 | val getRParam : string -> real | |
722 | val getBParam : string -> bool | |
723 | ||
724 | end = struct | |
725 | ||
726 | exception EOF | |
727 | ||
728 | val defaults = ref ([] : string list) | |
729 | ||
730 | (* ignore arg vector, remember defaults. *) | |
731 | fun initParam (argv, defl) = defaults := defl | |
732 | ||
733 | fun prompt items = ( | |
734 | TextIO.output(TextIO.stdOut, String.concat items); | |
735 | TextIO.flushOut TextIO.stdOut) | |
736 | ||
737 | structure SS = Substring | |
738 | ||
739 | (* export version prompts user for value. *) | |
740 | fun getParam name = let | |
741 | fun scanBind [] = NONE | |
742 | | scanBind (s::r) = let | |
743 | val (_, suffix) = SS.position name (SS.full s) | |
744 | in | |
745 | if (SS.isEmpty suffix) | |
746 | then scanBind r | |
747 | else SOME(SS.string(SS.triml (size name+1) suffix)) | |
748 | end | |
749 | fun get default = (case (TextIO.inputLine TextIO.stdIn) | |
750 | of NONE => raise EOF | |
751 | | SOME "\n" => default | |
752 | | SOME s => substring(s, 0, size s - 1) | |
753 | (* end case *)) | |
754 | in | |
755 | if (null (! defaults)) | |
756 | then raise Fail "getParam called before initParam" | |
757 | else (); | |
758 | case (scanBind (! defaults)) | |
759 | of (SOME s) => ( | |
760 | prompt ["enter ", name, " [", s, "]: "]; | |
761 | get s) | |
762 | | NONE => (prompt ["enter ", name, ": "]; get "") | |
763 | (* end case *) | |
764 | end | |
765 | ||
766 | local | |
767 | fun cvt scanFn = let | |
768 | fun cvt' name = let | |
769 | fun get () = (case getParam name of "" => get () | s => s) | |
770 | val param = get () | |
771 | in | |
772 | (valOf (scanFn param)) handle _ => (cvt' name) | |
773 | end | |
774 | in | |
775 | cvt' | |
776 | end | |
777 | in | |
778 | (* get integer parameter *) | |
779 | val getIParam = cvt Int.fromString | |
780 | (* get real parameter *) | |
781 | val getRParam = cvt Real.fromString | |
782 | (* get bool parameter *) | |
783 | val getBParam = cvt Bool.fromString | |
784 | end (* local *) | |
785 | ||
786 | end; (* GetParam *) | |
787 | ||
788 | (* rand-sig.sml | |
789 | * | |
790 | * COPYRIGHT (c) 1993 by AT&T Bell Laboratories. See COPYRIGHT file for details. | |
791 | * | |
792 | * Signature for a simple random number generator. | |
793 | * | |
794 | *) | |
795 | ||
796 | signature RAND = | |
797 | sig | |
798 | ||
799 | val randMin : real | |
800 | val randMax : real | |
801 | val random : real -> real | |
802 | (* Given seed, return value randMin <= v <= randMax | |
803 | * Iteratively using the value returned by random as the | |
804 | * next seed to random will produce a sequence of pseudo-random | |
805 | * numbers. | |
806 | *) | |
807 | ||
808 | val mkRandom : real -> unit -> real | |
809 | (* Given seed, return function generating a sequence of | |
810 | * random numbers randMin <= v <= randMax | |
811 | *) | |
812 | ||
813 | val norm : real -> real | |
814 | (* r -> r / (randMax + 1.0) *) | |
815 | ||
816 | val range : (int * int) -> real -> int | |
817 | (* Map v, randMin <= v <= randMax to integer range [i,j] | |
818 | * Exception - | |
819 | * BadArg if j < i | |
820 | *) | |
821 | ||
822 | end (* RAND *) | |
823 | (* rand.sml | |
824 | * | |
825 | * COPYRIGHT (c) 1991 by AT&T Bell Laboratories. See COPYRIGHT file for details | |
826 | * | |
827 | * Random number generator taken from Paulson, pp. 170-171. | |
828 | * Recommended by Stephen K. Park and Keith W. Miller, | |
829 | * Random number generators: good ones are hard to find, | |
830 | * CACM 31 (1988), 1192-1201 | |
831 | * | |
832 | * Note: The Random structure provides a better generator. | |
833 | *) | |
834 | ||
835 | structure Rand : RAND = | |
836 | struct | |
837 | ||
838 | (* real number version for systems with 46-bit mantissas *) | |
839 | val a = 16807.0 and m = 2147483647.0 | |
840 | ||
841 | val randMin = 1.0 | |
842 | val randMax = m - 1.0 | |
843 | ||
844 | fun random seed = let | |
845 | val t = a*seed | |
846 | in | |
847 | t - m * real(floor(t/m)) | |
848 | end | |
849 | ||
850 | fun mkRandom seed = let | |
851 | val seed = ref seed | |
852 | in | |
853 | fn () => (seed := random (!seed); !seed) | |
854 | end | |
855 | ||
856 | fun norm r = r / m | |
857 | ||
858 | fun range (i,j) = | |
859 | if j < i | |
860 | then raise Fail "Rand.range" | |
861 | else let | |
862 | val R = real(j - i + 1) | |
863 | in | |
864 | fn r => i + floor(R*(r/m)) | |
865 | end handle _ => let | |
866 | val ri = real i | |
867 | val R = (real j)-ri+1.0 | |
868 | in | |
869 | fn r => floor(ri + R*(r/m)) | |
870 | end | |
871 | ||
872 | end (* Rand *) | |
873 | ||
874 | (* main.sml | |
875 | * | |
876 | * COPYRIGHT (c) 1993, AT&T Bell Laboratories. | |
877 | * | |
878 | * This is the main driver for the Barnse-HutN-body code. | |
879 | *) | |
880 | ||
881 | functor Main (V : VECTOR) : sig | |
882 | ||
883 | structure S : SPACE | |
884 | structure V : VECTOR | |
885 | structure L : LOAD | |
886 | ||
887 | val srand : int -> unit | |
888 | (* reset the random number generator *) | |
889 | ||
890 | val testdata : int -> S.body list | |
891 | (* generate the Plummer model data *) | |
892 | ||
893 | val go : { | |
894 | output : {n2bcalc:int, nbccalc:int, nstep:int, selfint:int, tnow:real} | |
895 | -> unit, | |
896 | bodies : S.body list, tnow : real, tstop : real, | |
897 | dtime : real, eps : real, tol : real, | |
898 | rmin : real V.vec, rsize : real | |
899 | } -> unit | |
900 | ||
901 | val doit : unit -> unit | |
902 | ||
903 | end = struct | |
904 | ||
905 | structure V = V | |
906 | structure S = Space(V) | |
907 | structure L = Load(S) | |
908 | structure G = Grav(S) | |
909 | structure DataIO = DataIO(S) | |
910 | ||
911 | (* some math utilities *) | |
912 | (** NOTE: these are part of the Math structure in the new basis *) | |
913 | val pi = 3.14159265358979323846 | |
914 | fun pow(x, y) = | |
915 | if Real.==(y, 0.0) then 1.0 else Math.exp (y * Math.ln x) | |
916 | ||
917 | (* random numbers *) | |
918 | local | |
919 | val seed = ref 0.0 | |
920 | in | |
921 | fun srand s = (seed := real s) | |
922 | fun xrand (xl, xh) = let | |
923 | val r = Rand.random (! seed) | |
924 | in | |
925 | seed := r; | |
926 | xl + (((xh - xl) * r) / 2147483647.0) | |
927 | end | |
928 | end (* local *) | |
929 | ||
930 | (* default parameter values *) | |
931 | val defaults = [ | |
932 | (* file names for input/output *) | |
933 | "in=", (* snapshot of initial conditions *) | |
934 | "out=", (* stream of output snapshots *) | |
935 | ||
936 | (* params, used if no input specified, to make a Plummer Model*) | |
937 | "nbody=128", (* number of particles to generate *) | |
938 | "seed=123", (* random number generator seed *) | |
939 | ||
940 | (* params to control N-body integration *) | |
941 | "dtime=0.025", (* integration time-step *) | |
942 | "eps=0.05", (* usual potential softening *) | |
943 | "tol=1.0", (* cell subdivision tolerence *) | |
944 | "fcells=0.75", (* cell allocation parameter *) | |
945 | ||
946 | "tstop=2.0", (* time to stop integration *) | |
947 | "dtout=0.25", (* data-output interval *) | |
948 | ||
949 | "debug=false", (* turn on debugging messages *) | |
950 | "VERSION=1.0" (* JEB 06 March 1988 *) | |
951 | ] | |
952 | ||
953 | (* pick a random point on a sphere of specified radius. *) | |
954 | fun pickshell rad = let | |
955 | fun pickvec () = let | |
956 | val vec = V.tabulate (fn _ => xrand(~1.0, 1.0)) | |
957 | val rsq = V.dotvp(vec, vec) | |
958 | in | |
959 | if (rsq > 1.0) | |
960 | then pickvec () | |
961 | else V.mulvs (vec, rad / Math.sqrt(rsq)) | |
962 | end | |
963 | in | |
964 | pickvec () | |
965 | end | |
966 | ||
967 | (* generate Plummer model initial conditions for test runs, scaled | |
968 | * to units such that M = -4E = G = 1 (Henon, Hegge, etc). | |
969 | * See Aarseth, SJ, Henon, M, & Wielen, R (1974) Astr & Ap, 37, 183. | |
970 | *) | |
971 | fun testdata n = let | |
972 | val mfrac = 0.999 (* mass cut off at mfrac of total *) | |
973 | val rn = real n | |
974 | val rsc = (3.0 * pi) / 16.0 | |
975 | val vsc = Math.sqrt(1.0 / rsc) | |
976 | fun mkBodies (0, cmr, cmv, l) = let | |
977 | (* offset bodies by normalized cm coordinates. Also, reverse | |
978 | * the list to get the same order of bodies as in the C version. | |
979 | *) | |
980 | val cmr = V.divvs(cmr, rn) | |
981 | val cmv = V.divvs(cmv, rn) | |
982 | fun norm ([], l) = l | |
983 | | norm ((p as S.Body{pos, vel, ...}) :: r, l) = ( | |
984 | pos := V.subv(!pos, cmr); | |
985 | vel := V.subv(!vel, cmv); | |
986 | norm (r, p::l)) | |
987 | in | |
988 | norm (l, []) | |
989 | end | |
990 | | mkBodies (i, cmr, cmv, l) = let | |
991 | val r = 1.0 / Math.sqrt (pow(xrand(0.0, mfrac), ~2.0/3.0) - 1.0) | |
992 | val pos = pickshell (rsc * r) | |
993 | fun vN () = let (* von Neumann technique *) | |
994 | val x = xrand(0.0,1.0) | |
995 | val y = xrand(0.0,0.1) | |
996 | in | |
997 | if (y > x*x * (pow (1.0-x*x, 3.5))) then vN () else x | |
998 | end | |
999 | val v = ((Math.sqrt 2.0) * vN()) / pow(1.0 + r*r, 0.25) | |
1000 | val vel = pickshell (vsc * v) | |
1001 | val body = S.Body{ | |
1002 | mass = 1.0 / rn, | |
1003 | pos = ref pos, | |
1004 | vel = ref vel, | |
1005 | acc = ref V.zerov, | |
1006 | phi = ref 0.0 | |
1007 | } | |
1008 | in | |
1009 | mkBodies (i-1, V.addv(cmr, pos), V.addv(cmv, vel), body :: l) | |
1010 | end | |
1011 | in | |
1012 | mkBodies (n, V.zerov, V.zerov, []) | |
1013 | end (* testdata *) | |
1014 | ||
1015 | (* startup hierarchical N-body code. This either reads in or generates | |
1016 | * an initial set of bodies, and other parameters. | |
1017 | *) | |
1018 | ||
1019 | fun startrun argv = let | |
1020 | val _ = GetParam.initParam(argv, defaults) | |
1021 | val {nbody, bodies, tnow, headline} = (case (GetParam.getParam "in") | |
1022 | of "" => let | |
1023 | val nbody = GetParam.getIParam "nbody" | |
1024 | in | |
1025 | if (nbody < 1) | |
1026 | then raise Fail "startrun: absurd nbody" | |
1027 | else (); | |
1028 | srand (GetParam.getIParam "seed"); | |
1029 | { nbody = nbody, | |
1030 | bodies = testdata nbody, | |
1031 | tnow = 0.0, | |
1032 | headline = "Hack code: Plummer model" | |
1033 | } | |
1034 | end | |
1035 | | fname => DataIO.inputData fname | |
1036 | (* end case *)) | |
1037 | in | |
1038 | { nbody = nbody, | |
1039 | bodies = bodies, | |
1040 | headline = headline, | |
1041 | outfile = GetParam.getParam "out", | |
1042 | dtime = GetParam.getRParam "dtime", | |
1043 | eps = GetParam.getRParam "eps", | |
1044 | tol = GetParam.getRParam "tol", | |
1045 | tnow = tnow, | |
1046 | tstop = GetParam.getRParam "tstop", | |
1047 | dtout = GetParam.getRParam "dtout", | |
1048 | debug = GetParam.getBParam "debug", | |
1049 | rmin = V.tabulate (fn _ => ~2.0), | |
1050 | rsize = 4.0 | |
1051 | } | |
1052 | end | |
1053 | ||
1054 | (* advance N-body system one time-step. *) | |
1055 | fun stepSystem output {plist, dtime, eps, nstep, rmin, rsize, tnow, tol} = let | |
1056 | val dthf = 0.5 * dtime | |
1057 | val S.Space{rmin, rsize, root} = L.makeTree (plist, rmin, rsize) | |
1058 | (* recalculate accelaration *) | |
1059 | fun recalc ([], n2bcalc, nbccalc, selfint) = (n2bcalc, nbccalc, selfint) | |
1060 | | recalc (p::r, n2bcalc, nbccalc, selfint) = let | |
1061 | val S.Body{acc as ref acc1, vel, ...} = p | |
1062 | val {n2bterm, nbcterm, skipSelf} = G.hackGrav { | |
1063 | body = p, root = root, rsize = rsize, eps = eps, tol = tol | |
1064 | } | |
1065 | in | |
1066 | if (nstep > 0) | |
1067 | then (* use change in accel to make 2nd order *) | |
1068 | (* correction to vel. *) | |
1069 | vel := V.addv(!vel, V.mulvs(V.subv(!acc, acc1), dthf)) | |
1070 | else (); | |
1071 | recalc (r, n2bcalc+n2bterm, nbccalc+nbcterm, | |
1072 | if skipSelf then selfint else (selfint+1)) | |
1073 | end | |
1074 | (* advance bodies *) | |
1075 | fun advance (S.Body{pos, acc, vel, ...}) = let | |
1076 | val dvel = V.mulvs (!acc, dthf) | |
1077 | val vel1 = V.addv (!vel, dvel) | |
1078 | val dpos = V.mulvs (vel1, dtime) | |
1079 | in | |
1080 | pos := V.addv (!pos, dpos); | |
1081 | vel := V.addv (vel1, dvel) | |
1082 | end | |
1083 | val (n2bcalc, nbccalc, selfint) = recalc (plist, 0, 0, 0) | |
1084 | in | |
1085 | output {nstep=nstep, tnow=tnow, n2bcalc=n2bcalc, nbccalc=nbccalc, selfint=selfint}; | |
1086 | app advance plist; | |
1087 | (nstep+1, tnow + dtime) | |
1088 | end | |
1089 | ||
1090 | (* given an initial configuration, run the simulation *) | |
1091 | fun go { | |
1092 | output, bodies, tnow, tstop, | |
1093 | dtime, eps, tol, rsize, rmin | |
1094 | } = let | |
1095 | val step = stepSystem output | |
1096 | fun loop (nstep, tnow) = if (tnow < tstop + (0.1 * dtime)) | |
1097 | then loop (step { | |
1098 | plist = bodies, dtime = dtime, eps = eps, nstep = nstep, | |
1099 | rmin = rmin, rsize = rsize, tnow = tnow, tol = tol | |
1100 | }) | |
1101 | else () | |
1102 | in | |
1103 | loop (0, tnow) | |
1104 | end | |
1105 | ||
1106 | fun doit () = let | |
1107 | val { nbody, bodies, headline, outfile, | |
1108 | dtime, eps, tol, tnow, tstop, dtout, | |
1109 | debug, rsize, rmin | |
1110 | } = startrun [] | |
1111 | fun output {nstep, tnow, n2bcalc, nbccalc, selfint} = DataIO.output{ | |
1112 | bodies = bodies, nbody = nbody, | |
1113 | n2bcalc = n2bcalc, nbccalc = nbccalc, | |
1114 | selfint = selfint, tnow = tnow | |
1115 | } | |
1116 | in | |
1117 | DataIO.initOutput { | |
1118 | outfile = outfile, headline = headline, nbody = nbody, tnow = tnow, | |
1119 | dtime = dtime, eps = eps, tol = tol, dtout = dtout, tstop = tstop | |
1120 | }; | |
1121 | go { | |
1122 | output=output, bodies=bodies, tnow=tnow, tstop=tstop, | |
1123 | dtime=dtime, eps=eps, tol=tol, rsize=rsize, rmin=rmin | |
1124 | }; | |
1125 | DataIO.stopOutput() | |
1126 | end (* doit *) | |
1127 | ||
1128 | end; (* Main *) | |
1129 | (* vector3.sml | |
1130 | * | |
1131 | * COPYRIGHT (c) 1993, AT&T Bell Laboratories. | |
1132 | * | |
1133 | * 3 dimensional vector arithmetic. | |
1134 | *) | |
1135 | ||
1136 | structure Vector3 : VECTOR = | |
1137 | struct | |
1138 | ||
1139 | type 'a vec = {x : 'a, y : 'a, z : 'a} | |
1140 | type realvec = real vec | |
1141 | ||
1142 | val dim = 3 | |
1143 | ||
1144 | fun tabulate f = {x = f 0, y = f 1, z = f 2} | |
1145 | ||
1146 | val zerov = {x = 0.0, y = 0.0, z = 0.0} | |
1147 | fun equal({x, y, z}, {x=x1, y=y1, z=z1}) = | |
1148 | Real.==(x, x1) andalso Real.==(y, y1) andalso Real.==(z, z1) | |
1149 | ||
1150 | fun addv ({x=x1, y=y1, z=z1} : realvec, {x=x2, y=y2, z=z2}) = | |
1151 | {x=x1+x2, y=y1+y2, z=z1+z2} | |
1152 | ||
1153 | fun subv ({x=x1, y=y1, z=z1} : realvec, {x=x2, y=y2, z=z2}) = | |
1154 | {x=x1-x2, y=y1-y2, z=z1-z2} | |
1155 | ||
1156 | fun dotvp ({x=x1, y=y1, z=z1} : realvec, {x=x2, y=y2, z=z2}) = | |
1157 | x1*x2 + y1*y2 + z1*z2 | |
1158 | ||
1159 | fun crossvp ({x=x1, y=y1, z=z1} : realvec, {x=x2, y=y2, z=z2}) = | |
1160 | {x = y1*z2 - z1*y2, y = x1*z2 - z1*x2, z = x1*y2 - y1*x2} | |
1161 | ||
1162 | fun addvs ({x, y, z} : realvec, s) = {x=x+s, y=y+s, z=z+s} | |
1163 | ||
1164 | fun mulvs ({x, y, z} : realvec, s) = {x=x*s, y=y*s, z=z*s} | |
1165 | ||
1166 | fun divvs ({x, y, z} : realvec, s) = {x=x/s, y=y/s, z=z/s} | |
1167 | ||
1168 | fun mapv f {x, y, z} = {x = f x, y = f y, z = f z} | |
1169 | ||
1170 | fun map3v f ({x=x1, y=y1, z=z1}, {x=x2, y=y2, z=z2}, {x=x3, y=y3, z=z3}) = | |
1171 | {x = f(x1, x2, x3), y = f(y1, y2, y3), z = f(z1, z2, z3)} | |
1172 | ||
1173 | fun foldv f {x, y, z} init = f(z, f(y, f(x, init))) | |
1174 | ||
1175 | fun format {lp, rp, sep, cvt} {x, y, z} = String.concat[ | |
1176 | lp, cvt x, sep, cvt y, sep, cvt z, rp | |
1177 | ] | |
1178 | ||
1179 | fun explode {x, y, z} = [x, y, z] | |
1180 | ||
1181 | fun implode [x, y, z] = {x=x, y=y, z=z} | |
1182 | | implode _ = raise Fail "implode: bad dimension" | |
1183 | ||
1184 | type matrix = { | |
1185 | m00 : real, m01 : real, m02 : real, | |
1186 | m10 : real, m11 : real, m12 : real, | |
1187 | m20 : real, m21 : real, m22 : real | |
1188 | } | |
1189 | ||
1190 | val zerom = { | |
1191 | m00 = 0.0, m01 = 0.0, m02 = 0.0, | |
1192 | m10 = 0.0, m11 = 0.0, m12 = 0.0, | |
1193 | m20 = 0.0, m21 = 0.0, m22 = 0.0 | |
1194 | } | |
1195 | ||
1196 | fun addm (a : matrix, b : matrix) = { | |
1197 | m00=(#m00 a + #m00 b), m01=(#m01 a + #m01 b), m02=(#m02 a + #m02 b), | |
1198 | m10=(#m10 a + #m10 b), m11=(#m11 a + #m11 b), m12=(#m12 a + #m12 b), | |
1199 | m20=(#m20 a + #m20 b), m21=(#m21 a + #m21 b), m22=(#m22 a + #m22 b) | |
1200 | } | |
1201 | ||
1202 | fun outvp ({x=a0, y=a1, z=a2} : realvec, {x=b0, y=b1, z=b2}) = { | |
1203 | m00=(a0*b0), m01=(a0*b1), m02=(a0*b2), | |
1204 | m10=(a1*b0), m11=(a1*b1), m12=(a1*b2), | |
1205 | m20=(a2*b0), m21=(a2*b1), m22=(a2*b2) | |
1206 | } | |
1207 | ||
1208 | end (* VectMath *) | |
1209 | ||
1210 | signature BMARK = | |
1211 | sig | |
1212 | val doit : int -> unit | |
1213 | val testit : TextIO.outstream -> unit | |
1214 | end; | |
1215 | (* load file for bmark version *) | |
1216 | ||
1217 | (* | |
1218 | app use [ | |
1219 | "rand-sig.sml", | |
1220 | "rand.sml", | |
1221 | "vector-sig.sml", | |
1222 | "space.sml", | |
1223 | "load.sml", | |
1224 | "grav.sml", | |
1225 | "getparam.sml", | |
1226 | "data-io.sml", | |
1227 | "main.sml", | |
1228 | "vector3.sml" | |
1229 | ]; | |
1230 | *) | |
1231 | structure Main : BMARK = | |
1232 | struct | |
1233 | structure M3 = Main(Vector3); | |
1234 | ||
1235 | val name = "Barnes-Hut (3d)" | |
1236 | ||
1237 | fun testit strm = () | |
1238 | ||
1239 | fun doit n = ( | |
1240 | M3.srand 123; | |
1241 | M3.go { | |
1242 | output = fn _ => (), | |
1243 | bodies = M3.testdata n, | |
1244 | tnow = 0.0, tstop = 2.0, | |
1245 | dtime = 0.025, eps = 0.05, tol = 1.0, | |
1246 | rmin = M3.S.V.tabulate (fn _ => ~2.0), | |
1247 | rsize = 4.0 | |
1248 | }) | |
1249 | end; | |
1250 |