00001
00002
00003
00004
00005 #include "double_support.h"
00006 #include "double_star.h"
00007
00008 binary_type extract_binary_type_string(char *type) {
00009
00010 if (!strcmp("synchronized", type)) return Synchronized;
00011 else if (!strcmp("detached", type)) return Detached;
00012 else if(!strcmp("semi_detached", type)) return Semi_Detached;
00013 else if(!strcmp("contact", type)) return Contact;
00014 else if(!strcmp("common_envelope", type)) return Common_Envelope;
00015 else if(!strcmp("merged", type)) return Merged;
00016 else if(!strcmp("disrupted", type)) return Disrupted;
00017 else if(!strcmp("spiral_in", type)) return Spiral_In;
00018
00019 return Unknown_Binary_Type;
00020
00021 }
00022
00023 char* type_string(binary_type tpe) {
00024
00025 switch(tpe) {
00026 case Synchronized: return "synchronized";
00027 case Detached: return "detached";
00028 case Semi_Detached: return "semi_detached";
00029 case Contact: return "contact";
00030 case Common_Envelope: return "common_envelope";
00031 case Merged: return "merged";
00032 case Disrupted: return "disrupted";
00033 case Spiral_In: return "spiral_in";
00034 default: return "unknown_binary_type";
00035 }
00036 }
00037
00038
00039
00040
00041
00042 void double_hist::put_double_hist() {
00043
00044 cout << binary_age<< " " << semi <<" "<< eccentricity<<endl;
00045 }
00046
00047
00048
00049
00050
00051
00052 void double_init::read_element() {
00053
00054 cin >> mass_prim
00055 >> q
00056 >> semi
00057 >> eccentricity;
00058 }
00059
00060 void double_init::put_element() {
00061
00062 ofstream outfile("binev.data", ios::app|ios::out);
00063 if(!outfile) cerr << "error: couldn't create file binev.data"<<endl;
00064
00065 outfile << mass_prim << " "
00066 << q << " "
00067 << semi << " "
00068 << eccentricity << endl;
00069 }
00070
00071 void double_init::dump(ostream & s) {
00072
00073 s << " " << mass_prim
00074 << " " << mass_prim*q
00075 << " " << semi
00076 << " " << eccentricity
00077 << endl;
00078
00079 }
00080
00081 void double_init::dump(char* filename) {
00082
00083 ofstream s(filename, ios::app|ios::out);
00084 if(!s) cerr << "error: couldn't create file "<<filename<<endl;
00085
00086 s << " " << mass_prim
00087 << " " << mass_prim*q
00088 << " " << semi
00089 << " " << eccentricity
00090 << endl;
00091 }
00092
00093
00094 double_star * get_new_binary(double_init& init, const int id) {
00095
00096
00097
00098
00099
00100
00101 return NULL;
00102 }
00103
00104 void pptime(real time,
00105 ostream & s,
00106 char *t) {
00107
00108 if (time<1.e-3)
00109 ppperiod(365.25*time*1.e+6, s, t);
00110 else {
00111 s << " " << t << " = ";
00112 if (time < 1)
00113 s << time*1000 << " kyr";
00114 else if (time < 1000)
00115 s << time << " Myr";
00116 else if (time < 1.e+6)
00117 s << time/1000 << " Gyr";
00118 else if (time < 1.e+9)
00119 s << time/1.e+6 << " Tyr";
00120 else if (time < 1.e+12)
00121 s << time/1.e+9 << " Pyr";
00122 else
00123 s << time/1.e+12 << " Eyr";
00124 }
00125 }
00126
00127 void ppperiod(real period,
00128 ostream & s,
00129 char *p) {
00130
00131 if (period>1000*365.25)
00132 pptime(period*1.e-6/365.25, s, p);
00133 else {
00134 s << " " << p << " = ";
00135 if (period < 6.94e-4)
00136 s << period*1440 << " seconds";
00137 else if (period < 0.0417)
00138 s << period*1440 << " minutes";
00139 else if (period < 1)
00140 s << period*24 << " hours";
00141 else if (period < 365.25)
00142 s << period << " days";
00143 else
00144 s << period/365.25 << " yr";
00145 }
00146
00147 }
00148
00149
00150
00151
00152
00153
00154 void put_state(double_state d,
00155 ostream & s) {
00156
00157 switch (d.type) {
00158 case Merged:
00159 s << "{";
00160 d.primary.put_star_state(s);
00161 s << "}";
00162 break;
00163 case Disrupted:
00164 s << ")";
00165 d.primary.put_star_state(s);
00166 s << ", ";
00167 d.secondary.put_star_state(s);
00168 s << "(";
00169 break;
00170 default:
00171 if (d.primary.class_spec[Rl_filling])
00172 s << "[";
00173 else
00174 s << "(";
00175 d.primary.put_star_state(s);
00176 s << ", ";
00177 d.secondary.put_star_state(s);
00178 if (d.secondary.class_spec[Rl_filling])
00179 s << "]";
00180 else
00181 s << ")";
00182 real period = 2*PI*sqrt(pow(d.semi
00183 *cnsts.parameters(solar_radius), 3.)
00184 / (cnsts.physics(G)*cnsts.parameters(solar_mass)
00185 *(d.primary.mass+d.secondary.mass)))
00186 / cnsts.physics(days);
00187
00188 ppperiod(period, s, "Porb");
00189 }
00190 s << endl;
00191 }
00192
00193
00194
00195
00196
00197 void make_profile(int id, real start_time,
00198 double_profile& binary, double_init& init) {
00199
00200 cerr<<"double_support make_profile"<<endl;
00201 real dt = (init.end_time - start_time)/init.n_steps;
00202
00203 star_state primary, secondary;
00204 double_star * b = get_new_binary(init, id);
00205
00206 primary.init_star_state((star*)b->get_primary());
00207 secondary.init_star_state((star*)b->get_secondary());
00208 binary.init_double_profile((double_star*) b, primary, secondary);
00209
00210
00211
00212 for (int j = 0; j<init.n_steps; j++)
00213 b->evolve_element(start_time+dt*(j+1));
00214
00215
00216 if (primary.identity == b->get_primary()->get_identity()) {
00217 primary.make_star_state((star*)b->get_primary());
00218 secondary.make_star_state((star*)b->get_secondary());
00219
00220
00221 binary.enhance_double_profile((double_star*) b, primary, secondary);
00222 }
00223 else {
00224 secondary.make_star_state((star*)b->get_primary());
00225 primary.make_star_state((star*)b->get_secondary());
00226
00227
00228 binary.enhance_double_profile((double_star*) b, secondary, primary);
00229 }
00230
00231 }
00232
00233 void double_profile::init_double_profile(double_star* b,
00234 star_state& p, star_state& s) {
00235
00236 init.identity = b->get_identity();
00237 init.time = b->get_current_time();
00238 init.type = b->get_bin_type();
00239 init.total_mass = b->get_total_mass();
00240 init.primary = p;
00241 init.secondary = s;
00242 }
00243
00244 void double_profile::enhance_double_profile(double_star* b,
00245 star_state& p, star_state& s) {
00246
00247 final.identity = b->get_identity();
00248 final.time = b->get_current_time();
00249 final.type = b->get_bin_type();
00250 final.total_mass = b->get_total_mass();
00251 final.primary = p;
00252 final.secondary = s;
00253
00254 mdot = init.total_mass - final.total_mass;
00255 }
00256
00257 void double_profile::init_double_profile(double_star* b) {
00258
00259 init.identity = b->get_identity();
00260 init.time = b->get_current_time();
00261 init.type = b->get_bin_type();
00262 init.semi = b->get_semi();
00263 init.ecc = b->get_eccentricity();
00264 init.velocity = b->get_velocity();
00265 init.total_mass = b->get_total_mass();
00266
00267 star_state pa, sa;
00268 if(b->get_use_hdyn()) {
00269 pa.make_star_state(dynamic_cast(star*, b->get_primary()));
00270 sa.make_star_state(dynamic_cast(star*, b->get_secondary()));
00271 }
00272 else {
00273 pa.make_star_state(dynamic_cast(star*, b->get_initial_primary()));
00274 sa.make_star_state(dynamic_cast(star*, b->get_initial_secondary()));
00275 }
00276
00277 init.primary = pa;
00278 init.secondary = sa;
00279 }
00280
00281 void double_profile::enhance_double_profile(double_star* b) {
00282
00283 final.identity = b->get_identity();
00284 final.time = b->get_current_time();
00285 final.type = b->get_bin_type();
00286 final.semi = b->get_semi();
00287 final.ecc = b->get_eccentricity();
00288 final.velocity = b->get_velocity();
00289 final.total_mass = b->get_total_mass();
00290
00291 star_state pa, sa;
00292 if(b->get_use_hdyn()) {
00293 pa.make_star_state(dynamic_cast(star*, b->get_primary()));
00294 sa.make_star_state(dynamic_cast(star*, b->get_secondary()));
00295 }
00296 else {
00297 pa.make_star_state(dynamic_cast(star*, b->get_initial_primary()));
00298 sa.make_star_state(dynamic_cast(star*, b->get_initial_secondary()));
00299 }
00300
00301 final.primary = pa;
00302 final.secondary = sa;
00303 }
00304
00305 void put_profile(double_profile& d) {
00306
00307 cerr <<"Initial: " << endl;
00308 put_state(d.init, cerr);
00309 cerr << "final: " << endl;
00310 put_state(d.final, cerr);
00311
00312 }
00313
00314 #if 0
00315 double_star * triple_star(double_init& inner,
00316 double_init& outer, int id) {
00317
00318 single_star * p = new single_star();
00319 p->initialize(inner.mass_prim, id, inner.start_time);
00320 single_star * s = new single_star();
00321 s->initialize(inner.mass_prim*inner.q, id+1, inner.start_time);
00322
00323
00324 double_star * b = new double_star(p, s, inner);
00325
00326 single_star * t = new single_star();
00327 t->initialize(outer.mass_prim*outer.q, id+2, outer.start_time);
00328
00329
00330 double_star * triple = new double_star(b, t, outer);
00331
00332 return triple;
00333 }
00334 #endif
00335
00336 double_state make_state(double_star* b) {
00337
00338 double_state dbl;
00339 dbl.identity = b->get_identity();
00340 dbl.time = b->get_current_time();
00341
00342 dbl.type = b->obtain_binary_type();
00343 dbl.semi = b->get_semi();
00344 dbl.ecc = b->get_eccentricity();
00345 dbl.velocity = b->get_velocity();
00346 dbl.total_mass = b->get_total_mass();
00347
00348 star_state pa, sa;
00349 if(b->get_use_hdyn()) {
00350 pa.make_star_state(dynamic_cast(star*, b->get_primary()));
00351 sa.make_star_state(dynamic_cast(star*, b->get_secondary()));
00352 }
00353 else {
00354 pa.make_star_state(dynamic_cast(star*, b->get_initial_primary()));
00355 sa.make_star_state(dynamic_cast(star*, b->get_initial_secondary()));
00356 }
00357
00358 dbl.primary = pa;
00359 dbl.secondary = sa;
00360
00361 return dbl;
00362 }
00363