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_jdot
22 :
23 : use const_def, only: dp, pi, clight, standard_cgrav, rsun, convective_mixing
24 : use star_lib
25 : use star_def
26 : use math_lib
27 : use binary_def
28 :
29 : implicit none
30 :
31 : contains
32 :
33 0 : real(dp) function get_jdot(b)
34 : type (binary_info), pointer :: b
35 :
36 : integer :: ierr
37 :
38 : ! calculate jdot from gravitational wave radiation
39 0 : if (.not. b% do_jdot_gr) then
40 0 : b% jdot_gr = 0d0
41 0 : else if (.not. b% use_other_jdot_gr) then
42 0 : call default_jdot_gr(b% binary_id, ierr)
43 : else
44 0 : call b% other_jdot_gr(b% binary_id, ierr)
45 : end if
46 :
47 : ! calculate jdot for mass ejected from system
48 0 : if (.not. b% do_jdot_ml) then
49 0 : b% jdot_ml = 0d0
50 0 : else if (.not. b% use_other_jdot_ml) then
51 0 : call default_jdot_ml(b% binary_id, ierr)
52 : else
53 0 : call b% other_jdot_ml(b% binary_id, ierr)
54 : end if
55 :
56 : ! solve jdot due to L-S coupling
57 0 : if (.not. b% do_jdot_ls) then
58 0 : b% jdot_ls = 0d0
59 0 : else if (.not. b% use_other_jdot_ls) then
60 0 : call default_jdot_ls(b% binary_id, ierr)
61 : else
62 0 : call b% other_jdot_ls(b% binary_id, ierr)
63 : end if
64 :
65 : ! solve jdot due to "missing wind" (see binary_controls.defaults)
66 0 : if (.not. b% do_jdot_missing_wind) then
67 0 : b% jdot_missing_wind = 0d0
68 0 : else if (.not. b% use_other_jdot_missing_wind) then
69 0 : call default_jdot_missing_wind(b% binary_id, ierr)
70 : else
71 0 : call b% other_jdot_missing_wind(b% binary_id, ierr)
72 : end if
73 :
74 : ! calculate jdot from magnetic braking
75 0 : if (.not. b% do_jdot_mb) then
76 0 : b% jdot_mb = 0d0
77 0 : else if (.not. b% use_other_jdot_mb) then
78 0 : call default_jdot_mb(b% binary_id, ierr)
79 : else
80 0 : call b% other_jdot_mb(b% binary_id, ierr)
81 : end if
82 :
83 : ! calculate extra jdot
84 0 : if (.not. b% use_other_extra_jdot) then
85 0 : b% extra_jdot = 0
86 : else
87 0 : call b% other_extra_jdot(b% binary_id, ierr)
88 : end if
89 :
90 : get_jdot = (b% jdot_mb + b% jdot_gr + b% jdot_ml + b% jdot_missing_wind + &
91 0 : b% extra_jdot) * b% jdot_multiplier + b% jdot_ls
92 :
93 0 : end function get_jdot
94 :
95 0 : subroutine default_jdot_gr(binary_id, ierr)
96 : integer, intent(in) :: binary_id
97 : integer, intent(out) :: ierr
98 : type (binary_info), pointer :: b
99 0 : real(dp) :: bs4, clight5, cgrav3
100 : ierr = 0
101 0 : call binary_ptr(binary_id, b, ierr)
102 0 : if (ierr /= 0) then
103 0 : write(*,*) 'failed in binary_ptr'
104 0 : return
105 : end if
106 0 : bs4 = pow4(b% separation)
107 0 : clight5 = pow5(clight)
108 0 : cgrav3 = standard_cgrav*standard_cgrav*standard_cgrav
109 : b% jdot_gr = -32d0 * cgrav3 * b% m(b% a_i) * b% m(b% d_i) * (b% m(b% a_i) + b% m(b% d_i)) / &
110 0 : (5d0 * clight5 * bs4) * b% angular_momentum_j
111 : end subroutine default_jdot_gr
112 :
113 0 : subroutine default_jdot_ml(binary_id, ierr)
114 : integer, intent(in) :: binary_id
115 : integer, intent(out) :: ierr
116 : type (binary_info), pointer :: b
117 : ierr = 0
118 0 : call binary_ptr(binary_id, b, ierr)
119 0 : if (ierr /= 0) then
120 0 : write(*,*) 'failed in binary_ptr'
121 0 : return
122 : end if
123 : !mass lost from vicinity of donor
124 : b% jdot_ml = (b% mdot_system_transfer(b% d_i) + b% mdot_system_wind(b% d_i))*&
125 : pow2(b% m(b% a_i)/(b% m(b% a_i)+b% m(b% d_i))*b% separation)*2*pi/b% period *&
126 0 : sqrt(1 - pow2(b% eccentricity))
127 : !mass lost from vicinity of accretor
128 : b% jdot_ml = b% jdot_ml + (b% mdot_system_transfer(b% a_i) + b% mdot_system_wind(b% a_i))*&
129 : pow2(b% m(b% d_i)/(b% m(b% a_i)+b% m(b% d_i))*b% separation)*2*pi/b% period *&
130 0 : sqrt(1 - pow2(b% eccentricity))
131 : !mass lost from circumbinary coplanar toroid
132 : b% jdot_ml = b% jdot_ml + b% mdot_system_cct * b% mass_transfer_gamma * &
133 0 : sqrt(standard_cgrav * (b% m(1) + b% m(2)) * b% separation)
134 : end subroutine default_jdot_ml
135 :
136 0 : subroutine default_jdot_ls(binary_id, ierr)
137 : integer, intent(in) :: binary_id
138 : integer, intent(out) :: ierr
139 : type (binary_info), pointer :: b
140 0 : real(dp) :: delta_J
141 : ierr = 0
142 0 : call binary_ptr(binary_id, b, ierr)
143 0 : if (ierr /= 0) then
144 0 : write(*,*) 'failed in binary_ptr'
145 0 : return
146 : end if
147 0 : b% jdot_ls = 0
148 : ! ignore in first step, or if not doing rotation
149 0 : if (b% doing_first_model_of_run) &
150 : return
151 : ! bulk change in spin angular momentum takes tides into account
152 : delta_J = b% s_donor% total_angular_momentum_old - &
153 0 : b% s_donor% total_angular_momentum
154 : ! ignore angular momentum lost through winds
155 0 : if (b% s_donor% mstar_dot < 0) &
156 : delta_J = delta_J - b% s_donor% angular_momentum_removed * &
157 0 : abs(b% mdot_system_wind(b% d_i) / b% s_donor% mstar_dot)
158 0 : b% jdot_ls = b% jdot_ls + delta_J
159 :
160 : ! Repeat for accretor
161 0 : if (b% point_mass_i == 0) then
162 : delta_J = b% s_accretor% total_angular_momentum_old - &
163 0 : b% s_accretor% total_angular_momentum
164 0 : if (b% s_accretor% mstar_dot < 0) then
165 : ! all AM lost from the accretor is lost from the system
166 0 : delta_J = delta_J - b% s_accretor% angular_momentum_removed
167 : end if
168 0 : b% jdot_ls = b% jdot_ls + delta_J
169 0 : else if (b% model_twins_flag) then
170 0 : b% jdot_ls = b% jdot_ls + b% jdot_ls
171 : end if
172 :
173 0 : b% jdot_ls = b% jdot_ls / b% s_donor% dt
174 : end subroutine default_jdot_ls
175 :
176 0 : subroutine default_jdot_missing_wind(binary_id, ierr)
177 : integer, intent(in) :: binary_id
178 : integer, intent(out) :: ierr
179 : type (binary_info), pointer :: b
180 : type (star_info), pointer :: s
181 : ierr = 0
182 0 : call binary_ptr(binary_id, b, ierr)
183 0 : if (ierr /= 0) then
184 0 : write(*,*) 'failed in binary_ptr'
185 0 : return
186 : end if
187 0 : b% jdot_missing_wind = 0
188 0 : if (b% point_mass_i /= 0) return
189 :
190 0 : s => b% s_accretor
191 :
192 0 : if (s% mstar_dot < 0) then
193 0 : b% jdot_missing_wind = b% mtransfer_rate * b% fixed_xfer_fraction
194 : else
195 0 : b% jdot_missing_wind = b% mdot_system_wind(b% a_i)
196 : end if
197 0 : b% jdot_missing_wind = b% jdot_missing_wind * s% j_rot(1)
198 :
199 : end subroutine default_jdot_missing_wind
200 :
201 0 : subroutine check_jdot_mb_conditions(b, s, apply_jdot_mb, qconv_env)
202 : type (binary_info), pointer :: b
203 : type (star_info), pointer :: s
204 : logical, intent(out) :: apply_jdot_mb
205 : real(dp), intent(out) :: qconv_env
206 :
207 : real(dp) :: qconv_core
208 : integer :: k
209 :
210 : include 'formats'
211 :
212 : ! calculate how much of inner region is convective
213 0 : qconv_core = 0d0
214 0 : do k = s% nz, 1, -1
215 0 : if (s% q(k) > b% jdot_mb_qlim_for_check_rad_core .and. &
216 : (qconv_core == 0d0 .or. s% mixing_type(k) /= convective_mixing)) exit
217 0 : if (s% mixing_type(k) == convective_mixing) &
218 0 : qconv_core = qconv_core + s% dq(k)
219 : end do
220 :
221 : ! calculate how much of the envelope
222 0 : qconv_env = 0d0
223 0 : do k = 1, s% nz
224 0 : if (s% q(k) < b% jdot_mb_qlim_for_check_conv_env .and. &
225 : (qconv_env == 0d0 .or. s% mixing_type(k) /= convective_mixing)) exit
226 0 : if (s% mixing_type(k) == convective_mixing) &
227 0 : qconv_env = qconv_env + s% dq(k)
228 : end do
229 :
230 0 : apply_jdot_mb = .true.
231 0 : if (qconv_env < b% jdot_mb_min_qconv_env) then
232 0 : apply_jdot_mb = .false.
233 0 : return
234 : end if
235 :
236 0 : if (qconv_env > b% jdot_mb_max_qconv_env) then
237 0 : apply_jdot_mb = .false.
238 0 : return
239 : end if
240 :
241 0 : if (qconv_core > b% jdot_mb_max_qconv_core) then
242 0 : apply_jdot_mb = .false.
243 0 : return
244 : end if
245 :
246 : end subroutine check_jdot_mb_conditions
247 :
248 0 : subroutine default_jdot_mb(binary_id, ierr)
249 : integer, intent(in) :: binary_id
250 : integer, intent(out) :: ierr
251 : type (binary_info), pointer :: b
252 0 : real(dp) :: rsun4,two_pi_div_p3, qconv_env, jdot_scale
253 : logical :: apply_jdot_mb
254 : ierr = 0
255 0 : call binary_ptr(binary_id, b, ierr)
256 0 : if (ierr /= 0) then
257 0 : write(*,*) 'failed in binary_ptr'
258 0 : return
259 : end if
260 0 : b% jdot_mb = 0
261 0 : rsun4 = pow4(rsun)
262 0 : two_pi_div_p3 = (2.0d0*pi/b% period)*(2.0d0*pi/b% period)*(2.0d0*pi/b% period)
263 :
264 : ! use the formula from rappaport, verbunt, and joss. apj, 275, 713-731. 1983.
265 0 : call check_jdot_mb_conditions(b, b% s_donor, apply_jdot_mb, qconv_env)
266 0 : if (apply_jdot_mb .or. b% keep_mb_on) then
267 0 : jdot_scale = 1d0
268 0 : if (b% jdot_mb_scale_for_low_qconv_env) then
269 : !scale jdot for tiny convective envelope, from Podsiadlowski et al. (2002)
270 : !The Astrophysical Journal, Volume 565, Issue 2, pp. 1107-1133
271 0 : if (qconv_env > b% jdot_mb_mass_frac_for_scale) then
272 : jdot_scale = 1d0
273 : else
274 0 : jdot_scale = exp(-b% jdot_mb_mass_frac_for_scale/max(1d-99,qconv_env)+1)
275 : end if
276 : end if
277 : b% jdot_mb = -3.8d-30*b% m(b% d_i)*rsun4* &
278 : pow(min(b% r(b% d_i),b% rl(b% d_i))/rsun,b% magnetic_braking_gamma)* &
279 0 : two_pi_div_p3*jdot_scale
280 0 : write(*,*) "check jdot_scale", 1, jdot_scale, b% jdot_mb
281 0 : b% using_jdot_mb(b% d_i) = .true.
282 0 : if ((apply_jdot_mb .or. b% keep_mb_on) .and. .not. b% using_jdot_mb_old(b% d_i)) then
283 0 : write(*,*) 'turn on magnetic braking for star ', b% d_i
284 : end if
285 0 : else if (.not. (apply_jdot_mb .or. b% keep_mb_on) .and. b% using_jdot_mb_old(b% d_i)) then
286 : ! required mdot for the implicit scheme may drop drastically,
287 : ! so its necessary to increase change factor to avoid implicit
288 : ! scheme from getting stuck
289 0 : b% change_factor = b% max_change_factor
290 0 : b% using_jdot_mb(b% d_i) = .false.
291 0 : write(*,*) 'turn off magnetic braking for star ', b% d_i
292 : end if
293 :
294 0 : if (b% point_mass_i == 0 .and. b% include_accretor_mb) then
295 0 : call check_jdot_mb_conditions(b, b% s_accretor, apply_jdot_mb, qconv_env)
296 0 : if (apply_jdot_mb .or. b% keep_mb_on) then
297 0 : jdot_scale = 1d0
298 0 : if (b% jdot_mb_scale_for_low_qconv_env) then
299 : !scale jdot for tiny convective envelope, from Podsiadlowski et al. (2002)
300 : !The Astrophysical Journal, Volume 565, Issue 2, pp. 1107-1133
301 0 : if (qconv_env > b% jdot_mb_mass_frac_for_scale) then
302 : jdot_scale = 1d0
303 : else
304 0 : jdot_scale = exp(-b% jdot_mb_mass_frac_for_scale/max(1d-99,qconv_env)+1)
305 : end if
306 : end if
307 : b% jdot_mb = b% jdot_mb - &
308 : 3.8d-30*b% m(b% a_i)*rsun4* &
309 : pow(min(b% r(b% a_i),b% rl(b% a_i))/rsun,b% magnetic_braking_gamma)* &
310 0 : two_pi_div_p3*jdot_scale
311 0 : b% using_jdot_mb(b% a_i) = .true.
312 0 : if ((apply_jdot_mb .or. b% keep_mb_on) .and. .not. b% using_jdot_mb_old(b% a_i)) then
313 0 : write(*,*) 'turn on magnetic braking for star ', b% a_i
314 : end if
315 0 : else if (.not. (apply_jdot_mb .or. b% keep_mb_on) .and. b% using_jdot_mb_old(b% a_i)) then
316 : ! required mdot for the implicit scheme may drop drastically,
317 : ! so its necessary to increase change factor to avoid implicit
318 : ! scheme from getting stuck
319 0 : b% change_factor = b% max_change_factor
320 0 : b% using_jdot_mb(b% a_i) = .false.
321 0 : write(*,*) 'turn off magnetic braking for star ', b% a_i
322 : end if
323 0 : else if (b% model_twins_flag .and. b% include_accretor_mb) then
324 0 : b% jdot_mb = 2*b% jdot_mb
325 0 : b% using_jdot_mb(b% a_i) = .true.
326 : end if
327 :
328 : end subroutine default_jdot_mb
329 :
330 : end module binary_jdot
|
