make_aniso_king.C Source File

00001 
00002        //=======================================================//    _\|/_
00003       //  __  _____           ___                    ___       //      /|\ ~
00004      //  /      |      ^     |   \  |         ^     |   \     //          _\|/_
00005     //   \__    |     / \    |___/  |        / \    |___/    //            /|\ ~
00006    //       \   |    /___\   |  \   |       /___\   |   \   // _\|/_
00007   //     ___/   |   /     \  |   \  |____  /     \  |___/  //   /|\ ~
00008  //                                                       //            _\|/_
00009 //=======================================================//              /|\ ~
00010 
00045 
00046 // version 1:  Jul 1998     Steve McMillan
00047 //             Jan 1999     steve@zonker.drexel.edu
00048 //                          Physics Dept., Drexel University, Phila PA, USA
00049 //                          (Wrapper for D.C. Heggie's Fortran 77 program)
00050 //             Jun 1999     Steve McMillan: added non-point-mass fields
00051 //                          to this program and kira
00052 
00053 #include "dyn.h"
00054 
00055 #ifdef TOOLBOX
00056 
00057 #define  SEED_STRING_LENGTH  256
00058 char  tmp_string[SEED_STRING_LENGTH];
00059 
00060 #ifdef FORTRAN_TRAILING_UNDERSCORE
00061 #   define AKING aking_
00062 #   define RAND ran2_ 
00063 #else
00064 #   define AKING aking
00065 #   define RAND ran2
00066 #endif
00067 
00068 // Should really be 'extern "F77"'...
00069 
00070 extern "C" void AKING(real*, real*, real*, real*, real*, real*, real*,
00071                       real*, real*, int*, int*, real*, int*,
00072                       real*, real*, real*, real*);
00073 
00074 extern "C" real RAND(int &seed)
00075 {
00076     return randinter(0, 1);
00077 }
00078 
00079 local void mk_aniso_king(dyn * b, int n, real w0, real alpha1, real alpha3,
00080                          bool u_flag, int test, int seed)
00081 
00082 // Create an anisotropic King model using Heggie's Fortran code, and
00083 // optionally initialize an N-body system with G = 1, total mass = 1,
00084 // and length scale set by alpha1.
00085 
00086 {
00087     // These tests are redundant: they duplicate the tests at the
00088     // start of aking.
00089 
00090     if (alpha3 != 0 && alpha1 >= 0)
00091         err_exit("mk_aniso_king: must specify alpha1 < 0");
00092     if (w0 < 1)
00093         err_exit("mk_aniso_king: must specify w0 > 1");
00094     if (w0 > 12)
00095         err_exit("mk_aniso_king: must specify w0 < 12");
00096 
00097     // Compute the cluster density/velocity/potential profile
00098 
00099     real* m = new real[n];
00100     real* x = new real[n];
00101     real* y = new real[n];
00102     real* z = new real[n];
00103     real* vx = new real[n];
00104     real* vy = new real[n];
00105     real* vz = new real[n];
00106 
00107     // Enforce limit on seed required by internal random-number generator
00108     // in aking (ran2, from Fortran Numerical Recipes, 1ed).
00109 
00110     seed %= 150999;
00111 
00112     real potential, kinetic, tidal_energy;      // will be computed in aking;
00113                                                 // no point recomputing...
00114 
00115     real* coord = new real[3*n];                // workspace for aking...
00116                                                 // note that this duplicates
00117                                                 // the x, y, and z arrays
00118 
00119     AKING(m, x, y, z, vx, vy, vz,                       // DCH's Fortran
00120           &alpha1, &alpha3, &seed, &n, &w0, &test,      // function
00121           &potential, &kinetic, &tidal_energy, coord);
00122 
00123     if (b == NULL || n < 1) return;
00124 
00125     // Initialize the N-body system.
00126 
00127     sprintf(tmp_string,
00128     "       Anisotropic King model, w0 = %.2f, alpha1 = %f, alpha3 = %f",
00129             w0, alpha1, alpha3);
00130     b->log_comment(tmp_string);
00131 
00132     // Assign positions and velocities.
00133 
00134     int i = 0;
00135     for_all_daughters(dyn, b, bi) {
00136 
00137         bi->set_mass(m[i]);
00138         bi->set_pos(vector(x[i], y[i], z[i]));
00139         bi->set_vel(vector(vx[i], vy[i], vz[i]));
00140 
00141         i++;
00142     }
00143 
00144     delete m, x, y, z, vx, vy, vz;
00145 
00146     // System is in virial equilibrium in a consistent set of units
00147     // with G and total mass = 1 and length scale determined by the
00148     // fact that the cluster exactly fills its Roche lobe in the
00149     // external field.
00150 
00151     // Transform to center-of-mass coordinates.
00152     // Note: must recompute the kinetic energy and the tidal potential.
00153 
00154     b->to_com();        
00155 
00156     kinetic = get_kinetic_energy(b);
00157     tidal_energy = get_tidal_pot(b);
00158 
00159     // Definition of r_virial assumes GM = 1 and *does not include* the
00160     // tidal potential.
00161 
00162     // real r_virial = -0.5/(potential + tidal_energy);
00163     real r_virial = -0.5/potential;
00164 
00165     // Fake "Jacobi radius" is used to transmit tidal (alpha1) information
00166     // to kira -- r_jacobi actually is the Jacobi radius only for a 1/r
00167     // cluster potential (actually, not a bad approximation).  Assumes GM = 1.
00168 
00169     real r_jacobi = pow(-alpha1, -1.0/3);
00170 
00171     // Optionally scale to standard parameters.  Scaling is equivalent
00172     // to using "scale -s" as a tool.
00173 
00174     if (!u_flag && n > 1) {
00175 
00176         // Scaling is OK because cluster is Roche-lobe filling.
00177         // Function scale_virial adjusts kinetic.
00178 
00179         scale_virial(b, -0.5, potential, kinetic);  // potential + tidal_energy
00180 
00181         real energy = kinetic + potential + tidal_energy;
00182 
00183         // Scale_energy uses the specified energy to rescale velocities,
00184         // adjusts energy accordingly, and returns the factor by which
00185         // the positions were scaled.
00186 
00187         real fac = scale_energy(b, -0.25, energy);
00188 
00189         // Note: when recomputing energies for test purposes, we must
00190         // remember to rescale alpha1,3 to preserve Roche-lobe filling.
00191 
00192         alpha1 /= pow(fac,3);
00193         alpha3 /= pow(fac,3);
00194 
00195         sprintf(tmp_string,
00196                 "       Rescaled alpha1 = %f, alpha3 = %f", alpha1, alpha3);
00197         b->log_comment(tmp_string);
00198 
00199         // Recompute other quantities mainly for completeness.
00200 
00201         r_virial *= fac;                // should be 1
00202         kinetic /= fac;                 // should be 0.25
00203         potential /= fac;               // should be -0.5
00204         tidal_energy /= fac;
00205         r_jacobi *= fac;
00206 
00207         // Story output mimics mkking where possible.
00208 
00209         putrq(b->get_dyn_story(), "initial_total_energy", -0.25);
00210     }
00211 
00212     // Write essential model information to the root dyn story.
00213     // Story output mimics mkking where possible.
00214 
00215     putrq(b->get_log_story(), "initial_mass", 1.0);
00216     putrq(b->get_log_story(), "initial_rvirial", r_virial);
00217     putrq(b->get_log_story(), "initial_rtidal_over_rvirial",
00218           r_jacobi/r_virial, 10);
00219 
00220     cerr << " mk_aniso_king: "; PRL(r_jacobi/r_virial);
00221 
00222     // Note that, for this model, kira *must* use this Jacobi
00223     // radius to be consistent.
00224 
00225     // Additional information (indicator of anisotropic King model):
00226 
00227     putrq(b->get_log_story(), "alpha3_over_alpha1", alpha3/alpha1, 10);
00228 }
00229 
00230 #define OMEGA_2 3.e-4
00231 
00232 // Express Galactic parameters in "stellar" units (see Mihalas 1968):
00233 //
00234 //      G               =  1
00235 //      length unit     =  1 pc
00236 //      velocity unit   =  1 km/s
00237 //
00238 // ==>  time unit       =  0.978 Myr
00239 //      mass unit       =  232 Msun
00240 
00241 #define OORT_A  (0.0144)        // km/s/pc
00242 #define OORT_B  (-0.012)        // km/s/pc
00243 #define RHO_G   (0.11/232)      // (232 Msun)/pc^3
00244 
00245 #define OORT_ALPHA_RATIO \
00246         ((4*M_PI*RHO_G + 2*(OORT_A*OORT_A - OORT_B*OORT_B)) \
00247                          / (-4*OORT_A*(OORT_A-OORT_B)))
00248 
00249 main(int argc, char ** argv)
00250 {
00251     real w0;
00252     int  n = 0;
00253 
00254     real Oort_A = OORT_A;       // km/s/pc
00255     real Oort_B = OORT_B;       // km/s/pc
00256     real rho_G  = RHO_G;        // (232 Msun)/pc^3
00257 
00258     // Default tidal parameters are appropriate to a cluster in a
00259     // a circular orbit of radius D = 1000 length units around a
00260     // point-mass "galaxy" of mass Mg = 3.0e5
00261     //
00262     //      alpha1  = -3 G Mg / D^3  = -3 Omega^2  (Binney & Tremaine, p. 450)
00263     //      alpha3  =    G Mg / D^3  =    Omega^2
00264         
00265     real alpha1 = -3 * OMEGA_2;
00266     real alpha3 = OMEGA_2;
00267     real alpha3_over_alpha1 = alpha3/alpha1;
00268     bool a_flag = false;
00269     bool b_flag = false;
00270 
00271     int  input_seed, actual_seed;
00272 
00273     bool A_flag = false;
00274     bool B_flag = false;
00275     bool G_flag = false;
00276 
00277     bool c_flag = false;
00278     bool i_flag = false;
00279     bool n_flag = false;
00280     bool o_flag = false;
00281     bool s_flag = false; 
00282     bool w_flag = false;
00283     bool u_flag = false;
00284 
00285     int tidal_type = 1;
00286     bool F_flag = false;
00287 
00288     int test = 0;
00289     char  *comment;
00290 
00291     check_help();
00292 
00293     extern char *poptarg;
00294     int c;
00295     char* param_string = "A:a:B:b:c:F:G:in:os:T:uw:";
00296 
00297     while ((c = pgetopt(argc, argv, param_string)) != -1)
00298         switch(c) {
00299             case 'A': Oort_A = atof(poptarg);
00300                       A_flag = true;
00301                       break;
00302             case 'B': Oort_B = atof(poptarg);
00303                       B_flag = true;
00304                       break;
00305             case 'G': rho_G = atof(poptarg);
00306                       G_flag = true;
00307                       break;
00308             case 'a': alpha1 = atof(poptarg);
00309                       a_flag = true;
00310                       break;
00311             case 'b': alpha3_over_alpha1 = atof(poptarg);
00312                       b_flag = true;
00313                       F_flag = false;
00314                       break;
00315             case 'c': c_flag = true; 
00316                       comment = poptarg;
00317                       break;
00318             case 'F': tidal_type = atoi(poptarg);
00319                       b_flag = false;
00320                       F_flag = true;
00321                       break;
00322             case 'i': i_flag = true;
00323                       break;
00324             case 'n': n_flag = true;
00325                       n = atoi(poptarg);
00326                       break;
00327             case 'o': o_flag = true;
00328                       break;
00329             case 's': s_flag = true;
00330                       input_seed = atoi(poptarg);
00331                       break;
00332             case 'T': test = atoi(poptarg);
00333                       break;
00334             case 'u': u_flag = true;
00335                       break;
00336             case 'w': w_flag = true;
00337                       w0 = atof(poptarg);
00338                       break;
00339             case '?': params_to_usage(cerr, argv[0], param_string);
00340                       get_help();
00341         }
00342 
00343     if (!w_flag) {
00344         cerr << "mk_aniso_king: please specify the dimensionless depth";
00345         cerr << " with -w #\n";
00346         exit(1);
00347     }
00348 
00349 //    if (test == 0) {
00350 //        cerr << "mk_aniso_king: please specify the number # of";
00351 //        cerr << " particles with -n #\n";
00352 //        exit(1);
00353 //    }
00354 
00355     if (n < 0)
00356         err_exit("mk_aniso_king: n > 0 required");
00357 
00358     if (A_flag && B_flag && G_flag) {
00359 
00360       cerr << " mk_aniso_king: Custom tidal field with Oort constants" 
00361            << endl;
00362       PRI(4);PRC(Oort_A);PRC(Oort_B);PRL(rho_G);
00363       tidal_type = 4;
00364       F_flag = true;
00365 
00366     }
00367 
00368     if (F_flag) {
00369 
00370         if (tidal_type == 1)
00371 
00372             alpha3_over_alpha1 = -1.0/3;
00373 
00374         else if (tidal_type == 2)
00375 
00376             alpha3_over_alpha1 = -0.5;
00377 
00378         else if (tidal_type == 3) {
00379 
00380             // An odd way to do things...
00381 
00382             alpha1 = -4*OORT_A*(OORT_A-OORT_B);
00383             alpha3_over_alpha1 = OORT_ALPHA_RATIO;
00384             a_flag = true;
00385 
00386         } else if(tidal_type == 4) {
00387           
00388             alpha1 = -4*Oort_A*(Oort_A-Oort_B);
00389             alpha3_over_alpha1 = ((4*M_PI*rho_G 
00390                                    + 2*(Oort_A*Oort_A - Oort_B*Oort_B)) \
00391                                / (-4*Oort_A*(Oort_A-Oort_B)));
00392 
00393             a_flag = true;
00394         }
00395         else
00396             F_flag = false;
00397 
00398         if (F_flag) {
00399             cerr << " mk_aniso_king: tidal_type = " << tidal_type;
00400             if (tidal_type >= 3) cerr << ", alpha1 = " << alpha1;
00401             cerr << ", alpha3/alpha1 = " << alpha3_over_alpha1
00402                  << endl;
00403             b_flag = true;
00404         }
00405     }
00406 
00407     if (!b_flag && !F_flag)
00408         cerr << " mk_aniso_king: using default alpha3/alpha1 = "
00409              << alpha3_over_alpha1 << endl;
00410 
00411     alpha3 = alpha3_over_alpha1 * alpha1;
00412 
00413     if (!a_flag) {
00414         cerr << " mk_aniso_king: using default "; PR(alpha1);
00415         if (!b_flag) {
00416             cerr << ", "; PR(alpha3);
00417         }
00418         cerr << endl;
00419     }
00420 
00421     dyn *b = NULL;
00422 
00423     if (!n_flag) {
00424 
00425         b = get_dyn(cin);
00426         n = b->n_leaves();
00427     }
00428     else {
00429 
00430         b = new dyn();
00431         dyn *by, *bo;
00432 
00433         bo = new dyn();
00434         if (i_flag)
00435             bo->set_label(1);
00436         b->set_oldest_daughter(bo); 
00437         bo->set_parent(b);
00438 
00439         for (int i = 1; i < n; i++) {
00440             by = new dyn();
00441             if (i_flag)
00442                 by->set_label(i+1);
00443             by->set_parent(b);
00444             bo->set_younger_sister(by);
00445             by->set_elder_sister(bo);
00446             bo = by;
00447         }
00448     }
00449 
00450     if (c_flag)
00451         b->log_comment(comment);                // add comment to story
00452     b->log_history(argc, argv);
00453 
00454     if (s_flag == false) input_seed = 0;        // default
00455     actual_seed = srandinter(input_seed);
00456 
00457     if (o_flag)
00458         cerr << "mk_aniso_king: random seed = " << actual_seed << endl;
00459 
00460     sprintf(tmp_string,
00461             "       random number generator seed = %d",
00462             actual_seed);
00463     b->log_comment(tmp_string);
00464 
00465     mk_aniso_king(b, n, w0, alpha1, alpha3, u_flag, test, actual_seed);
00466 
00467     b->set_tidal_field(alpha1, alpha3);
00468     putiq(b->get_log_story(), "kira_tidal_field_type", tidal_type);
00469 
00470     if(tidal_type == 4) {
00471         if (A_flag) putrq(b->get_log_story(), "Oort_A_constant", Oort_A, 10);
00472         if (B_flag) putrq(b->get_log_story(), "Oort_B_constant", Oort_B, 10);
00473         if (G_flag) putrq(b->get_log_story(), "local_mass_density", rho_G, 10);
00474     }
00475 
00476     if (n_flag)
00477         put_node(cout, *b);
00478 }
00479 
00480 #endif
00481 
00482 /* end of: mk_aniso_king.C */