Line data Source code
1 1 : program test_colors
2 :
3 1 : use colors_def
4 : use colors_lib
5 : use math_lib
6 : use utils_lib, only: mesa_error, mkdir
7 :
8 : implicit none
9 :
10 1 : call do_test_colors
11 :
12 : contains
13 :
14 1 : subroutine do_test_colors
15 1 : use const_lib, only: const_init
16 :
17 : character(len=256) :: my_mesa_dir
18 : integer :: info
19 :
20 : logical, parameter :: do_one = .true.
21 : integer, parameter :: n_colors = 11
22 : ! logical, parameter :: do_one = .false.
23 :
24 1 : my_mesa_dir = '../..'
25 1 : call const_init(my_mesa_dir, info)
26 1 : if (info /= 0) then
27 0 : write (*, *) 'const_init failed'
28 0 : call mesa_error(__FILE__, __LINE__)
29 : end if
30 1 : call math_init()
31 :
32 2 : call colors_init(1, ['../data/lcb98cor.dat'], [n_colors], info)
33 :
34 1 : if (info /= 0) then
35 0 : write (*, *) 'colors_init failed during initialization'
36 0 : return
37 : end if
38 :
39 1 : if (do_one) then
40 : call do_one_colors
41 : else
42 : call create_plot_files
43 : end if
44 1 : call colors_shutdown()
45 :
46 : end subroutine do_test_colors
47 :
48 1 : subroutine do_one_colors
49 :
50 : integer, parameter :: num_results = 16
51 1 : real(dp) :: log_teff ! log10 of surface temp
52 1 : real(dp) :: log_l ! log10 of luminosity in solar units
53 1 : real(dp) :: mass ! mass in solar units
54 1 : real(dp) :: M_div_h ! [M/h], or as an approximation, log10[z/zsun]
55 1 : real(dp) :: boloMag
56 17 : real(dp), dimension(num_results) :: results
57 1 : real(dp) ::log_g, x
58 : integer :: info, i
59 : character(len=8) :: vname
60 : character(len=strlen), dimension(num_results):: colors_name
61 : logical, parameter :: doing_solar = .true.
62 : real(dp), dimension(num_results), parameter :: solar_expected_results = [ &
63 : 4.75d0, -0.11510d0, -0.14211d0, -0.61768d0, -0.36199d0, -0.68894d0, &
64 : -1.46926d0, -0.32695d0, -0.78032d0, -0.39024d0, 0.05223d0, &
65 : -0.10512d0, -0.33801d0, -0.44312d0, -0.44123d0, -0.43080d0]
66 :
67 : ! solar values
68 1 : log_teff = log10(5780d0)
69 1 : log_l = 0d0
70 1 : mass = 1d0
71 1 : M_div_h = 0d0
72 1 : log_g = log10(mass) + 4.0D0*log_Teff - log_L - 10.6071D0
73 :
74 1 : boloMag = get_abs_bolometric_mag(10**log_l)
75 :
76 : !Store answers in results array
77 :
78 : ! These are bololmetric correction differences NOT magnitude differences
79 : ! Thus they are -1*mag diff
80 1 : results(1) = boloMag
81 1 : results(2) = get_bc_by_name('V', log_Teff, log_g, M_div_h, info)
82 1 : results(3) = get_bc_by_name('U', log_Teff, log_g, M_div_h, info) - get_bc_by_name('B', log_Teff, log_g, M_div_h, info)
83 1 : results(4) = get_bc_by_name('B', log_Teff, log_g, M_div_h, info) - get_bc_by_name('V', log_Teff, log_g, M_div_h, info)
84 1 : results(5) = get_bc_by_name('V', log_Teff, log_g, M_div_h, info) - get_bc_by_name('R', log_Teff, log_g, M_div_h, info)
85 1 : results(6) = get_bc_by_name('V', log_Teff, log_g, M_div_h, info) - get_bc_by_name('I', log_Teff, log_g, M_div_h, info)
86 1 : results(7) = get_bc_by_name('V', log_Teff, log_g, M_div_h, info) - get_bc_by_name('K', log_Teff, log_g, M_div_h, info)
87 1 : results(8) = get_bc_by_name('R', log_Teff, log_g, M_div_h, info) - get_bc_by_name('I', log_Teff, log_g, M_div_h, info)
88 1 : results(9) = get_bc_by_name('I', log_Teff, log_g, M_div_h, info) - get_bc_by_name('K', log_Teff, log_g, M_div_h, info)
89 1 : results(10) = get_bc_by_name('J', log_Teff, log_g, M_div_h, info) - get_bc_by_name('H', log_Teff, log_g, M_div_h, info)
90 1 : results(11) = get_bc_by_name('H', log_Teff, log_g, M_div_h, info) - get_bc_by_name('K', log_Teff, log_g, M_div_h, info)
91 1 : results(12) = get_bc_by_name('K', log_Teff, log_g, M_div_h, info) - get_bc_by_name('L', log_Teff, log_g, M_div_h, info)
92 1 : results(13) = get_bc_by_name('J', log_Teff, log_g, M_div_h, info) - get_bc_by_name('K', log_Teff, log_g, M_div_h, info)
93 1 : results(14) = get_bc_by_name('J', log_Teff, log_g, M_div_h, info) - get_bc_by_name('L', log_Teff, log_g, M_div_h, info)
94 1 : results(15) = get_bc_by_name('J', log_Teff, log_g, M_div_h, info) - get_bc_by_name('Lprime', log_Teff, log_g, M_div_h, info)
95 1 : results(16) = get_bc_by_name('K', log_Teff, log_g, M_div_h, info) - get_bc_by_name('M', log_Teff, log_g, M_div_h, info)
96 :
97 1 : if (info /= 0) then
98 0 : call mesa_error(__FILE__, __LINE__, 'bad return from colors_get')
99 : end if
100 :
101 1 : write (*, '(A)')
102 1 : write (*, *) 'color magnitude results'
103 1 : write (*, '(A)')
104 1 : write (*, '(6a12)') 'teff', 'log_teff', 'log_l', 'mass', '[M_div_h]', 'log_g'
105 1 : write (*, '(i12,5f12.2)') floor(exp10(log_teff) + 0.5d0), log_teff, log_l, mass, M_div_h, log_g
106 1 : write (*, '(A)')
107 :
108 : !call get_all_bc_names(colors_name,total_num_colors,info)
109 : colors_name = ['bol ', 'bcv ', 'U-B ', 'B-V ', 'V-R ', 'V-I ', 'V-K ', 'R-I ', &
110 17 : 'I-K ', 'J-H ', 'H-K ', 'K-L ', 'J-K ', 'J-L ', 'J-Lp', 'K-M ']
111 :
112 17 : do i = 1, num_results
113 17 : write (*, '(9x,a8,f10.5)') colors_name(i), results(i)
114 : end do
115 1 : write (*, '(A)')
116 1 : vname = 'vcolors'
117 1 : write (*, '(9x,a8,f10.5)') vname, results(1) - results(2)
118 1 : write (*, '(A)')
119 :
120 : if (doing_solar) then
121 :
122 17 : do i = 1, num_results
123 17 : if (abs(results(i) - solar_expected_results(i)) > 0.02d0) then
124 0 : write (*, '(A)')
125 0 : write (*, *) "warning colors don't match solar expected values"
126 0 : write (*, *) "Color: ", trim(colors_name(i))
127 0 : write (*, *) "Got ", results(i), ' but expected ', solar_expected_results(i)
128 0 : write (*, '(A)')
129 0 : stop
130 : end if
131 : end do
132 :
133 1 : write (*, *) 'matches expected solar results'
134 1 : write (*, '(A)')
135 :
136 : end if
137 :
138 : !Check some extreme values
139 1 : x = get_bc_by_name('V', 10.d0, log_g, M_div_h, info)
140 1 : write (*, '(A,1pes40.16e3)') 'high logT', x
141 1 : x = get_bc_by_name('V', 0.d0, log_g, M_div_h, info)
142 1 : write (*, '(A,1pes40.16e3)') 'low logT', x
143 1 : x = get_bc_by_name('V', log_teff, -10.d0, M_div_h, info)
144 1 : write (*, '(A,1pes40.16e3)') 'low logg', x
145 1 : x = get_bc_by_name('V', log_teff, 10.d0, M_div_h, info)
146 1 : write (*, '(A,1pes40.16e3)') 'high logg', x
147 1 : x = get_bc_by_name('V', log_teff, -10.d0, 1.0d0, info)
148 1 : write (*, '(A,1pes40.16e3)') 'high m/h', x
149 :
150 1 : end subroutine do_one_colors
151 :
152 : subroutine create_plot_files
153 :
154 : character(len=256) :: fname, dir
155 : integer, parameter :: max_num_masses = 100
156 : integer, parameter :: num_results = 16
157 : real(dp), dimension(max_num_masses) :: mass, logl, logteff
158 : real(dp) :: read_junk, log_g, M_div_h, boloMag
159 : real(dp), dimension(num_results) :: results
160 : integer :: info, i, num_masses, io_unit, ios, iread_junk
161 :
162 : M_div_h = 0d0
163 : fname = 'zams_data/z02.log'
164 : io_unit = 40
165 : open (unit=io_unit, file=trim(fname), action='read', status='old', iostat=ios)
166 : if (ios /= 0) then
167 : write (*, *) 'failed to open the zams data'
168 : return
169 : end if
170 :
171 : num_masses = 0
172 : do i = 1, max_num_masses
173 : read (io_unit, fmt=*, iostat=ios) iread_junk, mass(i), logl(i), read_junk, logteff(i)
174 : if (ios /= 0) then
175 : num_masses = i - 1
176 : exit
177 : end if
178 : end do
179 : read_junk = read_junk; iread_junk = iread_junk ! to keep g95 quiet
180 :
181 : close (io_unit)
182 :
183 : dir = 'plot_data'
184 : call mkdir(dir)
185 : fname = trim(dir)//'/'//'colors.data'
186 : open (unit=io_unit, file=trim(fname))
187 :
188 : do i = 1, num_masses
189 :
190 : !call colors_get(logteff(i), logl(i), mass(i), M_div_h, results, log_g, info)
191 :
192 : boloMag = get_abs_bolometric_mag(10**logl(i))
193 :
194 : log_g = log10(mass(i)) + 4.0D0*logteff(i) - logl(i) - 10.6071D0
195 : !Store answers in results array
196 : results(1) = boloMag
197 : results(2) = get_bc_by_name('V', logteff(i), log_g, M_div_h, info)
198 : results(3) = get_bc_by_name('U', logteff(i), log_g, M_div_h, info) - get_bc_by_name('B', logteff(i), log_g, M_div_h, info)
199 : results(4) = get_bc_by_name('B', logteff(i), log_g, M_div_h, info) - get_bc_by_name('V', logteff(i), log_g, M_div_h, info)
200 : results(5) = get_bc_by_name('V', logteff(i), log_g, M_div_h, info) - get_bc_by_name('R', logteff(i), log_g, M_div_h, info)
201 : results(6) = get_bc_by_name('V', logteff(i), log_g, M_div_h, info) - get_bc_by_name('I', logteff(i), log_g, M_div_h, info)
202 : results(7) = get_bc_by_name('V', logteff(i), log_g, M_div_h, info) - get_bc_by_name('K', logteff(i), log_g, M_div_h, info)
203 : results(8) = get_bc_by_name('R', logteff(i), log_g, M_div_h, info) - get_bc_by_name('I', logteff(i), log_g, M_div_h, info)
204 : results(9) = get_bc_by_name('I', logteff(i), log_g, M_div_h, info) - get_bc_by_name('K', logteff(i), log_g, M_div_h, info)
205 : results(10) = get_bc_by_name('J', logteff(i), log_g, M_div_h, info) - get_bc_by_name('H', logteff(i), log_g, M_div_h, info)
206 : results(11) = get_bc_by_name('H', logteff(i), log_g, M_div_h, info) - get_bc_by_name('K', logteff(i), log_g, M_div_h, info)
207 : results(12) = get_bc_by_name('K', logteff(i), log_g, M_div_h, info) - get_bc_by_name('L', logteff(i), log_g, M_div_h, info)
208 : results(13) = get_bc_by_name('J', logteff(i), log_g, M_div_h, info) - get_bc_by_name('K', logteff(i), log_g, M_div_h, info)
209 : results(14) = get_bc_by_name('J', logteff(i), log_g, M_div_h, info) - get_bc_by_name('L', logteff(i), log_g, M_div_h, info)
210 : results(15) = &
211 : get_bc_by_name('J', logteff(i), log_g, M_div_h, info) - get_bc_by_name('Lprime', logteff(i), log_g, M_div_h, info)
212 : results(16) = get_bc_by_name('K', logteff(i), log_g, M_div_h, info) - get_bc_by_name('M', logteff(i), log_g, M_div_h, info)
213 :
214 : if (info == 0) write (io_unit, '(99f15.8)') mass(i), logl(i), logteff(i), log_g, results
215 :
216 : end do
217 :
218 : close (io_unit)
219 :
220 : write (*, *) 'finished creating plot files'
221 :
222 : end subroutine create_plot_files
223 :
224 : end program test_colors
|
