triple_star.C Source File

00001 //
00002 // triple_star.C
00003 //
00004 
00005 #include "double_star.h"
00006 #include "single_star.h"
00007 #include "main_sequence.h"
00008 
00009 
00010 //              Mass transfer utilities.
00011 /*
00012 real get_relative_mass() {
00013       return get_primary()->get_relative_mass() + get_secondary()->get_relative_mass();
00014    }
00015 */
00016 
00017 star* double_star::reduce_mass(const real mdot) {
00018 
00019 cerr <<"void double_star::reduce_mass(mdot="<<mdot<<")"<<endl;
00020 
00021         real envelope_mass = get_primary()->get_envelope_mass()
00022                           + get_secondary()->get_envelope_mass();
00023 
00024         if (envelope_mass<=0) 
00025 cerr <<"Really in trouble! envelope_mass = "<<envelope_mass<<endl;
00026 
00027 /*
00028         if (envelope_mass<=mdot) {
00029            get_primary()->reduce_mass(primary->get_envelope_mass());
00030            secondary->reduce_mass(secondary->get_envelope_mass());
00031            return;
00032         }
00033 */
00034 cerr<<"reduce_envelope of both components with the fraction: "
00035     <<mdot/envelope_mass<<endl;
00036 
00037 
00038       return this;
00039      }
00040 
00041 star* double_star::subtrac_mass_from_donor(const real dt, real& mdot) {
00042 
00043 cerr <<"real double_star::subtrac_mass_from_donor(dt= "<<dt<<")"<<endl;
00044 cerr<<"double_star not allowed to fill Roche-lobe"<<endl;
00045 exit(1);
00046 
00047 /*
00048         real mdot = relative_mass*dt/root->get_donor_timescale();
00049 //        mdot = mass_ratio_mdot_limit(mdot);
00050 
00051         real evelope_mass = primary->get_envelope_mass()
00052                           + secondary->get_envelope_mass();
00053 
00054 */
00055         //return 0;
00056         mdot = 0;
00057  
00058       return this;
00059      }
00060 
00061 real double_star::wind_velocity() {
00062     cerr<<"double_star::wind_velocity()"<<endl;
00063     real v_wind = 2.5*sqrt(cnsts.physics(G)*cnsts.parameters(solar_mass)*get_total_mass()
00064                            / (get_effective_radius()*cnsts.parameters(solar_radius)))
00065                   / cnsts.physics(km_per_s);
00066     return v_wind;
00067 
00068 }
00069 void double_star::stellar_wind(const real dt) {
00070 cerr<<"void double_star::stellar_wind(const real "<<dt<<")"<<endl;
00071       perform_wind_loss(dt);
00072 }
00073 
00074 
00075 real double_star::temperature() {
00076 //cerr<<"real double_star::temperature()"<<endl;
00077 
00078       real temp = get_primary()->temperature()
00079                 + get_secondary()->temperature();
00080 
00081       return temp;
00082    }
00083 
00084 real double_star::bolometric_correction() {
00085 cerr<<"real double_star::bolometric_correction()"<<endl;
00086 
00087       real bc = get_primary()->bolometric_correction()
00088               + get_secondary()->bolometric_correction();
00089 
00090       return bc;
00091    } 
00092 
00093 real double_star::mass_ratio_mdot_limit(real mdot) {
00094 cerr<<"real double_star::mass_ratio_mdot_limit(mdot="<<mdot<<")"<<endl;
00095 
00096         real accretor_mass = get_companion()->get_total_mass();
00097            if (accretor_mass<get_total_mass()) {
00098               real mdot_max = get_total_mass() - accretor_mass;
00099               if (mdot>mdot_max) mdot = mdot_max;
00100            }
00101 
00102         return mdot;
00103      }
00104 
00105 real double_star::kelvin_helmholds_timescale() {
00106 cerr<<"real double_star::kelvin_helmholds_timescale()"<<endl;
00107 
00108         return 1;
00109      }
00110 
00111 real double_star::nucleair_evolution_timescale() {
00112 cerr<<"real double_star::nucleair_evolution_timescale()"<<endl;
00113 
00114 
00115         return cnsts.physics(nucleair_efficiency);
00116      }
00117 
00118 real double_star::dynamic_timescale() {
00119 cerr<<"real double_star::dynamic_timescale()"<<endl;
00120 
00121         return 1;
00122      }
00123 
00124 
00125 real double_star::accretion_limit(const real mdot, const real dt) {
00126 cerr<<"real double_star::accretion_limit(mdot="<<mdot<<", dt="<<dt<<")"<<endl;
00127 
00128         real r_rl = get_binary()->roche_radius(this);
00129 cerr <<"rl "<<r_rl ;
00130         real prim_mdot = get_primary()->accretion_limit(mdot, dt);
00131         real sec_mdot  = get_secondary()->accretion_limit(mdot, dt);
00132 cerr<<"mdot_p="<<prim_mdot<<" mdot_s="<<sec_mdot<<endl;
00133         real accretion = prim_mdot + sec_mdot;
00134 
00135         return accretion;
00136 
00137 }
00138 /*
00139 //              Should be inplemented in base_element.h
00140 void double_star::adjust_accretor_radius(const real mdot, const real dt) {
00141 cerr<<"d double_star::adjust_accretor_radius(mdot="<<mdot
00142     <<", dt="<<dt<<")"<<endl;
00143 
00144       primary->adjust_accretor_radius(mdot, dt);
00145       secondary->adjust_accretor_radius(mdot, dt);
00146    }
00147 */
00148 
00149 real double_star::mass_transfer_timescale(mass_transfer_type &type) {
00150 cerr<<"real double_star::mass_transfer_timescale()"<<endl;
00151 
00152       real mtt_s = get_secondary()->mass_transfer_timescale(type);
00153       real mtt_p = get_primary()->mass_transfer_timescale(type);
00154 
00155       real mtt = 1./(1./mtt_p + 1./mtt_s);
00156       
00157       return mtt;
00158 }
00159 
00160 real double_star::zeta_adiabatic() {
00161     cerr<<"real double_star::zeta_adiabatic()"<<endl;
00162     return 0;
00163 }
00164 
00165 real double_star::zeta_thermal() {
00166     cerr<<"real double_star::zeta_thermal()"<<endl;
00167     return 0;
00168 }
00169 
00170 real double_star::add_mass_to_accretor(const real mdot) {
00171     cerr<<"real double_star::add_mass_to_accretor(mdot="<<mdot<<")"<<endl;
00172 
00173       if (bin_type==Merged || bin_type==Disrupted) 
00174          return 0;
00175 
00176       if (mdot<=0) {
00177          cerr << "double_star::add_mass_to_accretor(mdot=" << mdot << ")"<<endl;
00178          cerr << "mdot (" << mdot << ") smaller than zero!" << endl;
00179          cerr << "Action: No action!" << endl;
00180          return 0;
00181       }
00182       
00183       mass_transfer_type type = Unknown;
00184       real mtt_p = get_primary()->mass_transfer_timescale(type);
00185       real mtt_s = get_secondary()->mass_transfer_timescale(type);
00186       real mtt_b = mass_transfer_timescale(type);
00187 
00188       real mdot_p = mdot*mtt_b/mtt_p;
00189       real mdot_s = mdot*mtt_b/mtt_s;
00190 
00191       real M_old = get_total_mass();
00192       mdot_p = get_primary()->add_mass_to_accretor(mdot_p);
00193       mdot_s = get_secondary()->add_mass_to_accretor(mdot_s);
00194      
00195       real mdot_left = mdot - (mdot_p + mdot_s);
00196 
00197 /*
00198 //              Spiral in anyway (temporary treatment)
00199        real r_lobe = min(primary->get_effective_radius()/semi,
00200                          secondary->get_effective_radius()/semi);
00201        real a_spi = semi*(primary->get_total_mass()
00202                           /(primary->get_total_mass()+mdot))
00203                           /(1. + (2.*mdot_left
00204                    /      (cnsts.parameters(common_envelope_efficiency)
00205                    *       cnsts.parameters(envelope_binding_energy)
00206                    *       r_lobe*secondary->get_total_mass())));
00207        
00208 */
00209       //real beta = 3;
00210        real M_new = get_total_mass();
00211        semi *= pow(M_old/M_new,
00212                    2*cnsts.parameters(specific_angular_momentum_loss) + 1);
00213 
00214 //       if (a_spi <= get_primary()->get_effective_radius() 
00215 //                 + secondary->get_effective_radius()) 
00216        if(!stable() || semi<get_primary()->get_radius() + get_secondary()->get_radius())
00217           merge_elements(get_primary(), get_secondary());
00218        //else
00219        //   semi = a_spi;
00220 
00221       return mdot_p + mdot_s;
00222       
00223    }
00224 
00225 real double_star::add_mass_to_accretor(real mdot, const real dt) {
00226 cerr<<"real double_star::add_mass_to_accretor(mdot="<<mdot
00227     <<", dt="<<dt<<")"<<endl;
00228 
00229       if (bin_type==Merged || bin_type==Disrupted) 
00230          return 0;
00231 
00232       if (mdot<=0) {
00233          cerr << "double_star::add_mass_to_accretor(mdot=" << mdot << ")"<<endl;
00234          cerr << "mdot (mdot=" << mdot << ", dt="<<dt
00235               <<"): mdot smaller than zero!" << endl;
00236          cerr << "Action: No action!" << endl;
00237          return 0;
00238       }
00239 
00240       mass_transfer_type type = Unknown;
00241       real mtt_p = get_primary()->mass_transfer_timescale(type);
00242       real mtt_s = get_secondary()->mass_transfer_timescale(type);
00243       real mtt_b = mass_transfer_timescale(type);
00244 
00245       real mdot_p = mdot*mtt_b/mtt_p;
00246       real mdot_s = mdot*mtt_b/mtt_s;
00247 
00248       mdot_p = get_primary()->add_mass_to_accretor(mdot_p, dt);
00249       mdot_s = get_secondary()->add_mass_to_accretor(mdot_s, dt);
00250 
00251       real mdot_left = mdot - (mdot_p + mdot_s);
00252 
00253 //              Spiral in anyway (temporary treatment)
00254        real r_lobe = min(get_primary()->get_effective_radius()/semi,
00255                          get_secondary()->get_effective_radius()/semi);
00256        real a_spi = semi*(get_primary()->get_total_mass()
00257                           /(get_primary()->get_total_mass()+mdot))
00258                           /(1. + (2.*mdot_left
00259                    /      (cnsts.parameters(common_envelope_efficiency)
00260                    *       cnsts.parameters(envelope_binding_energy)
00261                    *       r_lobe*get_secondary()->get_total_mass())));
00262 
00263        if (a_spi <= get_primary()->get_effective_radius()
00264                  + get_secondary()->get_effective_radius())
00265           merge_elements(get_primary(), get_secondary());
00266        else
00267           semi = a_spi;
00268 
00269       return mdot_p + mdot_s;
00270 
00271 }
00272 real double_star::accrete_from_stellar_wind(const real mdot, const real dt) {
00273 
00274     // cerr<<"real double_star::accrete_from_stellar_wind(mdot="<<mdot
00275     //     <<", dt="<<dt<<")"<<endl;
00276 
00277     //          Just try to prevent error in star::accrete_from_stellar_wind()
00278 
00279       if(bin_type!=Merged && bin_type!=Disrupted) {
00280       real mdot_p = get_primary()->accrete_from_stellar_wind(mdot, dt);
00281       real mdot_s = mdot - mdot_p;
00282       if (mdot_s>0) 
00283          get_secondary()->accrete_from_stellar_wind(mdot_s, dt);
00284       }
00285 
00286       return 0;         // Added by Steve 8/4/97 to satisfy DEC C++
00287   }
00288 
00289 void double_star::set_rotational_velocity(const real v) {
00290     cerr<<"void double_star::set_rotational_velocity(v="<<v<<")"<<endl;
00291 }
00292 
00293 void double_star::set_previous_radius(const real r) {
00294     //cerr<<"void double_star::set_previous_radius(r="<<r<<")"<<endl;
00295 }
00296 
00297 void double_star::adjust_triple_after_wind_loss(star * donor,
00298                                                 const real md_dot,
00299                                                 const real dt) {
00300 
00301 cerr <<"void double_star::adjust_triple_after_wind_loss(d="
00302      <<donor<<", m.="<<md_dot<<" dt="<<dt<<"): "<<semi;
00303      if (bin_type!=Merged && bin_type!=Disrupted) {
00304         star* accretor = donor->get_companion();
00305 
00306      real M_old = get_total_mass() + md_dot;
00307      real old_donor_mass = donor->get_total_mass() + md_dot;
00308      real old_accretor_mass = accretor->get_total_mass();
00309 cerr<<"m's: "<<M_old<<" "<<old_donor_mass<<" "<<old_accretor_mass
00310              <<" "<<md_dot<<endl;
00311 
00312      real ma_dot = accretor->accrete_from_stellar_wind(md_dot, dt);
00313 
00314      real M_new = get_total_mass();
00315      real new_donor_mass = donor->get_total_mass();
00316      real new_accretor_mass = accretor->get_total_mass();
00317 cerr<<"M's: "<<M_new<<" "<<new_donor_mass<<" "<<new_accretor_mass
00318              <<" "<<md_dot<<endl;
00319 
00320 cerr<<"mdot: "<<md_dot<<" "<<ma_dot<<endl;
00321      if(md_dot>0 && ma_dot>=0) {
00322         real alpha = 1 - ma_dot/md_dot;
00323         semi *= pow(pow(new_donor_mass/old_donor_mass, 1-alpha)
00324              *  new_accretor_mass/old_accretor_mass, -2)
00325              *  M_old/M_new;
00326      }
00327      }
00328      else  {
00329 //        base_element * old_d = donor;
00330         donor->reduce_mass(md_dot);
00331 //        if (old_d==get_primary()) donor=primary;
00332 //        else donor=secondary;
00333      }
00334 cerr<<" a-->"<<semi<<endl;
00335 }
00336 
00337 
00338 // What is the gyration radius of a binary?
00339 real double_star::gyration_radius_sq() {
00340   cerr << "double_star::gyration_radius_sq()"<<endl;
00341   cerr << "Unknown: return 0"<<endl; 
00342 
00343   return 0; // Dummy
00344 }