makebinary.C Source File

00001 
00037 
00038 //   Steve McMillan, July 1996
00039 //   Modified by SPZ 21 April 1998
00040 
00041 #include "dyn.h"
00042 
00043 #ifdef TOOLBOX
00044 
00045 local void add_dynamics(dyn* cm, real ecc, real energy)
00046 {
00047     dyn* primary = cm->get_oldest_daughter();
00048     dyn* secondary = primary->get_younger_sister();
00049 
00050     real m_total = cm->get_mass();
00051     real m1 = primary->get_mass();
00052     real m2 = secondary->get_mass();
00053 
00054     // Set internal orbital elements:
00055 
00056     kepler k;
00057 
00058     real peri = 1; // Default value (unimportant unless ecc = 1).
00059     if (ecc == 1) peri = 0;
00060 
00061     // For now, binary phase is random.
00062 
00063     real mean_anomaly = randinter(-PI, PI);
00064 
00065     // Energies here are really binding energies (i.e. > 0 for a bound
00066     // orbit) per unit mass; kepler package expects energy < 0.
00067 
00068     energy = -energy;
00069 
00070     make_standard_kepler(k, 0, m_total, energy, ecc,
00071                          peri, mean_anomaly);
00072     //PRC(m_total);
00073     //PRC(energy);
00074     //PRL(ecc);
00075 
00076     set_random_orientation(k);
00077 
00078     // Set positions and velocities.
00079 
00080     primary->set_pos(-m2 * k.get_rel_pos() / m_total);
00081     primary->set_vel(-m2 * k.get_rel_vel() / m_total);
00082 
00083     secondary->set_pos(m1 * k.get_rel_pos() / m_total);
00084     secondary->set_vel(m1 * k.get_rel_vel() / m_total);
00085 }
00086 
00087 real zero_age_main_sequnece_radius(const real mass) {
00088 
00089     real alpha, beta, gamma, delta, kappa, lambda;
00090 
00091     real log_mass = log10(mass);
00092 
00093     if (mass > 1.334) {
00094 
00095         alpha = 0.1509 + 0.1709*log_mass;
00096         beta  = 0.06656 - 0.4805*log_mass;
00097         gamma = 0.007395 + 0.5083*log_mass;
00098         delta = (0.7388*pow(mass, 1.679) - 1.968*pow(mass, 2.887))
00099               / (1.0 - 1.821*pow(mass, 2.337));
00100     } 
00101     else {
00102       
00103         alpha = 0.08353 + 0.0565*log_mass;
00104         beta  = 0.01291 + 0.2226*log_mass;
00105         gamma = 0.1151 + 0.06267*log_mass;
00106         delta = pow(mass, 1.25) * (0.1148 + 0.8604*mass*mass)
00107               / (0.04651 + mass*mass);
00108     }
00109 
00110     real radius = delta;
00111 
00112     if(radius<0.1) {
00113 
00114         // assumed brown dwarf or planet.
00115         radius = 0.1;
00116     }
00117 
00118     return radius;
00119 }
00120 
00121 local real roche_radius(const real m1, const real m2) {
00122 
00123   real q = m1/m2;
00124   real q1_3 = pow(q, ONE_THIRD);
00125   real q2_3 = pow(q1_3, 2);   
00126   
00127   return 0.49*q2_3/(0.6*q2_3 + log(1 + q1_3));
00128 }
00129 
00130 
00131 local real minimum_semi_major_axis(dyn* b1, dyn* b2)
00132 {
00133     real ms_prim = b1->get_starbase()->conv_m_dyn_to_star(b1->get_mass());
00134     real ms_sec  = b2->get_starbase()->conv_m_dyn_to_star(b2->get_mass());
00135 
00136     // Use stellar mass as radius indicater.
00137     // mkbinary known little about stars
00138     real rs_prim = zero_age_main_sequnece_radius(ms_prim);
00139     real rs_sec = zero_age_main_sequnece_radius(ms_sec);
00140 
00141     // real rs_prim = b1->get_starbase()->conv_r_star_to_dyn(ms_prim);
00142     // real rs_sec  = b2->get_starbase()->conv_r_star_to_dyn(ms_sec);
00143   
00144     real sma_prim = rs_prim/roche_radius(ms_prim, ms_sec);
00145     real sma_sec = rs_sec/roche_radius(ms_sec, ms_prim);
00146 
00147     return max(sma_prim, sma_sec);
00148 }
00149 
00150 local bool dyn_present(dyn* bi) {
00151 
00152   bool dynamics_present = true;
00153 
00154   if(bi->get_pos()[0]==0 && bi->get_pos()[1]==0 && bi->get_pos()[2]==0 &&
00155      bi->get_vel()[0]==0 && bi->get_vel()[1]==0 && bi->get_vel()[2]==0)
00156     dynamics_present = false;
00157 
00158   return dynamics_present;
00159 }
00160 
00161 local void mkbinary(dyn* b, real lower, real upper,
00162                     int select, int option, real emax)
00163 {
00164     // Binary parameters are drawn from a 1/x distribution,
00165     // between the specified limits.
00166     // For now, use thermal eccentricities.
00167 
00168     for_all_daughters(dyn, b, bi)
00169         if (bi->get_oldest_daughter()) {
00170 
00171             dyn* primary = bi->get_oldest_daughter();
00172             dyn* secondary = primary->get_younger_sister();
00173             vector nul = 0;
00174             if(!dyn_present(primary) || !dyn_present(secondary)) {
00175 
00176             real m_total = bi->get_mass();
00177             real mu = primary->get_mass() * secondary->get_mass() / m_total;
00178 
00179             // Function select options:
00180             //
00181             // Choose binary parameters:  select = 1:  angular momentum
00182             //                            select = 2:  semi-major axis
00183             //                            select = 3:  energy
00184             //
00185             // 1. Select on angular momentum:
00186             //
00187             //   option = 1 ==> limits refer to angular momentum
00188             //   option = 2 ==> limits refer to angular momentum + detached
00189             //
00190             // 2. Select on semi-major axis:
00191             //
00192             //   option = 1 ==> limits refer to semi-major axis
00193             //   option = 2 ==> limits refer to semi-major axis + detached
00194             //   option = 3 ==> limits refer to pericenter & apocenter
00195             //                                                  + detached
00196             //
00197             // 3. Select on energy:
00198             //
00199             //   option = 1 ==> limits refer to binary energy.
00200             //   option = 2 ==> limits refer to binary energy per unit
00201             //                                                  reduced mass.
00202             //   option = 3 ==> limits refer to binary energy per unit
00203             //                                                  binary mass.
00204             //
00205             // However, add_dynamics expects energy per unit reduced mass,
00206             // so scale the energy in cases 1 and 3.
00207             //
00208             // Selection on angular momentum and semi-major axis are
00209             // determined iteratively using a do-loop.
00210             // Infinite looping is prohibited by an initial check.
00211             // Beware, however, such looping might be expensive/dangerous.
00212 
00213             real ecc, energy = 0;
00214 
00215             if (select == 1) {
00216 
00217                 real l_const, sma, angular_momentum;
00218 
00219                 if (option == 2) {
00220                     bool detached = false;
00221                     real a_min = minimum_semi_major_axis(primary, secondary);
00222                     a_min = b->get_starbase()->conv_r_star_to_dyn(a_min);
00223 
00224                     if (pow(upper, 2) <= a_min*m_total*(1-pow(ecc, 2))) {
00225                       PRC(upper);PRL(a_min*m_total*(1-pow(ecc, 2)));
00226                       err_exit("mkbinary: Illegal limits on angular momentum");
00227                     }
00228                     a_min = max(lower, a_min);
00229                     real frac = upper/lower;
00230 
00231                     do {
00232                         ecc = sqrt(randinter(0,emax));
00233                         l_const = log(upper) - log(a_min);
00234                         angular_momentum = a_min
00235                                          * pow(frac, randinter(0,1));
00236                         sma = pow(angular_momentum, 2)
00237                                 / (m_total*(1-pow(ecc, 2)));
00238                         if (sma*(1-ecc) > a_min) 
00239                           detached = true;
00240                     }
00241                     while(!detached);
00242                 }
00243                 else {
00244                     real frac = upper/lower;
00245                     ecc = sqrt(randinter(0,emax));
00246                     l_const = log(upper) - log(lower);
00247                     angular_momentum = lower*pow(frac, randinter(0,1));
00248                     sma = pow(angular_momentum, 2) / (m_total*(1-pow(ecc, 2)));
00249                 }
00250                 energy = 0.5 * m_total / sma;
00251             }
00252 
00253             else if (select == 2) {
00254 
00255                 real semi_major_axis, a_const;
00256 
00257                 if (option >= 2) {
00258 
00259                     bool detached = false;
00260                     real a_min = minimum_semi_major_axis(primary, secondary);
00261                     a_min = b->get_starbase()->conv_r_star_to_dyn(a_min);
00262 
00263                     //binding energy = 0.5 * m_total / sma;
00264                     real min_binding_energy = 10; // [kT]
00265                     real a_max = 0.5 * m_total/min_binding_energy;
00266                     if(upper>a_min) {
00267                       a_max = min(a_max, upper);
00268                     }
00269 
00270                     if(a_max<=a_min) {
00271                       PRC(a_max);PRL(a_min);
00272                       err_exit("mkbinary: Illegal limits on semi major axis");
00273                     }
00274                     a_min = max(lower, a_min);
00275 
00276                     if (option == 3) {// limits between peri- and apocenter.
00277                         real pericenter, apocenter;
00278                         do {
00279                             ecc = sqrt(randinter(0, emax));
00280                             pericenter = a_min*(1-ecc);
00281                             apocenter = a_max*(1+ecc);
00282                             a_const = log(a_max) - log(a_min);
00283                             semi_major_axis = a_min
00284                                             * exp(randinter(0., a_const));
00285                             if (semi_major_axis*(1-ecc) > a_min  &&
00286                                 (semi_major_axis > pericenter &&
00287                                  semi_major_axis < apocenter))
00288                                 detached = true;
00289                         }
00290                         while(!detached);
00291                     }
00292                     else {      // separation limited by semi-detached.
00293 
00294                         do {
00295                             ecc = sqrt(randinter(0, emax));
00296                             a_const = log(a_max) - log(a_min);
00297                             semi_major_axis = a_min
00298                                             * exp(randinter(0., a_const));
00299                             // if(semi_major_axis*(1-ecc)>a_min)
00300                             if(semi_major_axis*(1-pow(ecc, 2))>a_min)
00301                                 detached = true;
00302                         }
00303                         while(!detached);
00304                     }
00305                 }
00306                 else {
00307                     ecc = sqrt(randinter(0,emax));
00308                     a_const = log(upper) - log(lower);
00309                     semi_major_axis = lower*exp(randinter(0., a_const));
00310                     semi_major_axis = b->get_starbase()
00311                                        ->conv_r_star_to_dyn(semi_major_axis);
00312 
00313 
00314                 }
00315                 energy = 0.5 * m_total / semi_major_axis;
00316             }
00317 
00318             else {                              // default is select = 3
00319 
00320                 real frac = upper/lower;
00321                 ecc = sqrt(randinter(0,emax));
00322                 energy = lower * pow(frac, randinter(0,1));
00323                 if (option == 1)
00324                     energy /= mu;
00325                 else if (option == 3)
00326                     energy *= m_total / mu;
00327             }
00328 
00329             //PR(m_total);
00330             //PR(primary->get_mass());
00331             //PRL(secondary->get_mass());
00332             //PRL(energy);
00333             //real sma = 0.5 * (primary->get_mass() + secondary->get_mass())
00334             //               / energy;
00335             //real L = sqrt(m_total*sma*(1-pow(ecc, 2)));
00336             //PR(sma); PR(ecc); PRL(L);
00337 
00338             add_dynamics(bi, ecc, energy);
00339         }
00340         }
00341 }
00342 
00343 local void scale_limits(real& e1, real& e2, int option,
00344                         int N, real M, real K)
00345 {
00346     real scale = 1.0;
00347 
00348     // PRL(K);
00349 
00350     if (option == 3) {
00351 
00352         // Limits refer to kinetic energy per unit mass.  Scale to K/M.
00353 
00354         scale = K / M;
00355 
00356     } else if (option != 2) {
00357 
00358         // Limits refer to energy.  Scale to kT = (2/3) K/N.
00359 
00360         scale = K / (1.5*N);
00361         // PRL(scale);
00362     }
00363 
00364     e1 *= scale;
00365     e2 *= scale;
00366 }
00367 
00368 void main(int argc, char ** argv)
00369 {
00370   real lower = 1.0, upper = 1.0;   //Units depends on -o and -f
00371     real emax = 1;
00372     int function_select = 3;
00373     int option = 1;
00374     int random_seed = 0;
00375     char seedlog[64];
00376 
00377     check_help();
00378 
00379     extern char *poptarg;
00380     int c;
00381     char* param_string = "f:l:s:o:u:e:";
00382 
00383     while ((c = pgetopt(argc, argv, param_string)) != -1)
00384         switch(c) {
00385 
00386             case 'f': function_select = atoi(poptarg);
00387                       break;
00388             case 'l': lower = atof(poptarg);
00389                       break;
00390             case 'o': option= atoi(poptarg);
00391                       break;
00392             case 'e': emax = atof(poptarg);
00393                       break;
00394             case 's': random_seed = atoi(poptarg);
00395                       break;
00396             case 'u': upper = atof(poptarg);
00397                       break;
00398             case '?': params_to_usage(cerr, argv[0], param_string);
00399                       get_help();
00400                       exit(1);
00401         }
00402 
00403     if (lower < 0 || upper < 0 || lower > upper)
00404         err_exit("mkbinary: Illegal limits");
00405 
00406     dyn* b;
00407     b = get_dyn(cin);
00408 
00409     b->log_history(argc, argv);
00410 
00411     int actual_seed = srandinter(random_seed);
00412 
00413     sprintf(seedlog, "       random number generator seed = %d",
00414             actual_seed);
00415     b->log_comment(seedlog);
00416 
00417     // Look in the dyn story to see if a total system energy is known.
00418     // If it is, scale lower and upper accordingly.
00419 
00420     // Thus, in case function_select = 3, if the energy is known,
00421     // lower and upper are specified in "kT" units.  If the energy
00422     // is not known, lower and upper are absolute limits on the
00423     // binary kinetic energy.
00424 
00425     char* energy_string = "initial_total_energy";
00426 
00427     // Scale input to N-body units where appropirate
00428     // Possible options are:
00429     //   1: Angular momentum   [cgs]
00430     //   2: Semi-major axis    [Rsun]
00431     //   3: Binding energy     [N-body]
00432     if (function_select == 1) {
00433         if (b->get_starbase()->get_stellar_evolution_scaling()) {
00434             real scale = sqrt(b->get_starbase()->conv_m_star_to_dyn(1)
00435                               * b->get_starbase()->conv_r_star_to_dyn(1));
00436             lower *= scale;
00437             upper *= scale;
00438         }
00439         else
00440             cerr << "mkbinary: Using unscaled angular-momentum limits.\n";
00441     }
00442     else if (function_select == 2) {
00443         if (b->get_starbase()->conv_r_star_to_dyn(1)>0) {
00444           //        lower = b->get_starbase()->conv_r_star_to_dyn(lower);
00445           //        upper = b->get_starbase()->conv_r_star_to_dyn(upper);
00446           cerr << "mkbinary: Do not transform upper limit to dynamical units"
00447                << endl;
00448         }
00449         else
00450             cerr << "mkbinary: Using unscaled semi-major axis limits.\n";
00451     }
00452     else if (function_select == 3) {
00453         if (find_qmatch(b->get_dyn_story(), energy_string)) {
00454 
00455             // Use energy_string as an indicator that some energies
00456             // are known, but don't use its value as an estimator
00457             // of the kinetic energy, as it may include an external
00458             // field and also includes any CM motion of the system.
00459 
00460             // Recompute the top-level kinetic energy in the CM frame.
00461 
00462             vector com_pos, com_vel;
00463             compute_com(b, com_pos, com_vel);
00464 
00465             real kin = get_top_level_kinetic_energy(b)
00466                             - 0.5*b->get_mass()*square(com_vel);
00467 
00468             // PRL(get_top_level_kinetic_energy(b));
00469             // PRL(0.5*b->get_mass()*square(com_vel));
00470 
00471             scale_limits(lower, upper, option,
00472                          b->n_daughters(), b->get_mass(), kin);
00473 
00474         } else
00475             cerr << "mkbinary: Using unscaled energy limits.\n";
00476     }
00477 
00478     mkbinary(b, lower, upper, function_select, option, emax);
00479 
00480     put_dyn(cout, *b);
00481     rmtree(b);
00482 
00483 }
00484 #endif