Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 colors_lib
21 : use math_lib
22 : use colors_def
23 : use const_def, only : dp, strlen, mbolsun, loggsun, Teffsun
24 : ! library for calculating theoretical estimates of magnitudes and colors
25 : ! from Teff, L, M, and [M/H].
26 :
27 : ! Color-magnitude data shipped with MESA is from:
28 : ! Lejeune, Cuisinier, Buser (1998) A&AS 130, 65-75.
29 : ! However, you add your own bolometric corrections files for mesa to use
30 :
31 : ! The data interface for the library is defined in colors_def
32 : ! Th easiest way to get output is to add the columns to your history_columns.list file
33 :
34 : ! The preferred way for users (in a run_star_extras routine) for accessing the colors data is to
35 : ! call either get_by_by_name, get_abs_mag_by_name or get_abs_bolometric_mag. Other routines are there
36 : ! to hook into the rest of MESA.
37 :
38 : ! Routines get_bc will return the coefficients from interpolating over log Teff, log g, [M/H]
39 : ! even though the tables are defined as Teff, log g, [M/H]. get_abs_mag routines return
40 : ! data that's been turned into an absolute magnitude. A color can be computed by taking the difference between
41 : ! two get_bc or two get_abs_mag calls.
42 :
43 : ! Names for the filters should be unique across all data files (left to the user to enforce this).
44 : ! Name matching is performed in a case sensitive manner.
45 : ! The names themselves are not important as far as MESA is concerned, you can name each filter (including the
46 : ! ones MESA ships by defaults) by what ever name you want by editing the data file(s) and changing the names in the header.
47 : ! MESA does not rely on any particlaur band existing.
48 :
49 : implicit none
50 :
51 :
52 : contains ! the procedure interface for the library
53 : ! client programs should only call these routines.
54 :
55 :
56 2 : subroutine colors_init(num_files,fnames,num_colors,ierr)
57 : use mod_colors, only : do_colors_init
58 : integer, intent(in) :: num_files
59 : integer, dimension(:), intent(in) :: num_colors
60 : character(len=*), dimension(:), intent(in) :: fnames
61 : integer, intent(out) :: ierr
62 :
63 2 : ierr=0
64 :
65 4 : !$OMP critical (color_init)
66 2 : if (.not. color_is_initialized) then
67 2 : call do_colors_init(num_files,fnames,num_colors,ierr)
68 : end if
69 : !$OMP end critical (color_init)
70 :
71 2 : if(ierr/=0)THEN
72 0 : ierr=-1
73 0 : write(*,*) "colors_init failed"
74 : return
75 : end if
76 :
77 2 : end subroutine colors_init
78 :
79 :
80 2 : subroutine colors_shutdown ()
81 :
82 2 : use mod_colors, only : free_colors_all
83 :
84 2 : if (.not. color_is_initialized) return
85 :
86 2 : call free_colors_all()
87 :
88 2 : color_is_initialized = .FALSE.
89 :
90 2 : end subroutine colors_shutdown
91 :
92 :
93 34 : subroutine get_bcs_one(log_Teff, log_g, M_div_h, results, thead,n_colors, ierr)
94 2 : use mod_colors, only : Eval_Colors
95 : ! input
96 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
97 : real(dp), intent(in) :: M_div_H ! [M/H]
98 : ! output
99 : real(dp), dimension(:), intent(out) :: results
100 : real(dp), intent(in) :: log_g
101 : integer, intent(in) :: n_colors
102 : integer, intent(inout) :: ierr
103 : type (lgt_list),intent(inout), pointer :: thead
104 :
105 714 : results(:)=-99.9d0
106 : ierr=0
107 :
108 34 : if (.not. color_is_initialized) then
109 0 : ierr=-1
110 : return
111 : end if
112 :
113 34 : call Eval_Colors(log_Teff, log_g, M_div_h, results,thead,n_colors, ierr)
114 :
115 34 : end subroutine get_bcs_one
116 :
117 34 : real(dp) function get_bc_by_name(name,log_Teff,log_g, M_div_h, ierr)
118 : ! input
119 : character(len=*),intent(in) :: name
120 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
121 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
122 : real(dp), intent(in) :: M_div_h ! [M/H]
123 714 : real(dp), dimension(max_num_bcs_per_file) :: results
124 : type (lgt_list), pointer :: thead => null()
125 : integer, intent(inout) :: ierr
126 : integer :: i,j,n_colors
127 :
128 34 : get_bc_by_name=-99.9d0
129 34 : ierr=0
130 :
131 34 : if (.not. color_is_initialized) then
132 0 : ierr=-1
133 0 : return
134 : end if
135 :
136 34 : do i=1,num_thead
137 34 : thead=>thead_all(i)%thead
138 34 : n_colors=thead_all(i)%n_colors
139 183 : do j=1,n_colors
140 : if(trim(name)==trim(thead_all(i)%color_names(j)).or. &
141 : trim(name)=='bc_'//trim(thead_all(i)%color_names(j)).or. &
142 : trim(name)=='abs_mag_'//trim(thead_all(i)%color_names(j)).or. &
143 183 : trim(name)=='lum_band_'//trim(thead_all(i)%color_names(j)).or. &
144 : trim(name)=='log_lum_band_'//trim(thead_all(i)%color_names(j))&
145 0 : ) then
146 :
147 34 : call get_bcs_one(log_Teff,log_g, M_div_h, results,thead,n_colors, ierr)
148 68 : if(ierr/=0) return
149 :
150 34 : get_bc_by_name=results(j)
151 :
152 34 : return
153 : end if
154 : end do
155 :
156 : end do
157 :
158 :
159 34 : end function get_bc_by_name
160 :
161 0 : real(dp) function get_bc_by_id(id,log_Teff,log_g, M_div_h, ierr)
162 : ! input
163 : integer, intent(in) :: id
164 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
165 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
166 : real(dp), intent(in) :: M_div_h ! [M/H]
167 : integer, intent(inout) :: ierr
168 : character(len=strlen) :: name
169 :
170 0 : get_bc_by_id=-99.9d0
171 0 : ierr=0
172 :
173 0 : if (.not. color_is_initialized) then
174 0 : ierr=-1
175 0 : return
176 : end if
177 :
178 0 : name=get_bc_name_by_id(id,ierr)
179 0 : if(ierr/=0) return
180 :
181 0 : get_bc_by_id=get_bc_by_name(name,log_Teff,log_g, M_div_h, ierr)
182 :
183 0 : end function get_bc_by_id
184 :
185 0 : integer function get_bc_id_by_name(name,ierr)
186 : ! input
187 : character(len=*), intent(in) :: name
188 : integer, intent(inout) :: ierr
189 : integer :: i,j,k
190 :
191 0 : get_bc_id_by_name=-1
192 0 : ierr=0
193 :
194 0 : if (.not. color_is_initialized) then
195 0 : ierr=-1
196 0 : return
197 : end if
198 :
199 0 : k=0
200 0 : do i=1,num_thead
201 0 : do j=1,thead_all(i)%n_colors
202 0 : k=k+1
203 : if(trim(name)==trim(thead_all(i)%color_names(j)).or. &
204 0 : trim(name)=='bc_'//trim(thead_all(i)%color_names(j)).or. &
205 0 : trim(name)=='abs_mag_'//trim(thead_all(i)%color_names(j))) then
206 0 : get_bc_id_by_name=k
207 0 : return
208 : end if
209 : end do
210 : end do
211 :
212 : end function get_bc_id_by_name
213 :
214 0 : character(len=strlen) function get_bc_name_by_id(id,ierr)
215 : ! input
216 : integer, intent(in) :: id
217 : integer, intent(inout) :: ierr
218 : integer :: i,j,k
219 :
220 0 : get_bc_name_by_id=''
221 0 : ierr=0
222 :
223 0 : if (.not. color_is_initialized) then
224 0 : ierr=-1
225 0 : return
226 : end if
227 :
228 0 : k=1
229 0 : do i=1,num_thead
230 0 : do j=1,thead_all(i)%n_colors
231 0 : if(k==id) then
232 0 : get_bc_name_by_id=thead_all(i)%color_names(j)
233 0 : return
234 : end if
235 0 : k=k+1
236 : end do
237 : end do
238 :
239 : end function get_bc_name_by_id
240 :
241 1 : real(dp) function get_abs_bolometric_mag(lum)
242 : use const_def, only: dp
243 : real(dp), intent(in) :: lum ! Luminsoity in lsun units
244 :
245 1 : get_abs_bolometric_mag = mbolsun - 2.5d0*log10(lum)
246 :
247 1 : end function get_abs_bolometric_mag
248 :
249 0 : real(dp) function get_abs_mag_by_name(name,log_Teff,log_g, M_div_h,lum, ierr)
250 : ! input
251 : character(len=*) :: name
252 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
253 : real(dp), intent(in) :: M_div_h ! [M/H]
254 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
255 : real(dp), intent(in) :: lum ! Luminsoity in lsun units
256 : integer, intent(inout) :: ierr
257 :
258 0 : ierr=0
259 0 : get_abs_mag_by_name=-99.9d0
260 :
261 0 : if (.not. color_is_initialized) then
262 0 : ierr=-1
263 0 : return
264 : end if
265 :
266 : get_abs_mag_by_name=get_abs_bolometric_mag(lum)-&
267 0 : get_bc_by_name(name,log_Teff,log_g, M_div_h,ierr)
268 :
269 0 : end function get_abs_mag_by_name
270 :
271 0 : real(dp) function get_abs_mag_by_id(id,log_Teff,log_g, M_div_h,lum, ierr)
272 : ! input
273 : integer, intent(in) :: id
274 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
275 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
276 : real(dp), intent(in) :: M_div_h ! [M/H]
277 : real(dp), intent(in) :: lum ! Luminsoity in lsun units
278 : integer, intent(inout) :: ierr
279 : character(len=strlen) :: name
280 :
281 0 : ierr=0
282 0 : get_abs_mag_by_id=-99.9d0
283 :
284 0 : if (.not. color_is_initialized) then
285 0 : ierr=-1
286 0 : return
287 : end if
288 :
289 0 : name=get_bc_name_by_id(id,ierr)
290 0 : if(ierr/=0) return
291 :
292 0 : get_abs_mag_by_id=get_abs_mag_by_name(name,log_Teff,log_g, M_div_h,lum, ierr)
293 :
294 0 : end function get_abs_mag_by_id
295 :
296 0 : subroutine get_all_bc_names(names, ierr)
297 : character(len=strlen),dimension(:) :: names
298 : integer, intent(inout) :: ierr
299 : integer :: i,j,cnt
300 :
301 0 : names(:)=''
302 :
303 0 : if (.not. color_is_initialized) then
304 0 : ierr=-1
305 0 : return
306 : end if
307 :
308 0 : cnt=1
309 0 : do i=1,num_thead
310 0 : do j=1,thead_all(i)%n_colors
311 0 : names(cnt)=trim(thead_all(i)%color_names(j))
312 0 : cnt=cnt+1
313 : end do
314 : end do
315 :
316 : end subroutine get_all_bc_names
317 :
318 0 : subroutine get_bcs_all(log_Teff, log_g, M_div_h, results, ierr)
319 : ! input
320 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
321 : real(dp), intent(in) :: M_div_h ! [M/H]
322 : ! output
323 : real(dp),dimension(:), intent(out) :: results
324 : real(dp), intent(in) :: log_g
325 : integer, intent(inout) :: ierr
326 : type (lgt_list), pointer :: thead => null()
327 : integer :: i,iStart,iEnd
328 :
329 0 : ierr=0
330 0 : results(:)=-99.d0
331 :
332 0 : if (.not. color_is_initialized) then
333 0 : ierr=-1
334 0 : return
335 : end if
336 :
337 0 : do i=1,num_thead
338 0 : thead=>thead_all(i)%thead
339 0 : iStart=(i-1)*thead_all(i)%n_colors+1
340 0 : iEnd=i*thead_all(i)%n_colors
341 0 : call get_bcs_one(log_Teff, log_g, M_div_h, results(iStart:iEnd),thead,thead_all(i)%n_colors, ierr)
342 0 : if(ierr/=0) return
343 : end do
344 :
345 : end subroutine get_bcs_all
346 :
347 : !Returns in lsun units
348 0 : real(dp) function get_lum_band_by_name(name,log_Teff,log_g, M_div_h, lum, ierr)
349 : ! input
350 : character(len=*) :: name
351 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
352 : real(dp), intent(in) :: M_div_h ! [M/H]
353 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
354 : real(dp), intent(in) :: lum ! Total luminsoity in lsun units
355 0 : real(dp) :: solar_abs_mag, star_abs_mag
356 : integer, intent(inout) :: ierr
357 :
358 0 : ierr=0
359 0 : get_lum_band_by_name=-99.d0
360 :
361 0 : if (.not. color_is_initialized) then
362 0 : ierr=-1
363 0 : return
364 : end if
365 :
366 : ! Filter dependent terms
367 0 : solar_abs_mag=get_abs_mag_by_name(name, safe_log10(Teffsun), loggsun, 0.d0, 1.d0, ierr)
368 0 : if(ierr/=0) return
369 :
370 0 : star_abs_mag=get_abs_mag_by_name(name, log_Teff, log_g, M_div_h, lum, ierr)
371 0 : if(ierr/=0) return
372 :
373 0 : get_lum_band_by_name=exp10((star_abs_mag-solar_abs_mag)/(-2.5d0))
374 :
375 0 : end function get_lum_band_by_name
376 :
377 : !Returns in lsun units
378 0 : real(dp) function get_lum_band_by_id(id,log_Teff,log_g, M_div_h, lum, ierr)
379 : ! input
380 : integer, intent(in) :: id
381 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
382 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
383 : real(dp), intent(in) :: M_div_h ! [M/H]
384 : real(dp), intent(in) :: lum ! Total luminsoity in lsun units
385 0 : real(dp) :: solar_abs_mag, star_abs_mag
386 : integer, intent(inout) :: ierr
387 :
388 0 : ierr=0
389 0 : get_lum_band_by_id=-99.d0
390 :
391 0 : if (.not. color_is_initialized) then
392 0 : ierr=-1
393 0 : return
394 : end if
395 :
396 : ! Filter dependent terms
397 0 : solar_abs_mag=get_abs_mag_by_id(id, safe_log10(Teffsun), loggsun, 0.d0, 1.d0, ierr)
398 0 : if(ierr/=0) return
399 :
400 0 : star_abs_mag=get_abs_mag_by_id(id, log_Teff, log_g, M_div_h,lum, ierr)
401 0 : if(ierr/=0) return
402 :
403 0 : get_lum_band_by_id=exp10((star_abs_mag-solar_abs_mag)/(-2.5d0))
404 :
405 0 : end function get_lum_band_by_id
406 :
407 : end module colors_lib
408 :
|
