Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2022 The MESA Team & Matthias Fabry
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 pgbinary_orbit
21 :
22 : use binary_private_def
23 : use pgbinary_support
24 :
25 : implicit none
26 :
27 :
28 : contains
29 :
30 :
31 0 : subroutine Orbit_plot(id, device_id, ierr)
32 : integer, intent(in) :: id, device_id
33 : integer, intent(out) :: ierr
34 : type (binary_info), pointer :: b
35 : ierr = 0
36 0 : call get_binary_ptr(id, b, ierr)
37 0 : if (ierr /= 0) return
38 0 : call pgslct(device_id)
39 0 : call pgbbuf()
40 0 : call pgeras()
41 : call do_Orbit_plot(b, id, device_id, &
42 : b% pg% Orbit_xleft, b% pg% Orbit_xright, &
43 : b% pg% Orbit_ybot, b% pg% Orbit_ytop, .false., &
44 0 : b% pg% Orbit_title, b% pg% Orbit_txt_scale_factor, ierr)
45 0 : if (ierr /= 0) return
46 0 : call pgebuf()
47 : end subroutine Orbit_plot
48 :
49 :
50 0 : subroutine do_Orbit_plot(b, id, device_id, &
51 : winxmin, winxmax, winymin, winymax, subplot, title, txt_scale, ierr)
52 : type (binary_info), pointer :: b
53 : integer, intent(in) :: id, device_id
54 : real, intent(in) :: winxmin, winxmax, winymin, winymax, txt_scale
55 : logical, intent(in) :: subplot
56 : character (len = *), intent(in) :: title
57 : integer, intent(out) :: ierr
58 : call orbit_panel(b, device_id, &
59 0 : winxmin, winxmax, winymin, winymax, subplot, title, txt_scale, ierr)
60 0 : end subroutine do_Orbit_plot
61 :
62 :
63 0 : subroutine orbit_panel(b, device_id, &
64 : winxmin, winxmax, winymin, winymax, subplot, title, txt_scale, ierr)
65 :
66 : use num_lib, only : safe_root_with_guess
67 : use math_lib, only : pow
68 : use const_def, only : dp, pi
69 : use pgstar_colors
70 :
71 : type (binary_info), pointer :: b
72 : integer, intent(in) :: device_id
73 : real, intent(in) :: winxmin, winxmax, winymin, winymax, txt_scale
74 : logical, intent(in) :: subplot
75 : character (len = *), intent(in) :: title
76 : integer, intent(out) :: ierr
77 : integer :: i, icut
78 : integer, parameter :: num_points = 50
79 : ! cannot be dp's! PGPLOT doesn't know what to do with them...
80 : real, dimension(2 * num_points + 1) :: &
81 0 : thetas, r1s, r2s, x1s, x2s, y1s, y2s, rs, phis, &
82 0 : x1s_RL, x2s_RL, y1s_RL, y2s_RL
83 0 : real :: a1, a2, e, x1max, x2max, xmax
84 0 : integer, pointer :: ipar(:) ! (lipar)
85 0 : real(dp), pointer :: rpar(:) ! (lrpar)
86 0 : real(dp) :: cosp, q, this_psi, xl1
87 :
88 : include 'formats'
89 :
90 0 : ierr = 0
91 0 : call pgsave
92 0 : call pgsvp(winxmin, winxmax, winymin, winymax)
93 0 : if (.not. subplot) then
94 0 : call show_model_number_pgbinary(b)
95 0 : call show_age_pgbinary(b)
96 : end if
97 0 : call show_title_pgbinary(b, title)
98 0 : call pgunsa
99 :
100 0 : a1 = 1 / (1 + b% m(1) / b% m(2))
101 0 : a2 = 1 / (1 + b% m(2) / b% m(1))
102 0 : e = b% eccentricity
103 :
104 0 : !$OMP PARALLEL DO PRIVATE(i) SCHEDULE(dynamic,2)
105 : do i = 1, num_points
106 : thetas(i) = (i - 0.5) * pi / num_points
107 : r1s(i) = a1 * (1 - e**2) / (1 + e * cos(thetas(i)))
108 : r2s(i) = a2 / a1 * r1s(i)
109 :
110 : x1s(i) = -r1s(i) * cos(thetas(i)) ! minus to flip orbit
111 : x1s(2 * num_points - i + 1) = x1s(i)
112 : y1s(i) = -r1s(i) * sin(thetas(i))
113 : y1s(2 * num_points - i + 1) = -y1s(i) ! flip y for other side of orbit
114 : x2s(i) = r2s(i) * cos(thetas(i))
115 : x2s(2 * num_points - i + 1) = x2s(i)
116 : y2s(i) = r2s(i) * sin(thetas(i))
117 : y2s(2 * num_points - i + 1) = -y2s(i)
118 : end do
119 : !$OMP END PARALLEL DO
120 0 : x1s(2 * num_points + 1) = x1s(1)
121 0 : y1s(2 * num_points + 1) = y1s(1)
122 0 : x2s(2 * num_points + 1) = x2s(1)
123 0 : y2s(2 * num_points + 1) = y2s(1)
124 :
125 0 : x1max = maxval(abs(x1s))
126 0 : x2max = maxval(abs(x2s))
127 0 : xmax = max(x1max, x2max)
128 :
129 0 : q = b% m(2) / b% m(1)
130 0 : if (b% pg% Orbit_show_stars .and. abs(log10(q)) <= 2) then
131 0 : if (b% point_mass_i /= 1) then
132 0 : this_psi = Psi_fit(b% r(1) / b% separation, q)
133 0 : xl1 = xl1_fit(q)
134 0 : do i = 1, num_points
135 0 : phis(i) = (i - 0.5) * pi / num_points
136 0 : cosp = cos(phis(i))
137 : rs(i) = safe_root_with_guess(f, 1d-1, 1d-3, & ! function, guess, dx for bracket
138 : 1d-3, 1d0 - 1d-3, & ! left, right bracket
139 : roche(1d-3, cosp), roche(1d0 - 1d-3, cosp), & ! f(left, right bracket)
140 : 50, 100, 1d-6, 1d-8, & ! i_next, imax, x_tol, y_tol
141 : 0, rpar, 0, ipar, & ! func_params
142 0 : ierr)
143 :
144 0 : x1s_RL(i) = rs(i) * cosp
145 0 : y1s_RL(i) = rs(i) * sin(phis(i))
146 0 : y1s_RL(2 * num_points - i + 1) = -y1s_RL(i)
147 : end do
148 : icut = 1
149 0 : do while (x1s_RL(icut) > xl1) ! find most extreme x within cut
150 0 : icut = icut + 1
151 : end do
152 0 : do i = 1, icut - 1 ! move points beyond l1 to point close to xl1
153 0 : x1s_RL(i) = x1s_RL(icut)
154 0 : y1s_RL(i) = y1s_RL(icut)
155 : end do
156 0 : do i = 1, num_points ! displace the xs
157 0 : x1s_RL(i) = x1s_RL(i) - a1 * (1 - e)
158 0 : x1s_RL(2 * num_points - i + 1) = x1s_RL(i)
159 : end do
160 0 : x1s_RL(2 * num_points + 1) = x1s_RL(1) ! close contour
161 0 : y1s_RL(2 * num_points + 1) = y1s_RL(1)
162 0 : x1max = maxval(abs(x1s_RL))
163 0 : xmax = max(x1max, xmax)
164 : else
165 0 : x1s_RL = 0d0
166 0 : y1s_RL = 0d0
167 : end if
168 :
169 0 : if (b% point_mass_i /= 2) then
170 0 : q = 1d0 / q ! flip q for other star
171 0 : this_psi = Psi_fit(b% r(2) / b% separation, q)
172 0 : xl1 = xl1_fit(q)
173 0 : do i = 1, num_points
174 0 : phis(i) = (i - 0.5) * pi / num_points
175 0 : cosp = cos(phis(i))
176 : rs(i) = safe_root_with_guess(f, 1d-1, 1d-3, & ! function, guess, dx for bracket
177 : 1d-3, 1d0 - 1d-3, & ! left, right bracket
178 : roche(1d-3, cosp), roche(1d0 - 1d-3, cosp), & ! f(left, right bracket)
179 : 25, 50, 1d-4, 1d-6, & ! i_next, imax, x_tol, y_tol
180 : 0, rpar, 0, ipar, & ! func_params
181 0 : ierr)
182 :
183 0 : x2s_RL(i) = rs(i) * cosp
184 0 : y2s_RL(i) = rs(i) * sin(phis(i))
185 0 : y2s_RL(2 * num_points - i + 1) = -y2s_RL(i)
186 : end do
187 : icut = 1
188 0 : do while (x2s_RL(icut) > xl1) ! find most extreme x within cut
189 0 : icut = icut + 1
190 : end do
191 0 : do i = 1, icut - 1 ! move points beyond l1 to point close to xl1
192 0 : x2s_RL(i) = x2s_RL(icut)
193 0 : y2s_RL(i) = y2s_RL(icut)
194 : end do
195 0 : do i = 1, num_points ! displace the xs
196 0 : x2s_RL(i) = -(x2s_RL(i) - a2 * (1 - e)) ! flip x for 2nd star!
197 0 : x2s_RL(2 * num_points - i + 1) = x2s_RL(i)
198 : end do
199 0 : x2s_RL(2 * num_points + 1) = x2s_RL(1)
200 0 : y2s_RL(2 * num_points + 1) = y2s_RL(1)
201 0 : x2max = maxval(abs(x2s_RL))
202 0 : xmax = max(x2max, xmax)
203 : else
204 0 : x2s_RL = 0d0
205 0 : y2s_RL = 0d0
206 : end if
207 0 : else if (b% pg% Orbit_show_stars .and. abs(log10(q)) > 2) then
208 0 : write(*, 1) "pgbinary: Not plotting RL, q too extreme: abs(log(q)) = ", abs(log10(q))
209 : end if
210 :
211 0 : call pgsave
212 0 : call pgsci(clr_Foreground)
213 0 : call pgscf(1)
214 0 : call pgsch(txt_scale)
215 0 : call pgslw(b% pg% pgbinary_lw / 2)
216 0 : call pgwnad(-1.1 * xmax, 1.1 * xmax, -1.1 * xmax, 1.1 * xmax)
217 0 : call show_box_pgbinary(b, 'BCSTN', 'BCSTNMV')
218 0 : call show_xaxis_label_pgbinary(b, 'separation')
219 0 : call show_left_yaxis_label_pgbinary(b, 'separation')
220 :
221 0 : call pgsci(clr_Goldenrod)
222 0 : call pgline(2 * num_points + 1, x1s, y1s)
223 0 : call pgslw(1)
224 0 : call pgmtxt('T', -2.0, 0.05, 0.0, 'Star 1')
225 0 : call pgsci(clr_LightSkyBlue)
226 0 : call pgslw(b% pg% pgbinary_lw / 2)
227 0 : call pgline(2 * num_points + 1, x2s, y2s)
228 :
229 0 : if (b% pg% Orbit_show_stars .and. abs(log10(q)) <= 2) then
230 0 : call pgslw(int(2.0 * b% pg% pgbinary_lw / 3.0))
231 0 : call pgsfs(3)
232 0 : call pgshs(45.0, 0.33, 0.0)
233 0 : call pgsci(clr_Goldenrod)
234 0 : call pgline(2 * num_points + 1, x1s_RL, y1s_RL)
235 0 : call pgpoly(2 * num_points + 1, x1s_RL, y1s_RL)
236 0 : call pgsci(clr_LightSkyBlue)
237 0 : call pgline(2 * num_points + 1, x2s_RL, y2s_RL)
238 0 : call pgpoly(2 * num_points + 1, x2s_RL, y2s_RL)
239 : end if
240 :
241 0 : call pgslw(1)
242 0 : call pgmtxt('T', -2.0 - 1.3, 0.05, 0.0, 'Star 2')
243 :
244 0 : call pgsci(clr_Foreground)
245 0 : call pgpt1(0.0, 0.0, 5)
246 0 : call pgunsa
247 :
248 : contains
249 :
250 0 : real function xl1_fit(q)
251 : real(dp), intent(in) :: q
252 0 : real(dp) :: logq
253 :
254 0 : logq = log10(q)
255 0 : if (q > 1) logq = -logq
256 : xl1_fit = - 1.72452947 / pi * atan(logq * 0.21625699) + 0.5 &
257 : + 0.01559149 * logq &
258 0 : - 1.3924d-05 * logq * (logq + 1.5) * (logq + 4.0)
259 0 : if (q > 1) xl1_fit = 1 - xl1_fit
260 :
261 0 : end function xl1_fit
262 :
263 0 : real(dp) function Psi_fit(req, q)
264 : ! fit of Roche potential versus q = m_other / m_this and r_eq, the &
265 : ! dimensionless volume equivalent radius (== r / separation of the model)
266 : real(dp), intent(in) :: req, q
267 : Psi_fit = -1d0 / req - q &
268 : - 0.5533 * (1 + q) * req ** 2 &
269 : - 0.3642 * req ** 2 * (req ** 2 - 1) &
270 0 : - 1.8693 * req * (req - 0.1) * (req - 0.3) * (req - 0.7) * (req - 1.0414)
271 0 : end function Psi_fit
272 :
273 0 : real(dp) function roche(r, cosp)
274 : real(dp), intent(in) :: r, cosp
275 :
276 : roche = -1d0 / r &
277 : - q * (pow(1 - 2 * r * cosp + r**2, -0.5d0) - r * cosp) &
278 0 : - (1 + q) / 2 * r**2
279 0 : end function roche
280 :
281 0 : real(dp) function f(r, dfdx, lrpar, rpar, lipar, ipar, ierr)
282 : real(dp), intent(in) :: r
283 : integer, intent(in) :: lrpar, lipar
284 : real(dp), intent(out) :: dfdx
285 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
286 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
287 : integer, intent(out) :: ierr
288 :
289 0 : f = roche(r, cosp) - this_psi
290 : dfdx = 1d0 / r**2 &
291 : + q * (pow(1 - 2 * r * cosp + r**2, -1.5d0) * (r - cosp) + cosp) &
292 0 : - (1 + q) * r
293 0 : ierr = 0
294 0 : end function f
295 :
296 : end subroutine orbit_panel
297 :
298 :
299 : end module pgbinary_orbit
300 :
|
