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