00001
00018
00019 //............................................................................
00020 // version 1: July 1989 Piet Hut email: piet@iassns.bitnet
00021 // Institute for Advanced Study, Princeton, NJ, USA
00022 // version 2: Dec 1992 Piet Hut -- adopted to the new C++-based starlab
00023 //............................................................................
00024 // litt: S.J. Aarseth, M. Henon and R. Wielen (1974),
00025 // Astron. and Astrophys. 37, p. 183.
00026 //............................................................................
00027
00028 #include "dyn.h"
00029
00030 #ifdef TOOLBOX
00031
00032 #define MFRAC_DEFAULT 0.999 // radial cutoff
00033 #define RFRAC_DEFAULT 22.8042468 // corresponding radial cutoff
00034 // (not used right now)
00035 #define SEED_STRING_LENGTH 60
00036
00037 local void mkplummer(dyn * b, int n, real mfrac, real rfrac, bool u_flag)
00038 {
00039 int i;
00040
00041 real partmass; // equal mass of each particle
00042 real radius; // absolute value of position vector
00043 real velocity; // absolute value of velocity vector
00044 real theta, phi; // direction angles of above vectors
00045 real x, y; // for use in rejection technique
00046 real scalefactor; // for converting between different units
00047 real inv_scalefactor; // inverse scale factor
00048 real sqrt_scalefactor; // sqare root of scale factor
00049 real mrfrac; // m( rfrac )
00050 real m_min, m_max; // mass shell limits for quiet start
00051 dyn * bi;
00052
00053 b->set_mass(1.0);
00054 partmass = 1.0 / ((real) n);
00055
00056 // Calculating the coordinates is easiest in STRUCTURAL units;
00057 // conversion to VIRIAL units will be performed below.
00058 //
00059 // Recipe for scaling to the proper system of units:
00060 //
00061 // Since G = M = 1, if we switch from a coordinate system with
00062 // length unit r_old to a coordinate system with length unit r_new ,
00063 // the length units simply scale by a factor C = r_new / r_old .
00064 // Consequently, the coordinate values of physical quantities
00065 // such as positions should transform inversely to maintain the same
00066 // coordinate-invariant meaning. Similarly, the square of the velocities
00067 // should transform inversely proportional to the positions,
00068 // since GM = 1 (cf. a relation such as v*v = G*M/r ).
00069 //
00070 // To summarize: If
00071 // r_unit(new) = C * r_unit(old) ,
00072 // then
00073 // pos(new) = (1/C) * pos(old)
00074 // and
00075 // vel(new) = sqrt(C) * vel(old) .
00076
00077 scalefactor = 16.0 / (3.0 * PI);
00078 inv_scalefactor = 1.0 / scalefactor;
00079 sqrt_scalefactor = sqrt( scalefactor );
00080
00081 // Now convert rfrac into an equivalent mfrac, if necessary:
00082
00083 if (rfrac < VERY_LARGE_NUMBER) { // Note: the following powers can cause
00084 // problems on some compilers...
00085
00086 rfrac *= scalefactor; // Convert from VIRIAL to STRUCTURAL units
00087 mrfrac = rfrac*rfrac*rfrac / pow(1.0 + rfrac*rfrac, 1.5);
00088 if (mrfrac < mfrac)
00089 mfrac = mrfrac; // mfrac = min(mfrac, m(rfrac))
00090 }
00091
00092
00093 // Now construct the individual particles. Note that the "raw"
00094 // Plummer model has a virial radius of approximately 1.695 units
00095 // and is close to virial equilibrium. Include this factor so
00096 // that the unscaled model (-u) is already almost in standard
00097 // units.
00098
00099 real xfac = 1/1.695;
00100 real vfac = 1/sqrt(xfac);
00101
00102 for (i = 0, bi = b->get_oldest_daughter(); i < n;
00103 i++, bi = bi->get_younger_sister()) {
00104
00105 bi->set_mass(partmass);
00106
00107 // The position coordinates are determined by inverting the cumulative
00108 // mass-radius relation, with the cumulative mass drawn randomly from
00109 // [0, mfrac]; cf. Aarseth et al. (1974), eq. (A2).
00110
00111 m_min = (i * mfrac)/n;
00112 m_max = ((i+1) * mfrac)/n;
00113 real rrrr = randinter(m_min, m_max);
00114 radius = xfac / sqrt( pow (rrrr, -2.0/3.0) - 1.0);
00115
00116 // Note that this procedure arranges the particles approximately
00117 // in order of increasing distance fro the cluster center, which may
00118 // sometimes be desirable, and sometimes not. Extra option added
00119 // by Steve (7/99) to allow radii to be randomized (reshuffled).
00120
00121 theta = acos(randinter(-1.0, 1.0));
00122 phi = randinter(0.0, TWO_PI);
00123
00124 bi->set_pos(vector(radius * sin( theta ) * cos( phi ),
00125 radius * sin( theta ) * sin( phi ),
00126 radius * cos( theta )));
00127
00128 // The velocity coordinates are determined using von Neumann's
00129 // rejection technique, cf. Aarseth et al. (1974), eq. (A4,5).
00130 // First we take initial values for x, the ratio of velocity and
00131 // escape velocity (q in Aarseth et al.), and y, as a trick to
00132 // enter the body of the while loop.
00133
00134 x = 0.0;
00135 y = 0.1;
00136
00137 // Then we keep spinning the random number generator until we find a
00138 // pair of values (x,y), so that y < g(x) = x*x*pow( 1.0 - x*x, 3.5).
00139 // Whenever an y-value lies above the g(x) curve, the (x,y) pair is
00140 // discarded, and a new pair is selected. The value 0.1 is chosen as
00141 // a good upper limit for g(x) in [0,1] : 0.1 > max g(x) = 0.092 for
00142 // 0 < x < 1.
00143
00144 while (y > x*x*pow( 1.0 - x*x, 3.5)) {
00145 x = randinter(0.0,1.0);
00146 y = randinter(0.0,0.1);
00147 }
00148
00149 // If y < g(x), proceed to calculate the velocity components:
00150
00151 velocity = vfac* x * sqrt(2.0) * pow( 1.0 + radius*radius, -0.25);
00152 theta = acos(randinter(-1.0, 1.0));
00153 phi = randinter(0.0,TWO_PI);
00154
00155 bi->set_vel(vector(velocity * sin( theta ) * cos( phi ),
00156 velocity * sin( theta ) * sin( phi ),
00157 velocity * cos( theta )));
00158 }
00159
00160 // Transform to center-of-mass coordinates, and scale to standard
00161 // parameters, assuming a softening parameter of zero.
00162
00163 b->to_com();
00164
00165 putrq(b->get_log_story(), "initial_mass", 1.0);
00166
00167 if (!u_flag && n > 1) {
00168
00169 real kinetic, potential;
00170
00171 get_top_level_energies(b, 0.0, potential, kinetic);
00172 scale_virial(b, -0.5, potential, kinetic); // scales kinetic
00173 real energy = kinetic + potential;
00174 scale_energy(b, -0.25, energy); // scales energy
00175
00176 putrq(b->get_log_story(), "initial_rvirial", 1.0);
00177 putrq(b->get_dyn_story(), "initial_total_energy", -0.25);
00178 }
00179 }
00180
00181 local void swap(dyn* bi, dyn* bj)
00182 {
00183 if (bi != bj) {
00184
00185 real mass = bi->get_mass();
00186 vector pos = bi->get_pos();
00187 vector vel = bi->get_vel();
00188
00189 bi->set_mass(bj->get_mass());
00190 bi->set_pos(bj->get_pos());
00191 bi->set_vel(bj->get_vel());
00192
00193 bj->set_mass(mass);
00194 bj->set_pos(pos);
00195 bj->set_vel(vel);
00196 }
00197 }
00198
00199 // Original reshuffling code was O(N^2)! Modified by Steve, July 2000.
00200
00201 local void reshuffle_all(dyn* b, int n)
00202 {
00203 // Reorder the tree by swapping mass, pos, and vel of each node
00204 // with a randomly chosen node.
00205
00206 // First make a list of nodes.
00207
00208 int i = 0;
00209 dynptr *list = new dynptr[n];
00210 for_all_daughters(dyn, b, bi)
00211 list[i++] = bi;
00212
00213 // Then go down the list and randomize.
00214
00215 for (i = 0; i < n; i++)
00216 swap(list[i], list[(int)randinter(0, n)]);
00217
00218 delete [] list;
00219 }
00220
00221 //-----------------------------------------------------------------------------
00222 // main --
00223 // usage:
00224 // mkplummer -n# [options] ,
00225 //
00226 // where # is the number of bodies in the Nbody system.
00227 // options:
00228 // The following options are allowed:
00229 // seed:
00230 // -s # where # stands for an integer either in the range
00231 // [1,2147483647], or # = 0 which causes a seed to be
00232 // chosen as the value of the UNIX clock in seconds; this
00233 // guarantees that no two calls will give the same value for
00234 // the seed if they are more than 2.0 seconds apart.
00235 // cutoff:
00236 // -m # mass fraction of the (infinitely extended) Plummer model
00237 // -r # radius fraction of the (infinitely extended) Plummer model
00238 //
00239 // If the mass fraction and radius fraction are both 1.0 then
00240 // particles will be sprinkled in all over space. If mfrac < 1.0 or
00241 // rfrac < 1.0 then each particle is constrained to lie within
00242 // both the radial and (cumulative) mass bound.
00243 // For example, if rfrac( mfrac ) > rfrac then rfrac is the
00244 // limiting factor, but if rfrac( mfrac ) < rfrac then mfrac limits
00245 // the extent of the Plummer realization.
00246 //-----------------------------------------------------------------------------
00247
00248 main(int argc, char ** argv)
00249 {
00250 int i;
00251 int n;
00252 int input_seed, actual_seed;
00253 real mfrac;
00254 real rfrac;
00255
00256 bool c_flag = false;
00257 bool i_flag = false;
00258 bool m_flag = false;
00259 bool n_flag = false;
00260 bool o_flag = false;
00261 bool r_flag = false;
00262 bool R_flag = true;
00263 bool s_flag = false;
00264 bool u_flag = false;
00265
00266 char *comment;
00267 char seedlog[SEED_STRING_LENGTH];
00268
00269 check_help();
00270
00271 extern char *poptarg;
00272 int c;
00273 char* param_string = "c:im:n:os:r:Ru";
00274
00275 mfrac = rfrac = VERY_LARGE_NUMBER;
00276
00277 while ((c = pgetopt(argc, argv, param_string)) != -1)
00278 switch(c) {
00279
00280 case 'c': c_flag = TRUE;
00281 comment = poptarg;
00282 break;
00283 case 'i': i_flag = TRUE;
00284 break;
00285 case 'm': m_flag = TRUE;
00286 mfrac = atof(poptarg);
00287 break;
00288 case 'n': n_flag = TRUE;
00289 n = atoi(poptarg);
00290 break;
00291 case 'o': o_flag = TRUE;
00292 break;
00293 case 'r': r_flag = TRUE;
00294 rfrac = atof(poptarg);
00295 break;
00296 case 'R': R_flag = !R_flag;
00297 break;
00298 case 's': s_flag = TRUE;
00299 input_seed = atoi(poptarg);
00300 break;
00301 case 'u': u_flag = TRUE;
00302 break;
00303 case '?': params_to_usage(cerr, argv[0], param_string);
00304 get_help();
00305 exit(1);
00306 }
00307
00308 if (!n_flag) {
00309 cerr << "mkplummer: specify the number # of";
00310 cerr << " particles with -n#\n";
00311 exit(1);
00312 }
00313
00314 dyn *b, *by, *bo;
00315 b = new dyn();
00316
00317 if (n > 0) {
00318 bo = new dyn();
00319 if (i_flag)
00320 bo->set_label(1);
00321 b->set_oldest_daughter(bo);
00322 bo->set_parent(b);
00323 }
00324
00325 for (i = 1; i < n; i++) {
00326 by = new dyn();
00327 if (i_flag)
00328 by->set_label(i+1);
00329 by->set_parent(b);
00330 bo->set_younger_sister(by);
00331 by->set_elder_sister(bo);
00332 bo = by;
00333 }
00334
00335 if (c_flag == TRUE)
00336 b->log_comment(comment);
00337 b->log_history(argc, argv);
00338
00339 if (s_flag == FALSE)
00340 input_seed = 0; // default
00341 actual_seed = srandinter(input_seed);
00342
00343 if (o_flag) cerr << "mkplummer: random seed = " << actual_seed << endl;
00344
00345 sprintf(seedlog, " random number generator seed = %d",actual_seed);
00346 b->log_comment(seedlog);
00347
00348 if (m_flag == FALSE && r_flag == FALSE)
00349 mfrac = MFRAC_DEFAULT; // default
00350
00351 if (n != 0) {
00352 mkplummer(b, n, mfrac, rfrac, u_flag);
00353 if (R_flag) reshuffle_all(b, n);
00354 }
00355
00356 put_node(cout, *b);
00357 rmtree(b);
00358 }
00359
00360 #endif
00361
00362 // end of: mkplummer.C
1.2.6 written by Dimitri van Heesch,
© 1997-2001