makeplanetary.C Source File

00001 
00035 
00036 //   Steve McMillan, July 1996
00037 //   Modified by SPZ 21 April 1998
00038 
00039 #include "dyn.h"
00040 
00041 #ifdef TOOLBOX
00042 
00043 local void add_dynamics(dyn* cm, real ecc, real energy)
00044 {
00045     dyn* primary = cm->get_oldest_daughter();
00046     dyn* secondary = primary->get_younger_sister();
00047 
00048     real m_total = cm->get_mass();
00049     real m1 = primary->get_mass();
00050     real m2 = secondary->get_mass();
00051 
00052     // Set internal orbital elements:
00053 
00054     kepler k;
00055 
00056     real peri = 1; // Default value (unimportant unless ecc = 1).
00057     if (ecc == 1) peri = 0;
00058 
00059     // For now, binary phase is random.
00060 
00061     real mean_anomaly = randinter(-PI, PI);
00062 
00063     // Energies here are really binding energies (i.e. > 0 for a bound
00064     // orbit) per unit mass; kepler package expects energy < 0.
00065 
00066     energy = -energy;
00067 
00068     make_standard_kepler(k, 0, m_total, energy, ecc,
00069                          peri, mean_anomaly);
00070     //PRC(m_total);
00071     //PRC(energy);
00072     //PRL(ecc);
00073 
00074     set_random_orientation(k);
00075 
00076     // Set positions and velocities.
00077 
00078     primary->set_pos(-m2 * k.get_rel_pos() / m_total);
00079     primary->set_vel(-m2 * k.get_rel_vel() / m_total);
00080 
00081     secondary->set_pos(m1 * k.get_rel_pos() / m_total);
00082     secondary->set_vel(m1 * k.get_rel_vel() / m_total);
00083 }
00084 
00085 local real minimum_semi_major_axis(dyn* b1, dyn* b2)
00086 {
00087     real ms_prim = b1->get_starbase()->conv_m_dyn_to_star(b1->get_mass());
00088     real ms_sec  = b2->get_starbase()->conv_m_dyn_to_star( b2->get_mass());
00089     real rs_prim = b1->get_starbase()->conv_r_star_to_dyn(ms_prim);
00090     real rs_sec  = b2->get_starbase()->conv_r_star_to_dyn(ms_sec);
00091     real ms_tot  = ms_prim + ms_sec;
00092   
00093     real sma_prim = 2.2*rs_prim*pow(ms_tot/ms_prim, ONE_THIRD);
00094     real sma_sec  = 2.2*rs_sec*pow(ms_tot/ms_sec, ONE_THIRD);
00095     // a_min is Roche-lobe limited.
00096 
00097     return min(sma_prim, sma_sec);
00098 }
00099 
00100 local void add_secondary(node* original, real mass_ratio)
00101 {
00102     node* primary = new node;
00103     node* secondary = new node;
00104 
00105     // Add new links.
00106 
00107     original->set_oldest_daughter(primary);
00108 
00109     primary->set_parent(original);
00110     secondary->set_parent(original);
00111 
00112     primary->set_younger_sister(secondary);
00113     secondary->set_elder_sister(primary);
00114 
00115     // Set new masses.
00116 
00117     primary->set_mass(original->get_mass());
00118     secondary->set_mass(mass_ratio*original->get_mass());
00119     original->inc_mass(secondary->get_mass());
00120 
00121     // Naming convention:
00122 
00123     if (original->get_name() == NULL)
00124         if (original->get_index() >= 0) {
00125             char tmp[64];
00126             sprintf(tmp, "%d", original->get_index());
00127             original->set_name(tmp);
00128         }
00129 
00130     if (original->get_name() != NULL) {
00131 
00132         // Label components "a" and "b".
00133 
00134         primary->set_name(original->get_name());
00135         secondary->set_name(original->get_name());
00136         strcat(primary->get_name(), "a");
00137         strcat(secondary->get_name(), "b");
00138 
00139         // Make standard CM name.
00140 
00141         char tmp[256];
00142         sprintf(tmp, "(%s,%s)", primary->get_name(), secondary->get_name());
00143         original->set_name(tmp);
00144     }
00145 }
00146 
00147 local void mkbinary(dyn* b, real lower, real upper,
00148                     int select, int option, real emax)
00149 {
00150     // Binary parameters are drawn from a 1/x distribution,
00151     // between the specified limits.
00152     // For now, use thermal eccentricities.
00153 
00154     real frac = upper/lower;
00155 
00156     for_all_daughters(dyn, b, bi)
00157         if(!bi->is_parent()) {
00158             add_secondary(bi, 0.0001);
00159         }
00160 
00161     for_all_daughters(dyn, b, bi)
00162         if (bi->get_oldest_daughter()) {
00163 
00164             dyn* primary = bi->get_oldest_daughter();
00165             dyn* secondary = primary->get_younger_sister();
00166             real m_total = bi->get_mass();
00167             real mu = primary->get_mass() * secondary->get_mass() / m_total;
00168 
00169             // Function select options:
00170             //
00171             // Choose binary parameters:  select = 1:  angular momentum
00172             //                            select = 2:  semi-major axis
00173             //                            select = 3:  energy
00174             //
00175             // 1. Select on angular momentum:
00176             //
00177             //   option = 1 ==> limits refer to angular momentum
00178             //   option = 2 ==> limits refer to angular momentum + detached
00179             //
00180             // 2. Select on semi-major axis:
00181             //
00182             //   option = 1 ==> limits refer to semi-major axis
00183             //   option = 2 ==> limits refer to semi-major axis + detached
00184             //   option = 3 ==> limits refer to pericenter & apocenter
00185             //                                                  + detached
00186             //
00187             // 3. Select on energy:
00188             //
00189             //   option = 1 ==> limits refer to binary energy.
00190             //   option = 2 ==> limits refer to binary energy per unit
00191             //                                                  reduced mass.
00192             //   option = 3 ==> limits refer to binary energy per unit
00193             //                                                  binary mass.
00194             //
00195             // However, add_dynamics expects energy per unit reduced mass,
00196             // so scale the energy in cases 1 and 3.
00197             //
00198             // Selection on angular momentum and semi-major axis are
00199             // determined iteratively using a do-loop.
00200             // Infinite looping is prohibited by an initial check.
00201             // Beware, however, such looping might be expensive/dangerous.
00202 
00203             real ecc, energy = 0;
00204 
00205             if (select == 1) {
00206 
00207                 real l_const, sma, angular_momentum;
00208 
00209                 if (option == 2) {
00210                     bool detached = false;
00211                     real a_min = minimum_semi_major_axis(primary, secondary);
00212                     if (pow(upper, 2) <= a_min*m_total*(1-pow(ecc, 2)))
00213                         err_exit(
00214                             "mkbinary: Illegal limits on angular momentum");
00215                     do {
00216                         ecc = sqrt(randinter(0,emax));
00217                         l_const = log(upper) - log(lower);
00218                         angular_momentum = lower * pow(frac, randinter(0,1));
00219                         sma = pow(angular_momentum, 2)
00220                                 / (m_total*(1-pow(ecc, 2)));
00221                         if (sma*(1-pow(ecc, 2)) > a_min) detached = true;
00222                     }
00223                     while(!detached);
00224                 }
00225                 else {
00226                     ecc = sqrt(randinter(0,emax));
00227                     l_const = log(upper) - log(lower);
00228                     angular_momentum = lower * pow(frac, randinter(0,1));
00229                     sma = pow(angular_momentum, 2) / (m_total*(1-pow(ecc, 2)));
00230                 }
00231                 energy = 0.5 * m_total / sma;
00232             }
00233 
00234             else if (select == 2) {
00235 
00236                 real semi_major_axis, a_const;
00237 
00238                 if (option >= 2) {
00239                     bool detached = false;
00240                     real a_min = minimum_semi_major_axis(primary, secondary);
00241                     if(upper<=a_min)
00242                         err_exit(
00243                             "mkbinary: Illegal limits on angular momentum");
00244 
00245                     if (option == 3) {// limits between peri- and apocenter.
00246                         real pericenter, apocenter;
00247                         do {
00248                             ecc = sqrt(randinter(0,emax));
00249                             pericenter = lower*(1-ecc);
00250                             apocenter = upper*(1+ecc);
00251                             a_const = log(upper) - log(lower);
00252                             semi_major_axis = lower*exp(randinter(0., a_const));
00253                             if (semi_major_axis*(1-pow(ecc, 2)) > a_min  &&
00254                                 (semi_major_axis > pericenter &&
00255                                  semi_major_axis < apocenter))
00256                                 detached = true;
00257                         }
00258                         while(!detached);
00259                     }
00260                     else {      // separation limited by semi-detached.
00261 
00262                         // Assume for now that
00263                         //      stellar radius [Rsun] = mass [Msun].
00264                         // This is of course not correct, but for the moment
00265                         // good enough.  mkbinary does not know much about
00266                         // stars.
00267                         do {
00268                             ecc = sqrt(randinter(0,emax));
00269                             a_const = log(upper) - log(lower);
00270                             semi_major_axis = lower*exp(randinter(0., a_const));
00271                             if(semi_major_axis*(1-pow(ecc, 2))>a_min)
00272                                 detached = true;
00273                         }
00274                         while(!detached);
00275                     }
00276                 }
00277                 else {
00278                     ecc = sqrt(randinter(0,emax));
00279                     a_const = log(upper) - log(lower);
00280                     semi_major_axis = lower*exp(randinter(0., a_const));
00281                 }
00282                 energy = 0.5 * m_total / semi_major_axis;
00283             }
00284 
00285             else {                              // default is select = 3
00286 
00287                 ecc = sqrt(randinter(0,emax));
00288                 energy = lower * pow(frac, randinter(0,1));
00289                 if (option == 1)
00290                     energy /= mu;
00291                 else if (option == 3)
00292                     energy *= m_total / mu;
00293             }
00294 
00295             //PR(m_total);
00296             //PR(primary->get_mass());
00297             //PRL(secondary->get_mass());
00298             //PRL(energy);
00299             //real sma = 0.5 * (primary->get_mass() + secondary->get_mass())
00300             //               / energy;
00301             //real L = sqrt(m_total*sma*(1-pow(ecc, 2)));
00302             //PR(sma); PR(ecc); PRL(L);
00303 
00304             add_dynamics(bi, ecc, energy);
00305             add_dynamics(bi, ecc, 2*energy);
00306             add_dynamics(bi, ecc, 4*energy);
00307             add_dynamics(bi, ecc, 8*energy);
00308         }
00309 }
00310 
00311 local void scale_limits(real& e1, real& e2, int option,
00312                         int N, real M, real E)
00313 {
00314     real scale = 1.0;
00315 
00316 PRL(E);
00317     if (option == 3) {
00318 
00319         // Limits refer to kinetic energy per unit mass.  Scale to -E/M.
00320 
00321         scale = -E / M;
00322 
00323     } else if (option != 2) {
00324 
00325         // Limits refer to energy.  Scale to kT = -(2/3) E/N.
00326 
00327         scale = -E / (1.5*N);
00328 PRL(scale);
00329     }
00330 
00331     e1 *= scale;
00332     e2 *= scale;
00333 }
00334 
00335 void main(int argc, char ** argv)
00336 {
00337     real lower = 1.0, upper = 1.0;
00338     real emax = 1;
00339     int function_select = 3;
00340     int option = 1;
00341     int random_seed = 0;
00342     char seedlog[64];
00343 
00344     check_help();
00345 
00346     extern char *poptarg;
00347     int c;
00348     char* param_string = "f:l:s:o:u:e:";
00349 
00350     while ((c = pgetopt(argc, argv, param_string)) != -1)
00351         switch(c) {
00352 
00353             case 'f': function_select = atoi(poptarg);
00354                       break;
00355             case 'l': lower = atof(poptarg);
00356                       break;
00357             case 'o': option= atoi(poptarg);
00358                       break;
00359             case 'e': emax = atof(poptarg);
00360                       break;
00361             case 's': random_seed = atoi(poptarg);
00362                       break;
00363             case 'u': upper = atof(poptarg);
00364                       break;
00365             case '?': params_to_usage(cerr, argv[0], param_string);
00366                       get_help();
00367                       exit(1);
00368         }
00369 
00370     if (lower <= 0 || upper <= 0 || lower > upper)
00371         err_exit("mkbinary: Illegal limits");
00372 
00373     dyn* b;
00374     b = get_dyn(cin);
00375 
00376     b->log_history(argc, argv);
00377 
00378     int actual_seed = srandinter(random_seed);
00379 
00380     sprintf(seedlog, "       random number generator seed = %d",
00381             actual_seed);
00382     b->log_comment(seedlog);
00383 
00384     // Look in the dyn story to see if a total system energy is known.
00385     // If it is, scale lower and upper accordingly.
00386 
00387     // Thus, in case function_select = 3, if the energy is known,
00388     // lower and upper are specified in "kT" units.  If the energy
00389     // is not known, lower and upper are absolute limits on the
00390     // binary kinetic energy.
00391 
00392     char* energy_string = "initial_total_energy";
00393 
00394     if (function_select == 1) {
00395         if (b->get_starbase()->get_stellar_evolution_scaling()) {
00396             real scale = sqrt(b->get_starbase()->conv_m_star_to_dyn(1)
00397                               * b->get_starbase()->conv_r_star_to_dyn(1));
00398             lower *= scale;
00399             upper *= scale;
00400         }
00401         else
00402             cerr << "mkbinary: Using unscaled angular-momentum limits.\n";
00403     }
00404     else if (function_select == 2) {
00405         if (b->get_starbase()->conv_r_star_to_dyn(1)>0) {
00406             lower = b->get_starbase()->conv_r_star_to_dyn(lower);
00407             upper = b->get_starbase()->conv_r_star_to_dyn(upper);
00408         }
00409         else
00410             cerr << "mkbinary: Using unscaled semi-major axis limits.\n";
00411     }
00412     else if (function_select == 3) {
00413         if (find_qmatch(b->get_dyn_story(), energy_string))
00414             scale_limits(lower, upper, option,
00415                          b->n_daughters(), b->get_mass(),
00416                          getrq(b->get_dyn_story(), energy_string));
00417         else
00418             cerr << "mkbinary: Using unscaled energy limits.\n";
00419     }
00420 
00421     mkbinary(b, lower, upper, function_select, option, emax);
00422 
00423     put_dyn(cout, *b);
00424     rmtree(b);
00425 
00426 }
00427 #endif