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 : module binary_do_one_utils
21 :
22 : use binary_def
23 : use const_def, only: dp, pi, msun, rsun, secyer, standard_cgrav
24 : use math_lib
25 :
26 : implicit none
27 :
28 : contains
29 :
30 0 : subroutine write_binary_terminal_header(b)
31 : type (binary_info), pointer :: b
32 0 : if (b% model_number <= b% recent_binary_log_header) return
33 0 : if (b% just_wrote_binary_terminal_header) return
34 0 : b% recent_binary_log_header = b% model_number
35 0 : call do_show_binary_terminal_header(b)
36 0 : b% just_wrote_binary_terminal_header = .true.
37 : end subroutine write_binary_terminal_header
38 :
39 :
40 0 : subroutine do_show_binary_log_description(id, ierr)
41 : integer, intent(in) :: id
42 : integer, intent(out) :: ierr
43 : type (binary_info), pointer :: b
44 : ierr = 0
45 0 : call get_binary_ptr(id, b, ierr)
46 0 : if (ierr /= 0) return
47 0 : write(*,'(A)')
48 0 : write(*,'(a)') " The binary terminal output contains the following information"
49 0 : write(*,'(A)')
50 0 : write(*,'(a)') " 'step' is the number of steps since the start of the run,"
51 0 : write(*,'(a)') " 'lg_dt' is log10 timestep in years,"
52 0 : write(*,'(a)') " 'age_yr' is the simulated years since the start run,"
53 0 : write(*,'(a)') " 'M1+M2' is the total mass of the system (Msun),"
54 0 : write(*,'(a)') " 'M1' is the mass of the primary (Msun)"
55 0 : write(*,'(a)') " 'M2' is the mass of the secondary (Msun)"
56 0 : write(*,'(a)') " 'separ' is the semi-major axis of the orbit (Rsun),"
57 0 : write(*,'(a)') " 'R1' is the radius of the primary (Rsun)"
58 0 : write(*,'(a)') " 'R2' is the radius of the secondary (Rsun)"
59 0 : write(*,'(a)') " 'Porb' is the orbital period (days),"
60 0 : write(*,'(a)') " 'P1' is the rotation period of star 1 (days, zero if not modeling rotation),"
61 0 : write(*,'(a)') " 'P2' is the rotation period of star 2 (days, zero if not modeling rotation),"
62 0 : write(*,'(a)') " 'e' orbital eccentricity,"
63 0 : write(*,'(a)') " 'dot_e' time derivative of e (1/yr),"
64 0 : write(*,'(a)') " 'Eorb' orbital energy G*M1*M2/2*separation (ergs),"
65 0 : write(*,'(a)') " 'M2/M1' mass ratio,"
66 0 : write(*,'(a)') " 'vorb1' orbital velocity of star 1 (km/s),"
67 0 : write(*,'(a)') " 'vorb2' orbital velocity of star 2 (km/s),"
68 0 : write(*,'(a)') " 'pm_i' index of star evolved as point mass, zero if both stars are modeled,"
69 0 : write(*,'(a)') " 'RL1' Roche lobe radius of star 1 (Rsun),"
70 0 : write(*,'(a)') " 'Rl2' Roche lobe radius of star 2 (Rsun),"
71 0 : write(*,'(a)') " 'donor_i' index of star taken as donor,"
72 0 : write(*,'(a)') " 'RL_gap1' (R1-Rl1)/Rl1,"
73 0 : write(*,'(a)') " 'RL_gap2' (R2-Rl2)/Rl2,"
74 0 : write(*,'(a)') " 'dot_Mmt', mass transfer rate (Msun/yr),"
75 0 : write(*,'(a)') " 'dot_M1', time derivative for the mass of star 1 (Msun/yr),"
76 0 : write(*,'(a)') " 'dot_M2', time derivative for the mass of star 2 (Msun/yr),"
77 0 : write(*,'(a)') " 'eff', mass transfer efficiency, computed as -dot_M2/dot_M1 (zero if dot_M1=0),"
78 0 : write(*,'(a)') " 'dot_Medd', Eddington accretion rate (Msun/yr),"
79 0 : write(*,'(a)') " 'L_acc', accretion luminosity when accreting to a point mass (ergs/s),"
80 0 : write(*,'(a)') " 'Jorb', orbital angular momentum (g*cm^2/s)"
81 0 : write(*,'(a)') " 'spin1', spin angular momentum of star 1 (g*cm^2/s),"
82 0 : write(*,'(a)') " 'spin2', spin angular momentum of star 2 (g*cm^2/s),"
83 0 : write(*,'(a)') " 'dot_J', time derivative of Jorb (g*cm^2/s^2),"
84 0 : write(*,'(a)') " 'dot_Jgr', time derivative of Jorb due to gravitational waves (g*cm^2/s^2),"
85 0 : write(*,'(a)') " 'dot_Jml', time derivative of Jorb due to mass loss (g*cm^2/s^2),"
86 0 : write(*,'(a)') " 'dot_Jmb', time derivative of Jorb due to magnetic braking (g*cm^2/s^2),"
87 0 : write(*,'(a)') " 'dot_Jls', time derivative of Jorb due to spin-orbit coupling (g*cm^2/s^2),"
88 0 : write(*,'(a)') " 'rlo_iters', number of iterations for implicit calculation of mass transfer,"
89 0 : write(*,'(A)')
90 0 : write(*,'(a)') " All this and more can be saved in binary_history.data during the run."
91 : end subroutine do_show_binary_log_description
92 :
93 :
94 0 : subroutine do_show_binary_terminal_header(b)
95 : type (binary_info), pointer :: b
96 0 : call output_binary_terminal_header(b,terminal_iounit)
97 0 : if (b% extra_binary_terminal_iounit > 0) &
98 0 : call output_binary_terminal_header(b,b% extra_binary_terminal_iounit)
99 0 : end subroutine do_show_binary_terminal_header
100 :
101 :
102 0 : subroutine output_binary_terminal_header(b,io)
103 : type (binary_info), pointer :: b
104 : integer, intent(in) :: io
105 : write(io,'(a)') &
106 : '_______________________________________________________________________' // &
107 0 : '___________________________________________________________________________'
108 0 : write(io,'(A)')
109 : write(io,'(a)') &
110 : 'binary_step M1+M2 separ Porb e M2/M1' // &
111 0 : ' pm_i donor_i dot_Mmt eff Jorb dot_J dot_Jmb'
112 : write(io,'(a)') &
113 : ' lg_dt M1 R1 P1 dot_e vorb1' // &
114 0 : ' RL1 Rl_gap1 dot_M1 dot_Medd spin1 dot_Jgr dot_Jls'
115 : write(io,'(a)') &
116 : ' age_yr M2 R2 P2 Eorb vorb2' // &
117 0 : ' RL2 Rl_gap2 dot_M2 L_acc spin2 dot_Jml rlo_iters'
118 : write(io,'(a)') &
119 : '_______________________________________________________________________' // &
120 0 : '___________________________________________________________________________'
121 0 : write(io,'(A)')
122 :
123 0 : end subroutine output_binary_terminal_header
124 :
125 :
126 0 : subroutine do_binary_terminal_summary(b)
127 : type (binary_info), pointer :: b
128 0 : call output_binary_terminal_summary(b,terminal_iounit)
129 0 : if (b% extra_binary_terminal_iounit > 0) then
130 0 : call output_binary_terminal_summary(b,b% extra_binary_terminal_iounit)
131 0 : flush(b% extra_binary_terminal_iounit)
132 : end if
133 0 : end subroutine do_binary_terminal_summary
134 :
135 :
136 0 : subroutine output_binary_terminal_summary(b,io)
137 : type (binary_info), pointer :: b
138 : integer, intent(in) :: io
139 :
140 : real(dp) :: age, time_step, total_mass
141 0 : real(dp) :: Eorb, vorb1, vorb2, dot_M1, dot_M2, eff, dot_Medd, spin1, spin2, P1, P2
142 : integer :: model, ierr, rlo_iters
143 : character (len=90) :: fmt, fmt1, fmt2, fmt3, fmt4, fmt5, fmt6
144 :
145 : include 'formats'
146 :
147 0 : age = b% binary_age
148 0 : time_step = b% time_step
149 0 : model = b% model_number
150 0 : rlo_iters = b% num_tries
151 :
152 0 : total_mass = (b% m(1) + b% m(2))/Msun
153 :
154 0 : vorb1 = 2.0d0 * pi * b% m(2)/(b% m(1) + b% m(2)) * b% separation / b% period / 1.0d5
155 0 : vorb2 = 2.0d0 * pi * b% m(1)/(b% m(1) + b% m(2)) * b% separation / b% period / 1.0d5
156 0 : Eorb = -standard_cgrav * b% m(1) * b% m(2) / (2*b% separation)
157 0 : dot_M1 = b% component_mdot(1)/Msun*secyer
158 0 : if (abs(dot_M1) < 1d-99) dot_M1 = 0d0
159 0 : dot_M2 = b% component_mdot(2)/Msun*secyer
160 0 : if (abs(dot_M2) < 1d-99) dot_M2 = 0d0
161 0 : if (b% component_mdot(b% d_i) == 0d0) then
162 0 : eff = 1d0
163 : else
164 0 : eff = -b% component_mdot(b% a_i)/b% component_mdot(b% d_i)
165 : end if
166 0 : if (b% limit_retention_by_mdot_edd) then
167 0 : dot_Medd = b% mdot_edd/Msun*secyer
168 : else
169 0 : dot_Medd = 1d99
170 : end if
171 0 : if (b% point_mass_i /= 1) then
172 0 : spin1 = b% s1% total_angular_momentum
173 : else
174 0 : spin1 = 0d0
175 : end if
176 :
177 0 : if (b% point_mass_i /= 2) then
178 0 : spin2 = b% s2% total_angular_momentum
179 : else
180 0 : if (.not. b% model_twins_flag) then
181 0 : spin2 = 0d0
182 : else
183 0 : spin2 = b% s1% total_angular_momentum
184 : end if
185 : end if
186 :
187 0 : P1 = 0d0
188 0 : if (b% point_mass_i /= 1) then
189 0 : if (b% s1% rotation_flag) then
190 0 : if (abs(b% s1% omega_avg_surf) > 0) then
191 0 : P1 = 2*pi/(b% s1% omega_avg_surf*24d0*3600d0)
192 : end if
193 : end if
194 : end if
195 :
196 0 : P2 = 0d0
197 0 : if (b% point_mass_i /= 2)then
198 0 : if (b% s2% rotation_flag) then
199 0 : if (abs(b% s2% omega_avg_surf) > 0) then
200 0 : P2 = 2*pi/(b% s2% omega_avg_surf*24d0*3600d0)
201 : end if
202 : end if
203 0 : else if (b% model_twins_flag) then
204 0 : if (b% s1% rotation_flag) then
205 0 : if (abs(b% s1% omega_avg_surf) > 0) then
206 0 : P2 = 2*pi/(b% s1% omega_avg_surf*24d0*3600d0)
207 : end if
208 : end if
209 : end if
210 :
211 0 : ierr = 0
212 :
213 : !make format strings for first line of output
214 : !step and m1+m2
215 0 : if (total_mass >= 1d2) then
216 0 : fmt1 = '(a3,i8,1pe11.3,0p,'
217 : else
218 0 : fmt1 = '(a3,i8,f11.6,'
219 : end if
220 : !separation
221 0 : if (b% separation/Rsun >= 1d2) then
222 0 : fmt2 = '1pe11.3,0p,'
223 : else
224 0 : fmt2 = 'f11.6,'
225 : end if
226 : !Porb
227 0 : if (b% period/(3600d0*24d0) >= 1d2) then
228 0 : fmt3 = '1pe11.3,0p,'
229 : else
230 0 : fmt3 = 'f11.6,'
231 : end if
232 : !eccentricity
233 0 : if (b% eccentricity >= 1d-2) then
234 0 : fmt4 = 'f11.6,'
235 : else
236 0 : fmt4 = '1pe11.3,0p,'
237 : end if
238 : !mass ratio, pm_i, donor_index and dot_Mmt
239 0 : if (b% m(2)/b% m(1) >= 1d-2) then
240 0 : fmt5 = 'f11.6,2i11,1pe11.3,0p,'
241 : else
242 0 : fmt5 = '1pe11.3,0p,2i11,1pe11.3,0p,'
243 : end if
244 : !eff
245 0 : if (eff >= 1d-2) then
246 0 : fmt6 = 'f11.6,3(1pe11.3))'
247 : else
248 0 : fmt6 = '4(1pe11.3))'
249 : end if
250 :
251 :
252 0 : fmt = trim(fmt1) // trim(fmt2) // trim(fmt3) // trim(fmt4) // trim(fmt5) // trim(fmt6)
253 : write(io,fmt=fmt) &
254 0 : 'bin', &
255 0 : model, &
256 0 : total_mass, &
257 0 : b% separation / Rsun, &
258 0 : b% period / (3600d0*24d0), &
259 0 : b% eccentricity, &
260 0 : b% m(2)/b% m(1), &
261 0 : b% point_mass_i, &
262 0 : b% d_i, &
263 0 : b% step_mtransfer_rate/Msun*secyer, &
264 0 : eff, &
265 0 : b% angular_momentum_j, &
266 0 : b% jdot, &
267 0 : b% jdot_mb
268 :
269 : !make format strings for second line of output
270 : !lg_dt and M1
271 0 : if (b% m(1)/Msun >= 1d2) then
272 0 : fmt1 = '(f11.6,1pe11.3,0p,'
273 : else
274 0 : fmt1 = '(2f11.6,'
275 : end if
276 : !R1
277 0 : if (b% r(1)/Rsun >= 1d2) then
278 0 : fmt2 = '1pe11.3,0p,'
279 : else
280 0 : fmt2 = 'f11.6,'
281 : end if
282 : !P1, dot_e
283 0 : if (P1 >= 1d2) then
284 0 : fmt3 = '2(1pe11.3),0p'
285 : else
286 0 : fmt3 = 'f11.6,1pe11.3,0p'
287 : end if
288 : !vorb1
289 0 : if (vorb1 >= 1d-2 .and. vorb1 < 1d4) then
290 0 : fmt4 = 'f11.6,'
291 : else
292 0 : fmt4 = '1pe11.3,0p,'
293 : end if
294 : !RL1 and the rest
295 0 : if (b% rl(1)/Rsun >= 1d2) then
296 0 : fmt5 = '7(1pe11.3))'
297 : else
298 0 : fmt5 = 'f11.6,6(1pe11.3))'
299 : end if
300 :
301 0 : fmt = trim(fmt1) // trim(fmt2) // trim(fmt3) // trim(fmt4) // trim(fmt5)
302 : write(io,fmt=fmt) &
303 0 : safe_log10(time_step), &
304 0 : b% m(1)/Msun, &
305 0 : b% r(1)/Rsun, &
306 0 : P1, &
307 0 : b% edot, &
308 0 : vorb1, &
309 0 : b% rl(1)/Rsun, &
310 0 : b% rl_relative_gap(1), &
311 0 : dot_M1, &
312 0 : dot_Medd, &
313 0 : spin1, &
314 0 : b% jdot_gr, &
315 0 : b% jdot_ls
316 :
317 : !make format strings for third line of output
318 : !age_yr and M2
319 0 : if (b% m(2)/Msun >= 1d2) then
320 0 : fmt1 = '(1pe11.4,0p,1pe11.3,0p,'
321 : else
322 0 : fmt1 = '(1pe11.4,0p,f11.6,'
323 : end if
324 : !R2
325 0 : if (b% r(2)/Rsun >= 1d2) then
326 0 : fmt2 = '1pe11.3,0p,'
327 : else
328 0 : fmt2 = 'f11.6,'
329 : end if
330 : !P2, Eorb
331 0 : if (P2 >= 1d2) then
332 0 : fmt3 = '2(1pe11.3),0p'
333 : else
334 0 : fmt3 = 'f11.6,1pe11.3,0p'
335 : end if
336 : !vorb2
337 0 : if (vorb2 >= 1d-2 .and. vorb2 < 1d4) then
338 0 : fmt4 = 'f11.6,'
339 : else
340 0 : fmt4 = '1pe11.3,0p,'
341 : end if
342 : !RL2 and the rest
343 0 : if (b% rl(2)/Rsun >= 1d2) then
344 0 : fmt5 = '6(1pe11.3),i11)'
345 : else
346 0 : fmt5 = 'f11.6,5(1pe11.3),0p,i11)'
347 : end if
348 :
349 0 : fmt = trim(fmt1) // trim(fmt2) // trim(fmt3) // trim(fmt4) // trim(fmt5)
350 : write(io,fmt=fmt) &
351 0 : age, &
352 0 : b% m(2)/Msun, &
353 0 : b% r(2)/Rsun, &
354 0 : P2, &
355 0 : Eorb, &
356 0 : vorb2, &
357 0 : b% rl(2)/Rsun, &
358 0 : b% rl_relative_gap(2), &
359 0 : dot_M2, &
360 0 : b% accretion_luminosity, &
361 0 : spin2, &
362 0 : b% jdot_ml, &
363 0 : rlo_iters
364 :
365 0 : write(io,'(A)')
366 :
367 0 : b% just_wrote_binary_terminal_header = .false.
368 :
369 0 : end subroutine output_binary_terminal_summary
370 :
371 :
372 : end module binary_do_one_utils
373 :
|
