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_evolve
22 :
23 : use const_def, only: dp, pi, msun, rsun, secyer, secday, one_third, standard_cgrav
24 : use math_lib
25 : use star_lib
26 : use star_def
27 : use binary_def
28 : use binary_utils, only: eval_rlobe, set_angular_momentum_j, &
29 : set_separation_eccentricity, set_period_eccentricity
30 :
31 : implicit none
32 :
33 : contains
34 :
35 0 : subroutine binarydata_init(b, doing_restart)
36 : use utils_lib, only: is_bad
37 : type (binary_info), pointer :: b
38 : logical, intent(in) :: doing_restart
39 : integer :: finish_step_result
40 : integer :: i, ierr
41 0 : real(dp) :: r_isco, Z1, Z2
42 : include 'formats'
43 :
44 0 : b% evolve_both_stars = b% job% evolve_both_stars
45 0 : b% warn_binary_extra = b% job% warn_binary_extra
46 :
47 : ! Initialize arrays for phase dependent calculations
48 : allocate(b% theta_co(b% anomaly_steps), b% time_co(b% anomaly_steps), &
49 0 : b% mdot_donor_theta(b% anomaly_steps))
50 : allocate(b% edot_theta(b% anomaly_steps), b% e1(b% anomaly_steps), &
51 0 : b% e2(b% anomaly_steps), b% e3(b% anomaly_steps))
52 :
53 0 : if (.not. doing_restart) then
54 0 : if (.not. b% evolve_both_stars) then
55 0 : b% point_mass_i = 2
56 : else if ((b% job% change_initial_model_twins_flag .or. b% job% change_model_twins_flag) &
57 0 : .and. b% job% new_model_twins_flag) then
58 0 : b% point_mass_i = 2
59 : else
60 0 : b% point_mass_i = 0
61 : end if
62 :
63 0 : b% doing_first_model_of_run = .true.
64 0 : b% open_new_history_file = .true.
65 :
66 0 : b% s_donor => b% s1
67 0 : b% s_accretor => b% s2
68 0 : b% d_i = 1
69 0 : b% a_i = 2
70 :
71 0 : b% max_timestep = 1d99
72 0 : b% change_factor = b% initial_change_factor
73 :
74 0 : if (b% point_mass_i /= 1) then
75 0 : initial_mass(1) = b% s1% mstar / Msun
76 : else
77 0 : initial_mass(1) = b% m1
78 : end if
79 0 : if (b% point_mass_i /= 2) then
80 0 : initial_mass(2) = b% s2% mstar / Msun
81 : else
82 0 : if (.not. b% model_twins_flag) then
83 0 : initial_mass(2) = b% m2
84 : else
85 0 : initial_mass(2) = initial_mass(1)
86 : end if
87 : end if
88 :
89 0 : b% m(1) = initial_mass(1)*Msun
90 0 : b% m(2) = initial_mass(2)*Msun
91 0 : if (b% point_mass_i /= 1) then
92 0 : b% r(1) = Rsun*b% s1% photosphere_r
93 : else
94 0 : b% r(1) = 0
95 : end if
96 0 : if (b% point_mass_i /= 2) then
97 0 : b% r(2) = Rsun*b% s2% photosphere_r
98 : else
99 0 : if (.not. b% model_twins_flag) then
100 0 : b% r(2) = 0
101 : else
102 0 : b% r(2) = b% r(1)
103 : end if
104 : end if
105 :
106 0 : if (b% initial_period_in_days <= 0) then ! calculate from initial_separation_in_Rsuns
107 : call set_separation_eccentricity(b% binary_id, &
108 0 : b% initial_separation_in_Rsuns*Rsun, b% initial_eccentricity, ierr)
109 0 : if (ierr /= 0) then
110 0 : return
111 : end if
112 : else
113 : call set_period_eccentricity(b% binary_id, &
114 0 : b% initial_period_in_days*secday, b% initial_eccentricity, ierr)
115 0 : if (ierr /= 0) then
116 : return
117 : end if
118 : end if
119 :
120 : ! Set all parameters necessary for integration over the binary orbit
121 : ! 1) true anomaly = polar angle from periastron 0 -> 2pi
122 0 : do i = 1,b% anomaly_steps
123 0 : b% theta_co(i) = (i-1) * (2 * pi) / b% anomaly_steps
124 : end do
125 : ! 2) time between periastron and polar angle theta 0 -> 1 (fraction of the
126 : ! orbital period)
127 0 : do i = 1,b% anomaly_steps ! time between periastron and polar angle theta
128 : b% time_co(i) = ( 2 * atan( sqrt( (1-b% eccentricity)/(1 + b% eccentricity) ) * &
129 : tan(b% theta_co(i)/2d0) ) - b% eccentricity * &
130 : sqrt(1 - pow2(b% eccentricity)) * sin(b% theta_co(i)) / &
131 0 : (1 + b% eccentricity * cos(b% theta_co(i)) ) ) /2.d0 /pi
132 0 : if (i > b% anomaly_steps/2+1) then
133 0 : b% time_co(i) = b% time_co(i) + b% time_co(b% anomaly_steps/2+1) * 2
134 : end if
135 : end do
136 :
137 0 : if (is_bad(b% rl_relative_gap(1))) call mesa_error(__FILE__,__LINE__,'binarydata_init')
138 0 : if (is_bad(b% rl_relative_gap(2))) call mesa_error(__FILE__,__LINE__,'binarydata_init')
139 0 : b% using_jdot_mb(1) = .false.
140 0 : b% using_jdot_mb(2) = .false.
141 :
142 0 : if (b% point_mass_i /= 0) then
143 : ! this part is only relevant for BH accretors
144 0 : if (b% initial_bh_spin < 0d0) then
145 0 : b% initial_bh_spin = 0d0
146 0 : write(*,*) "initial_bh_spin is smaller than zero. It has been set to zero."
147 0 : else if (b% initial_bh_spin > 1d0) then
148 0 : b% initial_bh_spin = 1d0
149 0 : write(*,*) "initial_bh_spin is larger than one. It has been set to one."
150 : end if
151 : ! compute isco radius from eq. 2.21 of Bardeen et al. (1972), ApJ, 178, 347
152 : Z1 = 1d0 + pow(1d0 - pow2(b% initial_bh_spin),one_third) &
153 0 : * (pow(1d0 + b% initial_bh_spin,one_third) + pow(1d0 - b% initial_bh_spin,one_third))
154 0 : Z2 = sqrt(3d0*pow2(b% initial_bh_spin) + pow2(Z1))
155 0 : r_isco = 3d0 + Z2 - sqrt((3d0 - Z1)*(3d0 + Z1 + 2d0*Z2))
156 : ! compute equivalent mass at zero spin from eq. (3+1/2) (ie. the equation between (3) and (4))
157 : ! of Bardeen (1970), Nature, 226, 65, taking values with subscript zero to correspond to
158 : ! zero spin (r_isco = 6).
159 0 : b% eq_initial_bh_mass = b% m(b% point_mass_i) * sqrt(r_isco/6d0)
160 : end if
161 :
162 0 : write(*,'(A)')
163 0 : write(*,1) 'm2', b% m2
164 0 : write(*,1) 'm1', b% m1
165 0 : write(*,1) 'initial_period_in_days', b% initial_period_in_days
166 0 : write(*,1) 'initial_separation_in_Rsun', b% separation/Rsun
167 0 : write(*,1) 'jdot_multiplier', b% jdot_multiplier
168 0 : write(*,1) 'fr', b% fr
169 0 : write(*,'(A)')
170 :
171 0 : min_binary_period = b% period
172 0 : b% min_binary_separation = b% separation
173 0 : initial_binary_period = b% period
174 :
175 0 : b% ixtra(:) = 0
176 0 : b% xtra(:) = 0d0
177 0 : b% lxtra(:) = .false.
178 :
179 0 : b% CE_num1 = 0
180 0 : b% CE_num2 = 0
181 0 : b% CE_lambda1 = 0d0
182 0 : b% CE_lambda2 = 0d0
183 0 : b% CE_Ebind1 = 0d0
184 0 : b% CE_Ebind2 = 0d0
185 0 : b% mtransfer_rate = 0d0
186 :
187 0 : b% num_tries = 0
188 :
189 0 : finish_step_result = binary_finish_step(b)
190 : else
191 0 : if (b% d_i == b% point_mass_i) then
192 0 : write(*,*) "WARNING: restart has donor star set as point mass"
193 0 : write(*,*) "switching donor"
194 0 : b% d_i = b% a_i
195 0 : b% a_i = b% point_mass_i
196 : end if
197 0 : if (b% d_i == 1) then
198 0 : b% s_donor => b% s1
199 0 : b% s_accretor => b% s2
200 : else
201 0 : b% s_donor => b% s2
202 0 : b% s_accretor => b% s1
203 : end if
204 : end if
205 :
206 : end subroutine binarydata_init
207 :
208 0 : subroutine set_donor_star(b)
209 : type (binary_info), pointer :: b
210 : logical :: switch_donor
211 0 : real(dp) :: mdot_hi_temp
212 : include 'formats'
213 :
214 0 : switch_donor = .false.
215 :
216 0 : if (b% keep_donor_fixed .and. b% mdot_scheme /= "contact") return
217 :
218 : if (b% mdot_scheme == "roche_lobe" .and. &
219 0 : abs(b% mtransfer_rate/(Msun/secyer)) < b% mdot_limit_donor_switch .and. &
220 : b% rl_relative_gap_old(b% a_i) > b% rl_relative_gap_old(b% d_i)) then
221 : switch_donor = .true.
222 0 : else if (b% mtransfer_rate > 0d0) then
223 0 : switch_donor = .true.
224 0 : b% mtransfer_rate = - b% mtransfer_rate
225 0 : mdot_hi_temp = b% mdot_hi
226 0 : b% mdot_hi = - b% mdot_lo
227 0 : b% mdot_lo = - mdot_hi_temp
228 0 : if (.not. b% have_mdot_lo) then
229 0 : b% have_mdot_hi = .false.
230 : end if
231 0 : b% have_mdot_lo = .true.
232 0 : b% fixed_delta_mdot = b% fixed_delta_mdot / 2.0d0
233 : else if (b% mdot_scheme == "contact" .and. &
234 : b% rl_relative_gap_old(b% a_i) > b% rl_relative_gap_old(b% d_i) .and. &
235 : b% rl_relative_gap_old(b% a_i) < - b% implicit_scheme_tolerance .and. &
236 0 : b% rl_relative_gap_old(b% d_i) < - b% implicit_scheme_tolerance .and. &
237 : abs(b% mtransfer_rate/(Msun/secyer)) < b% mdot_limit_donor_switch) then
238 : switch_donor = .true.
239 : end if
240 :
241 : if (switch_donor) then
242 0 : if (b% report_rlo_solver_progress) write(*,*) "switching donor"
243 0 : if (b% d_i == 2) then
244 0 : b% d_i = 1
245 0 : b% a_i = 2
246 0 : b% s_donor => b% s1
247 0 : b% s_accretor => b% s2
248 : else
249 0 : b% d_i = 2
250 0 : b% a_i = 1
251 0 : b% s_donor => b% s2
252 0 : b% s_accretor => b% s1
253 : end if
254 : end if
255 : end subroutine set_donor_star
256 :
257 0 : integer function binary_evolve_step(b)
258 : use utils_lib, only: is_bad
259 : use binary_jdot, only: get_jdot
260 : use binary_edot, only: get_edot
261 : type(binary_info), pointer :: b
262 : integer :: i
263 :
264 : include 'formats'
265 :
266 : ! store the final mdots used for each star
267 : ! for a point mass accretor this is already set in the subroutine adjust_mdots of binary_mdot
268 0 : b% component_mdot(b% d_i) = b% s_donor% mstar_dot
269 0 : if (b% point_mass_i == 0) then
270 0 : b% component_mdot(b% a_i) = b% s_accretor% mstar_dot
271 0 : else if (b% model_twins_flag) then
272 0 : b% component_mdot(b% a_i) = b% component_mdot(b% d_i)
273 : end if
274 :
275 0 : b% m(b% d_i) = b% s_donor% mstar
276 0 : b% time_step = b% s_donor% time_step
277 0 : if (b% point_mass_i == 0) then
278 0 : b% m(b% a_i) = b% s_accretor% mstar
279 0 : else if (.not. b% model_twins_flag) then
280 0 : b% m(b% a_i) = b% m(b% a_i) + b% component_mdot(b% a_i)*b% s_donor% dt
281 : else
282 0 : b% m(b% a_i) = b% m(b% d_i)
283 : end if
284 :
285 0 : if (b% point_mass_i /= 1) then
286 0 : b% r(1) = Rsun*b% s1% photosphere_r
287 : else
288 0 : b% r(1) = 0
289 : end if
290 0 : if (b% point_mass_i /= 2) then
291 0 : b% r(2) = Rsun*b% s2% photosphere_r
292 0 : else if (.not. b% model_twins_flag) then
293 0 : b% r(2) = 0
294 : else
295 0 : b% r(2) = b% r(1)
296 : end if
297 :
298 : ! solve the winds in the system for jdot calculation,
299 : ! these don't include mass lost due to mass_transfer_efficiency < 1.0
300 : ! Since s% mstar_dot is not just mass loss, but includes the contribution from
301 : ! RLO mass transfer and wind mass transfer from the other component, these
302 : ! need to be removed to get the actual wind. Also, the fraction of the wind
303 : ! that is fransferred to the other component does not leave the system, and
304 : ! needs to be removed as well.
305 : b% mdot_system_wind(b% d_i) = b% s_donor% mstar_dot - b% mtransfer_rate &
306 0 : + b% mdot_wind_transfer(b% a_i) - b% mdot_wind_transfer(b% d_i)
307 0 : if (b% point_mass_i == 0) then
308 : b% mdot_system_wind(b% a_i) = b% s_accretor% mstar_dot &
309 : + b% mtransfer_rate * b% fixed_xfer_fraction + b% mdot_wind_transfer(b% d_i) &
310 0 : - b% mdot_wind_transfer(b% a_i)
311 : else
312 0 : if (.not. b% model_twins_flag) then
313 0 : b% mdot_system_wind(b% a_i) = 0.0d0
314 : else
315 0 : b% mdot_system_wind(b% a_i) = b% mdot_system_wind(b% d_i)
316 : end if
317 : end if
318 :
319 : ! get jdot and update orbital J
320 0 : b% jdot = get_jdot(b)
321 0 : b% angular_momentum_j = b% angular_momentum_j + b% jdot*b% time_step*secyer
322 :
323 0 : if (b% angular_momentum_j <= 0) then
324 0 : write(*,*) 'Retry due to negative angular_momentum_j', b% angular_momentum_j
325 0 : b% have_to_reduce_timestep_due_to_j = .true.
326 0 : binary_evolve_step = retry
327 0 : return
328 : end if
329 :
330 : ! update the eccentricity (ignore in first step)
331 0 : if (.not. b% doing_first_model_of_run) then
332 0 : b% eccentricity = b% eccentricity + get_edot(b) *b% time_step*secyer
333 0 : if (b% eccentricity < b% min_eccentricity) b% eccentricity = b% min_eccentricity
334 0 : if (b% eccentricity > b% max_eccentricity) b% eccentricity = b% max_eccentricity
335 : else
336 0 : b% edot_tidal = 0d0
337 0 : b% edot_enhance = 0d0
338 0 : b% extra_edot = 0d0
339 0 : b% edot = 0d0
340 : end if
341 :
342 : !use new eccentricity to calculate new time coordinate
343 0 : do i = 1,b% anomaly_steps ! time between periastron and polar angle theta
344 : b% time_co(i) = ( 2 * atan( sqrt( (1-b% eccentricity)/(1 + b% eccentricity) ) * &
345 : tan(b% theta_co(i)/2d0) ) - b% eccentricity * &
346 : sqrt(1 - pow2(b% eccentricity)) * sin(b% theta_co(i)) / &
347 0 : (1 + b% eccentricity * cos(b% theta_co(i)) ) ) /2.d0 /pi
348 0 : if (i > b% anomaly_steps/2+1) then
349 0 : b% time_co(i) = b% time_co(i) + b% time_co(b% anomaly_steps/2+1) * 2
350 : end if
351 : end do
352 :
353 : ! use the new j to calculate new separation
354 : b% separation = (pow2(b% angular_momentum_j/(b% m(1)*b% m(2)))) *&
355 0 : (b% m(1)+b% m(2)) / standard_cgrav * 1 / (1 - pow2(b% eccentricity))
356 0 : if (b% separation < b% min_binary_separation) &
357 0 : b% min_binary_separation = b% separation
358 :
359 : b% period = 2*pi*sqrt(pow3(b% separation)/&
360 0 : (standard_cgrav*(b% m(1)+b% m(2))))
361 0 : if (b% period < min_binary_period) min_binary_period = b% period
362 :
363 : ! use the new separation to calculate the new roche lobe radius
364 :
365 0 : b% rl(1) = eval_rlobe(b% m(1), b% m(2), b% separation)
366 0 : b% rl(2) = eval_rlobe(b% m(2), b% m(1), b% separation)
367 : b% rl_relative_gap(1) = (b% r(1) - b% rl(1) * (1 - b% eccentricity) ) / &
368 0 : b% rl(1) / (1 - b% eccentricity) ! gap < 0 means out of contact
369 : b% rl_relative_gap(2) = (b% r(2) - b% rl(2) * (1 - b% eccentricity) ) / &
370 0 : b% rl(2) / (1 - b% eccentricity) ! gap < 0 means out of contact
371 :
372 0 : if (is_bad(b% rl_relative_gap(1)) .or. is_bad(b% rl_relative_gap(2))) then
373 0 : write(*,*) "rl_relative_gap for each component", b% rl_relative_gap(1), b% rl_relative_gap(2)
374 0 : write(*,*) "rl for each component", b% rl(1), b% rl(2)
375 0 : write(*,*) "r for each component", b% r(1), b% r(2)
376 0 : write(*,*) "m for each component", b% m(1), b% m(2)
377 0 : write(*,*) "separation, angular momentum", b% separation, b% angular_momentum_j
378 0 : write(*,*) "jdot, jdot_mb, jdot_gr, jdot_ml:", b% jdot, b% jdot_mb, b% jdot_gr, b% jdot_ml
379 0 : write(*,*) "jdot_ls, jdot_missing_wind, extra_jdot:", b% jdot_ls, b% jdot_missing_wind, b% extra_jdot
380 0 : write(*,*) 'error solving rl_rel_gap'
381 0 : binary_evolve_step = retry
382 0 : return
383 : end if
384 :
385 0 : b% model_number = b% model_number + 1
386 0 : b% binary_age = b% binary_age + b% time_step
387 :
388 0 : binary_evolve_step = keep_going
389 :
390 0 : end function binary_evolve_step
391 :
392 0 : integer function binary_check_model(b)
393 0 : use binary_mdot, only: rlo_mdot, check_implicit_rlo
394 : use binary_irradiation
395 : type (binary_info), pointer :: b
396 :
397 : integer :: ierr, id
398 : logical :: implicit_rlo
399 0 : real(dp) :: new_mdot, q
400 :
401 :
402 : include 'formats'
403 :
404 0 : binary_check_model = retry
405 0 : ierr = 0
406 :
407 0 : implicit_rlo = (b% max_tries_to_achieve > 0 .and. b% implicit_scheme_tolerance > 0d0)
408 :
409 0 : binary_check_model = keep_going
410 :
411 0 : if (.not. b% ignore_rlof_flag) then
412 0 : if (implicit_rlo) then ! check agreement between new r and new rl
413 0 : if (.not. b% use_other_check_implicit_rlo) then
414 0 : binary_check_model = check_implicit_rlo(b% binary_id, new_mdot)
415 : else
416 0 : binary_check_model = b% other_check_implicit_rlo(b% binary_id, new_mdot)
417 : end if
418 0 : if (binary_check_model == keep_going) then
419 0 : b% donor_started_implicit_wind = .false.
420 : end if
421 : b% donor_started_implicit_wind = b% donor_started_implicit_wind .or. &
422 0 : b% s_donor% was_in_implicit_wind_limit
423 : else
424 0 : if (.not. b% use_other_rlo_mdot) then
425 0 : call rlo_mdot(b% binary_id, new_mdot, ierr) ! grams per second
426 0 : if (ierr /= 0) then
427 0 : write(*,*) 'failed in rlo_mdot'
428 0 : binary_check_model = retry
429 0 : return
430 : end if
431 : else
432 0 : call b% other_rlo_mdot(b% binary_id, new_mdot, ierr)
433 0 : if (ierr /= 0) then
434 0 : write(*,*) 'failed in other rlo_mdot'
435 0 : binary_check_model = retry
436 0 : return
437 : end if
438 : end if
439 0 : if (new_mdot > 0) then
440 0 : new_mdot = 0.0d0
441 0 : write(*,*) "WARNING: explicit computation of mass transfer results in accreting donor"
442 0 : write(*,*) "Not transferring mass"
443 : end if
444 : ! smooth out the changes in mdot
445 0 : new_mdot = b% cur_mdot_frac*b% mtransfer_rate + (1-b% cur_mdot_frac)*new_mdot
446 0 : if (-new_mdot/(Msun/secyer) > b% max_explicit_abs_mdot) new_mdot = -b% max_explicit_abs_mdot*Msun/secyer
447 : end if
448 0 : b% mtransfer_rate = new_mdot
449 : else
450 0 : b% mtransfer_rate = 0d0
451 : end if
452 0 : call adjust_irradiation(b)
453 :
454 0 : if (.not. b% CE_flag) then
455 : if ((b% point_mass_i == 0 .or. b% model_twins_flag) &
456 0 : .and. b% rl_relative_gap(b% a_i) >= 0.0d0) then
457 0 : if (b% rl_relative_gap(b% a_i) >= b% accretor_overflow_terminate) then
458 0 : binary_check_model = terminate
459 0 : b% s_donor% termination_code = t_xtra1
460 : termination_code_str(t_xtra1) = &
461 0 : "Terminate because accretor (r-rl)/rl > accretor_overflow_terminate"
462 : end if
463 : end if
464 : if (b% doing_first_model_of_run .and. b% terminate_if_initial_overflow &
465 0 : .and. (.not. b% ignore_rlof_flag .or. b% model_twins_flag)) then
466 : if (b% rl_relative_gap(b% d_i) >= 0.0d0 &
467 0 : .or. (b% point_mass_i == 0 .and. b% rl_relative_gap(b% a_i) >= 0.0d0)) then
468 0 : binary_check_model = terminate
469 0 : b% s_donor% termination_code = t_xtra1
470 : termination_code_str(t_xtra1) = &
471 0 : "Terminate because of overflowing initial model"
472 : end if
473 : end if
474 : if ((b% point_mass_i == 0 .or. b% model_twins_flag) &
475 0 : .and. b% terminate_if_L2_overflow) then
476 0 : if (b% m(1) > b% m(2)) then
477 0 : q = b% m(2) / b% m(1)
478 0 : id = 2
479 : else
480 0 : q = b% m(1) / b% m(2)
481 0 : id = 1
482 : end if
483 0 : if (b% rl_relative_gap(id) > 0.29858997d0*atan(1.83530121d0*pow(q,0.39661426d0))) then
484 0 : binary_check_model = terminate
485 0 : b% s_donor% termination_code = t_xtra1
486 : termination_code_str(t_xtra1) = &
487 0 : "Terminate because of L2 overflow"
488 : end if
489 : end if
490 :
491 0 : if (b% mtransfer_rate == -b% max_implicit_abs_mdot*Msun/secyer .and. &
492 : b% CE_begin_at_max_implicit_abs_mdot) then
493 0 : b% CE_flag = .true.
494 0 : b% CE_init = .false.
495 0 : binary_check_model = keep_going
496 : end if
497 : end if
498 :
499 0 : end function binary_check_model
500 :
501 0 : integer function binary_finish_step(b)
502 : type (binary_info), pointer :: b
503 :
504 0 : binary_finish_step = keep_going
505 : ! update change factor in case mtransfer_rate has changed
506 0 : if(.not. b% doing_first_model_of_run) then
507 : if(b% mtransfer_rate_old /= b% mtransfer_rate .and. &
508 0 : b% mtransfer_rate /= 0 .and. b% mtransfer_rate_old /= 0) then
509 0 : if(b% mtransfer_rate < b% mtransfer_rate_old) then
510 : b% change_factor = b% change_factor*(1d0-b% implicit_lambda) + b% implicit_lambda* &
511 0 : (1+b% change_factor_fraction*(b% mtransfer_rate/b% mtransfer_rate_old-1))
512 : else
513 : b% change_factor = b% change_factor*(1d0-b% implicit_lambda) + b% implicit_lambda* &
514 0 : (1+b% change_factor_fraction*(b% mtransfer_rate_old/b% mtransfer_rate-1))
515 : end if
516 0 : if(b% change_factor > b% max_change_factor) b% change_factor = b% max_change_factor
517 0 : if(b% change_factor < b% min_change_factor) b% change_factor = b% min_change_factor
518 : end if
519 : end if
520 :
521 : ! store all variables into "old"
522 :
523 0 : b% model_number_old = b% model_number
524 0 : b% binary_age_old = b% binary_age
525 0 : b% mtransfer_rate_old = b% mtransfer_rate
526 0 : b% angular_momentum_j_old = b% angular_momentum_j
527 0 : b% separation_old = b% separation
528 0 : b% eccentricity_old = b% eccentricity
529 0 : b% dt_old = b% dt
530 0 : b% env_old(1) = b% env(1)
531 0 : b% env_old(2) = b% env(2)
532 0 : b% period_old = b% period
533 0 : b% rl_relative_gap_old(1) = b% rl_relative_gap(1)
534 0 : b% rl_relative_gap_old(2) = b% rl_relative_gap(2)
535 0 : b% r_old(1) = b% r(1)
536 0 : b% r_old(2) = b% r(2)
537 0 : b% rl_old(1) = b% rl(1)
538 0 : b% rl_old(2) = b% rl(2)
539 0 : b% m_old(1) = b% m(1)
540 0 : b% m_old(2) = b% m(2)
541 0 : b% using_jdot_mb_old = b% using_jdot_mb
542 0 : b% max_timestep_old = b% max_timestep
543 0 : b% change_factor_old = b% change_factor
544 :
545 0 : b% d_i_old = b% d_i
546 0 : b% a_i_old = b% a_i
547 0 : b% point_mass_i_old = b% point_mass_i
548 :
549 0 : b% ignore_rlof_flag_old = b% ignore_rlof_flag
550 0 : b% model_twins_flag_old = b% model_twins_flag
551 :
552 0 : b% CE_flag_old = b% CE_flag
553 0 : b% CE_init_old = b% CE_init
554 :
555 0 : b% CE_num1_old = b% CE_num1
556 0 : b% CE_num2_old = b% CE_num2
557 0 : b% CE_lambda1_old = b% CE_lambda1
558 0 : b% CE_lambda2_old = b% CE_lambda2
559 0 : b% CE_Ebind1_old = b% CE_Ebind1
560 0 : b% CE_Ebind2_old = b% CE_Ebind2
561 0 : b% CE_years_detached_old = b% CE_years_detached
562 :
563 0 : b% dt_why_reason_old = b% dt_why_reason
564 :
565 : !set all xtra variables
566 0 : b% ixtra_old(:) = b% ixtra(:)
567 0 : b% xtra_old(:) = b% xtra(:)
568 0 : b% lxtra_old(:) = b% lxtra(:)
569 :
570 0 : end function binary_finish_step
571 :
572 0 : integer function binary_prepare_to_redo(b)
573 : type (binary_info), pointer :: b
574 :
575 0 : binary_prepare_to_redo = redo
576 0 : call binary_set_current_to_old(b)
577 :
578 0 : end function binary_prepare_to_redo
579 :
580 0 : integer function binary_prepare_to_retry(b)
581 : type (binary_info), pointer :: b
582 :
583 0 : b% num_tries = 0
584 : ! this call takes care of restoring variables
585 0 : call binary_set_current_to_old(b)
586 0 : binary_prepare_to_retry = retry
587 :
588 0 : end function binary_prepare_to_retry
589 :
590 0 : subroutine binary_set_current_to_old(b)
591 : type (binary_info), pointer :: b
592 : ! restore variables
593 0 : b% model_number = b% model_number_old
594 0 : b% binary_age = b% binary_age_old
595 : ! do not restore mtransfer_rate during implicit rlo
596 0 : if (b% num_tries == 0) b% mtransfer_rate = b% mtransfer_rate_old
597 0 : b% angular_momentum_j = b% angular_momentum_j_old
598 0 : b% separation = b% separation_old
599 0 : b% eccentricity = b% eccentricity_old
600 0 : b% dt = b% dt_old
601 0 : b% env(1) = b% env_old(1)
602 0 : b% env(2) = b% env_old(2)
603 0 : b% period = b% period_old
604 0 : b% rl_relative_gap(1) = b% rl_relative_gap_old(1)
605 0 : b% rl_relative_gap(2) = b% rl_relative_gap_old(2)
606 0 : b% r(1) = b% r_old(1)
607 0 : b% r(2) = b% r_old(2)
608 0 : b% rl(1) = b% rl_old(1)
609 0 : b% rl(2) = b% rl_old(2)
610 0 : b% m(1) = b% m_old(1)
611 0 : b% m(2) = b% m_old(2)
612 0 : b% using_jdot_mb = b% using_jdot_mb_old
613 0 : b% max_timestep = b% max_timestep_old
614 :
615 : ! the following need to be kept constant during implicit mdot
616 0 : if (b% num_tries == 0) then
617 0 : b% change_factor = b% change_factor_old
618 0 : b% d_i = b% d_i_old
619 0 : b% a_i = b% a_i_old
620 0 : b% point_mass_i = b% point_mass_i_old
621 0 : b% ignore_rlof_flag = b% ignore_rlof_flag_old
622 0 : b% model_twins_flag = b% model_twins_flag_old
623 0 : b% CE_flag = b% CE_flag_old
624 0 : b% CE_init = b% CE_init_old
625 : end if
626 :
627 0 : b% CE_num1 = b% CE_num1_old
628 0 : b% CE_num2 = b% CE_num2_old
629 0 : b% CE_lambda1 = b% CE_lambda1_old
630 0 : b% CE_lambda2 = b% CE_lambda2_old
631 0 : b% CE_Ebind1 = b% CE_Ebind1_old
632 0 : b% CE_Ebind2 = b% CE_Ebind2_old
633 0 : b% CE_years_detached = b% CE_years_detached_old
634 :
635 0 : b% dt_why_reason = b% dt_why_reason_old
636 :
637 : !set all xtra variables
638 0 : b% ixtra(:) = b% ixtra_old(:)
639 0 : b% xtra(:) = b% xtra_old(:)
640 0 : b% lxtra(:) = b% lxtra_old(:)
641 0 : end subroutine binary_set_current_to_old
642 :
643 0 : integer function binary_after_evolve(b)
644 : type (binary_info), pointer :: b
645 0 : binary_after_evolve = keep_going
646 :
647 : !take care of deallocating binary arrays here
648 0 : if (associated(b% theta_co)) then
649 0 : deallocate(b% theta_co)
650 0 : nullify(b% theta_co)
651 : end if
652 0 : if (associated(b% time_co)) then
653 0 : deallocate(b% time_co)
654 0 : nullify(b% time_co)
655 : end if
656 0 : if (associated(b% mdot_donor_theta)) then
657 0 : deallocate(b% mdot_donor_theta)
658 0 : nullify(b% mdot_donor_theta)
659 : end if
660 0 : if (associated(b% edot_theta)) then
661 0 : deallocate(b% edot_theta)
662 0 : nullify(b% edot_theta)
663 : end if
664 0 : if (associated(b% e1)) then
665 0 : deallocate(b% e1)
666 0 : nullify(b% e1)
667 : end if
668 0 : if (associated(b% e2)) then
669 0 : deallocate(b% e2)
670 0 : nullify(b% e2)
671 : end if
672 0 : if (associated(b% CE_m)) then
673 0 : deallocate(b% CE_m)
674 0 : nullify(b% CE_m)
675 : end if
676 0 : if (associated(b% CE_entropy)) then
677 0 : deallocate(b% CE_entropy)
678 0 : nullify(b% CE_entropy)
679 : end if
680 0 : if (associated(b% CE_U_in)) then
681 0 : deallocate(b% CE_U_in)
682 0 : nullify(b% CE_U_in)
683 : end if
684 0 : if (associated(b% CE_U_out)) then
685 0 : deallocate(b% CE_U_out)
686 0 : nullify(b% CE_U_out)
687 : end if
688 0 : if (associated(b% CE_Omega_in)) then
689 0 : deallocate(b% CE_Omega_in)
690 0 : nullify(b% CE_Omega_in)
691 : end if
692 0 : if (associated(b% CE_Omega_out)) then
693 0 : deallocate(b% CE_Omega_out)
694 0 : nullify(b% CE_Omega_out)
695 : end if
696 0 : end function binary_after_evolve
697 :
698 : end module binary_evolve
|
