thorne_zytkow.C Source File

00001 //
00002 // thorne_zytkow.C
00003 
00004 #include "thorne_zytkow.h"
00005 #include "main_sequence.h"
00006 
00007 // ANSI C++ first creates the base class before the dreived classes are
00008 // created.
00009 
00010 thorne_zytkow::thorne_zytkow(main_sequence & m) : single_star(m) {
00011 
00012       delete &m;
00013 
00014       // (SPZ: Removed 23March2000)
00015       //      birth_mass = 0;
00016       if (is_binary_component()) {
00017 
00018         birth_mass      = get_companion()->get_total_mass();
00019         if (get_companion()->remnant()) {
00020           magnetic_field  = get_companion()->get_magnetic_field(); 
00021           rotation_period = get_companion()->get_rotation_period();
00022         }
00023       }
00024       else if (birth_mass<=0) {
00025         birth_mass    = cnsts.parameters(kanonical_neutron_star_mass);
00026       }
00027 
00028       //(SPZ:23March2000)
00029       // This should have been done in main_sequence::merge_elementes
00030 //      real m_tot    = max(birth_mass, get_total_mass());
00031 //      core_mass     = birth_mass;
00032 //      envelope_mass = m_tot - core_mass;
00033 //      core_radius   = cnsts.parameters(kanonical_neutron_star_radius);
00034 
00035       wind_constant = envelope_mass;
00036       
00037 // (GN+SPZ May  4 1999) last update age is time of previous type change
00038       last_update_age = next_update_age;
00039 
00040       // Put star at start of super_giant stage.
00041       real t_ms = main_sequence_time();
00042       real t_giant = t_ms + hertzsprung_gap_time(t_ms)
00043                           + base_giant_branch_time(t_ms);
00044       real t_he = helium_giant_time(t_ms);
00045       relative_age = t_giant + t_he;
00046       adjust_next_update_age();
00047 
00048       instantaneous_element();
00049       update();
00050 
00051       post_constructor();
00052 } 
00053 
00054 #if 0
00055 void thorne_zytkow::adjust_initial_star() {
00056 
00057   if (relative_age<=0) {
00058     real t_ms = main_sequence_time();
00059     real t_giant = t_ms + hertzsprung_gap_time(t_ms)
00060       + base_giant_branch_time(t_ms);
00061     real t_he = helium_giant_time(t_ms);
00062     relative_age = max(t_giant + t_he, relative_age);
00063   }
00064 }
00065 #endif
00066 
00067 void thorne_zytkow::instantaneous_element() {
00068 
00069   real l_g = giant_luminosity();
00070   real t_ms = main_sequence_time();
00071   real t_gs = 0.15*t_ms;
00072   real t_b  = base_giant_time(t_ms);
00073 
00074   luminosity = l_g*pow(t_gs/(next_update_age
00075                              + t_b - relative_age), 1.17);
00076   luminosity = min(luminosity, maximum_luminosity());
00077 
00078   radius = (0.25*pow(luminosity, 0.4)
00079          + 0.8*pow(luminosity, 0.67))/pow(relative_mass, 0.27);
00080   radius = min(radius, 100.);
00081 
00082   // (SPZ+GN:  1 Aug 2000)
00083   // coming from previous type the effective readius should 
00084   // remain the same.
00085   //    effective_radius = radius;
00086   
00087 }
00088 
00089 // evolve a thorne_zytkow star upto time argument according to
00090 // the super_giant model discribed by Eggleton et al. 1989.
00091 void thorne_zytkow::evolve_element(const real end_time) {
00092 
00093       real dt = end_time - current_time;
00094       current_time = end_time;
00095       relative_age += dt;
00096 
00097       if (relative_age<=next_update_age) {
00098 
00099         // (neutron_star) core accretes from envelope.
00100          accrete_from_envelope(dt);
00101          if (core_mass>cnsts.parameters(maximum_neutron_star_mass)) {
00102            if (is_binary_component()) 
00103              get_binary()->dump("binev.data", false);
00104            else
00105              dump("binev.data", false);
00106            
00107             star_transformation_story(Black_Hole);
00108             new black_hole(*this);
00109             return;
00110          }
00111 
00112          real l_g = giant_luminosity();
00113          real t_ms = main_sequence_time();
00114          real t_gs = 0.15*t_ms;
00115          real t_b  = base_giant_time(t_ms);
00116 
00117          luminosity = l_g*pow(t_gs/(next_update_age
00118                     + t_b - relative_age), 1.17);
00119          luminosity = min(luminosity, maximum_luminosity());
00120          radius = (0.25*pow(luminosity, 0.4)
00121                 + 0.8*pow(luminosity, 0.67))/pow(relative_mass, 0.27);
00122 
00123          radius = min(radius, 100.);
00124       }
00125       else {
00126 
00127         if (is_binary_component()) 
00128           get_binary()->dump("binev.data", false);
00129         else
00130           dump("binev.data", false);
00131         
00132          // thorne_zytkow lifetime exceeded.
00133          if (core_mass>cnsts.parameters(maximum_neutron_star_mass)) {
00134             star_transformation_story(Black_Hole);
00135             new black_hole(*this);
00136             return;
00137          }
00138          else {
00139             star_transformation_story(Neutron_Star);
00140             new neutron_star(*this);
00141             return;
00142          }
00143       }
00144 
00145       update();
00146       stellar_wind(dt);
00147 }
00148 
00149 void thorne_zytkow::update() {
00150 
00151       detect_spectral_features();
00152 // (GN+SPZ May  4 1999) last_update_age now used as time of last type change
00153 //  last_update_age = relative_age;
00154       effective_radius = radius;
00155 
00156    }
00157 
00158 void thorne_zytkow::accrete_from_envelope(const real dt) {
00159 
00160       real mdot_edd = 1.5e-08*cnsts.parameters(solar_radius)*core_radius*dt;
00161       if (cnsts.parameters(hyper_critical))
00162           mdot_edd *= 1.e+8; 
00163 
00164       real mdot = min(mdot_edd, envelope_mass); 
00165 
00166       core_mass += mdot;
00167       envelope_mass -= mdot;
00168 }
00169 
00170 #if 0
00171 void thorne_zytkow::stellar_wind(const real dt) {
00172 
00173       real wind_mass = 5.5e-7*dt*radius*luminosity/get_total_mass();
00174       wind_mass = min(wind_mass, envelope_mass);
00175 
00176       if (is_binary_component())
00177          get_binary()->adjust_binary_after_wind_loss(this, wind_mass, dt);
00178       else
00179          reduce_donor_mass(wind_mass);
00180    }
00181 
00182 real thorne_zytkow::helium_core_mass() {
00183 
00184       real m_core = min(cnsts.parameters(kanonical_neutron_star_mass), get_total_mass());
00185 
00186       return m_core;
00187    }
00188 #endif
00189 
00190 //              Mass transfer utilities.
00191 real thorne_zytkow::add_mass_to_accretor(const real mdot) {
00192 
00193       adjust_accretor_age(mdot);
00194 
00195       envelope_mass += mdot;
00196       if (relative_mass<get_total_mass()) {
00197         relative_mass = get_total_mass();
00198         adjust_next_update_age();
00199       }
00200 
00201       set_spec_type(Accreting);
00202 
00203       return mdot;
00204 }
00205 
00206 real thorne_zytkow::add_mass_to_accretor(real mdot, const real dt) {
00207 
00208       mdot = accretion_limit(mdot, dt);
00209 
00210       adjust_accretor_age(mdot);
00211 
00212       envelope_mass += mdot;
00213       if (relative_mass<get_total_mass()) {
00214         relative_mass = get_total_mass();
00215         adjust_next_update_age();
00216       }
00217 
00218       adjust_accretor_radius(mdot, dt);
00219 
00220       set_spec_type(Accreting);
00221       
00222       return mdot;
00223    }
00224 
00225 star* thorne_zytkow::reduce_mass(const real mdot) {
00226 
00227     if (envelope_mass - mdot <=
00228         cnsts.parameters(neutron_star_accretion_limit)*core_mass) {
00229  
00230       envelope_mass = 0;
00231 
00232       if (is_binary_component()) 
00233         get_binary()->dump("binev.data", false);
00234       else
00235         dump("binev.data", false);
00236 
00237       // Envelope mass of neutron star cannot be supported by
00238       // gravitational pressure. Envelope is blown away by 
00239       // accretion luminosity of central compact object.
00240           
00241       if (core_mass>cnsts.parameters(maximum_neutron_star_mass)) {
00242         star_transformation_story(Black_Hole);
00243         return dynamic_cast(star*, new black_hole(*this));
00244       }
00245       else {
00246         star_transformation_story(Neutron_Star);
00247         return dynamic_cast(star*, new neutron_star(*this));
00248       }
00249     }
00250 
00251     envelope_mass -= mdot;
00252     return this;
00253 }
00254 
00255 star* thorne_zytkow::subtrac_mass_from_donor(const real dt, real& mdot) {
00256 
00257       real mdot_max = relative_mass*dt/get_binary()->get_donor_timescale();
00258       mdot = mass_ratio_mdot_limit(mdot_max);
00259 
00260       if (envelope_mass<mdot) {
00261 
00262         mdot = envelope_mass;
00263         envelope_mass = 0;
00264 
00265         if (is_binary_component()) 
00266           get_binary()->dump("binev.data", false);
00267         else
00268           dump("binev.data", false);
00269 
00270         if (core_mass>cnsts.parameters(maximum_neutron_star_mass)) {
00271           star_transformation_story(Black_Hole);
00272           return dynamic_cast(star*, new black_hole(*this));
00273         }
00274         else {
00275           star_transformation_story(Neutron_Star);
00276           return dynamic_cast(star*, new neutron_star(*this));
00277         }
00278       }
00279 
00280       envelope_mass -= mdot;
00281       return this;
00282 }
00283 
00284 void thorne_zytkow::adjust_accretor_age(const real mdot,
00285                                         const bool rejuvenate) {
00286 
00287       real m_tot_new = get_total_mass() + mdot;
00288       real m_rel_new = max(m_tot_new, relative_mass);
00289 
00290       real t_ms = main_sequence_time();
00291       real t_hg = hertzsprung_gap_time(t_ms);
00292       real t_bgb = base_giant_branch_time(t_ms);
00293       real t_he = helium_giant_time(t_ms);
00294       real t_super_old = t_ms + t_hg + t_bgb + t_he;
00295       real t_nuc = nucleair_evolution_time();
00296       real t_sg_old = t_nuc - t_super_old;
00297 
00298            t_ms = main_sequence_time(m_rel_new);
00299            t_hg = hertzsprung_gap_time(m_rel_new, t_ms);
00300            t_bgb = base_giant_branch_time(m_rel_new, t_ms);
00301            t_he = helium_giant_time(m_rel_new, t_ms);
00302       real t_super_new = t_ms + t_hg + t_bgb + t_he;
00303            t_nuc = nucleair_evolution_time(m_rel_new);
00304       real t_sg_new = t_nuc - t_super_new;
00305 
00306       real dtime = relative_age - t_super_old;
00307 
00308       relative_age = t_super_new
00309                    + dtime*(t_sg_new/t_sg_old);
00310 
00311       if (rejuvenate)
00312          relative_age *= rejuvenation_fraction(mdot/m_tot_new); 
00313 
00314 
00315       // next_update_age should not be reset here
00316       // next_update_age = t_nuc;
00317    }
00318 
00319 void thorne_zytkow::adjust_next_update_age() {
00320 
00321       next_update_age = nucleair_evolution_time();
00322    }
00323 
00324 #if 0
00325 real thorne_zytkow::stellar_radius(const real mass, const real age) {
00326 
00327       real t_nuc = nucleair_evolution_time(mass);
00328 
00329       real l_g = giant_luminosity();
00330       real t_ms = main_sequence_time();
00331       real t_gs = 0.15*t_ms;
00332       real t_b  = base_giant_time(t_ms);
00333 
00334       real l_agb = l_g*pow(t_gs/(t_nuc + t_b - age), 1.17);
00335       l_agb = min(l_agb, maximum_luminosity(mass));
00336 
00337       real r_agb = (0.25*pow(l_agb, 0.4)
00338              + 0.8*pow(l_agb, 0.67))/pow(mass, 0.27);
00339 
00340       return r_agb;
00341    }
00342 #endif
00343 
00344 real thorne_zytkow::gyration_radius_sq() {
00345 
00346      return cnsts.parameters(convective_star_gyration_radius_sq); 
00347 }