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_utils
22 :
23 : use const_def, only: dp, pi, one_third, standard_cgrav, msun
24 : use math_lib
25 : use star_lib
26 : use star_def
27 : use binary_def
28 :
29 : implicit none
30 :
31 : contains
32 :
33 0 : subroutine set_ignore_rlof_flag(binary_id, ignore_rlof_flag, ierr)
34 : integer, intent(in) :: binary_id
35 : logical, intent(in) :: ignore_rlof_flag
36 : integer, intent(out) :: ierr
37 :
38 : type (binary_info), pointer :: b
39 :
40 : ierr = 0
41 0 : call binary_ptr(binary_id, b, ierr)
42 0 : if (ierr /= 0) return
43 :
44 0 : b% ignore_rlof_flag = ignore_rlof_flag
45 : end subroutine set_ignore_rlof_flag
46 :
47 0 : subroutine set_point_mass_i(binary_id, point_mass_i, ierr)
48 : integer, intent(in) :: binary_id
49 : integer, intent(in) :: point_mass_i
50 : integer, intent(out) :: ierr
51 :
52 : type (binary_info), pointer :: b
53 :
54 : ierr = 0
55 0 : call binary_ptr(binary_id, b, ierr)
56 0 : if (ierr /= 0) return
57 :
58 0 : b% point_mass_i = point_mass_i
59 0 : if (point_mass_i == 1 .and. .not. b% have_star_2) then
60 0 : ierr = -1
61 0 : write(*,*) 'ERROR: setting point_mass_i=1 with have_star_2=.false.'
62 0 : else if (point_mass_i == 2 .and. .not. b% have_star_1) then
63 0 : ierr = -1
64 0 : write(*,*) 'ERROR: setting point_mass_i=2 with have_star_1=.false.'
65 : end if
66 :
67 0 : if (point_mass_i == 1 .and. b% d_i == 1) then
68 0 : b% d_i = 2
69 0 : b% a_i = 1
70 0 : b% s_donor => b% s2
71 0 : b% s_accretor => b% s1
72 0 : b% mtransfer_rate = 0d0
73 0 : b% change_factor = b% max_change_factor
74 0 : else if (point_mass_i == 2 .and. b% d_i == 2) then
75 0 : b% d_i = 1
76 0 : b% a_i = 2
77 0 : b% s_donor => b% s1
78 0 : b% s_accretor => b% s2
79 0 : b% mtransfer_rate = 0d0
80 0 : b% change_factor = b% max_change_factor
81 : end if
82 : end subroutine set_point_mass_i
83 :
84 0 : subroutine set_model_twins_flag(binary_id, model_twins_flag, ierr)
85 : integer, intent(in) :: binary_id
86 : logical, intent(in) :: model_twins_flag
87 : integer, intent(out) :: ierr
88 :
89 : type (binary_info), pointer :: b
90 :
91 : ierr = 0
92 0 : call binary_ptr(binary_id, b, ierr)
93 0 : if (ierr /= 0) return
94 :
95 0 : b% model_twins_flag = model_twins_flag
96 :
97 : ! also need to set ignore_rlog_flag to true
98 0 : if (model_twins_flag) then
99 0 : call set_ignore_rlof_flag(binary_id, .true., ierr)
100 0 : if (ierr /= 0) return
101 0 : call set_point_mass_i(binary_id, 2, ierr)
102 : end if
103 : end subroutine set_model_twins_flag
104 :
105 0 : subroutine set_m1(binary_id, m1, ierr)
106 : integer, intent(in) :: binary_id
107 : real(dp), intent(in) :: m1
108 : integer, intent(out) :: ierr
109 :
110 : type (binary_info), pointer :: b
111 :
112 : ierr = 0
113 0 : call binary_ptr(binary_id, b, ierr)
114 0 : if (ierr /= 0) return
115 :
116 0 : if (b% point_mass_i /= 1) then
117 0 : ierr = -1
118 0 : write(*,*) "ERROR: adjusting m1 with point_mass_i/=1 has no effect"
119 0 : return
120 : end if
121 0 : b% m(1) = m1*Msun
122 0 : if (b% model_twins_flag) b% m(2) = m1
123 0 : call set_separation_eccentricity(binary_id, b% separation, b% eccentricity, ierr)
124 : end subroutine set_m1
125 :
126 0 : subroutine set_m2(binary_id, m2, ierr)
127 : integer, intent(in) :: binary_id
128 : real(dp), intent(in) :: m2
129 : integer, intent(out) :: ierr
130 :
131 : type (binary_info), pointer :: b
132 :
133 : ierr = 0
134 0 : call binary_ptr(binary_id, b, ierr)
135 0 : if (ierr /= 0) return
136 :
137 0 : if (b% point_mass_i /= 2) then
138 0 : ierr = -1
139 0 : write(*,*) "ERROR: adjusting m2 with point_mass_i/=2 has no effect"
140 0 : return
141 : end if
142 0 : b% m(2) = m2*Msun
143 0 : call set_separation_eccentricity(binary_id, b% separation, b% eccentricity, ierr)
144 : end subroutine set_m2
145 :
146 0 : subroutine set_period_eccentricity(binary_id, period, eccentricity, ierr)
147 : integer, intent(in) :: binary_id
148 : real(dp), intent(in) :: period ! in seconds
149 : real(dp), intent(in) :: eccentricity
150 : type (binary_info), pointer :: b
151 : integer, intent(out) :: ierr
152 0 : call binary_ptr(binary_id,b,ierr)
153 0 : if (ierr /= 0) return
154 :
155 0 : b% eccentricity = eccentricity
156 0 : b% period = period
157 : b% separation = &
158 0 : pow((standard_cgrav*(b% m(1)+b% m(2)))*pow2(b% period/(2*pi)),one_third)
159 0 : call set_angular_momentum_j(binary_id)
160 :
161 : end subroutine set_period_eccentricity
162 :
163 0 : subroutine set_separation_eccentricity(binary_id, separation, eccentricity, ierr)
164 : integer, intent(in) :: binary_id
165 : real(dp), intent(in) :: separation ! in cm
166 : real(dp) :: eccentricity
167 : type (binary_info), pointer :: b
168 : integer, intent(out) :: ierr
169 0 : call binary_ptr(binary_id,b,ierr)
170 0 : if (ierr /= 0) return
171 :
172 0 : b% eccentricity = eccentricity
173 0 : b% separation = separation
174 : b% period = &
175 : (2*pi)*sqrt(b% separation*b% separation*b% separation/&
176 0 : (standard_cgrav*(b% m(1)+b% m(2))))
177 0 : call set_angular_momentum_j(binary_id)
178 :
179 : end subroutine set_separation_eccentricity
180 :
181 0 : subroutine set_angular_momentum_j(binary_id)
182 : ! Sets b% angular_momentum_j in terms of the masses, separation and eccentricity
183 : ! also sets the Roche lobe sizes and relative overflows
184 : integer, intent(in) :: binary_id
185 : type (binary_info), pointer :: b
186 : integer :: ierr
187 0 : call binary_ptr(binary_id,b,ierr)
188 0 : if (ierr /= 0) return
189 :
190 : b% angular_momentum_j = b% m(1) * b% m(2) * sqrt(standard_cgrav *&
191 0 : b% separation * (1 - pow2(b% eccentricity)) / (b% m(1) + b% m(2)) )
192 :
193 0 : b% rl(1) = eval_rlobe(b% m(1), b% m(2), b% separation)
194 0 : b% rl(2) = eval_rlobe(b% m(2), b% m(1), b% separation)
195 : b% rl_relative_gap(1) = (b% r(1) - b% rl(1) * (1 - b% eccentricity) ) / &
196 0 : b% rl(1) / (1 - b% eccentricity) ! gap < 0 means out of contact
197 : b% rl_relative_gap(2) = (b% r(2) - b% rl(2) * (1 - b% eccentricity) ) / &
198 0 : b% rl(2) / (1 - b% eccentricity) ! gap < 0 means out of contact
199 :
200 0 : b% ignore_hard_limits_this_step = .true.
201 :
202 : end subroutine set_angular_momentum_j
203 :
204 0 : real(dp) function eval_rlobe(m1, m2, a) result(rlobe)
205 : real(dp), intent(in) :: m1, m2, a
206 : real(dp) :: q
207 0 : q = pow(m1/m2,one_third)
208 : ! Roche lobe size for star of mass m1 with a
209 : ! companion of mass m2 at separation a, according to
210 : ! the approximation of Eggleton 1983, apj 268:368-369
211 0 : rlobe = a*0.49d0*q*q/(0.6d0*q*q + log1p(q))
212 0 : end function eval_rlobe
213 :
214 : end module binary_utils
215 :
|
