plummer generation files
authorPritish Jetley <pjetley2@illinois.edu>
Wed, 12 Oct 2011 22:12:40 +0000 (17:12 -0500)
committerPritish Jetley <pjetley2@illinois.edu>
Wed, 12 Oct 2011 22:12:40 +0000 (17:12 -0500)
gen_util.cpp [new file with mode: 0644]
plummer.cpp [new file with mode: 0644]

diff --git a/gen_util.cpp b/gen_util.cpp
new file mode 100644 (file)
index 0000000..2a88034
--- /dev/null
@@ -0,0 +1,65 @@
+/*************************************************************************/
+/*                                                                       */
+/*  Copyright (c) 1994 Stanford University                               */
+/*                                                                       */
+/*  All rights reserved.                                                 */
+/*                                                                       */
+/*  Permission is given to use, copy, and modify this software for any   */
+/*  non-commercial purpose as long as this copyright notice is not       */
+/*  removed.  All other uses, including redistribution in whole or in    */
+/*  part, are forbidden without prior written permission.                */
+/*                                                                       */
+/*  This software is provided with absolutely no warranty and no         */
+/*  support.                                                             */
+/*                                                                       */
+/*************************************************************************/
+
+#include <stdio.h>
+#include "stdinc.h"
+
+#include "define.h"
+
+#define HZ 60.0
+#define MULT 1103515245
+#define ADD 12345
+#define MASK (0x7FFFFFFF)
+#define TWOTO31 2147483648.0
+
+int A = 1;
+int B = 0;
+int randx = 1;
+int lastrand;   /* the last random number */
+
+/*
+ * XRAND: generate floating-point random number.
+ */
+
+Real prand();
+
+Real xrand(Real xl, Real xh)
+{
+   Real x;
+
+   return (xl + (xh - xl) * prand());
+}
+
+void pranset(int seed)
+{
+   int proc;
+  
+   A = 1;
+   B = 0;
+   randx = (A*seed+B) & MASK;
+   A = (MULT * A) & MASK;
+   B = (MULT*B + ADD) & MASK;
+}
+
+Real prand()
+/*
+       Return a random Real in [0, 1.0)
+*/
+{
+   lastrand = randx;
+   randx = (A*randx+B) & MASK;
+   return((Real)lastrand/TWOTO31);
+}
diff --git a/plummer.cpp b/plummer.cpp
new file mode 100644 (file)
index 0000000..4bb594b
--- /dev/null
@@ -0,0 +1,124 @@
+/* adapted from barnes benchmark in SPLASH2 */
+
+/*
+ * TESTDATA: generate Plummer model initial conditions for test runs,
+ * scaled to units such that M = -4E = G = 1 (Henon, Hegge, etc).
+ * See Aarseth, SJ, Henon, M, & Wielen, R (1974) Astr & Ap, 37, 183.
+ */
+
+#include "Vector3D.h"
+#include "Particle.h"
+
+#define MFRAC  0.999                /* mass cut off at MFRAC of total */
+
+Real xrand();
+Real pow();
+
+void testdata(int nbody)
+{
+   real rsc, vsc, rsq, r, v, x, y;
+   Vector3D<Real> cmr, cmv;
+   Particle *p;
+   int rejects = 0;
+   int k;
+   int halfnbody, i;
+   Real offset;
+   Particle *cp;
+   Real tmp;
+
+   Particle *bodytab = (bodyptr) G_MALLOC(nbody * sizeof(body));
+   if (bodytab == NULL) {
+      error("testdata: not enuf memory\n");
+   }
+   rsc = 9 * PI / 16;
+   vsc = sqrt(1.0 / rsc);
+
+   CLRV(cmr);
+   CLRV(cmv);
+
+   halfnbody = nbody / 2;
+   if (nbody % 2 != 0) halfnbody++;
+   for (p = bodytab; p < bodytab+halfnbody; p++) {
+      Type(p) = BODY;
+      Mass(p) = 1.0 / nbody;
+      Cost(p) = 1;
+
+      r = 1 / sqrt(pow(xrand(0.0, MFRAC), -2.0/3.0) - 1);
+      /*   reject radii greater than 10 */
+      while (r > 9.0) {
+        rejects++;
+        r = 1 / sqrt(pow(xrand(0.0, MFRAC), -2.0/3.0) - 1);
+      }        
+      pickshell(Pos(p), rsc * r);
+      ADDV(cmr, cmr, Pos(p));
+      do {
+        x = xrand(0.0, 1.0);
+        y = xrand(0.0, 0.1);
+
+      } while (y > x*x * pow(1 - x*x, 3.5));
+
+      v = sqrt(2.0) * x / pow(1 + r*r, 0.25);
+      pickshell(Vel(p), vsc * v);
+      ADDV(cmv, cmv, Vel(p));
+   }
+
+   offset = 4.0;
+
+   for (p = bodytab + halfnbody; p < bodytab+nbody; p++) {
+      Type(p) = BODY;
+      Mass(p) = 1.0 / nbody;
+      Cost(p) = 1;
+
+      cp = p - halfnbody;
+      for (i = 0; i < NDIM; i++){
+        Pos(p)[i] = Pos(cp)[i] + offset; 
+        ADDV(cmr, cmr, Pos(p));
+        Vel(p)[i] = Vel(cp)[i];
+        ADDV(cmv, cmv, Vel(p));
+      }
+   }
+
+   DIVVS(cmr, cmr, (real) nbody);
+   DIVVS(cmv, cmv, (real) nbody);
+
+   for (p = bodytab; p < bodytab+nbody; p++) {
+      SUBV(Pos(p), Pos(p), cmr);
+      SUBV(Vel(p), Vel(p), cmv);
+   }
+}
+
+/*
+ * PICKSHELL: pick a random point on a sphere of specified radius.
+ */
+
+pickshell(vec, rad)
+   real vec[];                     /* coordinate vector chosen */
+   real rad;                       /* radius of chosen point */
+{
+   register int k;
+   double rsq, xrand(), sqrt(), rsc;
+
+   do {
+      for (k = 0; k < NDIM; k++) {
+        vec[k] = xrand(-1.0, 1.0);
+      }
+      DOTVP(rsq, vec, vec);
+   } while (rsq > 1.0);
+
+   rsc = rad / sqrt(rsq);
+   MULVS(vec, vec, rsc);
+}
+
+
+int intpow(i,j)
+  int i,j;
+{   
+    int k;
+    int temp = 1;
+
+    for (k = 0; k < j; k++)
+        temp = temp*i;
+    return temp;
+}
+
+