nSeBa.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 "node.h" 
00096 #include "double_star.h"
00097 #include "main_sequence.h"
00098 //#include "dstar_to_dyn.h"
00099 
00100 #ifdef TOOLBOX
00101 
00102 local bool read_binary_params(ifstream& in, real &m_prim, 
00103                               real &m_sec, real &sma, real &ecc) {
00104 
00105   if(in.eof())
00106     return false;
00107 
00108   // reading from input file
00109   int id, time, tpp, tps; 
00110   real Rlp, Rls;
00111 
00112   in >>id >> time >> sma >> ecc 
00113      >> tpp >> m_prim >> Rlp >> tps >> m_sec >> Rls;
00114 
00115       PRC(m_prim);PRC(m_sec);PRC(sma);PRL(ecc);
00116 
00117   return true;
00118 }
00119 
00120 /*-----------------------------------------------------------------------------
00121  *  binev  --
00122  *-----------------------------------------------------------------------------
00123  */
00124 local bool  evolve_binary(node * bi,
00125                           real start_time, real end_time,
00126                           bool stop_at_merger_or_disruption,
00127                           bool stop_at_remnant_formation) { 
00128 
00129   double_star* ds = dynamic_cast(double_star*, 
00130                                  bi->get_starbase());
00131   
00132   //            Setup star from input data.
00133   real dt, time=start_time;
00134   ds->evolve_element(time);
00135   
00136   ds->dump("SeBa.data", true);
00137 
00138   if (!bi->is_root() &&
00139       bi->get_parent()->is_root()) 
00140 
00141     do {
00142 
00143       dt = ds->get_evolve_timestep() + cnsts.safety(minimum_timestep);
00144       time = min(time+dt, end_time);
00145 
00146       ds->evolve_element(time);
00147 
00148       if (stop_at_merger_or_disruption &&
00149           (ds->get_bin_type() == Merged || 
00150            ds->get_bin_type() == Disrupted))
00151         return false;
00152 
00153       if (stop_at_remnant_formation &&
00154          (ds->get_primary()->remnant() || ds->get_secondary()->remnant()))
00155         return false;
00156 
00157     }
00158     while (time<end_time);
00159 
00160   ds->dump("SeBa.data", true);
00161     
00162   return true;
00163 
00164 }
00165 
00166 void main(int argc, char ** argv) {
00167 
00168     bool e_flag = false;
00169     bool R_flag = false;
00170     bool F_flag = false;
00171     bool I_flag = false;
00172     bool P_flag = false;
00173     bool U_flag = false;
00174     bool G_flag = false;
00175 
00176     bool stop_at_merger_or_disruption = false;
00177     bool stop_at_remnant_formation = false;
00178     bool random_initialization = false;
00179     stellar_type type = Main_Sequence;
00180     binary_type bin_type = Detached;
00181     real binary_fraction = 1.0;
00182 
00183     int n_init = 0;
00184     int n = -1;
00185 
00186     real  m_prim;
00187     real  m_sec;
00188     real  sma;
00189     real  ecc;
00190 
00191     real m_tot = 1;
00192     real r_hm = 1;
00193     real t_hc = 1;
00194 
00195     char *mfc = new char[64];
00196     mass_function mf = mf_Power_Law;
00197     real m_min = 0.1;
00198     real m_max = 100;
00199     real m_exp = -2.35;
00200     char *qfc = new char[64];
00201     mass_ratio_distribution qf = Flat_q;
00202     real q_min = 0;
00203     real q_max = 1;
00204     real q_exp = 0;
00205     char *afc = new char[64];
00206     sma_distribution af = sma_Power_Law;
00207     real a_min = 0;
00208     real a_max = 1.e+6; 
00209     real a_exp = -1;                           
00210     char *efc = new char[64];
00211     ecc_distribution ef = Thermal_Distribution;
00212     real e_min = -1;    // allow detection of constant eccentricity
00213     real e_max = 1;
00214     real e_exp;
00215 
00216     real start_time = 0;
00217     real end_time   = 35;
00218 
00219     char* input_filename;
00220 
00221     int input_seed=0; 
00222     char seedlog[64];
00223     char paramlog[65];
00224 
00225     check_help();
00226 
00227     extern char *poptarg;
00228     int c;
00229     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:";
00230 
00231     while ((c = pgetopt(argc, argv, param_string)) != -1)
00232         switch(c) {
00233             case 'R': random_initialization = true;
00234                       break;
00235             case 'D': stop_at_merger_or_disruption = true;
00236                       break;
00237             case 'S': stop_at_remnant_formation = true;
00238                       break;
00239             case 'M': m_max = atof(poptarg);
00240                       break;
00241             case 'm': m_min = atof(poptarg);
00242                       break;
00243             case 'x': m_exp = atof(poptarg);
00244                       break;
00245             case 'F': F_flag = true;
00246                       strcpy(mfc, poptarg);
00247                       break;
00248             case 'f': mf = (mass_function)atoi(poptarg);
00249                       break;
00250             case 'A': a_max = atof(poptarg);
00251                       break;
00252             case 'a': a_min = atof(poptarg);
00253                       break;
00254             case 'y': a_exp = atof(poptarg);
00255                       break;
00256             case 'G': G_flag = true;
00257                       strcpy(afc, poptarg);
00258                       break;
00259             case 'g': af = (sma_distribution)atoi(poptarg);
00260                       break;
00261             case 'E': e_max = atof(poptarg);
00262                       break;
00263             case 'e': e_min = atof(poptarg);
00264                       break;
00265             case 'v': e_exp = atof(poptarg);
00266                       break;
00267             case 'U': U_flag = true;
00268                       strcpy(efc, poptarg);
00269                       break;
00270             case 'u': ef = (ecc_distribution)atoi(poptarg);
00271                       break;
00272             case 'Q': q_max = atof(poptarg);
00273                       break;
00274             case 'q': q_min = atof(poptarg);
00275                       break;
00276             case 'T': end_time = atof(poptarg);
00277                       break;
00278             case 'I': I_flag = true;
00279                       strcpy(input_filename, poptarg);
00280                       break;
00281             case 'w': q_exp = atof(poptarg);
00282                       break;
00283             case 'P': P_flag = true;
00284                       strcpy(qfc, poptarg);
00285                       break;
00286             case 'p': qf = (mass_ratio_distribution)atoi(poptarg);
00287                       break;
00288             case 'n': n = atoi(poptarg);
00289                       break;
00290             case 'N': n_init = atoi(poptarg);
00291                       break;
00292             case 's': input_seed = atoi(poptarg);
00293                       break;
00294             case '?': params_to_usage(cerr, argv[0], param_string);
00295                       exit(1);
00296         }
00297 
00298     int actual_seed = srandinter(input_seed);
00299     cerr << "random number generator seed = " << actual_seed << endl;
00300     sprintf(paramlog, 
00301    "         alpha  = %3.1f\n         lambda = %3.1f\n         beta   = %3.1f",
00302             cnsts.parameters(common_envelope_efficiency),
00303             cnsts.parameters(envelope_binding_energy),
00304             cnsts.parameters(specific_angular_momentum_loss));
00305 
00306 
00307     if (n <= 0) {
00308       hdyn *b;
00309 
00310       b = get_hdyn(cin);
00311       if (c_flag == TRUE)
00312         b->log_comment(comment);
00313       b->log_history(argc, argv);
00314       b->get_starbase()->set_stellar_evolution_scaling(m_tot, r_hm, t_hc);
00315 
00316       addstar(root, start_time, type);
00317 
00318       pp(root, cerr);
00319       cerr << endl;
00320 
00321       ds = new_double_star(the_binary, sma, ecc, start_time, 
00322                            i + n_init, bin_type);
00323 
00324     put_node(cout, *b);
00325 
00326     real t_out, dt = t_end/n;
00327     for (int i=0; i<n; i++) {
00328       t_out = dt * (i+1);
00329       
00330       for_all_daughters(node, b, bi) 
00331         bi->get_starbase()->evolve_element(t_out);
00332 
00333       put_node(cout, *b);
00334     }
00335 
00336 
00337     }
00338     else {
00339 
00340     if(F_flag)
00341         mf = extract_mass_function_type_string(mfc);
00342     delete mfc;
00343     if(G_flag)
00344         af = extract_semimajor_axis_distribution_type_string(afc);
00345     delete afc;
00346     if(U_flag)
00347         ef = extract_eccentricity_distribution_type_string(efc);
00348     delete efc;
00349     if(P_flag)
00350         qf = extract_mass_ratio_distribution_type_string(qfc);
00351     delete qfc;
00352 
00353     actual_seed = srandinter(input_seed);
00354     sprintf(seedlog, "       random number generator seed = %d",actual_seed);
00355 
00356     ifstream infile(input_filename);
00357     if(I_flag) {
00358       if (!infile) cerr << "error: couldn't create file "
00359                         << input_filename <<endl;
00360         cerr << "Reading input from file "<< input_filename <<endl;
00361     }
00362 
00363     node *root, *the_binary;
00364     double_star *ds;
00365 
00366     for (int i=0; i<n; i++) {
00367 
00368       if(I_flag) {
00369 
00370         if(read_binary_params(infile, m_prim, m_sec, sma, ecc)) 
00371           n=i+2;
00372         else
00373           break;
00374 
00375       }
00376       else if (random_initialization) 
00377         mkrandom_binary(m_min, m_max, mf, m_exp,
00378                         q_min, q_max, qf, q_exp,
00379                         a_min, a_max, af, a_exp,
00380                         e_min, e_max, ef, e_exp,
00381                         m_prim, m_sec, sma, ecc);
00382       else {
00383         m_prim = m_max;
00384         m_sec  = m_min;
00385         sma    = a_min;
00386         ecc    = e_min;
00387         n = 1;
00388       }
00389 
00390       //      PRC(m_prim);PRC(m_sec/m_prim);PRC(sma);PRL(ecc);
00391 
00392       // Create flat tree 
00393       root = mknode(1);
00394 
00395       root->log_history(argc, argv);
00396       root->log_comment(seedlog);
00397       root->log_comment(paramlog);
00398       root->print_log_story(cout);
00399       root->set_mass(1);
00400       root->get_starbase()->set_stellar_evolution_scaling(m_prim, r_hm, t_hc);
00401       
00402       the_binary = root->get_oldest_daughter();
00403       add_secondary(the_binary, m_sec/m_prim);
00404       
00405       addstar(root, start_time, type);
00406 
00407       pp(root, cerr);
00408       cerr << endl;
00409 
00410       ds = new_double_star(the_binary, sma, ecc, start_time, 
00411                            i + n_init, bin_type);
00412 
00413       //      PRL(((star*)the_binary->get_oldest_daughter()->get_starbase())->get_companion());
00414 
00415       //      PRL(the_binary);
00416       //      PRL(ds);
00417 
00418       ds->dump(cerr, false);
00419       ds->set_use_hdyn(false);
00420       ds->get_primary()->set_identity(0);
00421       ds->get_secondary()->set_identity(1);
00422 
00423       //      node *b      = root->get_oldest_daughter();
00424       //      starbase *s  = b->get_starbase();
00425       //      star *st     = dynamic_cast(star*, b->get_starbase());
00426    
00427       evolve_binary(the_binary, start_time, end_time, 
00428                     stop_at_merger_or_disruption, stop_at_remnant_formation);
00429 
00430       the_binary->get_starbase()->dump(cerr, false);
00431 
00432 //      if(the_binary
00433 //      detach_node_from_general_tree(node & n);
00434 
00435       //      put_node(cout, *root);
00436 
00437       rmtree(root);
00438       //      delete ds->get_primary();
00439       //      delete ds->get_secondary();
00440       //      delete ds;
00441       //      delete root->get_oldest_daughter()->get_younger_sister();
00442         
00443     }
00444   }
00445 }
00446 
00447 #endif