Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 Joris Vos & 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_edot
21 :
22 : use const_def, only: dp
23 : use star_lib
24 : use star_def
25 : use binary_def
26 : use binary_tides
27 : use math_lib
28 :
29 : implicit none
30 :
31 : contains
32 :
33 0 : real(dp) function get_edot(b) result(edot)
34 : type (binary_info), pointer :: b
35 : integer :: ierr
36 0 : ierr = 0
37 :
38 0 : if (b% eccentricity == 0d0) then
39 0 : b% edot_tidal = 0d0
40 0 : b% edot_enhance = 0d0
41 0 : b% extra_edot = 0d0
42 0 : b% edot = 0d0
43 : else
44 : ! tidal circularization
45 0 : if (b% do_tidal_circ) then
46 0 : if (.not. b% use_other_edot_tidal) then
47 0 : call edot_tidal(b% binary_id, ierr)
48 : else
49 0 : call b% other_edot_tidal(b% binary_id, ierr)
50 : end if
51 0 : if (ierr /= 0) then
52 0 : write(*,*) 'ERROR while calculating edot_tidal'
53 0 : return
54 : end if
55 : else
56 0 : b% edot_tidal = 0d0
57 : end if
58 :
59 0 : if (b% edot_tidal < -b% max_abs_edot_tidal) then
60 0 : b% edot_tidal = -b% max_abs_edot_tidal
61 : end if
62 :
63 : ! eccentricity enhancement
64 0 : if (b% use_eccentricity_enhancement) then
65 0 : if (.not. b% use_other_edot_enhance) then
66 0 : call edot_enhancement_Isotropic(b% binary_id, ierr)
67 : else
68 0 : call b% other_edot_enhance(b% binary_id, ierr)
69 : end if
70 0 : if (ierr /= 0) then
71 0 : write(*,*) 'ERROR while calculating edot_enhance'
72 0 : return
73 : end if
74 : else
75 0 : b% edot_enhance = 0d0
76 : end if
77 :
78 0 : if (b% edot_enhance > b% max_abs_edot_enhance) then
79 0 : b% edot_enhance = b% max_abs_edot_enhance
80 : end if
81 :
82 : ! user defined eccentricity changes
83 0 : if (b% use_other_extra_edot) then
84 0 : call b% other_extra_edot(b% binary_id, ierr)
85 0 : if (ierr /= 0) then
86 0 : write(*,*) 'ERROR while calculating extra_edot'
87 0 : return
88 : end if
89 : else
90 0 : b% extra_edot = 0d0
91 : end if
92 :
93 0 : b% edot = b% edot_tidal + b% edot_enhance + b% extra_edot
94 :
95 : end if
96 :
97 0 : edot = b% edot
98 :
99 0 : end function get_edot
100 :
101 : ! ==========================================
102 : ! edot TIDAL
103 : ! ==========================================
104 0 : subroutine edot_tidal(binary_id, ierr)
105 : integer, intent(in) :: binary_id
106 : integer, intent(out) :: ierr
107 : type (binary_info), pointer :: b
108 :
109 : ierr = 0
110 0 : call binary_ptr(binary_id, b, ierr)
111 0 : if (ierr /= 0) then
112 0 : write(*,*) 'failed in binary_ptr'
113 0 : return
114 : end if
115 :
116 0 : b% edot_tidal = 0d0
117 :
118 0 : if (b% point_mass_i /= 1) then
119 0 : if (b% circ_type_1 == "Hut_conv") then
120 0 : b% edot_tidal = edot_tidal_Hut(b, b% s1, .true., ierr)
121 0 : if(ierr/=0) return
122 0 : else if (b% circ_type_1 == "Hut_rad") then
123 0 : b% edot_tidal = edot_tidal_Hut(b, b% s1, .false., ierr)
124 0 : if(ierr/=0) return
125 : else
126 0 : write(*,*) "Unrecognized circ_type_1", b% circ_type_1
127 : end if
128 : end if
129 0 : if (b% point_mass_i /= 2) then
130 0 : if (b% circ_type_2 == "Hut_conv") then
131 0 : b% edot_tidal = b% edot_tidal + edot_tidal_Hut(b, b% s2, .true., ierr)
132 0 : if(ierr/=0) return
133 0 : else if (b% circ_type_2 == "Hut_rad") then
134 0 : b% edot_tidal = b% edot_tidal + edot_tidal_Hut(b, b% s2, .false., ierr)
135 0 : if(ierr/=0) return
136 : else
137 0 : write(*,*) "Unrecognized circ_type_2", b% circ_type_2
138 : end if
139 : end if
140 :
141 0 : if (b% model_twins_flag) then
142 0 : b% edot_tidal = b% edot_tidal + b% edot_tidal
143 : end if
144 : end subroutine edot_tidal
145 :
146 0 : real(dp) function edot_tidal_Hut(b, s, has_convective_envelope, ierr) result(edot_tidal)
147 : type (binary_info), pointer :: b
148 : type (star_info), pointer :: s
149 : logical, intent(in) :: has_convective_envelope
150 0 : real(dp) :: m, porb, r_phot, osep, qratio, omega_s, omega_sync, kdivt
151 : integer, intent(out) :: ierr
152 :
153 0 : edot_tidal = 0d0
154 :
155 0 : porb = b% period
156 0 : omega_sync = 2d0*pi/b% period
157 0 : omega_s = s% omega_avg_surf
158 0 : osep = b% separation
159 :
160 0 : qratio = b% m(b% a_i) / b% m(b% d_i)
161 0 : if (is_donor(b, s)) then
162 0 : m = b% m(b% d_i)
163 0 : r_phot = b% r(b% d_i)
164 : else
165 0 : qratio = 1.0d0/qratio
166 0 : m = b% m(b% a_i)
167 0 : r_phot = b% r(b% a_i)
168 : end if
169 :
170 : ! eq. (10) of Hut, P. 1981, A&A, 99, 126
171 : edot_tidal = -27.0d0*qratio*(1+qratio)*pow8(r_phot/osep) &
172 0 : * b% eccentricity*pow(1-pow2(b% eccentricity),-6.5d0)*b% Ftid_1
173 : ! add multiplication by (k/T), eq. (29) of Hurley et al. 2002
174 0 : kdivt = k_div_T(b, s, has_convective_envelope, ierr)
175 0 : if(ierr/=0) return
176 0 : edot_tidal = edot_tidal * kdivt
177 : ! add terms dependant on omega
178 : edot_tidal = edot_tidal*(f3(b% eccentricity) - &
179 : 11d0/18d0 * omega_s / omega_sync * f4(b% eccentricity) * &
180 0 : pow(1-pow2(b% eccentricity),1.5d0))
181 :
182 0 : end function edot_tidal_Hut
183 :
184 : ! ==========================================
185 : ! Edot MASS LOSS
186 : ! ==========================================
187 :
188 0 : subroutine edot_enhancement_Isotropic(binary_id, ierr)
189 : integer, intent(in) :: binary_id
190 : integer, intent(out) :: ierr
191 : type (binary_info), pointer :: b
192 : integer :: i
193 0 : real(dp) :: de, Mtot, costh
194 :
195 : ierr = 0
196 0 : call binary_ptr(binary_id, b, ierr)
197 0 : if (ierr /= 0) then
198 0 : write(*,*) 'failed in binary_ptr'
199 0 : return
200 : end if
201 :
202 0 : b% edot_enhance = 0d0
203 :
204 : ! cos_cr isn't vectorized, so we have to do this in a loop
205 0 : do i = 1, b% anomaly_steps
206 0 : costh = cos(b% theta_co(i))
207 :
208 0 : b% e1(i) = b% eccentricity + costh
209 0 : b% e2(i) = 2d0*costh + b% eccentricity*(1d0 + costh*costh)
210 0 : b% e3(i) = b% eccentricity*(1d0 - costh*costh) ! = b% eccentricity*sin(b% theta_co)**2
211 : end do
212 :
213 : ! xfer = min(b% wind_xfer_fraction, b% xfer_fraction)
214 0 : Mtot = b% m(1) + b% m(2) ! total mass in gr
215 :
216 0 : b% edot_theta = - b% mdot_donor_theta / Mtot * b% e1 !- &
217 : ! b% mdot_donor_theta * xfer / b% m(b% a_i) * (b% m(b% d_i) / Mtot * &
218 : ! ((b% m(b% a_i)**2 / b% m(b% d_i)**2 - 1 ) * e2 - e3 ))
219 :
220 : ! integrate to get total eccentricity enhancement
221 : de = 0d0
222 0 : do i = 2, b% anomaly_steps ! trapezoidal integration
223 0 : de = de + 0.5d0 * (b% edot_theta(i-1) + b% edot_theta(i)) * (b% time_co(i) - b% time_co(i-1))
224 : end do
225 :
226 0 : b% edot_enhance = de
227 :
228 : end subroutine edot_enhancement_Isotropic
229 :
230 :
231 : end module binary_edot
|
