Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 Pablo Marchant & The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 :
21 : module binary_mdot
22 :
23 : use const_def, only: dp, pi, pi4, clight, one_third, kerg, mp
24 : use math_lib
25 : use star_lib
26 : use star_def
27 : use binary_def
28 : use binary_wind
29 : use binary_ce
30 : use utils_lib, only: mesa_error
31 :
32 : implicit none
33 :
34 : contains
35 :
36 0 : integer function check_implicit_rlo(binary_id, new_mdot)
37 : integer, intent(in) :: binary_id
38 : real(dp), intent(out) :: new_mdot
39 :
40 : type (binary_info), pointer :: b
41 : type (star_info), pointer :: s
42 0 : real(dp) :: function_to_solve, explicit_mdot, q, slope_contact
43 : integer :: ierr
44 : logical :: use_sum, detached
45 : character (len=90) :: rlo_result
46 :
47 : include 'formats'
48 : ierr = 0
49 0 : call binary_ptr(binary_id, b, ierr)
50 0 : if (ierr /= 0) then
51 0 : write(*,*) 'failed in binary_ptr'
52 0 : return
53 : end if
54 0 : s => b% s_donor
55 0 : use_sum = .false.
56 0 : detached = .false.
57 :
58 : ! NOTE: keep in mind that for mass loss, mdot is negative.
59 : ! b% mtransfer_rate will be considered valid if function_to_solve = 0
60 : ! within the tolerance given by b% implicit_scheme_tolerance, i.e.
61 : ! |function_to_solve| < b% implicit_scheme_tolerance.
62 : !
63 : ! For the roche_lobe scheme, function_to_solve is set such that solutions
64 : ! will satisfy 0 > (r-rl)/rl > - b% implicit_scheme_tolerance. If
65 : ! (r-rl)/rl < 0
66 : ! and the new mass transfer rate that was attempted satisfies
67 : ! |new_mdot| < b% roche_min_mdot*Msun/secyer
68 : ! then the system is considered to be detached and mass loss is turned off.
69 : !
70 : ! For other schemes, function_to_solve is chosen as the difference between
71 : ! b% mtransfer_rate and the explicit transfer rate, divided by the
72 : ! explicit transfer rate.
73 :
74 0 : check_implicit_rlo = keep_going
75 0 : new_mdot = b% mtransfer_rate
76 0 : b% num_tries = b% num_tries + 1
77 0 : rlo_result = "missing message"
78 :
79 0 : if(b% report_rlo_solver_progress .and. b% num_tries == 1) then
80 : write(*,'(a)') &
81 : ' binary_step rlo_iter di ch_fact curr_mdot' // &
82 0 : ' next_mdot mdot_up mdot_low rl_gap1 rl_gap2 f result'
83 : write(*,'(a)') &
84 : ' __________________________________________' // &
85 0 : '______________________________________________________________________________________'
86 : end if
87 :
88 0 : if (b% use_other_implicit_function_to_solve) then
89 : call b% other_implicit_function_to_solve(b% binary_id, &
90 0 : function_to_solve, use_sum, detached, ierr)
91 0 : if (ierr /= 0) then
92 0 : write(*,*) 'failed in other_implicit_function_to_solve'
93 0 : check_implicit_rlo = retry
94 0 : return
95 : end if
96 0 : if (detached) then
97 0 : if (b% report_rlo_solver_progress) then
98 0 : rlo_result = 'OK (detached)'
99 0 : call report_rlo_iter
100 : end if
101 : end if
102 0 : else if (b% mdot_scheme == "roche_lobe" .and. .not. b% use_other_rlo_mdot) then
103 : function_to_solve = (b% rl_relative_gap(b% d_i) &
104 0 : + b% implicit_scheme_tolerance/2.0d0) * 2.0d0
105 :
106 0 : if (function_to_solve < 0 .and. abs(b% mtransfer_rate) == 0) then
107 0 : if (b% report_rlo_solver_progress) then
108 0 : rlo_result = 'OK (detached)'
109 0 : call report_rlo_iter
110 : end if
111 0 : detached = .true.
112 : end if
113 0 : else if (b% mdot_scheme == "contact" .and. .not. b% use_other_rlo_mdot .and. .not. b% CE_flag) then
114 0 : if (b% point_mass_i /= 0) then
115 0 : new_mdot = 0d0
116 0 : write(*,*) "WARNING: contact scheme requires evolve_both_stars=.true."
117 0 : write(*,*) "Not transferring mass"
118 0 : return
119 : end if
120 0 : q = b% m(b% a_i) / b% m(b% d_i)
121 0 : slope_contact = pow(q, -0.52d0)
122 : ! If accretor is overflowing its Roche lobe, then the contact scheme needs to be used.
123 : ! Otherwise, if accretor radius is (within tolerance) below the equipotential
124 : ! of the donor, or donor is below tolerance for detachment, then use regular roche_lobe scheme.
125 0 : if (b% rl_relative_gap(b% a_i) < 0 .and. &
126 : (b% rl_relative_gap(b% d_i)*slope_contact - b% rl_relative_gap(b% a_i) &
127 : > b% implicit_scheme_tolerance .or. &
128 : b% rl_relative_gap(b% d_i) < - b% implicit_scheme_tolerance)) then
129 :
130 : function_to_solve = (b% rl_relative_gap(b% d_i) &
131 0 : + b% implicit_scheme_tolerance/2.0d0) * 2.0d0
132 :
133 0 : if (function_to_solve < 0 .and. abs(b% mtransfer_rate) == 0) then
134 0 : if (b% report_rlo_solver_progress) then
135 0 : rlo_result = 'OK (detached)'
136 0 : call report_rlo_iter
137 : end if
138 0 : detached = .true.
139 : end if
140 : else
141 0 : if (q < 1d0) then
142 0 : function_to_solve = (b% rl_relative_gap(b% d_i)*slope_contact - b% rl_relative_gap(b% a_i))
143 : else
144 0 : function_to_solve = (b% rl_relative_gap(b% d_i) - b% rl_relative_gap(b% a_i)/slope_contact)
145 : end if
146 0 : use_sum = .true.
147 : end if
148 : else
149 0 : if (b% CE_flag) then
150 0 : call CE_rlo_mdot(b% binary_id, explicit_mdot, ierr)
151 0 : if (ierr /= 0) then
152 0 : write(*,*) 'failed in CE_rlo_mdot'
153 0 : check_implicit_rlo = retry
154 0 : return
155 : end if
156 0 : else if (.not. b% use_other_rlo_mdot) then
157 0 : call rlo_mdot(b% binary_id, explicit_mdot, ierr)
158 0 : if (ierr /= 0) then
159 0 : write(*,*) 'failed in rlo_mdot'
160 0 : check_implicit_rlo = retry
161 0 : return
162 : end if
163 : else
164 0 : call b% other_rlo_mdot(b% binary_id, explicit_mdot, ierr)
165 0 : if (ierr /= 0) then
166 0 : write(*,*) 'failed in other_rlo_mdot'
167 0 : check_implicit_rlo = retry
168 0 : return
169 : end if
170 : end if
171 0 : if (explicit_mdot > 0d0) then
172 0 : new_mdot = 0d0
173 0 : write(*,*) "WARNING: explicit computation of mass transfer results in accreting donor"
174 0 : write(*,*) "Not transferring mass"
175 0 : return
176 : end if
177 0 : function_to_solve = 0d0
178 0 : if(explicit_mdot /= 0d0) then
179 0 : function_to_solve = (explicit_mdot - new_mdot) / explicit_mdot
180 0 : else if (new_mdot /= 0d0) then
181 0 : function_to_solve = (explicit_mdot - new_mdot) / new_mdot
182 : end if
183 0 : if (abs(explicit_mdot) <= b% min_mdot_for_implicit*Msun/secyer .and. &
184 : abs(new_mdot) <= b% min_mdot_for_implicit*Msun/secyer) then
185 0 : new_mdot = explicit_mdot
186 0 : if (b% report_rlo_solver_progress) then
187 0 : rlo_result = 'OK (explicit)'
188 0 : call report_rlo_iter
189 : end if
190 0 : return
191 : end if
192 : end if
193 :
194 0 : if (detached) return
195 :
196 0 : if (abs(function_to_solve) <= b% implicit_scheme_tolerance) then
197 0 : if (b% report_rlo_solver_progress) then
198 0 : rlo_result = 'OK (in tol)'
199 0 : call report_rlo_iter
200 0 : write(*,'(A)')
201 : end if
202 0 : return
203 : end if
204 :
205 0 : if (b% num_tries > b% max_tries_to_achieve) then
206 0 : check_implicit_rlo = retry
207 0 : if (b% report_rlo_solver_progress) then
208 0 : rlo_result = 'retry (>max tries)'
209 0 : call report_rlo_iter
210 : end if
211 0 : return
212 : end if
213 :
214 0 : if (b% num_tries == 1) then
215 0 : b% have_mdot_lo = .false.
216 0 : b% have_mdot_hi = .false.
217 0 : b% mdot_lo = 0
218 0 : b% mdot_hi = 0
219 0 : b% fixed_delta_mdot = b% mtransfer_rate * (1 - b% change_factor)
220 : end if
221 :
222 0 : new_mdot = pick_mdot_for_implicit_rlo(b, function_to_solve, b% mtransfer_rate, use_sum, ierr)
223 :
224 : !if this iteration is done using the maximum mass transfer rate,
225 : !and the mass transfer rate increases in the following iteration,
226 : !then stop at this point and accept the step
227 : !NOTE: both b% mtransfer_rate and new_mdot are negative
228 : if (b% mtransfer_rate == -b% max_implicit_abs_mdot*Msun/secyer &
229 0 : .and. b% mtransfer_rate > new_mdot) then
230 0 : if (b% report_rlo_solver_progress) then
231 0 : rlo_result = 'OK (max mdot)'
232 0 : call report_rlo_iter
233 : end if
234 0 : new_mdot = - b% max_implicit_abs_mdot*Msun/secyer
235 0 : return
236 : end if
237 0 : if (-new_mdot > b% max_implicit_abs_mdot*Msun/secyer) then
238 : !limit to maximum mdot, following iteration will be made
239 : !using this value
240 0 : new_mdot = - b% max_implicit_abs_mdot*Msun/secyer
241 : end if
242 :
243 0 : if (ierr /= 0) then
244 0 : check_implicit_rlo = retry
245 0 : if (b% report_rlo_solver_progress) then
246 0 : rlo_result = 'retry (ierr)'
247 0 : call report_rlo_iter
248 : end if
249 0 : return
250 : end if
251 :
252 0 : if (-new_mdot < b% roche_min_mdot*Msun/secyer .and. function_to_solve < 0 .and. &
253 : (b% mdot_scheme == "roche_lobe" .or. (b% mdot_scheme == "contact" .and. &
254 : .not. use_sum))) then
255 0 : check_implicit_rlo = keep_going
256 0 : new_mdot = 0d0
257 0 : if (b% report_rlo_solver_progress) then
258 0 : rlo_result = 'OK (detachment)'
259 0 : call report_rlo_iter
260 : end if
261 0 : return
262 : end if
263 :
264 0 : if (b% have_mdot_hi .and. b% have_mdot_lo) then
265 0 : if (abs(b% mdot_hi - b% mdot_lo) < &
266 : b% implicit_scheme_tiny_factor*min(abs(b% mdot_hi),abs(b% mdot_lo))) then
267 0 : if (b% report_rlo_solver_progress) then
268 0 : rlo_result = 'OK (tiny change)'
269 0 : call report_rlo_iter
270 : end if
271 0 : new_mdot = b% mtransfer_rate
272 0 : check_implicit_rlo = keep_going
273 0 : return
274 : end if
275 : end if
276 :
277 : ! boost change_factor every num_tries_for_increase_change_factor tries, to avoid
278 : ! implicit scheme from becoming stuck with sudden changes.
279 : if (b% num_tries_for_increase_change_factor > 0 .and. b% change_factor_increase > 1d0 &
280 0 : .and. .not. (b% have_mdot_lo .and. b% have_mdot_hi)) then
281 0 : if(mod(b% num_tries, b% num_tries_for_increase_change_factor) == 0) then
282 : b% change_factor = min(b% max_change_factor, &
283 0 : b% change_factor_increase * b% change_factor)
284 : end if
285 : end if
286 0 : if (b% report_rlo_solver_progress) then
287 0 : rlo_result = 'redo'
288 0 : call report_rlo_iter
289 : end if
290 :
291 0 : check_implicit_rlo = redo
292 :
293 : contains
294 :
295 0 : subroutine report_rlo_iter
296 0 : real(dp) :: mdot_hi, mdot_lo
297 0 : if (.not. b% have_mdot_hi) then
298 0 : mdot_hi = huge(mdot_hi)
299 : else
300 0 : mdot_hi = -b% mdot_hi/Msun*secyer
301 : end if
302 0 : if (.not. b% have_mdot_lo) then
303 0 : mdot_lo = -huge(mdot_lo)
304 : else
305 0 : mdot_lo = -b% mdot_lo/Msun*secyer
306 : end if
307 : write(*,'(i13,i9,i3,f8.4,7(1pe11.3,0p),a20)') &
308 0 : b% model_number, &
309 0 : b% num_tries, &
310 0 : b% d_i, &
311 0 : b% change_factor, &
312 0 : -b% mtransfer_rate/Msun*secyer, &
313 0 : -new_mdot/Msun*secyer, &
314 0 : mdot_hi, &
315 0 : mdot_lo, &
316 0 : b% rl_relative_gap(1), &
317 0 : b% rl_relative_gap(2), &
318 0 : function_to_solve, &
319 0 : trim(rlo_result)
320 0 : end subroutine report_rlo_iter
321 :
322 : end function check_implicit_rlo
323 :
324 0 : real(dp) function pick_mdot_for_implicit_rlo( &
325 0 : b, new_function_to_solve, mdot_current, use_sum, ierr) result(mdot_next)
326 : use num_lib, only: find0_quadratic, find0
327 : type(binary_info), pointer :: b
328 : real(dp), intent(in) :: new_function_to_solve, mdot_current
329 : logical, intent(in) :: use_sum
330 : integer, intent(out) :: ierr
331 :
332 0 : real(dp) :: starting_mdot, current_change_factor
333 : logical :: do_cubic
334 : include 'formats'
335 :
336 : ! NOTE: keep in mind that for mass loss, mdot is negative
337 :
338 :
339 0 : starting_mdot = -b% starting_mdot*Msun/secyer
340 0 : current_change_factor = pow(b% change_factor, b% num_tries+1)
341 :
342 0 : if (.not.(b% solver_type == "cubic" .or. &
343 : b% solver_type == "bisect" .or. &
344 : b% solver_type == "both")) then
345 0 : write(*,*) "ERROR: unrecognized b% solver_type", trim(b% solver_type)
346 0 : write(*,*) "setting b% solver_type to 'both'"
347 0 : b% solver_type = "both"
348 : end if
349 :
350 0 : if (b% have_mdot_lo .and. b% have_mdot_hi) then
351 0 : ierr = 0
352 : do_cubic = (mod(b% num_tries,2)==0 .and. b% solver_type == "both") .or. &
353 0 : b% solver_type == "cubic"
354 : if (do_cubic) then
355 : mdot_next = find0_quadratic( &
356 : b% mdot_lo, b% implicit_function_lo, &
357 : mdot_current, new_function_to_solve, &
358 0 : b% mdot_hi, b% implicit_function_hi, ierr)
359 0 : if (ierr /= 0) then
360 : mdot_next = find0(b% mdot_lo, b% implicit_function_lo, &
361 0 : b% mdot_hi, b% implicit_function_hi)
362 0 : ierr = 0
363 : end if
364 : end if
365 0 : if (new_function_to_solve >= 0) then
366 0 : b% mdot_lo = mdot_current
367 0 : b% implicit_function_lo = new_function_to_solve
368 : else
369 0 : b% mdot_hi = mdot_current
370 0 : b% implicit_function_hi = new_function_to_solve
371 : end if
372 0 : if (.not. do_cubic) then
373 0 : mdot_next = (b% mdot_hi+b% mdot_lo)/2.0d0
374 : end if
375 0 : else if (b% have_mdot_lo) then ! don't have mdot_hi
376 0 : if (new_function_to_solve < 0) then
377 0 : b% mdot_hi = mdot_current
378 0 : b% implicit_function_hi = new_function_to_solve
379 0 : b% have_mdot_hi = .true.
380 0 : mdot_next = (b% mdot_hi+b% mdot_lo)/2.0d0
381 : !mdot_next = find0(b% mdot_lo, b% implicit_function_lo, &
382 : ! b% mdot_hi, b% implicit_function_hi)
383 : else ! still too low
384 0 : b% mdot_lo = mdot_current
385 0 : b% implicit_function_lo = new_function_to_solve
386 0 : mdot_next = b% mdot_lo*b% change_factor
387 : end if
388 0 : else if (b% have_mdot_hi) then ! don't have mdot_lo
389 0 : if (new_function_to_solve >= 0) then
390 0 : b% mdot_lo = mdot_current
391 0 : b% implicit_function_lo = new_function_to_solve
392 0 : b% have_mdot_lo = .true.
393 0 : mdot_next = (b% mdot_hi+b% mdot_lo)/2.0d0
394 : !mdot_next = find0(b% mdot_lo, b% implicit_function_lo, &
395 : ! b% mdot_hi, b% implicit_function_hi)
396 : else ! mdot still too high
397 0 : b% mdot_hi = mdot_current
398 0 : b% implicit_function_hi = new_function_to_solve
399 0 : if (use_sum .and. mod(b% num_tries, 2) == 0) then
400 0 : mdot_next = b% mdot_hi + b% fixed_delta_mdot
401 : else
402 0 : mdot_next = b% mdot_hi/b% change_factor
403 : end if
404 : end if
405 : else ! don't have either
406 0 : if (mdot_current > starting_mdot .and. (.not. abs(mdot_current) > 0)) then ! recall that both are negative
407 : mdot_next = starting_mdot
408 0 : else if (new_function_to_solve >= 0) then
409 0 : b% mdot_lo = mdot_current
410 0 : b% implicit_function_lo = new_function_to_solve
411 0 : b% have_mdot_lo = .true.
412 0 : mdot_next = b% mdot_lo*b% change_factor
413 : else
414 0 : b% mdot_hi = mdot_current
415 0 : b% implicit_function_hi = new_function_to_solve
416 0 : b% have_mdot_hi = .true.
417 0 : if (use_sum .and. mod(b% num_tries, 2) == 0) then
418 0 : mdot_next = b% mdot_hi + b% fixed_delta_mdot
419 : else
420 0 : mdot_next = b% mdot_hi/b% change_factor
421 : end if
422 : end if
423 : end if
424 :
425 0 : end function pick_mdot_for_implicit_rlo
426 :
427 :
428 0 : subroutine eval_mdot_edd(binary_id, mdot_edd, mdot_edd_eta, ierr)
429 0 : use utils_lib, only: is_bad
430 :
431 : integer, intent(in) :: binary_id
432 : real(dp), intent(out) :: mdot_edd ! eddington accretion rate
433 : real(dp), intent(out) :: mdot_edd_eta ! fraction of rest mass energy released as radiation
434 : integer, intent(out) :: ierr
435 : type(binary_info), pointer :: b
436 : include 'formats'
437 :
438 : ierr = 0
439 0 : call binary_ptr(binary_id, b, ierr)
440 0 : if (ierr /= 0) then
441 0 : write(*,*) 'failed in binary_ptr'
442 0 : return
443 : end if
444 :
445 0 : if (b% point_mass_i == 0) then
446 0 : if (b% limit_retention_by_mdot_edd) then
447 0 : write(*,*) "Default mdot_edd calculation cannot be used when evolving both stars"
448 0 : write(*,*) "Maybe you want to set limit_retention_by_mdot_edd=.false. in binary_controls?"
449 0 : write(*,*) "Setting mdot_edd to zero"
450 : end if
451 0 : mdot_edd = 0
452 0 : mdot_edd_eta = 0
453 0 : return
454 : end if
455 :
456 0 : if (b% use_this_for_mdot_edd_eta > 0) then
457 0 : mdot_edd_eta = b% use_this_for_mdot_edd_eta
458 : else
459 : ! eg., eq. (6) of Podsiadlowski, Rappaport & Han 2003, MNRAS, 341, 385
460 : mdot_edd_eta = 1d0 &
461 0 : - sqrt(1d0 - pow2(min(b% m(b% a_i),sqrt(6d0)*b% eq_initial_bh_mass)/(3d0*b% eq_initial_bh_mass)))
462 : end if
463 :
464 0 : if (b% use_this_for_mdot_edd > 0) then
465 0 : mdot_edd = b% use_this_for_mdot_edd*(Msun/secyer)
466 : else
467 : ! eg., eq. (9) of Podsiadlowski, Rappaport & Han 2003, MNRAS, 341, 385
468 0 : if (.not. b% use_es_opacity_for_mdot_edd) then
469 : mdot_edd = pi4*standard_cgrav*b% m(b% a_i) &
470 0 : /(clight*b% s_donor% opacity(1)*mdot_edd_eta)
471 : else
472 : mdot_edd = pi4*standard_cgrav*b% m(b% a_i)&
473 0 : /(clight*0.2d0*(1d0+b% s_donor% surface_h1)*mdot_edd_eta)
474 : end if
475 : end if
476 :
477 0 : if (is_bad(mdot_edd_eta) .or. mdot_edd_eta<0) then
478 0 : write(*,*) "ERROR while computing mdot_edd_eta"
479 0 : ierr = -1
480 0 : write(*,*) "mdot_edd_eta, b% m(b% a_i), b% eq_initial_bh_mass", &
481 0 : mdot_edd_eta, b% m(b% a_i), b% eq_initial_bh_mass
482 0 : stop
483 : end if
484 :
485 0 : if (is_bad(mdot_edd) .or. mdot_edd<0) then
486 0 : write(*,*) "ERROR while computing mdot_edd"
487 0 : ierr = -1
488 0 : write(*,*) "mdot_edd, b% m(b% a_i), b% s_donor% opacity(1), b% s_donor% surface_h1", &
489 0 : mdot_edd, b% m(b% a_i), b% s_donor% opacity(1), b% s_donor% surface_h1
490 0 : stop
491 : end if
492 :
493 : end subroutine eval_mdot_edd
494 :
495 0 : subroutine adjust_mdots(b)
496 : use binary_wind, only: eval_wind_xfer_fractions
497 : type (binary_info), pointer :: b
498 :
499 : real(dp) :: actual_mtransfer_rate
500 : integer :: ierr
501 :
502 0 : actual_mtransfer_rate = 0d0
503 :
504 0 : if (b% use_other_adjust_mdots) then
505 0 : call b% other_adjust_mdots(b% binary_id, ierr)
506 0 : if (ierr /= 0) then
507 0 : write(*,*) "Error in other_adjust_mdots"
508 0 : stop
509 : end if
510 : return
511 : end if
512 :
513 : b% fixed_xfer_fraction = 1 - b% mass_transfer_alpha - b% mass_transfer_beta - &
514 0 : b% mass_transfer_delta
515 :
516 0 : if (.not. b% use_other_mdot_edd) then
517 0 : call eval_mdot_edd(b% binary_id, b% mdot_edd, b% mdot_edd_eta, ierr)
518 : else
519 0 : call b% other_mdot_edd(b% binary_id, b% mdot_edd, b% mdot_edd_eta, ierr)
520 : end if
521 :
522 : ! Add tidal enhancement of wind
523 0 : call Tout_enhance_wind(b, b% s_donor)
524 0 : if (b% point_mass_i == 0) then
525 : ! do not repeat if using the implicit wind
526 0 : if (.not. (b% num_tries >0 .and. b% s_accretor% was_in_implicit_wind_limit)) &
527 0 : call Tout_enhance_wind(b, b% s_accretor)
528 : end if
529 :
530 : ! solve wind mass transfer
531 : ! b% mdot_wind_transfer(b% d_i) is a negative number that gives the
532 : ! amount of mass transferred by unit time from the donor to the
533 : ! accretor.
534 0 : call eval_wind_xfer_fractions(b% binary_id, ierr)
535 0 : if (ierr/=0) then
536 0 : write(*,*) "Error in eval_wind_xfer_fractions"
537 0 : return
538 : end if
539 : b% mdot_wind_transfer(b% d_i) = b% s_donor% mstar_dot * &
540 0 : b% wind_xfer_fraction(b% d_i)
541 0 : if (b% point_mass_i == 0) then
542 : b% mdot_wind_transfer(b% a_i) = b% s_accretor% mstar_dot * &
543 0 : b% wind_xfer_fraction(b% a_i)
544 : else
545 0 : b% mdot_wind_transfer(b% a_i) = 0d0
546 : end if
547 :
548 : ! Set mdot for the donor
549 : b% s_donor% mstar_dot = b% s_donor% mstar_dot + b% mtransfer_rate - &
550 0 : b% mdot_wind_transfer(b% a_i)
551 :
552 : ! Set mdot for the accretor
553 0 : if (b% point_mass_i == 0 .and. .not. b% CE_flag) then
554 : ! do not repeat if using the implicit wind
555 0 : if (.not. (b% num_tries >0 .and. b% s_accretor% was_in_implicit_wind_limit)) then
556 0 : b% accretion_mode = 0
557 0 : b% acc_am_div_kep_am = 0.0d0
558 : b% s_accretor% mstar_dot = b% s_accretor% mstar_dot - &
559 0 : b% mtransfer_rate*b% fixed_xfer_fraction - b% mdot_wind_transfer(b% d_i)
560 :
561 : !set angular momentum accretion as described in A.3.3 of de Mink et al. 2013
562 0 : if (b% do_j_accretion) then
563 0 : if (.not. b% use_other_accreted_material_j) then
564 0 : call eval_accreted_material_j(b% binary_id, ierr)
565 : else
566 0 : call b% other_accreted_material_j(b% binary_id, ierr)
567 : end if
568 0 : if (ierr /= 0) then
569 0 : write(*,*) 'error in accreted_material_j'
570 0 : return
571 : end if
572 : end if
573 : end if
574 :
575 0 : b% accretion_luminosity = 0d0 !only set for point mass
576 :
577 0 : else if (.not. b% CE_flag) then
578 : ! accretor is a point mass
579 0 : if (.not. b% model_twins_flag) then
580 : !combine wind and RLOF mass transfer
581 0 : actual_mtransfer_rate = b% mtransfer_rate*b% fixed_xfer_fraction+b% mdot_wind_transfer(b% d_i) !defined negative
582 0 : b% component_mdot(b% a_i) = -actual_mtransfer_rate
583 : ! restrict accretion to the Eddington limit
584 0 : if (b% limit_retention_by_mdot_edd .and. b% component_mdot(b% a_i) > b% mdot_edd) then
585 0 : b% component_mdot(b% a_i) = b% mdot_edd ! remove all accretion above the edd limit
586 : end if
587 : b% accretion_luminosity = &
588 0 : b% mdot_edd_eta*b% component_mdot(b% a_i)*clight*clight
589 : ! remove rest mass radiated away
590 0 : if (b% use_radiation_corrected_transfer_rate) then
591 0 : b% component_mdot(b% a_i) = (1 - b% mdot_edd_eta) * b% component_mdot(b% a_i)
592 : end if
593 : end if
594 : else
595 : !doing CE, just be sure to set mdot for a point mass to zero
596 0 : b % accretion_luminosity = 0d0
597 0 : if (b% point_mass_i /= 0) then
598 0 : b% component_mdot(b% a_i) = 0d0
599 : end if
600 : end if
601 :
602 : ! mdot_system_transfer is mass lost from the vicinity of each star
603 : ! due to inefficient rlof mass transfer, mdot_system_cct is mass lost
604 : ! from a circumbinary coplanar toroid.
605 0 : if (b% mtransfer_rate+b% mdot_wind_transfer(b% d_i) >= 0 .or. b% CE_flag) then
606 0 : b% mdot_system_transfer(b% d_i) = 0d0
607 0 : b% mdot_system_transfer(b% a_i) = 0d0
608 0 : b% mdot_system_cct = 0d0
609 : else
610 0 : b% mdot_system_transfer(b% d_i) = b% mtransfer_rate * b% mass_transfer_alpha
611 0 : b% mdot_system_cct = b% mtransfer_rate * b% mass_transfer_delta
612 0 : if (b% point_mass_i == 0 .or. b% model_twins_flag) then
613 0 : b% mdot_system_transfer(b% a_i) = b% mtransfer_rate * b% mass_transfer_beta
614 : else
615 : ! do not compute mass lost from the vicinity using just mass_transfer_beta, as
616 : ! mass transfer can be stopped also by going past the Eddington limit
617 : b% mdot_system_transfer(b% a_i) = (actual_mtransfer_rate + b% component_mdot(b% a_i)) &
618 0 : + b% mtransfer_rate * b% mass_transfer_beta
619 : end if
620 : end if
621 :
622 0 : end subroutine adjust_mdots
623 :
624 0 : subroutine rlo_mdot(binary_id, mdot, ierr) ! Adapted from a routine kindly provided by Anastasios Fragkos
625 : integer, intent(in) :: binary_id
626 : real(dp), intent(out) :: mdot
627 : integer, intent(out) :: ierr
628 : type (binary_info), pointer :: b
629 :
630 : include 'formats'
631 :
632 : ierr = 0
633 0 : call binary_ptr(binary_id, b, ierr)
634 0 : if (ierr /= 0) then
635 0 : write(*,*) 'failed in binary_ptr'
636 0 : return
637 : end if
638 :
639 0 : mdot = 0d0
640 :
641 0 : if (b% mdot_scheme == "roche_lobe") then
642 0 : write(*,*) "mdot_scheme = roche_lobe not applicable for explicit scheme"
643 0 : write(*,*) "Not transferring mass"
644 0 : mdot = 0
645 0 : return
646 0 : else if (b% mdot_scheme /= "Ritter" .and. b% mdot_scheme /= "Kolb" .and. b% mdot_scheme /= "Arras") then
647 0 : write(*,*) "mdot_scheme = " , b% mdot_scheme , " not recognized"
648 0 : write(*,*) "Not transferring mass"
649 0 : mdot = 0
650 0 : return
651 : end if
652 :
653 0 : if (b% mdot_scheme == "Kolb" .and. b% eccentricity <= 0.0d0) then
654 0 : call get_info_for_ritter(b)
655 0 : mdot = b% mdot_thin
656 0 : call get_info_for_kolb(b)
657 0 : mdot = mdot + b% mdot_thick
658 :
659 0 : else if (b% mdot_scheme == "Kolb" .and. b% eccentricity > 0.0d0) then
660 0 : call get_info_for_ritter_eccentric(b)
661 0 : mdot = b% mdot_thin
662 0 : call get_info_for_kolb_eccentric(b)
663 0 : mdot = mdot + b% mdot_thick
664 :
665 0 : else if (b% mdot_scheme == "Ritter" .and. b% eccentricity <= 0.0d0) then
666 0 : call get_info_for_ritter(b)
667 0 : mdot = b% mdot_thin
668 :
669 0 : else if (b% mdot_scheme == "Ritter" .and. b% eccentricity > 0.0d0) then
670 0 : call get_info_for_ritter_eccentric(b)
671 0 : mdot = b% mdot_thin
672 :
673 : end if
674 :
675 0 : if (b% mdot_scheme == "Arras") then
676 0 : if (b% eccentricity > 0d0) &
677 0 : write(*,*) "mdot_scheme = Arras is not properly implemented for e>0"
678 0 : call get_info_for_arras(b)
679 0 : mdot = b% mdot_thin
680 :
681 : end if
682 :
683 0 : end subroutine rlo_mdot
684 :
685 0 : subroutine get_info_for_ritter(b)
686 : type(binary_info), pointer :: b
687 0 : real(dp) :: F1, q, rho, p, grav, hp, v_th, rl3, q_temp
688 : include 'formats'
689 :
690 : !--------------------- Optically thin MT rate -----------------------------------------------
691 : ! As described in H. Ritter 1988, A&A 202,93-100 and U. Kolb and H. Ritter 1990, A&A 236,385-392
692 :
693 0 : rho = b% s_donor% rho(1) ! density at surface in g/cm^3
694 0 : p = b% s_donor% Peos(1) ! pressure at surface in dynes/cm^2
695 0 : grav = standard_cgrav*b% m(b% d_i)/pow2(b% r(b% d_i)) ! local gravitational acceleration
696 0 : hp = p/(grav*rho) ! pressure scale height
697 0 : v_th = sqrt(kerg * b% s_donor% T(1) / (mp * b% s_donor% mu(1)))
698 :
699 0 : q = b% m(b% a_i)/b% m(b% d_i) ! Mass ratio, as defined in Ritter 1988
700 : ! (Kolb & Ritter 1990 use the opposite!)
701 : ! consider range of validity for F1, do not extrapolate! Eq. A9 of Ritter 1988
702 0 : q_temp = min(max(q,0.5d0),10d0)
703 0 : F1 = (1.23d0 + 0.5D0* log10(q_temp))
704 0 : rl3 = (b% rl(b% d_i))*(b% rl(b% d_i))*(b% rl(b% d_i))
705 : b% mdot_thin0 = (2.0D0*pi/exp(0.5d0)) * v_th*v_th*v_th * &
706 0 : rl3/(standard_cgrav*b% m(b% d_i)) * rho * F1
707 : !Once again, do not extrapolate! Eq. (7) of Ritter 1988
708 0 : q_temp = min(max(q,0.04d0),20d0)
709 0 : if (q_temp < 1.0d0) then
710 0 : b% ritter_h = hp/( 0.954D0 + 0.025D0*log10(q_temp) - 0.038D0*pow2(log10(q_temp)) )
711 : else
712 0 : b% ritter_h = hp/( 0.954D0 + 0.039D0*log10(q_temp) + 0.114D0*pow2(log10(q_temp)) )
713 : end if
714 :
715 0 : b% ritter_exponent = (b% r(b% d_i)-b% rl(b% d_i))/b% ritter_h
716 :
717 0 : if (b% mdot_scheme == "Kolb") then
718 0 : if (b% ritter_exponent > 0) then
719 0 : b% mdot_thin = -b% mdot_thin0
720 : else
721 0 : b% mdot_thin = -b% mdot_thin0 * exp(b% ritter_exponent)
722 : end if
723 : else
724 0 : b% mdot_thin = -b% mdot_thin0 * exp(b% ritter_exponent)
725 : end if
726 :
727 0 : end subroutine get_info_for_ritter
728 :
729 0 : real(dp) function calculate_kolb_mdot_thick(b, indexR, rl_d) result(mdot_thick)
730 : real(dp), intent(in) :: rl_d
731 : integer, intent(in) :: indexR
732 0 : real(dp) :: F1, F3, G1, d_P, q, q_temp
733 : integer :: i
734 : type(binary_info), pointer :: b
735 : include 'formats'
736 :
737 : !--------------------- Optically thin MT rate -----------------------------------------------
738 : ! As described in Kolb and H. Ritter 1990, A&A 236,385-392
739 :
740 : ! compute integral in Eq. (A17 of Kolb & Ritter 1990)
741 0 : mdot_thick = 0d0
742 0 : do i=1,indexR-1
743 0 : G1 = b% s_donor% gamma1(i)
744 0 : F3 = sqrt(G1) * pow(2d0/(G1+1d0), (G1+1d0)/(2d0*G1-2d0))
745 : mdot_thick = mdot_thick + F3*sqrt(kerg * b% s_donor% T(i) / &
746 0 : (mp * b% s_donor% mu(i)))*(b% s_donor% Peos(i+1)-b% s_donor% Peos(i))
747 : end do
748 : ! only take a fraction of d_P for last cell
749 0 : G1 = b% s_donor% gamma1(i)
750 0 : F3 = sqrt(G1) * pow(2d0/(G1+1d0), (G1+1d0)/(2d0*G1-2d0))
751 : d_P = (b% s_donor% r(indexR) - rl_d) / &
752 0 : (b% s_donor% r(indexR) - b% s_donor% r(indexR+1)) * (b% s_donor% Peos(i+1)-b% s_donor% Peos(i))
753 0 : mdot_thick = mdot_thick + F3*sqrt(kerg * b% s_donor% T(i) / (mp*b% s_donor% mu(i)))*d_P
754 :
755 0 : q = b% m(b% a_i)/b% m(b% d_i) ! Mass ratio, as defined in Ritter 1988
756 : ! (Kolb & Ritter 1990 use the opposite!)
757 : ! consider range of validity for F1, do not extrapolate! Eq. A9 of Ritter 1988
758 0 : q_temp = min(max(q,0.5d0),10d0)
759 0 : F1 = (1.23d0 + 0.5D0* log10(q_temp))
760 0 : mdot_thick = -2.0D0*pi*F1*rl_d*rl_d*rl_d/(standard_cgrav*b% m(b% d_i))*mdot_thick
761 :
762 0 : end function calculate_kolb_mdot_thick
763 :
764 0 : subroutine get_info_for_kolb(b)
765 : type(binary_info), pointer :: b
766 : integer :: i, indexR
767 : include 'formats'
768 :
769 : !--------------------- Optically thick MT rate -----------------------------------------------
770 : ! As described in H. Ritter 1988, A&A 202,93-100 and U. Kolb and H. Ritter 1990, A&A 236,385-392
771 :
772 : ! First we need to find how deep inside the star the Roche lobe reaches. In other words the mesh point of the star at which R=R_RL
773 0 : b% mdot_thick = 0d0
774 0 : indexR=-1
775 0 : if(b% r(b% d_i)-b% rl(b% d_i) > 0.0d0) then
776 : i=1
777 0 : do while (b% s_donor% r(i) > b% rl(b% d_i))
778 0 : i=i+1
779 : end do
780 :
781 0 : if (i == 1) then
782 : b% mdot_thick = 0d0
783 : else
784 0 : b% mdot_thick = calculate_kolb_mdot_thick(b, i-1, b% rl(b% d_i))
785 : end if
786 : end if
787 :
788 0 : end subroutine get_info_for_kolb
789 :
790 0 : subroutine get_info_for_arras(b)
791 : type(binary_info), pointer :: b
792 0 : real(dp) :: q, rho, p, grav, hp, v_th
793 0 : real(dp) :: area,Asl,G,ma,md,mfac1,mfac2,my_mdot_thin,my_ritter_exponent,Omega,&
794 0 : phi,phiL1,q13,rfac,rhoL1,rv,rvL1,sep
795 : include 'formats'
796 :
797 : !--------------------- Optically thin MT rate -----------------------------------------------
798 : ! Ritter 1988 but with better fits for the various formulas that work at extreme q
799 :
800 0 : rho = b% s_donor% rho(1) ! density at surface in g/cm^3
801 0 : p = b% s_donor% Peos(1) ! pressure at surface in dynes/cm^2
802 0 : grav = standard_cgrav*b% m(b% d_i)/pow2(b% r(b% d_i)) ! local gravitational acceleration
803 0 : hp = p/(grav*rho) ! pressure scale height
804 0 : v_th = sqrt(kerg * b% s_donor% T(1) / (mp * b% s_donor% mu(1)))
805 :
806 0 : q = b% m(b% a_i) / b% m(b% d_i)
807 0 : G = standard_cgrav
808 0 : md = b% m(b% d_i)
809 0 : ma = b% m(b% a_i)
810 0 : sep = b% separation
811 0 : Omega = 2.d0*pi / b% period
812 0 : rvL1 = b% rl(b% d_i)
813 0 : rv = b% r(b% d_i)
814 0 : mfac1 = 1d0 + ma/md ! (md+ma)/md
815 : ! mfac2 = ( (md+ma)**2 + 3d0*ma*(md+ma) + 9d0*ma**2 ) / md**2
816 0 : mfac2 = 1d0 + 5d0*ma/md + 13d0*ma*ma/(md*md)
817 0 : rfac=rvL1/sep
818 : phiL1 = - G*md/rvL1 &
819 0 : * ( 1.d0 + mfac1*pow3(rfac)/3.d0 + 4.d0*mfac2*pow6(rfac)/45.d0 )
820 0 : rfac=rv/sep
821 : phi = - G*md/rv &
822 0 : * ( 1.d0 + mfac1*pow3(rfac)/3.d0 + 4.d0*mfac2*pow6(rfac)/45.d0 )
823 0 : my_ritter_exponent = - (phiL1-phi)/(v_th*v_th)
824 0 : rhoL1 = rho/sqrt(exp(1.d0)) * exp( my_ritter_exponent )
825 0 : q13=pow(q,one_third)
826 0 : Asl = 4.d0 + 4.16d0/(-0.96d0 + q13 + 1.d0/q13)
827 0 : area = 2.d0 * pi * pow2(v_th/Omega) / sqrt( Asl*(Asl-1d0) )
828 0 : my_mdot_thin = - rhoL1 * v_th * area
829 0 : b% mdot_thin = my_mdot_thin
830 :
831 0 : end subroutine get_info_for_arras
832 :
833 0 : subroutine get_info_for_ritter_eccentric(b)
834 : type(binary_info), pointer :: b
835 : integer :: i
836 0 : real(dp) :: F1, q, q_temp, rho, p, grav, hp, v_th, dm
837 0 : real(dp), DIMENSION(b% anomaly_steps):: mdot0, mdot, Erit, rl_d
838 : include 'formats'
839 :
840 : ! Optically thin MT rate adapted for eccentric orbits
841 : ! As described in H. Ritter 1988, A&A 202,93-100 and U. Kolb and H. Ritter 1990, A&A 236,385-392
842 :
843 0 : rho = b% s_donor% rho(1) ! density at surface in g/cm^3
844 0 : p = b% s_donor% Peos(1) ! pressure at surface in dynes/cm^2
845 0 : grav = standard_cgrav*b% m(b% d_i)/pow2(b% r(b% d_i)) ! local gravitational acceleration
846 0 : hp = p/(grav*rho) ! pressure scale height
847 0 : v_th = sqrt(kerg * b% s_donor% T(1) / (mp * b% s_donor% mu(1))) ! kerg = Boltzmann's constant
848 :
849 : ! phase dependant RL radius
850 0 : do i = 1, b% anomaly_steps
851 : rl_d(i) = b% rl(b% d_i) * (1d0 - pow2(b% eccentricity)) / &
852 0 : (1 + b% eccentricity * cos(b% theta_co(i)) )
853 : end do
854 :
855 0 : q = b% m(b% a_i)/b% m(b% d_i) ! Mass ratio, as defined in Ritter 1988
856 : ! (Kolb & Ritter 1990 use the opposite!)
857 0 : q_temp = min(max(q,0.5d0),10d0)
858 0 : F1 = (1.23d0 + 0.5D0* log10(q_temp))
859 :
860 : mdot0 = (2.0D0*pi/exp(0.5d0)) * pow3(v_th) * rl_d*rl_d*rl_d / &
861 0 : (standard_cgrav*b% m(b% d_i)) * rho * F1
862 :
863 0 : q_temp = min(max(q,0.04d0),20d0)
864 0 : if (q_temp < 1.0d0) then
865 0 : b% ritter_h = hp/( 0.954D0 + 0.025D0*log10(q_temp) - 0.038D0*pow2(log10(q_temp)) )
866 : else
867 0 : b% ritter_h = hp/( 0.954D0 + 0.039D0*log10(q_temp) + 0.114D0*pow2(log10(q_temp)) )
868 : end if
869 :
870 0 : Erit = (b% r(b% d_i)- rl_d) / b% ritter_h
871 :
872 0 : if (b% mdot_scheme == "Kolb") then
873 0 : do i = 1,b% anomaly_steps
874 0 : if (Erit(i) > 0) then
875 0 : mdot(i) = -1 * mdot0(i)
876 : else
877 0 : mdot(i) = -1 * mdot0(i) * exp(Erit(i))
878 : end if
879 : end do
880 : else
881 0 : do i = 1,b% anomaly_steps
882 0 : mdot(i) = -1 * mdot0(i) * exp(Erit(i))
883 : end do
884 : end if
885 :
886 0 : b% mdot_donor_theta = mdot
887 :
888 : !integrate to get total massloss
889 : dm = 0d0
890 0 : do i = 2,b% anomaly_steps ! trapezoidal integration
891 0 : dm = dm + 0.5d0 * (mdot(i-1) + mdot(i)) * (b% time_co(i) - b% time_co(i-1))
892 : end do
893 :
894 0 : b% mdot_thin = dm
895 :
896 0 : end subroutine get_info_for_ritter_eccentric
897 :
898 0 : subroutine get_info_for_kolb_eccentric(b)
899 : type(binary_info), pointer :: b
900 0 : real(dp) :: e, dm
901 : integer :: i, j
902 0 : real(dp), DIMENSION(b% anomaly_steps):: rl_d_i, mdot_thick_i
903 : include 'formats'
904 :
905 : ! Optically thick MT rate adapted for eccentric orbits
906 : ! As described in H. Ritter 1988, A&A 202,93-100 and U. Kolb and H. Ritter 1990, A&A 236,385-392
907 :
908 0 : b% mdot_thick = 0d0
909 0 : e = b% eccentricity
910 :
911 : ! If the radius of the donor is smaller as the smallest RL radius,
912 : ! there is only atmospheric RLOF, thus return.
913 0 : if ( b% r(b% d_i) < b% rl(b% d_i) * (1-e*e)/(1+e) ) then
914 : return
915 : end if
916 :
917 : ! For each point in the orbit calculate mdot_thick
918 0 : do i = 1,b% anomaly_steps
919 : ! phase dependent RL radius
920 : rl_d_i(i) = b% rl(b% d_i) * (1d0 - e*e) / &
921 0 : (1 + e*cos(b% theta_co(i)) )
922 :
923 : ! find how deep in the star we are
924 0 : j=1
925 0 : do while (b% s_donor% r(j) > rl_d_i(i))
926 0 : j=j+1
927 : end do
928 :
929 : ! calculate mdot_thick
930 0 : if (j == 1) then
931 0 : mdot_thick_i(i) = 0d0
932 : else
933 0 : mdot_thick_i(i) = calculate_kolb_mdot_thick(b, j-1, rl_d_i(i))
934 : end if
935 : end do
936 :
937 0 : b% mdot_donor_theta = b% mdot_donor_theta + mdot_thick_i
938 :
939 : ! Integrate mdot_thick over the orbit
940 : dm = 0d0
941 0 : do i = 2,b% anomaly_steps ! trapezoidal integration
942 : dm = dm + 0.5d0 * (mdot_thick_i(i-1) + mdot_thick_i(i)) * &
943 0 : (b% time_co(i) - b% time_co(i-1))
944 : end do
945 :
946 0 : b% mdot_thick = dm
947 :
948 0 : end subroutine get_info_for_kolb_eccentric
949 :
950 0 : subroutine eval_accreted_material_j(binary_id, ierr)
951 : integer, intent(in) :: binary_id
952 : integer, intent(out) :: ierr
953 : type(binary_info), pointer :: b
954 0 : real(dp) :: qratio, min_r
955 : logical, parameter :: dbg = .false.
956 : include 'formats'
957 :
958 : ierr = 0
959 0 : call binary_ptr(binary_id, b, ierr)
960 0 : if (ierr /= 0) then
961 0 : write(*,*) 'failed in binary_ptr'
962 0 : return
963 : end if
964 0 : qratio = b% m(b% a_i) / b% m(b% d_i)
965 0 : qratio = min(max(qratio,0.0667d0),15d0)
966 0 : min_r = 0.0425d0*b% separation*pow(qratio+qratio*qratio, 0.25d0)
967 :
968 : !TODO: MUST USE EQUATORIAL RADIUS
969 : if (dbg) write(*,*) "radius, impact_radius, separation: ", &
970 : b% r(b% a_i), min_r/rsun, b% separation/rsun
971 0 : if (b% r(b% a_i) < min_r) then
972 0 : b% accretion_mode = 2
973 : b% s_accretor% accreted_material_j = &
974 0 : sqrt(standard_cgrav * b% m(b% a_i) * b% r(b% a_i))
975 : else
976 0 : b% accretion_mode = 1
977 : b% s_accretor% accreted_material_j = &
978 0 : sqrt(standard_cgrav * b% m(b% a_i) * 1.7d0*min_r)
979 : end if
980 : b% acc_am_div_kep_am = b% s_accretor% accreted_material_j / &
981 0 : sqrt(standard_cgrav * b% m(b% a_i) * b% r(b% a_i))
982 :
983 : !TODO: when using wind mass transfer donor star can end up
984 : ! with positive mdot, need to properly set jdot in that case
985 :
986 : end subroutine eval_accreted_material_j
987 :
988 0 : subroutine set_accretion_composition(b, acc_index)
989 : use chem_def, only: chem_isos
990 : type (binary_info), pointer :: b
991 : integer, intent(in) :: acc_index ! index of star that gains mass
992 :
993 : integer :: j
994 :
995 0 : if (acc_index == b% a_i) then
996 : !set accreted material composition
997 0 : b% s_accretor% num_accretion_species = b% s_donor% species
998 :
999 0 : if(b% s_donor% species > size(b% s_accretor% accretion_species_id,dim=1)) then
1000 0 : call mesa_error(__FILE__,__LINE__,'Nuclear network is too large for accretor, increase max_num_accretion_species')
1001 : end if
1002 :
1003 0 : do j = 1, b% s_donor% species
1004 0 : b% s_accretor% accretion_species_id(j) = chem_isos% name(b% s_donor% chem_id(j))
1005 0 : b% s_accretor% accretion_species_xa(j) = b% s_donor% xa_removed(j)
1006 : end do
1007 : else
1008 : ! also for the donor to account for wind mass transfer
1009 0 : b% s_donor% num_accretion_species = b% s_accretor% species
1010 :
1011 0 : if(b% s_accretor% species > size(b% s_donor% accretion_species_id,dim=1)) then
1012 0 : call mesa_error(__FILE__,__LINE__,'Nuclear network is too large for donor, increase max_num_accretion_species')
1013 : end if
1014 :
1015 0 : do j = 1, b% s_accretor% species
1016 0 : b% s_donor% accretion_species_id(j) = chem_isos% name(b% s_accretor% chem_id(j))
1017 0 : b% s_donor% accretion_species_xa(j) = b% s_accretor% xa_removed(j)
1018 : end do
1019 : end if
1020 :
1021 0 : end subroutine set_accretion_composition
1022 :
1023 : end module binary_mdot
|
