makepowerlaw.C Source File

00001 
00044 
00045 #include "dyn.h"
00046 
00047 #ifdef TOOLBOX
00048 
00049 local void  makepowerlaw(dyn * root, int n,
00050                          real A, real a, real x,
00051                          real cutoff, real mass,
00052                          real r_min, real r_max)
00053 {
00054     real xi = 1/x, pmass = 1./n, armax = 0, vc;
00055     root->set_mass(1);
00056 
00057     if (a > 0) {
00058         armax = pow(a/r_max, x);
00059         vc =sqrt(3*A*pow(a, x-1)/(3-x));        // 3-D vrms at r = a
00060     } else
00061         vc = sqrt(3*A*pow(cutoff, x-1)/(3-x));  // 3-D vrms at r = cutoff
00062 
00063     for_all_daughters(dyn, root, bi) {
00064 
00065         bi->set_mass(pmass);
00066 
00067         // Radii are distributed between 0 and r_max, roughly uniformly
00068         // in r^x.  Power-law is OK if a << r_max; otherwise, better to
00069         // use a "core-halo" approximation to m(r):
00070         //
00071         //      m(r)  ~  (r/a)^3 (a/r_max)^x    (r < a)
00072         //               (r/r_max)^x            (r > a)
00073         //
00074         // Use this distribution even inside the cutoff radius.
00075 
00076         real radius = 0;
00077         while (radius < r_min) {        // lower limit for convenience only
00078             if (randinter(0, 1) < armax)
00079                 radius = a*pow(randinter(0, 1), 1./3);
00080             else
00081                 radius = r_max*pow(randinter(armax, 1), xi);
00082         }
00083 
00084         real costheta;
00085 
00086         // Details of orbits within the cutoff radius are to be determined.
00087         // For now, the orbits form a disk close to the central point mass.
00088 
00089         if (radius > cutoff)
00090 
00091             // Spherically symmetric.
00092 
00093             costheta = randinter(-1, 1);
00094 
00095         else {
00096 
00097             // Progressively force orbits into the x-y plane as r --> 0.
00098 
00099             real lim = radius/cutoff;
00100             costheta = randinter(-lim, lim);
00101         }
00102 
00103         real sintheta = 1 - costheta*costheta;
00104         if (sintheta > 0) sintheta = sqrt(sintheta);
00105         real phi = randinter(0.0, TWO_PI);
00106 
00107         bi->set_pos(radius*vector(sintheta * cos(phi),
00108                                   sintheta * sin(phi),
00109                                   costheta));
00110 
00111         // Choose velocities to ensure equilibrium in the external field.
00112         // Modify the velocity dispersion inside the cutoff to take the
00113         // central potential into account.
00114 
00115         real vrms;
00116 
00117         if (a > 0) {
00118 
00119             // Note:  scale > 0 ==> cutoff = 0.
00120 
00121             vrms = vc;
00122             if (radius > a) vrms *= pow(radius/a, (x-1)/2);
00123 
00124         } else {
00125 
00126             if (radius > cutoff) 
00127                 vrms = sqrt(3*A*pow(radius, x-1)/(3-x));
00128             else
00129                 // vrms = sqrt(mass/radius);
00130                 vrms = max(sqrt(mass/radius), vc);
00131         }
00132 
00133         // Play with orbital orientations.  At present, orbits become more
00134         // nearly planar and circular as we approach the center.
00135 
00136         if (radius > cutoff) {
00137 
00138             // Isotropic.
00139 
00140             bi->set_vel(vrms*vector(randinter(-1,1),
00141                                     randinter(-1,1),
00142                                     randinter(-1,1)));
00143         } else {
00144 
00145             // Some handy unit vectors.
00146 
00147             vector rhat = bi->get_pos()/radius;
00148             vector zhat = vector(0,0,1);
00149 
00150             vector xyhat = rhat^zhat;           // transverse unit vector
00151             xyhat /= abs(xyhat);                // in the x-y plane
00152 
00153             real lim = sqrt(radius/cutoff);
00154             real rfrac = randinter(-lim, lim);
00155             real zfrac = randinter(-lim, lim);
00156             real tfrac = sqrt(max(0.0, 1-rfrac*rfrac-zfrac*zfrac));
00157 
00158             bi->set_vel(vrms*(rfrac*rhat+tfrac*xyhat+zfrac*zhat));
00159         }
00160     }
00161 
00162     putrq(root->get_log_story(), "initial_mass", 1.0);
00163 }
00164 
00165 #define  SEED_STRING_LENGTH  60
00166 
00167 main(int argc, char ** argv) {
00168     int  i;
00169     int  n;
00170     int  input_seed, actual_seed;
00171     int  a_flag = false;
00172     int  c_flag = false;
00173     int  e_flag = false;
00174     int  i_flag = false;
00175     int  m_flag = false;
00176     int  n_flag = false;
00177     int  o_flag = false;
00178     int  s_flag = false;
00179 
00180     char  *comment;
00181     char  seedlog[SEED_STRING_LENGTH];
00182 
00183     real coeff = 1, scale = 1, exponent = 1;
00184     vector center = 0;
00185 
00186     real cutoff = 0, mass = 0, softening = 0;
00187 
00188     real r_max = 10, r_min = 0.1;
00189 
00190     check_help();
00191 
00192     extern char *poptarg;
00193     int c;
00194     char* param_string = "A:a:b:c:e:il:M:m:n:oR:r:s:x:";
00195 
00196     while ((c = pgetopt(argc, argv, param_string)) != -1)
00197         switch(c) {
00198             case 'A': coeff = atof(poptarg);
00199                       break;
00200             case 'a':
00201             case 'R':
00202                       a_flag = true;
00203                       scale = atof(poptarg);
00204                       break;
00205             case 'b': cutoff = atof(poptarg);
00206                       break;
00207             case 'c': c_flag = true;
00208                       comment = poptarg;
00209                       break;
00210             case 'e': e_flag = true;
00211                       softening = atof(poptarg);
00212                       break;
00213             case 'i': i_flag = true;
00214                       break;
00215             case 'l': r_min = atof(poptarg);
00216                       break;
00217             case 'M':
00218             case 'm': m_flag = true;
00219                       mass = atof(poptarg);
00220                       break;
00221             case 'n': n_flag = true;
00222                       n = atoi(poptarg);
00223                       break;
00224             case 'o': o_flag = true;
00225                       break;
00226             case 'r': r_max = atof(poptarg);
00227                       break;
00228             case 's': s_flag = true;
00229                       input_seed = atoi(poptarg);
00230                       break;
00231             case 'x': exponent = atof(poptarg);
00232                       break;
00233             case '?': params_to_usage(cerr, argv[0], param_string);
00234                       get_help();
00235                       exit(1);
00236         }            
00237     
00238     if (!n_flag) {
00239         cerr << "makepowerlaw: must specify the number # of";
00240         cerr << " particles with -n#" << endl;
00241         exit(1);
00242     }
00243     
00244     if (n < 1) {
00245         cerr << "makepowerlaw: n < 1 not allowed" << endl;
00246         exit(1);
00247     }
00248 
00249     if (exponent <= 0 || exponent >= 3) {
00250         cerr << "makepowerlaw: 0 < x < 3 required" << endl;
00251         exit(1);
00252     }
00253 
00254     if (cutoff > 0) {
00255         if (!m_flag) mass = coeff*pow(cutoff, exponent);
00256         if (!e_flag) softening = cutoff/1000;
00257         if (a_flag)
00258             cerr << "makepowerlaw: scale = " << scale
00259                  << " inconsistent with cutoff = " << cutoff
00260                  << "; setting scale = 0" << endl;
00261         scale = 0;
00262     }
00263 
00264     // Create a linked list of (labeled) nodes.
00265 
00266     dyn *b, *by, *bo;
00267 
00268     b = new dyn();
00269     if (i_flag) b->set_label("root");
00270 
00271     bo = new dyn();
00272     if (i_flag) bo->set_label(1);
00273     b->set_oldest_daughter(bo);
00274     bo->set_parent(b);
00275 
00276     for (i = 1; i < n; i++) {
00277         by = new dyn();
00278         if (i_flag) by->set_label(i+1);
00279         bo->set_younger_sister(by);
00280         by->set_elder_sister(bo);
00281         bo = by;
00282     }
00283 
00284     // Add comments, etc.
00285 
00286     if (c_flag) b->log_comment(comment);
00287 
00288     b->log_history(argc, argv);
00289 
00290     if (!s_flag) input_seed = 0;
00291     actual_seed = srandinter(input_seed);
00292 
00293     if (o_flag) cerr << "makepowerlaw: random seed = " << actual_seed << endl;
00294 
00295     sprintf(seedlog, "       random number generator seed = %d",actual_seed);
00296     b->log_comment(seedlog);
00297 
00298     // Create dyn data.
00299 
00300     makepowerlaw(b, n, coeff, scale, exponent, cutoff, mass, r_min, r_max);
00301 
00302     // Flag actions for use by kira.
00303 
00304     putrq(b->get_log_story(), "kira_pl_coeff", coeff);
00305     putrq(b->get_log_story(), "kira_pl_exponent", exponent);
00306     putrq(b->get_log_story(), "kira_pl_scale", scale);
00307     putvq(b->get_log_story(), "kira_pl_center", center);
00308 
00309     putrq(b->get_log_story(), "kira_pl_cutoff", cutoff);
00310     putrq(b->get_log_story(), "kira_pl_mass", mass);
00311     putrq(b->get_log_story(), "kira_pl_softening", softening);
00312 
00313     putiq(b->get_log_story(), "ignore_internal", 1);
00314     putrq(b->get_log_story(), "r_reflect", r_max);
00315 
00316     put_node(cout, *b);
00317 }
00318 
00319 #endif
00320 
00321 /* end of: makepowerlaw.c */