SeBa.C Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 
00014 
00015 
00016 
00017 
00018 
00019 
00020 
00021 
00022 
00023 
00024 
00025 
00026 
00027 
00028 
00029 
00030 
00031 
00032 
00033 
00034 
00035 
00036 
00037 
00038 
00039 
00040 
00041 
00042 
00043 
00044 
00045 
00046 
00047 
00048 
00049 
00050 
00051 
00052 
00053 
00054 
00055 
00056 
00057 
00058 
00059 
00060 
00061 
00062 
00063 
00064 
00065 
00066 
00067 
00068 
00069 
00070 
00071 
00072 
00073 
00074 
00075 
00076 
00077 
00078 
00079 
00080 
00081 
00082 
00083 
00084 //   Note:  libnode.a is referenced for the routines which produce the 
00085 //          mass function
00086 //
00087 //      version 1.0     Simon Portegies Zwart, Utrecht, 1992
00088 //                      -First version with class structure
00089 //      version 2.0     Simon Portegies Zwart, Utrecht, 1994
00090 //                      -Coupling to starlab
00091 //      version 3.0     Simon Portegies Zwart, Amsterdam, June 1997
00092 //      version 3.3     Simon Portegies Zwart, Cambridge, March 1999
00093 //
00094 
00095 #include "dyn.h" 
00096 #include "double_star.h"
00097 #include "main_sequence.h"
00098 //#include "dstar_to_dyn.h"
00099 //#include "seba.h"
00100 
00101 
00102 #ifdef TOOLBOX
00103 
00104 local bool read_binary_params(ifstream& in, real &m_prim, 
00105                               real &m_sec, real &sma, real &ecc) {
00106 
00107   if(in.eof())
00108     return false;
00109 
00110   // reading from input file
00111   int id, tpp, tps; 
00112   real time, Rlp, Rls;
00113 
00114   in >>id >> time >> sma >> ecc 
00115      >> tpp >> m_prim >> Rlp >> tps >> m_sec >> Rls;
00116 
00117       PRC(m_prim);PRC(m_sec);PRC(sma);PRL(ecc);
00118 
00119   return true;
00120 }
00121 
00122 /*-----------------------------------------------------------------------------
00123  *  binev  --
00124  *-----------------------------------------------------------------------------
00125  */
00126 local bool  evolve_binary(dyn * bi,
00127                           real start_time, real end_time,
00128                           bool stop_at_merger_or_disruption,
00129                           bool stop_at_remnant_formation) { 
00130 
00131   double_star* ds = dynamic_cast(double_star*, 
00132                                  bi->get_starbase());
00133   
00134   //            Setup star from input data.
00135   real dt, time=start_time;
00136   ds->evolve_element(time);
00137   
00138   ds->dump("SeBa.data", true);
00139 
00140   if (!bi->is_root() &&
00141       bi->get_parent()->is_root()) 
00142 
00143     do {
00144 
00145       dt = ds->get_evolve_timestep() + cnsts.safety(minimum_timestep);
00146       time = min(time+dt, end_time);
00147 
00148       ds->evolve_element(time);
00149 
00150       if (stop_at_merger_or_disruption &&
00151           (ds->get_bin_type() == Merged || 
00152            ds->get_bin_type() == Disrupted))
00153         return false;
00154 
00155       if (stop_at_remnant_formation &&
00156          (ds->get_primary()->remnant() || ds->get_secondary()->remnant()))
00157         return false;
00158 
00159     }
00160     while (time<end_time);
00161 
00162   ds->dump("SeBa.data", true);
00163     
00164   return true;
00165 
00166 }
00167 
00168 void main(int argc, char ** argv) {
00169 
00170     bool e_flag = false;
00171     bool R_flag = false;
00172     bool F_flag = false;
00173     bool I_flag = false;
00174     bool P_flag = false;
00175     bool U_flag = false;
00176     bool G_flag = false;
00177 
00178     bool stop_at_merger_or_disruption = false;
00179     bool stop_at_remnant_formation = false;
00180     bool random_initialization = false;
00181     stellar_type type = Main_Sequence;
00182     binary_type bin_type = Detached;
00183     real binary_fraction = 1.0;
00184 
00185     int n_init = 0;
00186     int n = 1;
00187 
00188     real  m_prim;
00189     real  m_sec;
00190     real  sma;
00191     real  ecc;
00192 
00193     real m_tot = 1;
00194     real r_hm = 1;
00195     real t_hc = 1;
00196 
00197     char *mfc = new char[64];
00198     mass_function mf = mf_Power_Law;
00199     real m_min = 0.1;
00200     real m_max = 100;
00201     real m_exp = -2.35;
00202     char *qfc = new char[64];
00203     mass_ratio_distribution qf = Flat_q;
00204     real q_min = 0;
00205     real q_max = 1;
00206     real q_exp = 0;
00207     char *afc = new char[64];
00208     sma_distribution af = sma_Power_Law;
00209     real a_min = 0;
00210     real a_max = 1.e+6; 
00211     real a_exp = -1;                           
00212     char *efc = new char[64];
00213     ecc_distribution ef = Thermal_Distribution;
00214     real e_min = -1;    // allow detection of constant eccentricity
00215     real e_max = 1;
00216     real e_exp;
00217 
00218     real start_time = 0;
00219     real end_time   = 35;
00220 
00221     char* input_filename;
00222 
00223     int input_seed=0; 
00224     char seedlog[64];
00225     char paramlog[90];
00226 
00227     check_help();
00228 
00229 //    seba_counters* new_seba_counters = new seba_counters;
00230 
00231     extern char *poptarg;
00232     int c;
00233     char* param_string = "n:N:RDSM:m:x:F:f:A:a:y:G:g:E:e:v:U:u:Q:q:T:I:w:P:p:n:s:";
00234 
00235     while ((c = pgetopt(argc, argv, param_string)) != -1)
00236         switch(c) {
00237             case 'R': random_initialization = true;
00238                       break;
00239             case 'D': stop_at_merger_or_disruption = true;
00240                       break;
00241             case 'S': stop_at_remnant_formation = true;
00242                       break;
00243             case 'M': m_max = atof(poptarg);
00244                       break;
00245             case 'm': m_min = atof(poptarg);
00246                       break;
00247             case 'x': m_exp = atof(poptarg);
00248                       break;
00249             case 'F': F_flag = true;
00250                       strcpy(mfc, poptarg);
00251                       break;
00252             case 'f': mf = (mass_function)atoi(poptarg);
00253                       break;
00254             case 'A': a_max = atof(poptarg);
00255                       break;
00256             case 'a': a_min = atof(poptarg);
00257                       break;
00258             case 'y': a_exp = atof(poptarg);
00259                       break;
00260             case 'G': G_flag = true;
00261                       strcpy(afc, poptarg);
00262                       break;
00263             case 'g': af = (sma_distribution)atoi(poptarg);
00264                       break;
00265             case 'E': e_max = atof(poptarg);
00266                       break;
00267             case 'e': e_min = atof(poptarg);
00268                       break;
00269             case 'v': e_exp = atof(poptarg);
00270                       break;
00271             case 'U': U_flag = true;
00272                       strcpy(efc, poptarg);
00273                       break;
00274             case 'u': ef = (ecc_distribution)atoi(poptarg);
00275                       break;
00276             case 'Q': q_max = atof(poptarg);
00277                       break;
00278             case 'q': q_min = atof(poptarg);
00279                       break;
00280             case 'T': end_time = atof(poptarg);
00281                       break;
00282             case 'I': I_flag = true;
00283                       input_filename = poptarg;
00284                       break;
00285             case 'w': q_exp = atof(poptarg);
00286                       break;
00287             case 'P': P_flag = true;
00288                       strcpy(qfc, poptarg);
00289                       break;
00290             case 'p': qf = (mass_ratio_distribution)atoi(poptarg);
00291                       break;
00292             case 'n': n = atoi(poptarg);
00293                       break;
00294             case 'N': n_init = atoi(poptarg);
00295                       break;
00296             case 's': input_seed = atoi(poptarg);
00297                       break;
00298             case '?': params_to_usage(cerr, argv[0], param_string);
00299                       exit(1);
00300         }
00301 
00302     int actual_seed = srandinter(input_seed);
00303     cerr << "random number generator seed = " << actual_seed << endl;
00304     sprintf(paramlog, 
00305 "         alpha  = %3.1f\n         lambda = %3.1f\n         beta   = %3.1f\n         gamma  = %4.2f",
00306             cnsts.parameters(common_envelope_efficiency),
00307             cnsts.parameters(envelope_binding_energy),
00308             cnsts.parameters(specific_angular_momentum_loss),
00309             cnsts.parameters(dynamic_mass_transfer_gamma));
00310 
00311     if (n <= 0) err_exit("mknodes: N > 0 required!");
00312 
00313     if(F_flag)
00314         mf = extract_mass_function_type_string(mfc);
00315     delete mfc;
00316     if(G_flag)
00317         af = extract_semimajor_axis_distribution_type_string(afc);
00318     delete afc;
00319     if(U_flag)
00320         ef = extract_eccentricity_distribution_type_string(efc);
00321     delete efc;
00322     if(P_flag)
00323         qf = extract_mass_ratio_distribution_type_string(qfc);
00324     delete qfc;
00325 
00326     actual_seed = srandinter(input_seed);
00327     sprintf(seedlog, "       random number generator seed = %d",actual_seed);
00328 
00329     ifstream infile(input_filename, ios::in);
00330     if(I_flag) {
00331       if (!infile) cerr << "error: couldn't read file "
00332                         << input_filename <<endl;
00333         cerr << "Reading input from file "<< input_filename <<endl;
00334     }
00335 
00336     dyn *root, *the_binary;
00337     double_star *ds;
00338 
00339     // Create flat tree 
00340     root = mkdyn(1);
00341     root->log_history(argc, argv);
00342     root->log_comment(seedlog);
00343     root->log_comment(paramlog);
00344     root->print_log_story(cout);
00345 
00346     print_initial_binary_distributions(m_min, m_max, mf, m_exp,
00347                                        q_min, q_max, qf, q_exp,
00348                                        a_min, a_max, af, a_exp,
00349                                        e_min, e_max, ef, e_exp);
00350 
00351     for (int i=0; i<n; i++) {
00352 
00353       if(I_flag) {
00354 
00355         if(read_binary_params(infile, m_prim, m_sec, sma, ecc)) 
00356           n=i+2;
00357         else
00358           break;
00359 
00360       }
00361       else if (random_initialization) 
00362 
00363         mkrandom_binary(m_min, m_max, mf, m_exp,
00364                         q_min, q_max, qf, q_exp,
00365                         a_min, a_max, af, a_exp,
00366                         e_min, e_max, ef, e_exp,
00367                         m_prim, m_sec, sma, ecc);
00368       else {
00369         m_prim = m_max;
00370         m_sec  = m_min;
00371         sma    = a_min;
00372         ecc    = e_min;
00373         n = 1;
00374       }
00375 
00376 //      PRC(m_prim);PRC(m_sec/m_prim);PRC(sma);PRL(ecc);
00377 
00378       // Create flat tree
00379 //      rmtree(root, false);
00380       root = mkdyn(1);
00381       root->set_mass(1);
00382       root->get_starbase()->set_stellar_evolution_scaling(m_prim, r_hm, t_hc);
00383       the_binary = root->get_oldest_daughter();
00384       add_secondary(the_binary, m_sec/m_prim);
00385 
00386       addstar(root, start_time, type);
00387 //      root->get_starbase()->set_seba_counters(new_seba_counters);
00388 
00389 //      pp(root, cerr);
00390 //      cerr << endl;
00391 
00392       ds = new_double_star(the_binary, sma, ecc, start_time, 
00393                            i + n_init, bin_type);
00394 
00395 //            PRL(((star*)the_binary->get_oldest_daughter()->get_starbase())->get_companion());
00396 
00397 //            PRL(the_binary);
00398 //      PRL(ds);
00399 
00400 //      ds->dump(cerr, false);
00401 
00402       ds->set_use_hdyn(false);
00403       ds->get_primary()->set_identity(0);
00404       ds->get_secondary()->set_identity(1);
00405 
00406       //      node *b      = root->get_oldest_daughter();
00407       //      starbase *s  = b->get_starbase();
00408       //      star *st     = dynamic_cast(star*, b->get_starbase());
00409    
00410       evolve_binary(the_binary, start_time, end_time, 
00411                     stop_at_merger_or_disruption, stop_at_remnant_formation);
00412 
00413 //      the_binary->get_starbase()->dump(cerr, false);
00414 
00415 //      ((star*)the_binary->get_starbase())->set_seba_counters(new_seba_counters);
00416 
00417 //      if(the_binary
00418 //      detach_node_from_general_tree(node & n);
00419 
00420       //      put_node(cout, *root);
00421 
00422 //      rmtree(root);
00423 
00424       //delete root->get_starbase()->get_seba_counters();
00425 
00426       //      delete ds->get_primary();
00427       //      delete ds->get_secondary();
00428       //      delete ds;
00429       //      delete root->get_oldest_daughter()->get_younger_sister();
00430         
00431     }
00432 
00433     root->log_history(argc, argv);
00434     root->log_comment(seedlog);
00435     root->log_comment(paramlog);
00436     root->print_log_story(cout);
00437 //    rmtree(root, false);
00438 }
00439 
00440 #endif