Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2022 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 : !!% This module contains code to interpolate gravity darkening coefficients.
21 : !!% The coefficients scale surface-averaged Teff and Luminosity as a function
22 : !!% of the orientation of the observer w.r.t. the rotation axis of the star
23 : !!% and the ratio of the surface angular velocity to the Keplerian angular
24 : !!% velocity. The orientation angle is denoted "inclination" and the angular
25 : !!% velocity is denoted "omega". Permissible ranges are
26 : !!%
27 : !!% 0 (no rotation) <= omega <= 1 (critical rotation)
28 : !!% 0 (equator) <= inclination in radians <= pi/2 (pole)
29 : !!%
30 : !!%
31 : !!% The coefficients are obtained via 2D interpolation using the interp_2d
32 : !!% module in tables. The tables were computed using the code at
33 : !!%
34 : !!% https://github.com/aarondotter/GDit
35 : !!%
36 : !!% where the reader can also find documentation on how the calculations are
37 : !!% performed. The table included in mesa/data/star_data consists of a square
38 : !!% array with 21 values each of omega and inclination, equally spaced between
39 : !!% the limits quoted above. The interpolated results fall with roughly 0.001
40 : !!% of exact values calculated using the code referenced above over the full
41 : !!% range.
42 :
43 : module gravity_darkening
44 :
45 : use interp_2d_lib_db
46 : use utils_lib
47 :
48 : implicit none
49 :
50 : logical :: GD_initialized = .false.
51 :
52 : integer, parameter :: ndim = 21
53 : integer, parameter :: size_per_point = 4
54 :
55 : real(dp) :: omega_grid(ndim), inclination_grid(ndim)
56 : real(dp), target :: C_T_ary(size_per_point * ndim * ndim)
57 : real(dp), target :: C_L_ary(size_per_point * ndim * ndim)
58 : real(dp), pointer :: C_T(:,:,:), C_T1(:), C_L(:,:,:), C_L1(:)
59 : character(len=256) :: coefficient_filename
60 :
61 : integer :: ibcxmin, ibcxmax, ibcymin, ibcymax
62 : real(dp) :: bcxmin(ndim), bcxmax(ndim), bcymin(ndim), bcymax(ndim)
63 :
64 : private
65 : public :: gravity_darkening_Teff_coeff, gravity_darkening_L_coeff
66 :
67 : contains
68 :
69 0 : subroutine GD_init(ierr)
70 : use const_def, only: mesa_data_dir
71 : integer, intent(out) :: ierr
72 0 : real(dp) :: dummy(4)
73 : integer :: io, i, j, ilinx, iliny
74 :
75 0 : C_T1 => C_T_ary
76 0 : C_T(1:size_per_point,1:ndim,1:ndim) => C_T_ary(1:size_per_point*ndim*ndim)
77 :
78 0 : C_L1 => C_L_ary
79 0 : C_L(1:size_per_point,1:ndim,1:ndim) => C_L_ary(1:size_per_point*ndim*ndim)
80 :
81 0 : ibcxmin = 0; bcxmin = 0
82 0 : ibcxmax = 0; bcxmax = 0
83 0 : ibcymin = 0; bcymin = 0
84 0 : ibcymax = 0; bcymax = 0
85 :
86 0 : coefficient_filename = trim(mesa_data_dir) // '/star_data/gravity_darkening_coefficients.data'
87 :
88 0 : open(newunit=io,file=trim(coefficient_filename),status='old',action='read')
89 0 : read(io,*) !skip header
90 0 : do i=1,ndim
91 0 : do j=1,ndim
92 0 : read(io,*) dummy(1:4)
93 0 : if(i==1) inclination_grid(j) = dummy(2)
94 0 : if(j==1) omega_grid(i) = dummy(1)
95 0 : C_T(1,i,j) = dummy(3)
96 0 : C_L(1,i,j) = dummy(4)
97 : end do
98 : end do
99 0 : close(io)
100 :
101 : ! construct interpolant for C_T
102 : call interp_mkbicub_db( omega_grid, ndim, inclination_grid, ndim, C_T1, ndim, &
103 : ibcxmin, bcxmin, ibcxmax, bcxmax, ibcymin, bcymin, ibcymax, bcymax, &
104 0 : ilinx, iliny, ierr)
105 :
106 : !construct interpolant for C_L
107 : call interp_mkbicub_db( omega_grid, ndim, inclination_grid, ndim, C_L1, ndim, &
108 : ibcxmin, bcxmin, ibcxmax, bcxmax, ibcymin, bcymin, ibcymax, bcymax, &
109 0 : ilinx, iliny, ierr)
110 :
111 0 : if(ierr==0) GD_initialized = .true.
112 0 : end subroutine GD_init
113 :
114 :
115 0 : function GD_coeff(omega,inclination,C1) result(coeff)
116 : use const_def, only: pi
117 : real(dp), intent(in) :: omega, inclination
118 : real(dp), pointer :: C1(:)
119 0 : real(dp) :: coeff, coeff_eval(6), safe_omega, safe_incl
120 : integer, parameter :: ict(6)=[1,0,0,0,0,0]
121 : integer :: ierr
122 0 : if(.not.GD_initialized) call GD_init(ierr)
123 : !ensure that omega and inclination are within table bounds
124 0 : safe_omega = min(max(omega,0.0d0),1.0d0)
125 0 : safe_incl = min(max(inclination,0.0d0), 0.5d0*pi)
126 : call interp_evbicub_db(safe_omega, safe_incl, omega_grid, ndim, inclination_grid, ndim, &
127 0 : 1, 1, C1, ndim, ict, coeff_eval, ierr)
128 0 : if(ierr==0)then
129 0 : coeff = coeff_eval(1)
130 : else
131 : coeff = 1.0d0
132 : end if
133 0 : end function GD_coeff
134 :
135 :
136 0 : function gravity_darkening_Teff_coeff(omega,inclination) result(Teff_coeff)
137 : real(dp), intent(in) :: omega, inclination
138 : real(dp) :: Teff_coeff
139 0 : Teff_coeff = GD_coeff(omega,inclination,C_T1)
140 0 : end function gravity_darkening_Teff_coeff
141 :
142 :
143 0 : function gravity_darkening_L_coeff(omega,inclination) result(L_coeff)
144 : real(dp), intent(in) :: omega, inclination
145 : real(dp) :: L_coeff
146 0 : L_coeff = GD_coeff(omega,inclination,C_L1)
147 0 : end function gravity_darkening_L_coeff
148 :
149 : end module gravity_darkening
|
