00001
00002
00003
00004
00005
00006
00007 #include "node.h"
00008 #include "double_star.h"
00009 #include "main_sequence.h"
00010
00011
00012 #ifdef TOOLBOX
00013
00014
00015
00016
00017
00018 local bool evolve_binary(node * bi,
00019 real start_time, real end_time) {
00020
00021 double_star* ds = dynamic_cast(double_star*,
00022 bi->get_starbase());
00023
00024
00025 real dt, time=start_time;
00026 ds->evolve_element(time);
00027
00028 ds->dump("init.dat", true);
00029
00030 if (!bi->is_root() &&
00031 bi->get_parent()->is_root())
00032
00033 do {
00034
00035
00036 dt = 100;
00037 time = min(time+dt, end_time);
00038
00039 PRC(bi->get_starbase()->get_evolve_timestep());PRC(time);
00040 PRC(end_time);PRL(ds->get_bin_type());
00041 PRL(((double_star*)bi->get_starbase())->get_donor_timescale());
00042
00043 ds->evolve_element(time);
00044
00045 if (ds->get_bin_type() == Merged ||
00046 ds->get_bin_type() == Disrupted)
00047 return false;
00048
00049 if (ds->get_primary()->remnant() || ds->get_secondary()->remnant())
00050 return false;
00051
00052 }
00053 while (time<end_time);
00054
00055 return true;
00056
00057 }
00058
00059 void main(int argc, char ** argv) {
00060
00061 bool e_flag = false;
00062 bool R_flag = true;
00063
00064 bool reandom_initialization = false;
00065 int n = 1;
00066 int n_init = 0;
00067 int id=1;
00068 stellar_type type = Main_Sequence;
00069 binary_type bin_type = Detached;
00070 real binary_fraction = 1.0;
00071
00072 real m_tot = 1;
00073 real r_hm = 1;
00074 real t_hc = 1;
00075
00076 real m_min = 0.5;
00077 real m_max = 100;
00078 real lm_min = log10(m_min);
00079 real lm_max = log10(m_max);
00080 real a_min = 10;
00081 real a_max = 1.e+4;
00082 real la_min = log10(a_min);
00083 real la_max = log10(a_max);
00084 real q_min = 0;
00085 real q_max = 1;
00086 real sma = 138;
00087 real ecc = 0.0;
00088 real m_prim = 13.1;
00089 real m_sec = 9.8;
00090
00091 real start_time = 0;
00092 real end_time = 35;
00093
00094 int input_seed=0, actual_seed;
00095 char seedlog[64];
00096 extern char *poptarg;
00097 int c;
00098 char* param_string = "A:a:e:M:m:N:n:Q;q:Rs:T:";
00099
00100 while ((c = pgetopt(argc, argv, param_string)) != -1)
00101 switch(c) {
00102 case 'A': a_max = atof(poptarg);
00103 la_max = log10(a_max);
00104 break;
00105 case 'a': a_min = atof(poptarg);
00106 la_min = log10(a_min);
00107 break;
00108 case 'e': e_flag = true;
00109 break;
00110 case 'M': m_max = atof(poptarg);
00111 lm_max = log10(m_max);
00112 break;
00113 case 'm': m_min = atof(poptarg);
00114 lm_min = log10(m_min);
00115 break;
00116 case 'n': n = atoi(poptarg);
00117 break;
00118 case 'N': n_init = atoi(poptarg);
00119 break;
00120 case 'Q': q_max = atof(poptarg);
00121 break;
00122 case 'q': q_min = atof(poptarg);
00123 break;
00124 case 'R': R_flag = false;
00125 break;
00126 case 'T': end_time = atof(poptarg);
00127 break;
00128 case 's': input_seed = atoi(poptarg);
00129 break;
00130 case '?': params_to_usage(cerr, argv[0], param_string);
00131 exit(1);
00132 }
00133
00134 if (n <= 0) err_exit("mknodes: N > 0 required!");
00135
00136 actual_seed = srandinter(input_seed);
00137 sprintf(seedlog, " random number generator seed = %d",actual_seed);
00138
00139
00140 node *root = mknode(n);
00141 root->log_history(argc, argv);
00142 root->log_comment(seedlog);
00143 root->get_starbase()->set_stellar_evolution_scaling(m_tot, r_hm, t_hc);
00144
00145 node* the_binary = root->get_oldest_daughter();
00146
00147 if (!R_flag) {
00148
00149 m_prim = m_max;
00150 m_sec = m_min;
00151 sma = a_min;
00152 ecc = 0;
00153
00154 the_binary->set_mass(m_prim);
00155
00156
00157 add_secondary(the_binary, m_sec);
00158 addstar(root, start_time, type);
00159 double_star* ds = new_double_star(the_binary, sma, ecc,
00160 start_time, n_init, bin_type);
00161 ds->set_use_hdyn(false);
00162 ds->get_primary()->set_identity(0);
00163 ds->get_secondary()->set_identity(1);
00164
00165 node *b = root->get_oldest_daughter();
00166 starbase *s = b->get_starbase();
00167 star *st = dynamic_cast(star*, b->get_starbase());
00168
00169 evolve_binary(the_binary, 0, end_time);
00170
00171 }
00172 else
00173 for (int i=0; i<n; i++) {
00174
00175
00176
00177
00178 mkrandom_binary(lm_min, lm_max,
00179 2.35,
00180 la_min, la_max,
00181 false,
00182 m_prim, m_sec, sma, ecc);
00183
00184 the_binary->set_mass(m_prim);
00185
00186
00187 add_secondary(the_binary, m_sec);
00188 addstar(dynamic_cast(node*, root), start_time, type);
00189 double_star* ds = new_double_star(the_binary, sma, ecc,
00190 start_time, i + n_init, bin_type);
00191
00192 ds->set_use_hdyn(false);
00193 ds->get_primary()->set_identity(0);
00194 ds->get_secondary()->set_identity(1);
00195
00196 node *b = root->get_oldest_daughter();
00197 starbase *s = b->get_starbase();
00198 star *st = dynamic_cast(star*, b->get_starbase());
00199
00200 evolve_binary(the_binary, 0, end_time);
00201
00202 delete ds->get_primary();
00203 delete ds->get_secondary();
00204 delete ds;
00205
00206 }
00207 }
00208
00209 #endif