Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2025 Matthias Fabry & The MESA Team
4 : !
5 : ! MESA is free software; you can use it and/or modify
6 : ! it under the combined terms and restrictions of the MESA MANIFESTO
7 : ! and the GNU General Library Public License as published
8 : ! by the Free Software Foundation; either version 2 of the License,
9 : ! or (at your option) any later version.
10 : !
11 : ! You should have received a copy of the MESA MANIFESTO along with
12 : ! this software; if not, it is available at the mesa website:
13 : ! http://mesa.sourceforge.net/
14 : !
15 : ! MESA is distributed in the hope that it will be useful,
16 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
17 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
18 : ! See the GNU Library General Public License for more details.
19 : !
20 : ! You should have received a copy of the GNU Library General Public License
21 : ! along with this software; if not, write to the Free Software
22 : ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 : !
24 : ! ***********************************************************************
25 :
26 : module binary_roche_deformation
27 :
28 : ! computes the stellar structure correction factors fp/ft and the specific
29 : ! moment of inertia i_rot assuming shellularity in the Roche potential of
30 : ! a binary star. Follows Fabry, Marchant and Sana, 2022, A&A 661, A123.
31 :
32 : use ieee_arithmetic, only: ieee_is_nan
33 : use interp_2d_lib_db
34 : use auto_diff
35 : use star_def
36 : use const_def
37 : use binary_def
38 : use binary_utils, only : eval_rlobe
39 :
40 :
41 : implicit none
42 :
43 : real(dp), parameter :: nudge = 1d-4
44 : real(dp), pointer :: xvals(:), yvals(:), fpfunc1d(:), ftfunc1d(:), &
45 : irotfunc1d(:), otherrfunc1d(:), afunc1d(:), psifunc1d(:)
46 : logical :: inter_ok = .false., dbg = .false.
47 : integer :: num_xpts, num_ypts
48 :
49 : contains
50 :
51 0 : subroutine build_roche_interpolators
52 0 : real(dp) :: xtest, ytest, testval
53 : integer :: ierr
54 : character(len=strlen) :: upstairs
55 : include 'formats'
56 0 : if (.not. inter_ok) then
57 0 : upstairs = trim(mesa_data_dir) // '/roche_data/' ! where fp/ft data lives
58 0 : if (dbg) write(*, 1) 'starting interpolator setup'
59 : call setup_interpolator(trim(upstairs) // 'fp_data.txt', xvals, num_xpts, yvals, &
60 0 : num_ypts, fpfunc1d, ierr)
61 : call setup_interpolator(trim(upstairs) // 'ft_data.txt', xvals, num_xpts, yvals, &
62 0 : num_ypts, ftfunc1d, ierr)
63 : call setup_interpolator(trim(upstairs) // 'irot_data.txt', xvals, num_xpts, &
64 0 : yvals, num_ypts, irotfunc1d, ierr)
65 : call setup_interpolator(trim(upstairs) // 'psi_data.txt', xvals, num_xpts, &
66 0 : yvals, num_ypts, psifunc1d, ierr)
67 :
68 0 : if (dbg) then
69 0 : xtest = -0.5
70 0 : ytest = 1.35
71 0 : write(*, 11) 'grid size', num_xpts, num_ypts
72 0 : write(*, 1) 'setup interpolators succesful,'
73 :
74 : call interp_evbipm_db(xtest, ytest, xvals, num_xpts, yvals, num_ypts,&
75 0 : fpfunc1d, num_xpts, testval, ierr)
76 0 : write(*, 1) 'fp test gave should be close to 0.6', testval
77 : call interp_evbipm_db(xtest, ytest, xvals, num_xpts, yvals, num_ypts,&
78 0 : ftfunc1d, num_xpts, testval, ierr)
79 0 : write(*, 1) 'ft test gave should be close to 0.8', testval
80 : call interp_evbipm_db(xtest, ytest, xvals, num_xpts, yvals, num_ypts,&
81 0 : irotfunc1d, num_xpts, testval, ierr)
82 0 : write(*, 1) 'irot test gave should be close to 1.4', testval
83 : end if
84 0 : inter_ok = .true.
85 : end if
86 :
87 : contains
88 : ! interpolator data reading
89 0 : subroutine setup_interpolator(filename, xs, num_xs, ys, num_ys, func1d, ierr)
90 : integer, intent(out) :: ierr, num_xs, num_ys
91 : character(len = *) :: filename
92 : integer :: k, iounit
93 : real(dp), pointer, intent(out) :: xs(:), ys(:), func1d(:)
94 0 : real(dp), pointer :: func(:,:,:)
95 :
96 : include 'formats'
97 :
98 0 : if (dbg) then
99 0 : write(*, '(a100)') 'loading ' // filename
100 : end if
101 : ! open data to interpolate
102 : open(newunit = iounit, file = trim(filename), status = 'old', action = 'read',&
103 0 : iostat = ierr)
104 :
105 0 : read(iounit, *, iostat = ierr) num_xs
106 0 : allocate(xs(num_xs))
107 0 : do k = 1, num_xs
108 0 : read(iounit, *, iostat = ierr) xs(k)
109 : end do
110 :
111 0 : read(iounit, *, iostat = ierr) num_ys
112 0 : allocate(ys(num_ys))
113 0 : do k = 1, num_ys
114 0 : read(iounit, *, iostat = ierr) ys(k)
115 : end do
116 :
117 : ! create a 1d array with all the data, point func to it
118 0 : allocate(func1d(4 * num_xs * num_ys))
119 0 : func(1:4, 1:num_xs, 1:num_ys) => func1d(1:4 * num_xs * num_ys)
120 0 : do k = 1, num_xs
121 0 : read(iounit, *, iostat = ierr) func(1, k, :)
122 : end do
123 :
124 0 : if (ierr /= 0) then
125 0 : close(iounit)
126 : end if
127 : ! create interpolator
128 0 : call interp_mkbipm_db(xs, num_xs, ys, num_ys, func1d, num_xs, ierr)
129 0 : end subroutine setup_interpolator
130 :
131 : end subroutine build_roche_interpolators
132 :
133 0 : real(dp) function eval_fp(lq, ar, ierr) result(fp)
134 : ! evaluates fp of the equipotential shell with fractional equivalent radius
135 : ! ar = r/rl and lq = log10(m(other)/m(this)), no units (dimensionless)
136 : real(dp), intent(in) :: lq
137 : real(dp) :: ar
138 : integer, intent(out) :: ierr
139 :
140 : include 'formats'
141 :
142 0 : if (ar >= yvals(num_ypts)) ar = yvals(num_ypts) - 2 * nudge
143 : call interp_evbipm_db(lq, ar, xvals, num_xpts, yvals, num_ypts,&
144 0 : fpfunc1d, num_xpts, fp, ierr)
145 0 : if (ierr /= 0) write(*, 1) "error in eval fp", ar, fp
146 0 : end function eval_fp
147 :
148 0 : real(dp) function eval_ft(lq, ar, ierr) result(ft)
149 : ! evaluates ft of the equipotential shell with fractional equivalent radius
150 : ! ar = r/rl and lq = log10(m(other)/m(this)), no units (dimensionless)
151 : real(dp), intent(in) :: lq
152 : real(dp) :: ar
153 : integer, intent(out) :: ierr
154 :
155 : include 'formats'
156 :
157 0 : if (ar >= yvals(num_ypts)) ar = yvals(num_ypts) - 2 * nudge
158 : call interp_evbipm_db(lq, ar, xvals, num_xpts, yvals, num_ypts,&
159 0 : ftfunc1d, num_xpts, ft, ierr)
160 0 : if (ierr /= 0) write(*, 1) "error in eval ft", ar, ft
161 :
162 0 : end function eval_ft
163 :
164 0 : real(dp) function eval_irot(lq, ar, ierr) result(irot)
165 : ! evaluates moment of inertia irot divided by the spherical equivalent moi (2/3 r_\psi^2)
166 : ! of the equipotential shell with fractional equivalent radius
167 : ! ar = r/rl and lq = log10(m(other)/m(this)), no units (dimensionless)
168 : real(dp), intent(in) :: lq
169 : real(dp) :: ar
170 : integer, intent(out) :: ierr
171 :
172 : include 'formats'
173 :
174 0 : if (ar >= yvals(num_ypts)) ar = yvals(num_ypts) - 2 * nudge
175 : call interp_evbipm_db(lq, ar, xvals, num_xpts, yvals, num_ypts,&
176 0 : irotfunc1d, num_xpts, irot, ierr)
177 0 : if (ierr /= 0) write(*, 1) "error in eval irot", ar, irot
178 :
179 0 : end function eval_irot
180 :
181 0 : real(dp) function eval_psi(lq, ar, ierr) result(psi)
182 : ! evaluates the dimensionless Roche potential of a shell of a given fractional equivalent radius ar = r/rl and
183 : ! lq = log10(m(other)/m(this)); result has no units. Scaling to dimensionfull units should be done by:
184 : ! Psi_dim_full = (Psi_dim_less - q^2/(2(1+q)) * G * M_this / separation, with q = m_other / m_this
185 : ! note there was an error in Fabry et al. (2022) regarding this scaling.
186 :
187 : real(dp), intent(in) :: lq
188 : real(dp) :: ar
189 : integer, intent(out) :: ierr
190 :
191 : include 'formats'
192 :
193 0 : if (ar >= yvals(num_ypts)) ar = yvals(num_ypts) - 2 * nudge
194 : call interp_evbipm_db(lq, ar, xvals, num_xpts, yvals, num_ypts,&
195 0 : psifunc1d, num_xpts, psi, ierr)
196 0 : if (ierr /= 0) write(*, 1) "error in eval psi", ar, psi
197 0 : end function eval_psi
198 :
199 0 : subroutine roche_psi(id, r, psi, ierr)
200 : integer, intent(in) :: id
201 : real(dp), intent(in) :: r
202 : integer, intent(out) :: ierr
203 : real(dp), intent(out) :: psi
204 :
205 : type (star_info), pointer :: s
206 : type (binary_info), pointer :: b
207 : integer :: j, this_star, other_star
208 0 : real(dp) :: m1, m2, a, r_roche, lq
209 :
210 : ierr = 0
211 0 : call star_ptr(id, s, ierr)
212 0 : if (ierr /= 0) return
213 0 : call binary_ptr(s% binary_id, b, ierr)
214 0 : if (ierr /= 0) return
215 :
216 0 : m1 = b% m(this_star)
217 0 : m2 = b% m(other_star)
218 0 : a = b% separation
219 0 : r_roche = b% rl(this_star)
220 : ! if masses or sep are not yet set at the binary level, use these
221 0 : if (m1 <= 0) m1 = b% m1 * Msun
222 0 : if (m2 <= 0) m2 = b% m2 * Msun
223 0 : if (a <= 0) a = pow(standard_cgrav * (m1 + m2) * &
224 0 : pow((b% initial_period_in_days) * 86400, 2) / (4 * pi2), one_third)
225 0 : if (r_roche <= 0) r_roche = eval_rlobe(m1, m2, a)
226 0 : lq = log10(m2 / m1)
227 :
228 0 : psi = eval_psi(lq, r / r_roche, ierr)
229 0 : if (ieee_is_nan(psi)) then
230 0 : psi = -9d99 ! let's put negative a lot for regions outside what we can interpolate
231 : ! this likely happens in the core where r ~= 0, and so psi == -infty
232 : end if
233 :
234 : end subroutine roche_psi
235 :
236 0 : subroutine roche_fp_ft(id, r, fp, ft, r_polar, r_equatorial, report_ierr, ierr)
237 : integer, intent(in) :: id
238 : real(dp), intent(in) :: r
239 : logical, intent(in) :: report_ierr
240 : real(dp), intent(inout) :: r_polar, r_equatorial
241 : type (auto_diff_real_star_order1), intent(out) :: fp, ft
242 : integer, intent(out) :: ierr
243 :
244 : type (star_info), pointer :: s
245 : type (binary_info), pointer :: b
246 : integer :: j, this_star, other_star
247 0 : real(dp) :: r_roche, lq, m1, m2, a, ar
248 :
249 : include 'formats'
250 :
251 0 : if (.not. inter_ok) then
252 0 : write(*, 1) "interpolators not setup, should happen in startup"
253 0 : stop
254 : end if
255 :
256 : ierr = 0
257 0 : call star_ptr(id, s, ierr)
258 0 : if (ierr /= 0) return
259 0 : call binary_ptr(s% binary_id, b, ierr)
260 0 : if (ierr /= 0) return
261 :
262 : this_star = 0
263 : other_star = 0
264 0 : call assign_stars(id, this_star, other_star, ierr)
265 0 : if (ierr /= 0) return
266 :
267 0 : m1 = b% m(this_star)
268 0 : m2 = b% m(other_star)
269 0 : a = b% separation
270 0 : r_roche = b% rl(this_star)
271 : ! if masses or sep are not yet set at the binary level, use these
272 0 : if (m1 <= 0) m1 = b% m1 * Msun
273 0 : if (m2 <= 0) m2 = b% m2 * Msun
274 0 : if (a <= 0) a = pow(standard_cgrav * (m1 + m2) * &
275 0 : pow((b% initial_period_in_days) * 86400, 2) / (4 * pi2), one_third)
276 0 : if (r_roche <= 0) r_roche = eval_rlobe(m1, m2, a)
277 0 : lq = log10(m2 / m1)
278 :
279 0 : ar = r / r_roche
280 : ! set values
281 0 : fp = eval_fp(lq, ar, ierr)
282 0 : ft = eval_ft(lq, ar, ierr)
283 : ! set log derivatives
284 0 : fp% d1Array(i_lnR_00) = (eval_fp(lq, ar + nudge, ierr) - fp% val) / nudge * ar
285 0 : ft% d1Array(i_lnR_00) = (eval_ft(lq, ar + nudge, ierr) - ft% val) / nudge * ar
286 : ! if (fp(j) == 0d0 .or. fp(j) == 0d0) write(*, *) j, ar, fp(j), ft(j), ierr ! debug
287 : ! fix these to the current radius, they're only used for some wind mass loss enhancement
288 0 : r_polar = r
289 0 : r_equatorial = r
290 : end subroutine roche_fp_ft
291 :
292 0 : subroutine roche_irot(id, r00, i_rot)
293 : integer, intent(in) :: id
294 : real(dp), intent(in) :: r00
295 : type (auto_diff_real_star_order1), intent(out) :: i_rot
296 : type (star_info), pointer :: s
297 : type (binary_info), pointer :: b
298 : integer :: ierr, j, this_star, other_star
299 0 : real(dp) :: r_roche, a, m1, m2, ar, lq, eval, d_eval_d_ar
300 :
301 : ierr = 0
302 0 : call star_ptr(id, s, ierr)
303 0 : if (ierr /= 0) return
304 0 : call binary_ptr(s% binary_id, b, ierr)
305 0 : if (ierr /= 0) return
306 :
307 : this_star = 0
308 : other_star = 0
309 0 : call assign_stars(id, this_star, other_star, ierr)
310 0 : if (ierr /= 0) return
311 :
312 0 : m1 = b% m(this_star)
313 0 : m2 = b% m(other_star)
314 0 : a = b% separation
315 0 : r_roche = b% rl(this_star)
316 0 : if (m1 <= 0d0) m1 = b% m1 * Msun
317 0 : if (m2 <= 0d0) m2 = b% m2 * Msun
318 0 : if (a <= 0d0) a = pow(standard_cgrav * (m1 + m2) * &
319 0 : pow((b% initial_period_in_days) * secday, 2) / (4d0 * pi2), one_third)
320 0 : if (r_roche <= 0d0) r_roche = eval_rlobe(m1, m2, a)
321 : !
322 0 : lq = log10(m2 / m1)
323 0 : i_rot = 0d0
324 :
325 0 : if (inter_ok) then
326 0 : ar = r00 / r_roche
327 : ! set value
328 0 : eval = eval_irot(lq, ar, ierr)
329 0 : d_eval_d_ar = (eval_irot(lq, ar + nudge, ierr) - eval) / nudge
330 : ! scale value with spherical moi (see eval_irot)
331 0 : i_rot = eval * (two_thirds * r00 * r00)
332 : ! set radius derivative
333 0 : i_rot% d1Array(i_lnR_00) = two_thirds * r00 * r00 * (ar * d_eval_d_ar + 2 * eval)
334 : ! write(*, *) r00, r_roche, i_rot, w_div_w_crit_roche ! debug
335 : end if
336 : end subroutine roche_irot
337 :
338 0 : subroutine assign_stars(id, this_star, other_star, ierr)
339 : ! determine which star is which in the binary
340 : integer, intent(in) :: id
341 : integer, intent(out) :: this_star, other_star, ierr
342 : type (binary_info), pointer :: b
343 : type (star_info), pointer :: s
344 :
345 : ierr = 0
346 0 : call star_ptr(id, s, ierr)
347 0 : call binary_ptr(s% binary_id, b, ierr)
348 :
349 0 : if (ierr /= 0) return
350 :
351 0 : if (b% s1% id == s% id .and. b% have_star_1) then
352 0 : this_star = 1
353 0 : other_star = 2
354 0 : else if (b% s2% id == s% id .and. b% have_star_2) then
355 0 : this_star = 2
356 0 : other_star = 1
357 : else
358 0 : ierr = 1
359 : end if
360 :
361 : end subroutine assign_stars
362 :
363 : ! determines by what fraction the tidal deformation corrections should be used
364 : ! eg fp = f_switch * fp_tidal + (1-f_switch) * fp_single
365 : ! it uses the synchronicity parameter to estimate this. When the shell is quite synchronous,
366 : ! you probably want to use the tidal corrections, so f_switch -> 1, when not synchronous, the single
367 : ! star rotation deformation is likely more accurate, so f_switch -> 0 in that case.
368 0 : subroutine synchronicity(id, k, omega_in, f_switch, ierr)
369 : integer, intent(in) :: id, k
370 : real(dp), intent(in) :: omega_in
371 : integer, intent(out) :: ierr
372 : real(dp), intent(out) :: f_switch
373 :
374 : type (star_info), pointer :: s
375 : type (binary_info), pointer :: b
376 0 : real(dp) :: omega, omega_sync, f_sync, p
377 :
378 : include 'formats'
379 :
380 0 : call star_ptr(id, s, ierr)
381 0 : call binary_ptr(s% binary_id, b, ierr)
382 0 : if (ierr /= 0) return
383 :
384 0 : if (b% use_other_tidal_deformation_switch_function) then
385 0 : call b% other_tidal_deformation_switch_function(id, k, omega_in, f_switch, ierr)
386 0 : return
387 : end if
388 :
389 0 : p = b% period
390 0 : if (p <= 0d0) p = b% initial_period_in_days * secday
391 0 : omega_sync = 2*pi/p
392 0 : if (ieee_is_nan(omega_in)) then
393 0 : omega = 0d0
394 0 : else if (omega_in < 0d0) then
395 0 : omega = abs(omega_in)
396 : else
397 0 : omega = omega_in
398 : end if
399 0 : f_sync = omega / omega_sync
400 0 : f_sync = min(f_sync, 1d0 / f_sync) ! we could be super or sub synchronous
401 :
402 0 : f_switch = 1d0 / (1d0 + exp(-b% f_sync_switch_width * (f_sync - b% f_sync_switch_from_rot_defor)))
403 : ! apply limiting values if f_switch is far up (or down) the sigmoid
404 0 : if (f_switch < b% f_sync_switch_lim) then
405 0 : f_switch = 0d0
406 0 : else if (1d0 - f_switch < b% f_sync_switch_lim) then
407 0 : f_switch = 1d0
408 : end if
409 :
410 0 : if (ieee_is_nan(f_switch) .and. k == 1) then
411 0 : write(*, 1) "error in synchronicity", f_switch, f_sync, omega, omega_sync, p
412 0 : ierr = 1
413 : end if
414 :
415 : end subroutine synchronicity
416 : end module binary_roche_deformation
|
