[mesa-users] A new parametrization of the surface effect
Warrick Ball
wball at astro.physik.uni-goettingen.de
Wed Aug 6 05:37:05 EDT 2014
Hi all,
Our recently accepted paper on surface effects appeared on arXiv this
morning:
"A new correction of stellar oscillation frequencies
for near-surface effects"
http://arxiv.org/abs/1408.0986
In it, we argue for the formula frequency^3/(mode inertia) as a
good-fitting functional form for the surface effect. (i.e. the systematic
difference between modelled and observed oscillation frequencies.) An
even better fit is found when there's an extra term proportional to
frequency^(-1)/(mode inertia), but, though it fits the difference with the
Sun quite well, the extra term doesn't really improve the fit for the
distant star we tried, HD 52265. (The same goes for other stars that I've
been fitting between submission and acceptance of the paper.)
We made use of MESA models, and the fit was included by modifying the
surface correction subroutines in astero_support.f. In the spirit of the
open sourceness of MESA, I've attached the modified versions that I wrote,
as well as the tweaked makefile. The .cube file corresponds to the
"cubic" surface term (eqn 3 in the paper), and the .both file to the
"combined" term (eqn 4). In short, I simply replaced the get_freq_corr
subroutine with my own routine, which calculates the best-fitting
coefficients for the surface term (in a least squares sense). I put the
modified astero_support.f files in the src/ directory, hence the
modification to the makefile.
Watch out though: I also modified the chi2 function because I didn't want
a weighted mean of reduced chi-squareds from different sources, which is
the default. If you use these files, you may wish to remove the changes
to the chi2 subroutines.
I'm no master coder, and my snippet is ugly, brutal, minimally-functional
code. When I have a bit of time, I intend to add an option for which
correction to use. e.g. to make an option in &astero_search_controls
something like
which_surface_correction = 'none' ! default
! options are:
! 'none' No frequency correction
! 'power_law' Correction of Kjeldsen et al. (2008)
! 'cubic' Cubic correction of Ball & Gizon (2014, eqn 3)
! 'combined' Combined correction of Ball & Gizon (2014, eq 4)
so I don't have to choose the correction by using differently compiled
binaries. I'll update the list if/when I get that far. I'm also
interested in implementing other new options that might be in the works.
All comments welcome!
Cheers,
Warrick
------------
Warrick Ball
Postdoc, Institut für Astrophysik Göttingen
wball at astro.physik.uni-goettingen.de
+49 (0) 551 39 5069
-------------- next part --------------
! ***********************************************************************
!
! Copyright (C) 2013 Bill Paxton
!
! this file is part of mesa.
!
! mesa is free software; you can redistribute it and/or modify
! it under the terms of the gnu general library public license as published
! by the free software foundation; either version 2 of the license, or
! (at your option) any later version.
!
! mesa is distributed in the hope that it will be useful,
! but without any warranty; without even the implied warranty of
! merchantability or fitness for a particular purpose. see the
! gnu library general public license for more details.
!
! you should have received a copy of the gnu library general public license
! along with this software; if not, write to the free software
! foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa
!
! ***********************************************************************
module astero_support
use astero_data
use star_lib
use star_def
use const_def
use crlibm_lib
use utils_lib
implicit none
contains
subroutine check_search_controls(ierr)
integer, intent(out) :: ierr
integer :: i
include 'formats'
ierr = 0
do i=2,nl0
if (l0_obs(i) <= l0_obs(i-1)) then
write(*,3) 'l0_obs values out of order', i-1, i, l0_obs(i-1), l0_obs(i)
ierr = -1
end if
end do
do i=2,nl1
if (l1_obs(i) <= l1_obs(i-1)) then
write(*,3) 'l1_obs values out of order', i-1, i, l1_obs(i-1), l1_obs(i)
ierr = -1
end if
end do
do i=2,nl2
if (l2_obs(i) <= l2_obs(i-1)) then
write(*,3) 'l2_obs values out of order', i-1, i, l2_obs(i-1), l2_obs(i)
ierr = -1
end if
end do
do i=2,nl3
if (l3_obs(i) <= l3_obs(i-1)) then
write(*,3) 'l3_obs values out of order', i-1, i, l3_obs(i-1), l3_obs(i)
ierr = -1
end if
end do
if (ierr /= 0) &
write(*,1) 'please put frequency values in ascending order'
end subroutine check_search_controls
subroutine get_one_el_info( &
s, l, nu1, nu2, iscan, i1, i2, store_model, code, ierr)
use num_lib, only: qsort
use adipls_support
use gyre_support
type (star_info), pointer :: s
integer, intent(in) :: l, iscan, i1, i2
real(dp), intent(in) :: nu1, nu2
logical, intent(in) :: store_model
character (len=*), intent(in) :: code
integer, intent(out) :: ierr
real(dp) :: nu_obs, dist_j, nu, dist, min_dist, min_freq, &
R, G, M, sig_fac, b, sum_1, sum_2, sum_3
integer :: min_dist_j, min_order, n, cnt
integer :: nsel, itrsig, nsig
real(dp) :: els1, dels, sig1, sig2, dfsig
integer :: num_l0_terms, k, i, j
integer, pointer :: index(:)
character (len=256) :: gyre_file
include 'formats'
ierr = 0
if (code == 'gyre') then
if (.not. gyre_is_enabled) then
ierr = -1
write(*,*)
write(*,'(a)') 'gyre is not currently enabled in your configuration of mesa.'
write(*,'(a)') 'check that your utils/makefile_header has USE_GYRE = YES'
write(*,*)
return
end if
if (l==0) then
gyre_file = gyre_l0_input_file
else if (l==1) then
gyre_file = gyre_l1_input_file
else if (l==2) then
gyre_file = gyre_l2_input_file
else if (l==3) then
gyre_file = gyre_l3_input_file
else
ierr = -1
write(*,*)
write(*,2) 'l must be in range 0 to 3', l
write(*,*)
return
end if
num_results = 0
call do_gyre_get_modes(s, l, gyre_file, store_model, ierr)
if (ierr /= 0) then
write(*,*) 'failed in do_gyre_get_modes'
stop 'get_one_el_info'
end if
else if (code == 'adipls') then
R = Rsun*s% photosphere_r
G = standard_cgrav
M = s% m_grav(1)
sig_fac = (2*pi)*(2*pi)*R*R*R/(G*M)
b = correction_b
! set controls for adipls
nsel = 0
dels = 1
els1 = dble(l)
itrsig = 1
sig1 = sig_fac*(nu1*1d-6)*(nu1*1d-6)
sig2 = sig_fac*(nu2*1d-6)*(nu2*1d-6)
dfsig = sig_fac*delta_nu_model*delta_nu_model
nsig = 2
call set_adipls_controls( &
l, nsel, els1, dels, itrsig, iscan, sig1, sig2, dfsig, nsig, &
adipls_irotkr, adipls_nprtkr, adipls_igm1kr, adipls_npgmkr)
num_results = 0
call run_adipls( &
s, .false., store_model, &
add_center_point, keep_surface_point, add_atmosphere, &
do_redistribute_mesh, ierr)
if (ierr /= 0) then
write(*,*) 'failed in run_adipls'
stop 'get_one_el_info'
end if
else
write(*,'(a)') 'invalid oscillation_code: ' // trim(oscillation_code)
ierr = -1
return
end if
! sort results by increasing frequency
allocate(index(num_results), stat=ierr)
if (ierr /= 0) then
stop 'failed in allocate before calling qsort'
end if
call qsort(index, num_results, cyclic_freq)
if (l == 0) then
call set_to_closest(l0_obs, l0_freq, l0_inertia, l0_order, ierr)
else if (l == 1) then
call set_to_closest(l1_obs, l1_freq, l1_inertia, l1_order, ierr)
else if (l == 2) then
call set_to_closest(l2_obs, l2_freq, l2_inertia, l2_order, ierr)
else if (l == 3) then
call set_to_closest(l3_obs, l3_freq, l3_inertia, l3_order, ierr)
else
stop 'bad value for l in get_one_el_info'
end if
if (ierr /= 0) then
return
end if
if (l == 0 .and. correction_factor > 0 .and. nl0 > 0 .and. &
delta_nu > 0 .and. nu_max > 0 .and. avg_nu_obs > 0) then
! calculate surface correction info
cnt = 0
sum_1 = 0
do i=1,nl0
if (l0_obs(i) < 0) cycle
cnt = cnt + 1
sum_1 = sum_1 + l0_freq(i)
end do
if (cnt == 0) return
avg_nu_model = sum_1/cnt
sum_1 = 0
sum_2 = 0
sum_3 = 0
do i=1,nl0
if (l0_obs(i) < 0) cycle
sum_1 = sum_1 + &
(l0_freq(i) - avg_nu_model)*(l0_n_obs(i) - avg_radial_n)
sum_2 = sum_2 + pow2(l0_n_obs(i) - avg_radial_n)
sum_3 = sum_3 + pow_cr(l0_obs(i)/nu_max,b)
end do
if (sum_2 == 0 .or. sum_3 == 0) return
delta_nu_model = sum_1/sum_2
correction_r = & ! K08 eqn 6
(b-1)/(b*avg_nu_model/avg_nu_obs - delta_nu_model/delta_nu)
if (correction_r <= 0) return
correction_a = & ! K08 eqn 10
min(0d0, avg_nu_obs - correction_r*avg_nu_model)*nl0/sum_3
a_div_r = correction_a/correction_r
end if
deallocate(index)
contains
subroutine set_to_closest(l_obs, l_freq, l_inertia, l_order, ierr)
real(dp), intent(in) :: l_obs(:)
real(dp), intent(out) :: l_freq(:), l_inertia(:)
integer, intent(out) :: l_order(:), ierr
integer :: i, j, jprev
jprev = 0
ierr = 0
do i = i1, i2
j = find_closest(l_obs(i),jprev)
if (j <= 0) then
l_freq(i) = 0
l_inertia(i) = 0
l_order(i) = 0
ierr = -1
else
l_freq(i) = cyclic_freq(j)
l_inertia(i) = inertia(j)
l_order(i) = order(j)
jprev = j
end if
end do
end subroutine set_to_closest
integer function find_closest(nu,jprev) ! find closest model frequency
real(dp), intent(in) :: nu
integer, intent(in) :: jprev
min_dist = 1d99; min_dist_j = -1
do j = jprev+1, num_results
if (el(j) /= l) cycle
dist = abs(cyclic_freq(j) - nu)
if (min_dist_j <= 0 .or. dist < min_dist) then
min_dist = dist; min_dist_j = j
end if
if (cyclic_freq(j) > nu) exit
end do
find_closest = min_dist_j
end function find_closest
end subroutine get_one_el_info
subroutine get_frequency_ratios( &
init, nl0, l0, nl1, l1, n, l0_first, l1_first, r01, r10)
logical, intent(in) :: init
integer, intent(in) :: nl0, nl1
real(dp), intent(in) :: l0(:), l1(:)
integer, intent(out) :: n, l0_first, l1_first
real(dp), intent(out) :: r01(:), r10(:)
integer :: l0_seq_n, l0_last, l1_seq_n, l1_last, i, i0, i1
real(dp) :: d01, d10, sd01, sd10, dnu, sdnu
logical :: dbg
include 'formats'
dbg = .false.
n = 0
if (nl1 <= 0) return
call get_max_sequence(nl0, l0, l0_first, l0_seq_n)
l0_last = l0_first + l0_seq_n - 1
if (dbg) write(*,4) 'l0_first l0_last l0_seq_n', l0_first, l0_last, l0_seq_n
call get_max_sequence(nl1, l1, l1_first, l1_seq_n)
l1_last = l1_first + l1_seq_n - 1
if (dbg) write(*,4) 'l1_first l1_last l1_seq_n', l1_first, l1_last, l1_seq_n
do ! trim high end of l0 until < last l1
if (l0(l0_last) < l1(l1_last)) exit
l0_last = l0_last - 1
if (l0_last < l0_first) then ! no overlap
if (dbg) then
write(*,*) 'l0_last < l0_first', l0_last, l0_first
end if
return
end if
end do
if (dbg) write(*,2) 'l0_last after trim', l0_last
do ! trim low end of l1 until > first l0
if (l1(l1_first) > l0(l0_first)) exit
l1_first = l1_first + 1
if (l1_first > l1_last) then ! no overlap
return
end if
end do
if (dbg) write(*,2) 'l1_first after trim', l1_first
do ! trim low end of l0 until only 1 < 1st l1
if (l0_first == l0_last) exit
if (l0(l0_first+1) >= l1(l1_first)) exit
l0_first = l0_first + 1
end do
if (dbg) write(*,2) 'l0_first after trim', l0_first
do ! trim high end of l1 until only 1 > last l0
if (l1_last == l1_first) exit
if (l1(l1_last-1) <= l0(l0_last)) exit
l1_last = l1_last - 1
end do
if (dbg) write(*,2) 'l1_last after trim', l1_last
l0_seq_n = l0_last - l0_first + 1
l1_seq_n = l1_last - l1_first + 1
n = l0_seq_n - 2
if (dbg) write(*,2) 'n', n
if (l0_seq_n /= l1_seq_n .or. n < 1) then
return
end if
do i = 1, n
i0 = i + l0_first
i1 = i + l1_first
d01 = (l0(i0-1) - 4*l1(i1-1) + 6*l0(i0) - 4*l1(i1) + l0(i0+1))/8d0
r01(i) = d01/(l1(i1) - l1(i1-1))
d10 = -(l1(i1-1) - 4*l0(i0) + 6*l1(i1) - 4*l0(i0+1) + l1(i1+1))/8d0
r10(i) = d10/(l0(i0+1) - l0(i0))
if (.not. init) cycle
! set ratio sigmas
sd01 = sqrt(pow2(l0_obs_sigma(i0-1)) + pow2(4*l1_obs_sigma(i1-1)) + &
pow2(6*l0_obs_sigma(i0)) + pow2(4*l1_obs_sigma(i1)) + pow2(l0_obs_sigma(i0+1)))/8d0
dnu = l1(i1) - l1(i1-1)
sdnu = sqrt(pow2(l1_obs_sigma(i1)) + pow2(l1_obs_sigma(i1-1)))
sigmas_r01(i) = sqrt(pow2(sd01/dnu) + pow2(sdnu*d01/(dnu*dnu)))
sd10 = sqrt(pow2(l1_obs_sigma(i1-1)) + pow2(4*l0_obs_sigma(i0)) + &
pow2(6*l1_obs_sigma(i1)) + pow2(4*l0_obs_sigma(i0+1)) + pow2(l1_obs_sigma(i1+1)))/8d0
dnu = l0(i0+1) - l0(i0)
sdnu = sqrt(pow2(l0_obs_sigma(i0+1)) + pow2(l0_obs_sigma(i0)))
sigmas_r10(i) = sqrt(pow2(sd10/dnu) + pow2(sdnu*d10/(dnu*dnu)) )
if (trace_chi2_seismo_ratios_info) then
write(*,'(a30,i4,99f16.6)') 'r01 r10 sigmas_r01 sigmas_r10', &
i, r01(i), r10(i), sigmas_r01(i), sigmas_r10(i)
end if
end do
end subroutine get_frequency_ratios
subroutine get_r02_frequency_ratios(init, nl0, l0, nl1, l1, nl2, l2, r02)
logical, intent(in) :: init
integer, intent(in) :: nl0, nl1, nl2
real(dp), intent(in) :: l0(:), l1(:), l2(:)
real(dp), intent(out) :: r02(:)
integer :: i, i0, i1, i2, jmin, j
real(dp) :: d02, sd02, dnu, sdnu, df, f0, f2, fmin, fmax, dfmin
logical :: dbg
include 'formats'
dbg = .false.
if (init) then ! set i2_for_r02
do i = 1, ratios_n
i0 = i + ratios_l0_first
i1 = i + ratios_l1_first
dnu = l1(i1) - l1(i1-1)
df = 0.25*dnu
f0 = l0(i0)
fmin = f0 - df
fmax = f0 + df
dfmin = 1d99
jmin = 0
do j = 1, nl2
f2 = l2(j)
!if (i==1) write(*,2) 'f2', j, fmin, f0, f2, fmax, abs(f2 - f0), dfmin
if (f2 <= fmax .and. f2 >= fmin .and. &
abs(f2 - f0) < dfmin) then
dfmin = abs(f2 - f0)
jmin = j
end if
if (f2 > f0) exit
end do
!if (.true.) write(*,3) 'i2_for_r02', i, jmin, fmin, f0, fmax, dnu
i2_for_r02(i) = jmin
sigmas_r02(i) = 0d0
end do
end if
!write(*,2) 'ratios_n', ratios_n
!stop
do i = 1, ratios_n
if ((.not. init) .and. sigmas_r02(i) == 0d0) cycle
i2 = i2_for_r02(i)
if (i2 == 0) cycle
i0 = i + ratios_l0_first
i1 = i + ratios_l1_first
d02 = l0(i0) - l2(i2)
dnu = l1(i1) - l1(i1-1)
r02(i) = d02/dnu
if (.not. init) cycle
! set ratio sigmas
sd02 = sqrt(pow2(l0_obs_sigma(i0)) + pow2(l2_obs_sigma(i2)))
sdnu = sqrt(pow2(l1_obs_sigma(i1)) + pow2(l1_obs_sigma(i1-1)))
sigmas_r02(i) = sqrt(pow2(sd02/dnu) + pow2(sdnu*d02/(dnu*dnu)))
if (trace_chi2_seismo_ratios_info) then
write(*,'(a30,i4,99f16.6)') 'r02 sigmas_r02', &
i, r02(i), sigmas_r02(i)
end if
end do
end subroutine get_r02_frequency_ratios
real(dp) function interpolate_ratio_r010( &
freq, first, model_freqs, model_ratios) result(ratio)
real(dp), intent(in) :: freq
integer, intent(in) :: first
real(dp), intent(in), dimension(:) :: model_freqs, model_ratios
integer :: i, j
real(dp) :: alfa, beta
ratio = 0
if (ratios_n == 0) return
j = 1 + first
if (freq <= model_freqs(j)) then
ratio = model_ratios(j)
return
end if
j = ratios_n + first
if (freq >= model_freqs(j)) then
ratio = model_ratios(j)
return
end if
do i=2,ratios_n
j = i+first
if (freq < model_freqs(j)) then
alfa = (freq - model_freqs(j-1))/(model_freqs(j) - model_freqs(j-1))
beta = 1d0 - alfa
ratio = alfa*model_ratios(j) + beta*model_ratios(j-1)
return
end if
end do
end function interpolate_ratio_r010
real(dp) function interpolate_ratio_r02( &
freq, model_freqs, model_ratios) result(ratio)
real(dp), intent(in) :: freq
real(dp), intent(in), dimension(:) :: model_freqs, model_ratios
integer :: i, i_lo, i_hi
real(dp) :: alfa, beta
ratio = 0
i_lo = 0
do i=1,nl0
if (sigmas_r02(i) == 0) cycle
i_lo = i
exit
end do
if (i_lo == 0) return
if (freq <= model_freqs(i_lo)) then
ratio = model_ratios(i_lo)
return
end if
i_hi = i_lo
do i=i_lo+1,nl0
if (sigmas_r02(i) == 0) cycle
i_hi = i
if (freq > model_freqs(i)) cycle
alfa = (freq - model_freqs(i_lo))/ &
(model_freqs(i_hi) - model_freqs(i_lo))
beta = 1d0 - alfa
ratio = alfa*model_ratios(i_hi) + beta*model_ratios(i_lo)
return
end do
ratio = model_ratios(i_hi)
end function interpolate_ratio_r02
subroutine get_max_sequence(nl, l_obs, max_seq_i, max_seq_n)
integer, intent(in) :: nl
real(dp), intent(in) :: l_obs(:)
integer, intent(out) :: max_seq_i, max_seq_n
integer :: i, j, seq_i, seq_n
max_seq_i = 0
max_seq_n = 0
seq_i = 0
seq_n = 0
do
i = seq_i + seq_n + 1 ! start of next sequence
if (i >= nl) exit
seq_i = i
seq_n = 1
do j = seq_i, nl-1 ! j is in series; try to add j+1
if (l_obs(j+1) - l_obs(j) > 1.5*delta_nu) then ! end of series
if (seq_n > max_seq_n) then
max_seq_i = seq_i
max_seq_n = seq_n
end if
exit
end if
seq_n = seq_n + 1
end do
end do
if (seq_n > max_seq_n) then
max_seq_i = seq_i
max_seq_n = seq_n
end if
end subroutine get_max_sequence
subroutine init_obs_data(ierr)
integer, intent(out) :: ierr
integer :: i, cnt, norders
integer, dimension(max_nl0) :: orders
real(dp) :: sum_1, sum_2, sum_3, range, nmax
real(dp) :: x, y, isig2, sum_xy, sum_x, sum_y, sum_x2, sum_isig2, d
logical, parameter :: dbg = .false.
include 'formats'
ierr = 0
!call test_get_frequency_ratios
if (nl0 <= 0) return
sigmas_r02 = 0d0
if (chi2_seismo_r_010_fraction > 0) then
call get_frequency_ratios( &
.true., nl0, l0_obs, nl1, l1_obs, &
ratios_n, ratios_l0_first, ratios_l1_first, &
ratios_r01, ratios_r10)
end if
if (chi2_seismo_r_02_fraction > 0) then
call get_r02_frequency_ratios( &
.true., nl0, l0_obs, nl1, l1_obs, nl2, l2_obs, ratios_r02)
end if
if (delta_nu <= 0 .and. nl0 > 1 .and. l0_n_obs(1) > 0) then
sum_xy = 0
sum_x = 0
sum_y = 0
sum_x2 = 0
sum_isig2 = 0
do i=1,nl0
isig2 = 1d0/pow2(l0_obs_sigma(i))
x = dble(l0_n_obs(i))
y = l0_obs(i)
sum_xy = sum_xy + x*y*isig2
sum_x = sum_x + x*isig2
sum_y = sum_y + y*isig2
sum_x2 = sum_x2 + x*x*isig2
sum_isig2 = sum_isig2 + isig2
end do
d = sum_isig2*sum_x2 - sum_x*sum_x
delta_nu = (sum_isig2*sum_xy - sum_x*sum_y)/d
if (delta_nu_sigma <= 0) delta_nu_sigma = sqrt(sum_isig2/d)
end if
if (correction_factor <= 0) return
if (l0_n_obs(1) <= 0) then
if (delta_nu <= 0) then
write(*,*) 'must supply value for delta_nu'
ierr = -1
return
end if
! set l0_n_obs(i) to order of l0_obs(i)
range = l0_obs(nl0) - l0_obs(1)
norders = int(range/delta_nu + 0.5d0) + 1
nmax = (nu_max/delta_nu)*(delta_nu_sun/nu_max_sun)*22.6 - 1.6
l0_n_obs(1) = int(nmax - (norders-1)/2)
if (dbg) write(*,3) 'l0_n_obs(i)', 1, l0_n_obs(1), l0_obs(1)
do i=2,norders
l0_n_obs(i) = l0_n_obs(1) + &
int((l0_obs(i) - l0_obs(1))/delta_nu + 0.5)
if (dbg) write(*,3) 'l0_n_obs(i)', i, l0_n_obs(i), l0_obs(i)
end do
if (dbg) then
write(*,1) 'range', range
write(*,2) 'norders', norders
write(*,1) 'nmax', nmax
write(*,2) '(norders+1)/2', (norders+1)/2
write(*,2) 'l0_n_obs(1)', l0_n_obs(1)
write(*,*)
!stop
end if
end if
cnt = 0
sum_1 = 0
sum_2 = 0
do i=1,nl0
if (l0_obs(i) < 0) cycle
cnt = cnt + 1
sum_1 = sum_1 + l0_obs(i)
sum_2 = sum_2 + l0_n_obs(i)
end do
avg_nu_obs = sum_1/cnt
avg_radial_n = sum_2/cnt
if (dbg) then
write(*,1) 'avg_nu_obs', avg_nu_obs
write(*,1) 'avg_radial_n', avg_radial_n
write(*,2) 'cnt', cnt
write(*,1) 'sum_1', sum_1
write(*,1) 'sum_2', sum_2
write(*,*)
stop 'init_obs_data'
end if
end subroutine init_obs_data
real(dp) function interpolate_l0_inertia(freq) result(inertia)
real(dp), intent(in) :: freq
integer :: i
real(dp) :: alfa, beta
inertia = 0
if (nl0 == 0) return
if (freq <= l0_freq(1)) then
inertia = l0_inertia(1)
return
end if
if (freq >= l0_freq(nl0)) then
inertia = l0_inertia(nl0)
return
end if
do i=2,nl0
if (freq < l0_freq(i)) then
alfa = (freq - l0_freq(i-1))/(l0_freq(i) - l0_freq(i-1))
beta = 1d0 - alfa
inertia = alfa*l0_inertia(i) + beta*l0_inertia(i-1)
return
end if
end do
end function interpolate_l0_inertia
subroutine get_l0_freq_corr( &
a_div_r, b, nu_max, correction_factor, check_obs, &
nl0, l0_obs, l0_freq, l0_freq_corr, l0_inertia)
real(dp), intent(in) :: a_div_r, b, nu_max, correction_factor
logical, intent(in) :: check_obs ! if false, then l0_obs is not used
integer, intent(in) :: nl0
real(dp), intent(in), dimension(:) :: &
l0_obs, l0_freq, l0_inertia
real(dp), intent(out) :: l0_freq_corr(:)
integer :: i
real(dp) :: Qnl
do i = 1, nl0
if (check_obs) then
if (l0_obs(i) < 0) cycle
end if
Qnl = 1
l0_freq_corr(i) = l0_freq(i)
if (b > 0 .and. correction_factor > 0 .and. l0_freq_corr(i) > 0) &
l0_freq_corr(i) = l0_freq_corr(i) + &
correction_factor*(a_div_r/Qnl)*pow_cr(l0_freq(i)/nu_max,b)
end do
end subroutine get_l0_freq_corr
subroutine get_l1_freq_corr( &
a_div_r, b, nu_max, correction_factor, check_obs, &
nl1, l1_obs, l1_freq, l1_freq_corr, l1_inertia, l0_inertia)
real(dp), intent(in) :: a_div_r, b, nu_max, correction_factor
logical, intent(in) :: check_obs ! if false, then l1_obs is not used
integer, intent(in) :: nl1
real(dp), intent(in), dimension(:) :: &
l1_obs, l1_freq, l1_inertia, l0_inertia
real(dp), intent(out) :: l1_freq_corr(:)
integer :: i
real(dp) :: Qnl, interp_l0_inertia
include 'formats'
do i = 1, nl1
if (check_obs) then
if (l1_obs(i) < 0) cycle
end if
l1_freq_corr(i) = l1_freq(i)
if (b > 0 .and. correction_factor > 0 .and. l1_freq_corr(i) > 0 .and. nl0 > 0) then
interp_l0_inertia = interpolate_l0_inertia(l1_freq(i))
!write(*,2) 'l1_freq_corr: l0_inertia interp prev', i, interp_l0_inertia, &
! (l0_inertia(min(nl0,i)) + l0_inertia(min(nl0,i+1)))/2
Qnl = l1_inertia(i)/interp_l0_inertia
l1_freq_corr(i) = l1_freq_corr(i) + &
correction_factor*(a_div_r/Qnl)*pow_cr(l1_freq(i)/nu_max,b)
end if
end do
end subroutine get_l1_freq_corr
subroutine get_l2_freq_corr( &
a_div_r, b, nu_max, correction_factor, check_obs, &
nl2, l2_obs, l2_freq, l2_freq_corr, l2_inertia, l0_inertia)
real(dp), intent(in) :: a_div_r, b, nu_max, correction_factor
logical, intent(in) :: check_obs ! if false, then l2_obs is not used
integer, intent(in) :: nl2
real(dp), intent(in), dimension(:) :: &
l2_obs, l2_freq, l2_inertia, l0_inertia
real(dp), intent(out) :: l2_freq_corr(:)
integer :: i
real(dp) :: Qnl, interp_l0_inertia
include 'formats'
do i = 1, nl2
if (check_obs) then
if (l2_obs(i) < 0) cycle
end if
l2_freq_corr(i) = l2_freq(i)
if (b > 0 .and. correction_factor > 0 .and. l2_freq_corr(i) > 0 .and. nl0 > 0) then
interp_l0_inertia = interpolate_l0_inertia(l2_freq(i))
Qnl = l2_inertia(i)/interp_l0_inertia
!write(*,2) 'l2_freq_corr: l0_inertia interp prev', i, interp_l0_inertia, &
! l0_inertia(min(nl0,i+1))
l2_freq_corr(i) = l2_freq_corr(i) + &
correction_factor*(a_div_r/Qnl)*pow_cr(l2_freq(i)/nu_max,b)
end if
end do
end subroutine get_l2_freq_corr
subroutine get_l3_freq_corr( &
a_div_r, b, nu_max, correction_factor, check_obs, &
nl3, l3_obs, l3_freq, l3_freq_corr, l3_inertia, l0_inertia)
real(dp), intent(in) :: a_div_r, b, nu_max, correction_factor
logical, intent(in) :: check_obs ! if false, then l3_obs is not used
integer, intent(in) :: nl3
real(dp), intent(in), dimension(:) :: &
l3_obs, l3_freq, l3_inertia, l0_inertia
real(dp), intent(out) :: l3_freq_corr(:)
integer :: i
real(dp) :: Qnl, interp_l0_inertia
include 'formats'
do i = 1, nl3
if (check_obs) then
if (l3_obs(i) < 0) cycle
end if
l3_freq_corr(i) = l3_freq(i)
if (b > 0 .and. correction_factor > 0 .and. l3_freq_corr(i) > 0 .and. nl0 > 0) then
interp_l0_inertia = interpolate_l0_inertia(l3_freq(i))
Qnl = l3_inertia(i)/interp_l0_inertia
l3_freq_corr(i) = l3_freq_corr(i) + &
correction_factor*(a_div_r/Qnl)*pow_cr(l3_freq(i)/nu_max,b)
end if
end do
end subroutine get_l3_freq_corr
subroutine get_freq_corr
real(dp), allocatable :: l_freq(:), l_obs(:), &
l_obs_sigma(:), y(:), X(:)
real(dp) :: XtX, Xty, z
real(dp) :: detXtX
integer :: i, j, k, N
logical :: only_radial
! fetch all frequencies into one array
if (l1_freq(1) > 0d0) then
only_radial = .false.
else
only_radial = .true.
end if
if (only_radial) then
N = nl0
else
N = nl0 + nl1 + nl2 + nl3
end if
allocate(l_freq(N), l_obs(N), l_obs_sigma(N), y(N))
allocate(X(N))
N = 0
do i = 1, nl0
l_freq(N+i) = l0_freq(i)
l_obs(N+i) = l0_obs(i)
l_obs_sigma(N+i) = l0_obs_sigma(i)
X(N+i) = l0_freq(i)**3/l0_inertia(i)/l0_obs_sigma(i)
end do
if (.not. only_radial) then
N = nl0
do i = 1, nl1
l_freq(N+i) = l1_freq(i)
l_obs(N+i) = l1_obs(i)
l_obs_sigma(N+i) = l1_obs_sigma(i)
X(N+i) = l1_freq(i)**3/l1_inertia(i)/l1_obs_sigma(i)
end do
N = nl0 + nl1
do i = 1, nl2
l_freq(N+i) = l2_freq(i)
l_obs(N+i) = l2_obs(i)
l_obs_sigma(N+i) = l2_obs_sigma(i)
X(N+i) = l2_freq(i)**3/l2_inertia(i)/l2_obs_sigma(i)
end do
N = nl0 + nl1 + nl2
do i = 1, nl3
l_freq(N+i) = l3_freq(i)
l_obs(N+i) = l3_obs(i)
l_obs_sigma(N+i) = l3_obs_sigma(i)
X(N+i) = l3_freq(i)**3/l3_inertia(i)/l3_obs_sigma(i)
end do
end if
if (only_radial) then
N = nl0
else
N = nl0 + nl1 + nl2 + nl3
end if
y = (l_obs - l_freq)/l_obs_sigma
! write(*,*) y(nl0), y(nl1), y(nl2), y(nl3)
! write(*,*) X(nl0,1), X(nl1,1), X(nl2,1), X(nl3,1)
! create matrices
XtX = sum(X*X)
Xty = sum(X*y)
! write(*,*) 'XtX ='
! write(*,*) XtX(1,1), XtX(1,2)
! write(*,*) XtX(2,1), XtX(2,2)
! evaluate matrix equation
z = Xty/XtX
write(*,*) 'a_3 =', z
l0_freq_corr = l0_freq + z*l0_freq**3/l0_inertia
l1_freq_corr = l1_freq + z*l1_freq**3/l1_inertia
l2_freq_corr = l2_freq + z*l2_freq**3/l2_inertia
l3_freq_corr = l3_freq + z*l3_freq**3/l3_inertia
deallocate(l_freq, l_obs, l_obs_sigma, y)
deallocate(X)
! call get_l0_freq_corr( &
! a_div_r, correction_b, nu_max, correction_factor, .true., &
! nl0, l0_obs, l0_freq, l0_freq_corr, l0_inertia)
! call get_l1_freq_corr( &
! a_div_r, correction_b, nu_max, correction_factor, .true., &
! nl1, l1_obs, l1_freq, l1_freq_corr, l1_inertia, l0_inertia)
! call get_l2_freq_corr( &
! a_div_r, correction_b, nu_max, correction_factor, .true., &
! nl2, l2_obs, l2_freq, l2_freq_corr, l2_inertia, l0_inertia)
! call get_l3_freq_corr( &
! a_div_r, correction_b, nu_max, correction_factor, .true., &
! nl3, l3_obs, l3_freq, l3_freq_corr, l3_inertia, l0_inertia)
end subroutine get_freq_corr
! chi2 = chi2_seismo*chi2_seismo_fraction &
! + chi2_spectro*(1 - chi2_seismo_fraction)
real(dp) function get_chi2(s, max_el, trace_okay, ierr)
type (star_info), pointer :: s
integer, intent(in) :: max_el
logical, intent(in) :: trace_okay
integer, intent(out) :: ierr
integer :: i, n, chi2N1, chi2N2
real(dp) :: chi2term, Teff, logL, chi2sum1, chi2sum2, frac, &
model_r01, model_r10, model_r02
! calculate chi^2 following Brandao et al, 2011, eqn 11
include 'formats'
ierr = 0
chi2sum1 = 0
chi2N1 = 0
chi2_r_010_ratios = -1
chi2_r_02_ratios = -1
chi2_frequencies = 0
if (chi2_seismo_freq_fraction > 0) then
if (trace_okay .and. trace_chi2_seismo_frequencies_info) &
write(*,'(a30,a6,99(a20))') &
'chi2term l0', 'model number', 'i', 'chi2term', &
'l0_freq_corr(i)', 'l0_obs(i)', 'l0_obs_sigma(i)'
do i = 1, nl0
if (l0_obs(i) < 0) cycle
chi2term = &
pow2((l0_freq_corr(i) - l0_obs(i))/l0_obs_sigma(i))
if (trace_okay .and. trace_chi2_seismo_frequencies_info) &
write(*,'(a30,2i6,99(1pe20.10))') 'chi2term l0', &
s% model_number, i, chi2term, &
l0_freq_corr(i), l0_obs(i), l0_obs_sigma(i)
chi2sum1 = chi2sum1 + chi2term
chi2N1 = chi2N1 + 1
end do
if (max_el >= 1) then
do i = 1, nl1
if (l1_obs(i) < 0) cycle
chi2term = &
pow2((l1_freq_corr(i) - l1_obs(i))/l1_obs_sigma(i))
if (trace_okay .and. trace_chi2_seismo_frequencies_info) &
write(*,'(a30,2i6,99(1pe20.10))') 'chi2term l1', &
s% model_number, i, chi2term, &
l1_freq_corr(i), l1_obs(i), l1_obs_sigma(i)
chi2sum1 = chi2sum1 + chi2term
chi2N1 = chi2N1 + 1
end do
end if
if (max_el >= 2) then
do i = 1, nl2
if (l2_obs(i) < 0) cycle
chi2term = &
pow2((l2_freq_corr(i) - l2_obs(i))/l2_obs_sigma(i))
if (trace_okay .and. trace_chi2_seismo_frequencies_info) &
write(*,'(a30,2i6,99(1pe20.10))') 'chi2term l2', &
s% model_number, i, chi2term, &
l2_freq_corr(i), l2_obs(i), l2_obs_sigma(i)
chi2sum1 = chi2sum1 + chi2term
chi2N1 = chi2N1 + 1
end do
end if
if (max_el >= 3) then
do i = 1, nl3
if (l3_obs(i) < 0) cycle
chi2term = &
pow2((l3_freq_corr(i) - l3_obs(i))/l3_obs_sigma(i))
if (trace_okay .and. trace_chi2_seismo_frequencies_info) &
write(*,'(a30,2i6,99(1pe20.10))') 'chi2term l3', &
s% model_number, i, chi2term, &
l3_freq_corr(i), l3_obs(i), l3_obs_sigma(i)
chi2sum1 = chi2sum1 + chi2term
chi2N1 = chi2N1 + 1
end do
end if
num_chi2_seismo_terms = chi2N1
chi2_frequencies = chi2sum1! /max(1,chi2N1)
end if
if (chi2_seismo_r_010_fraction > 0 .and. max_el >= 1) then
if (ratios_n == 0) then
write(*,*) 'ERROR: chi2_seismo_r_010_fraction > 0 but cannot evaluate r_010'
ierr = -1
return
end if
chi2sum1 = 0
do i=1,ratios_n
model_r01 = interpolate_ratio_r010( &
l0_obs(i + ratios_l0_first), ratios_l0_first, l0_freq, model_ratios_r01)
if (trace_okay .and. trace_chi2_seismo_ratios_info) &
write(*,2) 'r01 obs, model, interp, model - interp', &
i, ratios_r01(i), model_ratios_r01(i + ratios_l0_first), model_r01, &
model_ratios_r01(i + ratios_l0_first) - model_r01
model_r10 = interpolate_ratio_r010( &
l1_obs(i + ratios_l1_first), ratios_l1_first, l1_freq, model_ratios_r10)
if (trace_okay .and. trace_chi2_seismo_ratios_info) &
write(*,2) 'r10 obs, model, interp, model - interp', &
i, ratios_r10(i), model_ratios_r10(i + ratios_l1_first), model_r10, &
model_ratios_r10(i + ratios_l1_first) - model_r10
chi2term = &
pow2((model_r01 - ratios_r01(i))/sigmas_r01(i)) + &
pow2((model_r10 - ratios_r10(i))/sigmas_r10(i))
chi2sum1 = chi2sum1 + chi2term
if (trace_okay .and. trace_chi2_seismo_ratios_info) &
write(*,2) 'chi2 ratios terms r01 r10', i, chi2term, &
pow2((model_r01 - ratios_r01(i))/sigmas_r01(i)), &
pow2((model_r10 - ratios_r10(i))/sigmas_r10(i)), &
model_r01, model_r10
end do
n = 2*ratios_n
chi2_r_010_ratios = chi2sum1! /max(1,n)
end if
if (chi2_seismo_r_02_fraction > 0 .and. max_el >= 2) then
chi2sum1 = 0
n = 0
do i=2,nl0
if (sigmas_r02(i) == 0d0) cycle
model_r02 = interpolate_ratio_r02( &
l0_obs(i + ratios_l0_first), l0_freq, model_ratios_r02)
if (trace_okay .and. trace_chi2_seismo_ratios_info) &
write(*,2) 'r02 obs, model, interp, model - interp', &
i, ratios_r02(i), model_ratios_r02(i), model_r02, &
model_ratios_r02(i) - model_r10
chi2sum1 = chi2sum1 + &
pow2((model_r02 - ratios_r02(i))/sigmas_r02(i))
n = n+1
end do
if (n == 0) then
write(*,*) 'ERROR: chi2_seismo_r_02_fraction > 0 but cannot evaluate r_02'
ierr = -1
return
end if
chi2_r_02_ratios = chi2sum1! /max(1,n)
end if
chi2_seismo = &
chi2_seismo_r_010_fraction*chi2_r_010_ratios + &
chi2_seismo_r_02_fraction*chi2_r_02_ratios + &
chi2_seismo_freq_fraction*chi2_frequencies + &
chi2_seismo_delta_nu_fraction*chi2_delta_nu + &
chi2_seismo_nu_max_fraction*chi2_nu_max
chi2sum2 = 0
chi2N2 = 0
if (Teff_sigma > 0 .and. include_Teff_in_chi2_spectro) then
Teff = s% Teff
chi2term = pow2((Teff - Teff_target)/Teff_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term Teff', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (logL_sigma > 0 .and. include_logL_in_chi2_spectro) then
logL = s% log_surface_luminosity
chi2term = pow2((logL - logL_target)/logL_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term logL', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (logg_sigma > 0 .and. include_logg_in_chi2_spectro) then
chi2term = pow2((logg - logg_target)/logg_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term logg', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (FeH_sigma > 0 .and. include_FeH_in_chi2_spectro) then
chi2term = pow2((FeH - FeH_target)/FeH_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term FeH', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (logR_sigma > 0 .and. include_logR_in_chi2_spectro) then
chi2term = pow2((logR - logR_target)/logR_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term logR', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (age_sigma > 0 .and. include_age_in_chi2_spectro) then
chi2term = pow2((s% star_age - age_target)/age_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term age', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (surface_Z_div_X_sigma > 0 .and. include_surface_Z_div_X_in_chi2_spectro) then
chi2term = pow2((surface_Z_div_X - surface_Z_div_X_target)/surface_Z_div_X_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term surface_Z_div_X', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (surface_He_sigma > 0 .and. include_surface_He_in_chi2_spectro) then
chi2term = pow2((surface_He - surface_He_target)/surface_He_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term surface_He', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (Rcz_sigma > 0 .and. include_Rcz_in_chi2_spectro) then
chi2term = pow2((Rcz - Rcz_target)/Rcz_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term Rcz', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (csound_rms_sigma > 0 .and. include_csound_rms_in_chi2_spectro) then
chi2term = pow2((csound_rms - csound_rms_target)/csound_rms_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term csound_rms', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (my_var1_sigma > 0 .and. include_my_var1_in_chi2_spectro) then
chi2term = pow2((my_var1 - my_var1_target)/my_var1_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term ' // trim(my_var1_name), s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (my_var2_sigma > 0 .and. include_my_var2_in_chi2_spectro) then
chi2term = pow2((my_var2 - my_var2_target)/my_var2_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term ' // trim(my_var2_name), s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (my_var3_sigma > 0 .and. include_my_var3_in_chi2_spectro) then
chi2term = pow2((my_var3 - my_var3_target)/my_var3_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term ' // trim(my_var3_name), s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
num_chi2_spectro_terms = chi2N2
chi2_spectro = chi2sum2! /max(1,chi2N2)
frac = chi2_seismo_fraction
chi2 = frac*chi2_seismo + (1-frac)*chi2_spectro
get_chi2 = chi2
if (chi2_seismo_fraction < 0 .or. chi2_seismo_fraction > 1) then
write(*,1) 'chi2_seismo_fraction', chi2_seismo_fraction
stop
end if
!if (is_bad_num(chi2)) stop 'get_chi2'
end function get_chi2
end module astero_support
-------------- next part --------------
! ***********************************************************************
!
! Copyright (C) 2013 Bill Paxton
!
! this file is part of mesa.
!
! mesa is free software; you can redistribute it and/or modify
! it under the terms of the gnu general library public license as published
! by the free software foundation; either version 2 of the license, or
! (at your option) any later version.
!
! mesa is distributed in the hope that it will be useful,
! but without any warranty; without even the implied warranty of
! merchantability or fitness for a particular purpose. see the
! gnu library general public license for more details.
!
! you should have received a copy of the gnu library general public license
! along with this software; if not, write to the free software
! foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa
!
! ***********************************************************************
module astero_support
use astero_data
use star_lib
use star_def
use const_def
use crlibm_lib
use utils_lib
implicit none
contains
subroutine check_search_controls(ierr)
integer, intent(out) :: ierr
integer :: i
include 'formats'
ierr = 0
do i=2,nl0
if (l0_obs(i) <= l0_obs(i-1)) then
write(*,3) 'l0_obs values out of order', i-1, i, l0_obs(i-1), l0_obs(i)
ierr = -1
end if
end do
do i=2,nl1
if (l1_obs(i) <= l1_obs(i-1)) then
write(*,3) 'l1_obs values out of order', i-1, i, l1_obs(i-1), l1_obs(i)
ierr = -1
end if
end do
do i=2,nl2
if (l2_obs(i) <= l2_obs(i-1)) then
write(*,3) 'l2_obs values out of order', i-1, i, l2_obs(i-1), l2_obs(i)
ierr = -1
end if
end do
do i=2,nl3
if (l3_obs(i) <= l3_obs(i-1)) then
write(*,3) 'l3_obs values out of order', i-1, i, l3_obs(i-1), l3_obs(i)
ierr = -1
end if
end do
if (ierr /= 0) &
write(*,1) 'please put frequency values in ascending order'
end subroutine check_search_controls
subroutine get_one_el_info( &
s, l, nu1, nu2, iscan, i1, i2, store_model, code, ierr)
use num_lib, only: qsort
use adipls_support
use gyre_support
type (star_info), pointer :: s
integer, intent(in) :: l, iscan, i1, i2
real(dp), intent(in) :: nu1, nu2
logical, intent(in) :: store_model
character (len=*), intent(in) :: code
integer, intent(out) :: ierr
real(dp) :: nu_obs, dist_j, nu, dist, min_dist, min_freq, &
R, G, M, sig_fac, b, sum_1, sum_2, sum_3
integer :: min_dist_j, min_order, n, cnt
integer :: nsel, itrsig, nsig
real(dp) :: els1, dels, sig1, sig2, dfsig
integer :: num_l0_terms, k, i, j
integer, pointer :: index(:)
character (len=256) :: gyre_file
include 'formats'
ierr = 0
if (code == 'gyre') then
if (.not. gyre_is_enabled) then
ierr = -1
write(*,*)
write(*,'(a)') 'gyre is not currently enabled in your configuration of mesa.'
write(*,'(a)') 'check that your utils/makefile_header has USE_GYRE = YES'
write(*,*)
return
end if
if (l==0) then
gyre_file = gyre_l0_input_file
else if (l==1) then
gyre_file = gyre_l1_input_file
else if (l==2) then
gyre_file = gyre_l2_input_file
else if (l==3) then
gyre_file = gyre_l3_input_file
else
ierr = -1
write(*,*)
write(*,2) 'l must be in range 0 to 3', l
write(*,*)
return
end if
num_results = 0
call do_gyre_get_modes(s, l, gyre_file, store_model, ierr)
if (ierr /= 0) then
write(*,*) 'failed in do_gyre_get_modes'
stop 'get_one_el_info'
end if
else if (code == 'adipls') then
R = Rsun*s% photosphere_r
G = standard_cgrav
M = s% m_grav(1)
sig_fac = (2*pi)*(2*pi)*R*R*R/(G*M)
b = correction_b
! set controls for adipls
nsel = 0
dels = 1
els1 = dble(l)
itrsig = 1
sig1 = sig_fac*(nu1*1d-6)*(nu1*1d-6)
sig2 = sig_fac*(nu2*1d-6)*(nu2*1d-6)
dfsig = sig_fac*delta_nu_model*delta_nu_model
nsig = 2
call set_adipls_controls( &
l, nsel, els1, dels, itrsig, iscan, sig1, sig2, dfsig, nsig, &
adipls_irotkr, adipls_nprtkr, adipls_igm1kr, adipls_npgmkr)
num_results = 0
call run_adipls( &
s, .false., store_model, &
add_center_point, keep_surface_point, add_atmosphere, &
do_redistribute_mesh, ierr)
if (ierr /= 0) then
write(*,*) 'failed in run_adipls'
stop 'get_one_el_info'
end if
else
write(*,'(a)') 'invalid oscillation_code: ' // trim(oscillation_code)
ierr = -1
return
end if
! sort results by increasing frequency
allocate(index(num_results), stat=ierr)
if (ierr /= 0) then
stop 'failed in allocate before calling qsort'
end if
call qsort(index, num_results, cyclic_freq)
if (l == 0) then
call set_to_closest(l0_obs, l0_freq, l0_inertia, l0_order, ierr)
else if (l == 1) then
call set_to_closest(l1_obs, l1_freq, l1_inertia, l1_order, ierr)
else if (l == 2) then
call set_to_closest(l2_obs, l2_freq, l2_inertia, l2_order, ierr)
else if (l == 3) then
call set_to_closest(l3_obs, l3_freq, l3_inertia, l3_order, ierr)
else
stop 'bad value for l in get_one_el_info'
end if
if (ierr /= 0) then
return
end if
if (l == 0 .and. correction_factor > 0 .and. nl0 > 0 .and. &
delta_nu > 0 .and. nu_max > 0 .and. avg_nu_obs > 0) then
! calculate surface correction info
cnt = 0
sum_1 = 0
do i=1,nl0
if (l0_obs(i) < 0) cycle
cnt = cnt + 1
sum_1 = sum_1 + l0_freq(i)
end do
if (cnt == 0) return
avg_nu_model = sum_1/cnt
sum_1 = 0
sum_2 = 0
sum_3 = 0
do i=1,nl0
if (l0_obs(i) < 0) cycle
sum_1 = sum_1 + &
(l0_freq(i) - avg_nu_model)*(l0_n_obs(i) - avg_radial_n)
sum_2 = sum_2 + pow2(l0_n_obs(i) - avg_radial_n)
sum_3 = sum_3 + pow_cr(l0_obs(i)/nu_max,b)
end do
if (sum_2 == 0 .or. sum_3 == 0) return
delta_nu_model = sum_1/sum_2
correction_r = & ! K08 eqn 6
(b-1)/(b*avg_nu_model/avg_nu_obs - delta_nu_model/delta_nu)
if (correction_r <= 0) return
correction_a = & ! K08 eqn 10
min(0d0, avg_nu_obs - correction_r*avg_nu_model)*nl0/sum_3
a_div_r = correction_a/correction_r
end if
deallocate(index)
contains
subroutine set_to_closest(l_obs, l_freq, l_inertia, l_order, ierr)
real(dp), intent(in) :: l_obs(:)
real(dp), intent(out) :: l_freq(:), l_inertia(:)
integer, intent(out) :: l_order(:), ierr
integer :: i, j, jprev
jprev = 0
ierr = 0
do i = i1, i2
j = find_closest(l_obs(i),jprev)
if (j <= 0) then
l_freq(i) = 0
l_inertia(i) = 0
l_order(i) = 0
ierr = -1
else
l_freq(i) = cyclic_freq(j)
l_inertia(i) = inertia(j)
l_order(i) = order(j)
jprev = j
end if
end do
end subroutine set_to_closest
integer function find_closest(nu,jprev) ! find closest model frequency
real(dp), intent(in) :: nu
integer, intent(in) :: jprev
min_dist = 1d99; min_dist_j = -1
do j = jprev+1, num_results
if (el(j) /= l) cycle
dist = abs(cyclic_freq(j) - nu)
if (min_dist_j <= 0 .or. dist < min_dist) then
min_dist = dist; min_dist_j = j
end if
if (cyclic_freq(j) > nu) exit
end do
find_closest = min_dist_j
end function find_closest
end subroutine get_one_el_info
subroutine get_frequency_ratios( &
init, nl0, l0, nl1, l1, n, l0_first, l1_first, r01, r10)
logical, intent(in) :: init
integer, intent(in) :: nl0, nl1
real(dp), intent(in) :: l0(:), l1(:)
integer, intent(out) :: n, l0_first, l1_first
real(dp), intent(out) :: r01(:), r10(:)
integer :: l0_seq_n, l0_last, l1_seq_n, l1_last, i, i0, i1
real(dp) :: d01, d10, sd01, sd10, dnu, sdnu
logical :: dbg
include 'formats'
dbg = .false.
n = 0
if (nl1 <= 0) return
call get_max_sequence(nl0, l0, l0_first, l0_seq_n)
l0_last = l0_first + l0_seq_n - 1
if (dbg) write(*,4) 'l0_first l0_last l0_seq_n', l0_first, l0_last, l0_seq_n
call get_max_sequence(nl1, l1, l1_first, l1_seq_n)
l1_last = l1_first + l1_seq_n - 1
if (dbg) write(*,4) 'l1_first l1_last l1_seq_n', l1_first, l1_last, l1_seq_n
do ! trim high end of l0 until < last l1
if (l0(l0_last) < l1(l1_last)) exit
l0_last = l0_last - 1
if (l0_last < l0_first) then ! no overlap
if (dbg) then
write(*,*) 'l0_last < l0_first', l0_last, l0_first
end if
return
end if
end do
if (dbg) write(*,2) 'l0_last after trim', l0_last
do ! trim low end of l1 until > first l0
if (l1(l1_first) > l0(l0_first)) exit
l1_first = l1_first + 1
if (l1_first > l1_last) then ! no overlap
return
end if
end do
if (dbg) write(*,2) 'l1_first after trim', l1_first
do ! trim low end of l0 until only 1 < 1st l1
if (l0_first == l0_last) exit
if (l0(l0_first+1) >= l1(l1_first)) exit
l0_first = l0_first + 1
end do
if (dbg) write(*,2) 'l0_first after trim', l0_first
do ! trim high end of l1 until only 1 > last l0
if (l1_last == l1_first) exit
if (l1(l1_last-1) <= l0(l0_last)) exit
l1_last = l1_last - 1
end do
if (dbg) write(*,2) 'l1_last after trim', l1_last
l0_seq_n = l0_last - l0_first + 1
l1_seq_n = l1_last - l1_first + 1
n = l0_seq_n - 2
if (dbg) write(*,2) 'n', n
if (l0_seq_n /= l1_seq_n .or. n < 1) then
return
end if
do i = 1, n
i0 = i + l0_first
i1 = i + l1_first
d01 = (l0(i0-1) - 4*l1(i1-1) + 6*l0(i0) - 4*l1(i1) + l0(i0+1))/8d0
r01(i) = d01/(l1(i1) - l1(i1-1))
d10 = -(l1(i1-1) - 4*l0(i0) + 6*l1(i1) - 4*l0(i0+1) + l1(i1+1))/8d0
r10(i) = d10/(l0(i0+1) - l0(i0))
if (.not. init) cycle
! set ratio sigmas
sd01 = sqrt(pow2(l0_obs_sigma(i0-1)) + pow2(4*l1_obs_sigma(i1-1)) + &
pow2(6*l0_obs_sigma(i0)) + pow2(4*l1_obs_sigma(i1)) + pow2(l0_obs_sigma(i0+1)))/8d0
dnu = l1(i1) - l1(i1-1)
sdnu = sqrt(pow2(l1_obs_sigma(i1)) + pow2(l1_obs_sigma(i1-1)))
sigmas_r01(i) = sqrt(pow2(sd01/dnu) + pow2(sdnu*d01/(dnu*dnu)))
sd10 = sqrt(pow2(l1_obs_sigma(i1-1)) + pow2(4*l0_obs_sigma(i0)) + &
pow2(6*l1_obs_sigma(i1)) + pow2(4*l0_obs_sigma(i0+1)) + pow2(l1_obs_sigma(i1+1)))/8d0
dnu = l0(i0+1) - l0(i0)
sdnu = sqrt(pow2(l0_obs_sigma(i0+1)) + pow2(l0_obs_sigma(i0)))
sigmas_r10(i) = sqrt(pow2(sd10/dnu) + pow2(sdnu*d10/(dnu*dnu)) )
if (trace_chi2_seismo_ratios_info) then
write(*,'(a30,i4,99f16.6)') 'r01 r10 sigmas_r01 sigmas_r10', &
i, r01(i), r10(i), sigmas_r01(i), sigmas_r10(i)
end if
end do
end subroutine get_frequency_ratios
subroutine get_r02_frequency_ratios(init, nl0, l0, nl1, l1, nl2, l2, r02)
logical, intent(in) :: init
integer, intent(in) :: nl0, nl1, nl2
real(dp), intent(in) :: l0(:), l1(:), l2(:)
real(dp), intent(out) :: r02(:)
integer :: i, i0, i1, i2, jmin, j
real(dp) :: d02, sd02, dnu, sdnu, df, f0, f2, fmin, fmax, dfmin
logical :: dbg
include 'formats'
dbg = .false.
if (init) then ! set i2_for_r02
do i = 1, ratios_n
i0 = i + ratios_l0_first
i1 = i + ratios_l1_first
dnu = l1(i1) - l1(i1-1)
df = 0.25*dnu
f0 = l0(i0)
fmin = f0 - df
fmax = f0 + df
dfmin = 1d99
jmin = 0
do j = 1, nl2
f2 = l2(j)
!if (i==1) write(*,2) 'f2', j, fmin, f0, f2, fmax, abs(f2 - f0), dfmin
if (f2 <= fmax .and. f2 >= fmin .and. &
abs(f2 - f0) < dfmin) then
dfmin = abs(f2 - f0)
jmin = j
end if
if (f2 > f0) exit
end do
!if (.true.) write(*,3) 'i2_for_r02', i, jmin, fmin, f0, fmax, dnu
i2_for_r02(i) = jmin
sigmas_r02(i) = 0d0
end do
end if
!write(*,2) 'ratios_n', ratios_n
!stop
do i = 1, ratios_n
if ((.not. init) .and. sigmas_r02(i) == 0d0) cycle
i2 = i2_for_r02(i)
if (i2 == 0) cycle
i0 = i + ratios_l0_first
i1 = i + ratios_l1_first
d02 = l0(i0) - l2(i2)
dnu = l1(i1) - l1(i1-1)
r02(i) = d02/dnu
if (.not. init) cycle
! set ratio sigmas
sd02 = sqrt(pow2(l0_obs_sigma(i0)) + pow2(l2_obs_sigma(i2)))
sdnu = sqrt(pow2(l1_obs_sigma(i1)) + pow2(l1_obs_sigma(i1-1)))
sigmas_r02(i) = sqrt(pow2(sd02/dnu) + pow2(sdnu*d02/(dnu*dnu)))
if (trace_chi2_seismo_ratios_info) then
write(*,'(a30,i4,99f16.6)') 'r02 sigmas_r02', &
i, r02(i), sigmas_r02(i)
end if
end do
end subroutine get_r02_frequency_ratios
real(dp) function interpolate_ratio_r010( &
freq, first, model_freqs, model_ratios) result(ratio)
real(dp), intent(in) :: freq
integer, intent(in) :: first
real(dp), intent(in), dimension(:) :: model_freqs, model_ratios
integer :: i, j
real(dp) :: alfa, beta
ratio = 0
if (ratios_n == 0) return
j = 1 + first
if (freq <= model_freqs(j)) then
ratio = model_ratios(j)
return
end if
j = ratios_n + first
if (freq >= model_freqs(j)) then
ratio = model_ratios(j)
return
end if
do i=2,ratios_n
j = i+first
if (freq < model_freqs(j)) then
alfa = (freq - model_freqs(j-1))/(model_freqs(j) - model_freqs(j-1))
beta = 1d0 - alfa
ratio = alfa*model_ratios(j) + beta*model_ratios(j-1)
return
end if
end do
end function interpolate_ratio_r010
real(dp) function interpolate_ratio_r02( &
freq, model_freqs, model_ratios) result(ratio)
real(dp), intent(in) :: freq
real(dp), intent(in), dimension(:) :: model_freqs, model_ratios
integer :: i, i_lo, i_hi
real(dp) :: alfa, beta
ratio = 0
i_lo = 0
do i=1,nl0
if (sigmas_r02(i) == 0) cycle
i_lo = i
exit
end do
if (i_lo == 0) return
if (freq <= model_freqs(i_lo)) then
ratio = model_ratios(i_lo)
return
end if
i_hi = i_lo
do i=i_lo+1,nl0
if (sigmas_r02(i) == 0) cycle
i_hi = i
if (freq > model_freqs(i)) cycle
alfa = (freq - model_freqs(i_lo))/ &
(model_freqs(i_hi) - model_freqs(i_lo))
beta = 1d0 - alfa
ratio = alfa*model_ratios(i_hi) + beta*model_ratios(i_lo)
return
end do
ratio = model_ratios(i_hi)
end function interpolate_ratio_r02
subroutine get_max_sequence(nl, l_obs, max_seq_i, max_seq_n)
integer, intent(in) :: nl
real(dp), intent(in) :: l_obs(:)
integer, intent(out) :: max_seq_i, max_seq_n
integer :: i, j, seq_i, seq_n
max_seq_i = 0
max_seq_n = 0
seq_i = 0
seq_n = 0
do
i = seq_i + seq_n + 1 ! start of next sequence
if (i >= nl) exit
seq_i = i
seq_n = 1
do j = seq_i, nl-1 ! j is in series; try to add j+1
if (l_obs(j+1) - l_obs(j) > 1.5*delta_nu) then ! end of series
if (seq_n > max_seq_n) then
max_seq_i = seq_i
max_seq_n = seq_n
end if
exit
end if
seq_n = seq_n + 1
end do
end do
if (seq_n > max_seq_n) then
max_seq_i = seq_i
max_seq_n = seq_n
end if
end subroutine get_max_sequence
subroutine init_obs_data(ierr)
integer, intent(out) :: ierr
integer :: i, cnt, norders
integer, dimension(max_nl0) :: orders
real(dp) :: sum_1, sum_2, sum_3, range, nmax
real(dp) :: x, y, isig2, sum_xy, sum_x, sum_y, sum_x2, sum_isig2, d
logical, parameter :: dbg = .false.
include 'formats'
ierr = 0
!call test_get_frequency_ratios
if (nl0 <= 0) return
sigmas_r02 = 0d0
if (chi2_seismo_r_010_fraction > 0) then
call get_frequency_ratios( &
.true., nl0, l0_obs, nl1, l1_obs, &
ratios_n, ratios_l0_first, ratios_l1_first, &
ratios_r01, ratios_r10)
end if
if (chi2_seismo_r_02_fraction > 0) then
call get_r02_frequency_ratios( &
.true., nl0, l0_obs, nl1, l1_obs, nl2, l2_obs, ratios_r02)
end if
if (delta_nu <= 0 .and. nl0 > 1 .and. l0_n_obs(1) > 0) then
sum_xy = 0
sum_x = 0
sum_y = 0
sum_x2 = 0
sum_isig2 = 0
do i=1,nl0
isig2 = 1d0/pow2(l0_obs_sigma(i))
x = dble(l0_n_obs(i))
y = l0_obs(i)
sum_xy = sum_xy + x*y*isig2
sum_x = sum_x + x*isig2
sum_y = sum_y + y*isig2
sum_x2 = sum_x2 + x*x*isig2
sum_isig2 = sum_isig2 + isig2
end do
d = sum_isig2*sum_x2 - sum_x*sum_x
delta_nu = (sum_isig2*sum_xy - sum_x*sum_y)/d
if (delta_nu_sigma <= 0) delta_nu_sigma = sqrt(sum_isig2/d)
end if
if (correction_factor <= 0) return
if (l0_n_obs(1) <= 0) then
if (delta_nu <= 0) then
write(*,*) 'must supply value for delta_nu'
ierr = -1
return
end if
! set l0_n_obs(i) to order of l0_obs(i)
range = l0_obs(nl0) - l0_obs(1)
norders = int(range/delta_nu + 0.5d0) + 1
nmax = (nu_max/delta_nu)*(delta_nu_sun/nu_max_sun)*22.6 - 1.6
l0_n_obs(1) = int(nmax - (norders-1)/2)
if (dbg) write(*,3) 'l0_n_obs(i)', 1, l0_n_obs(1), l0_obs(1)
do i=2,norders
l0_n_obs(i) = l0_n_obs(1) + &
int((l0_obs(i) - l0_obs(1))/delta_nu + 0.5)
if (dbg) write(*,3) 'l0_n_obs(i)', i, l0_n_obs(i), l0_obs(i)
end do
if (dbg) then
write(*,1) 'range', range
write(*,2) 'norders', norders
write(*,1) 'nmax', nmax
write(*,2) '(norders+1)/2', (norders+1)/2
write(*,2) 'l0_n_obs(1)', l0_n_obs(1)
write(*,*)
!stop
end if
end if
cnt = 0
sum_1 = 0
sum_2 = 0
do i=1,nl0
if (l0_obs(i) < 0) cycle
cnt = cnt + 1
sum_1 = sum_1 + l0_obs(i)
sum_2 = sum_2 + l0_n_obs(i)
end do
avg_nu_obs = sum_1/cnt
avg_radial_n = sum_2/cnt
if (dbg) then
write(*,1) 'avg_nu_obs', avg_nu_obs
write(*,1) 'avg_radial_n', avg_radial_n
write(*,2) 'cnt', cnt
write(*,1) 'sum_1', sum_1
write(*,1) 'sum_2', sum_2
write(*,*)
stop 'init_obs_data'
end if
end subroutine init_obs_data
real(dp) function interpolate_l0_inertia(freq) result(inertia)
real(dp), intent(in) :: freq
integer :: i
real(dp) :: alfa, beta
inertia = 0
if (nl0 == 0) return
if (freq <= l0_freq(1)) then
inertia = l0_inertia(1)
return
end if
if (freq >= l0_freq(nl0)) then
inertia = l0_inertia(nl0)
return
end if
do i=2,nl0
if (freq < l0_freq(i)) then
alfa = (freq - l0_freq(i-1))/(l0_freq(i) - l0_freq(i-1))
beta = 1d0 - alfa
inertia = alfa*l0_inertia(i) + beta*l0_inertia(i-1)
return
end if
end do
end function interpolate_l0_inertia
subroutine get_l0_freq_corr( &
a_div_r, b, nu_max, correction_factor, check_obs, &
nl0, l0_obs, l0_freq, l0_freq_corr, l0_inertia)
real(dp), intent(in) :: a_div_r, b, nu_max, correction_factor
logical, intent(in) :: check_obs ! if false, then l0_obs is not used
integer, intent(in) :: nl0
real(dp), intent(in), dimension(:) :: &
l0_obs, l0_freq, l0_inertia
real(dp), intent(out) :: l0_freq_corr(:)
integer :: i
real(dp) :: Qnl
do i = 1, nl0
if (check_obs) then
if (l0_obs(i) < 0) cycle
end if
Qnl = 1
l0_freq_corr(i) = l0_freq(i)
if (b > 0 .and. correction_factor > 0 .and. l0_freq_corr(i) > 0) &
l0_freq_corr(i) = l0_freq_corr(i) + &
correction_factor*(a_div_r/Qnl)*pow_cr(l0_freq(i)/nu_max,b)
end do
end subroutine get_l0_freq_corr
subroutine get_l1_freq_corr( &
a_div_r, b, nu_max, correction_factor, check_obs, &
nl1, l1_obs, l1_freq, l1_freq_corr, l1_inertia, l0_inertia)
real(dp), intent(in) :: a_div_r, b, nu_max, correction_factor
logical, intent(in) :: check_obs ! if false, then l1_obs is not used
integer, intent(in) :: nl1
real(dp), intent(in), dimension(:) :: &
l1_obs, l1_freq, l1_inertia, l0_inertia
real(dp), intent(out) :: l1_freq_corr(:)
integer :: i
real(dp) :: Qnl, interp_l0_inertia
include 'formats'
do i = 1, nl1
if (check_obs) then
if (l1_obs(i) < 0) cycle
end if
l1_freq_corr(i) = l1_freq(i)
if (b > 0 .and. correction_factor > 0 .and. l1_freq_corr(i) > 0 .and. nl0 > 0) then
interp_l0_inertia = interpolate_l0_inertia(l1_freq(i))
!write(*,2) 'l1_freq_corr: l0_inertia interp prev', i, interp_l0_inertia, &
! (l0_inertia(min(nl0,i)) + l0_inertia(min(nl0,i+1)))/2
Qnl = l1_inertia(i)/interp_l0_inertia
l1_freq_corr(i) = l1_freq_corr(i) + &
correction_factor*(a_div_r/Qnl)*pow_cr(l1_freq(i)/nu_max,b)
end if
end do
end subroutine get_l1_freq_corr
subroutine get_l2_freq_corr( &
a_div_r, b, nu_max, correction_factor, check_obs, &
nl2, l2_obs, l2_freq, l2_freq_corr, l2_inertia, l0_inertia)
real(dp), intent(in) :: a_div_r, b, nu_max, correction_factor
logical, intent(in) :: check_obs ! if false, then l2_obs is not used
integer, intent(in) :: nl2
real(dp), intent(in), dimension(:) :: &
l2_obs, l2_freq, l2_inertia, l0_inertia
real(dp), intent(out) :: l2_freq_corr(:)
integer :: i
real(dp) :: Qnl, interp_l0_inertia
include 'formats'
do i = 1, nl2
if (check_obs) then
if (l2_obs(i) < 0) cycle
end if
l2_freq_corr(i) = l2_freq(i)
if (b > 0 .and. correction_factor > 0 .and. l2_freq_corr(i) > 0 .and. nl0 > 0) then
interp_l0_inertia = interpolate_l0_inertia(l2_freq(i))
Qnl = l2_inertia(i)/interp_l0_inertia
!write(*,2) 'l2_freq_corr: l0_inertia interp prev', i, interp_l0_inertia, &
! l0_inertia(min(nl0,i+1))
l2_freq_corr(i) = l2_freq_corr(i) + &
correction_factor*(a_div_r/Qnl)*pow_cr(l2_freq(i)/nu_max,b)
end if
end do
end subroutine get_l2_freq_corr
subroutine get_l3_freq_corr( &
a_div_r, b, nu_max, correction_factor, check_obs, &
nl3, l3_obs, l3_freq, l3_freq_corr, l3_inertia, l0_inertia)
real(dp), intent(in) :: a_div_r, b, nu_max, correction_factor
logical, intent(in) :: check_obs ! if false, then l3_obs is not used
integer, intent(in) :: nl3
real(dp), intent(in), dimension(:) :: &
l3_obs, l3_freq, l3_inertia, l0_inertia
real(dp), intent(out) :: l3_freq_corr(:)
integer :: i
real(dp) :: Qnl, interp_l0_inertia
include 'formats'
do i = 1, nl3
if (check_obs) then
if (l3_obs(i) < 0) cycle
end if
l3_freq_corr(i) = l3_freq(i)
if (b > 0 .and. correction_factor > 0 .and. l3_freq_corr(i) > 0 .and. nl0 > 0) then
interp_l0_inertia = interpolate_l0_inertia(l3_freq(i))
Qnl = l3_inertia(i)/interp_l0_inertia
l3_freq_corr(i) = l3_freq_corr(i) + &
correction_factor*(a_div_r/Qnl)*pow_cr(l3_freq(i)/nu_max,b)
end if
end do
end subroutine get_l3_freq_corr
subroutine get_freq_corr
real(dp), allocatable :: l_freq(:), l_obs(:), &
l_obs_sigma(:), y(:), X(:,:)
real(dp) :: XtX(2,2), XtXi(2,2), Xty(2), z(2)
real(dp) :: detXtX
integer :: i, j, k, N
logical :: only_radial
! fetch all frequencies into one array
if (l1_freq(1) > 0d0) then
only_radial = .false.
else
only_radial = .true.
end if
if (only_radial) then
N = nl0
else
N = nl0 + nl1 + nl2 + nl3
end if
allocate(l_freq(N), l_obs(N), l_obs_sigma(N), y(N))
allocate(X(N,2))
N = 0
do i = 1, nl0
l_freq(N+i) = l0_freq(i)
l_obs(N+i) = l0_obs(i)
l_obs_sigma(N+i) = l0_obs_sigma(i)
X(N+i, 1) = l0_freq(i)**(-1)/l0_inertia(i)/l0_obs_sigma(i)
X(N+i, 2) = l0_freq(i)**3/l0_inertia(i)/l0_obs_sigma(i)
end do
if (.not. only_radial) then
N = nl0
do i = 1, nl1
l_freq(N+i) = l1_freq(i)
l_obs(N+i) = l1_obs(i)
l_obs_sigma(N+i) = l1_obs_sigma(i)
X(N+i, 1) = l1_freq(i)**(-1)/l1_inertia(i)/l1_obs_sigma(i)
X(N+i, 2) = l1_freq(i)**3/l1_inertia(i)/l1_obs_sigma(i)
end do
N = nl0 + nl1
do i = 1, nl2
l_freq(N+i) = l2_freq(i)
l_obs(N+i) = l2_obs(i)
l_obs_sigma(N+i) = l2_obs_sigma(i)
X(N+i, 1) = l2_freq(i)**(-1)/l2_inertia(i)/l2_obs_sigma(i)
X(N+i, 2) = l2_freq(i)**3/l2_inertia(i)/l2_obs_sigma(i)
end do
N = nl0 + nl1 + nl2
do i = 1, nl3
l_freq(N+i) = l3_freq(i)
l_obs(N+i) = l3_obs(i)
l_obs_sigma(N+i) = l3_obs_sigma(i)
X(N+i, 1) = l3_freq(i)**(-1)/l3_inertia(i)/l3_obs_sigma(i)
X(N+i, 2) = l3_freq(i)**3/l3_inertia(i)/l3_obs_sigma(i)
end do
end if
if (only_radial) then
N = nl0
else
N = nl0 + nl1 + nl2 + nl3
end if
y = (l_obs - l_freq)/l_obs_sigma
! write(*,*) y(nl0), y(nl1), y(nl2), y(nl3)
! write(*,*) X(nl0,1), X(nl1,1), X(nl2,1), X(nl3,1)
! create matrices
XtX(1,1) = 0d0
XtX(1,2) = 0d0
XtX(2,1) = 0d0
XtX(2,2) = 0d0
Xty(1) = 0d0
Xty(2) = 0d0
do i = 1, N
XtX(1,1) = XtX(1,1) + X(i, 1)*X(i, 1)
XtX(1,2) = XtX(1,2) + X(i, 1)*X(i, 2)
XtX(2,2) = XtX(2,2) + X(i, 2)*X(i, 2)
Xty(1) = Xty(1) + X(i, 1)*y(i)
Xty(2) = Xty(2) + X(i, 2)*y(i)
end do
XtX(2,1) = XtX(1,2) ! (X^T.X) is symmetric
! write(*,*) 'XtX ='
! write(*,*) XtX(1,1), XtX(1,2)
! write(*,*) XtX(2,1), XtX(2,2)
! invert XtX
XtXi(1,1) = XtX(2,2)
XtXi(2,2) = XtX(1,1)
XtXi(1,2) = -XtX(1,2)
XtXi(2,1) = -XtX(2,1)
detXtX = XtX(1,1)*XtX(2,2) - XtX(1,2)*XtX(2,1)
write(*,*) 'det(XtX) =', detXtX
XtXi = XtXi/detXtX
! write(*,*) 'XtXi ='
! write(*,*) XtXi(1,1), XtX(1,2)
! write(*,*) XtXi(2,1), XtX(2,2)
! evaluate matrix equation
z(1) = XtXi(1, 1)*Xty(1) + XtXi(1,2)*Xty(2)
z(2) = XtXi(2, 1)*Xty(1) + XtXi(2,2)*Xty(2)
write(*,*) 'a_-1, a_3 =', z
l0_freq_corr = l0_freq + z(1)*l0_freq**(-1)/l0_inertia &
+ z(2)*l0_freq**3/l0_inertia
l1_freq_corr = l1_freq + z(1)*l1_freq**(-1)/l1_inertia &
+ z(2)*l1_freq**3/l1_inertia
l2_freq_corr = l2_freq + z(1)*l2_freq**(-1)/l2_inertia &
+ z(2)*l2_freq**3/l2_inertia
l3_freq_corr = l3_freq + z(1)*l3_freq**(-1)/l3_inertia &
+ z(2)*l3_freq**3/l3_inertia
deallocate(l_freq, l_obs, l_obs_sigma, y)
deallocate(X)
! call get_l0_freq_corr( &
! a_div_r, correction_b, nu_max, correction_factor, .true., &
! nl0, l0_obs, l0_freq, l0_freq_corr, l0_inertia)
! call get_l1_freq_corr( &
! a_div_r, correction_b, nu_max, correction_factor, .true., &
! nl1, l1_obs, l1_freq, l1_freq_corr, l1_inertia, l0_inertia)
! call get_l2_freq_corr( &
! a_div_r, correction_b, nu_max, correction_factor, .true., &
! nl2, l2_obs, l2_freq, l2_freq_corr, l2_inertia, l0_inertia)
! call get_l3_freq_corr( &
! a_div_r, correction_b, nu_max, correction_factor, .true., &
! nl3, l3_obs, l3_freq, l3_freq_corr, l3_inertia, l0_inertia)
end subroutine get_freq_corr
! chi2 = chi2_seismo*chi2_seismo_fraction &
! + chi2_spectro*(1 - chi2_seismo_fraction)
real(dp) function get_chi2(s, max_el, trace_okay, ierr)
type (star_info), pointer :: s
integer, intent(in) :: max_el
logical, intent(in) :: trace_okay
integer, intent(out) :: ierr
integer :: i, n, chi2N1, chi2N2
real(dp) :: chi2term, Teff, logL, chi2sum1, chi2sum2, frac, &
model_r01, model_r10, model_r02
! calculate chi^2 following Brandao et al, 2011, eqn 11
include 'formats'
ierr = 0
chi2sum1 = 0
chi2N1 = 0
chi2_r_010_ratios = -1
chi2_r_02_ratios = -1
chi2_frequencies = 0
if (chi2_seismo_freq_fraction > 0) then
if (trace_okay .and. trace_chi2_seismo_frequencies_info) &
write(*,'(a30,a6,99(a20))') &
'chi2term l0', 'model number', 'i', 'chi2term', &
'l0_freq_corr(i)', 'l0_obs(i)', 'l0_obs_sigma(i)'
do i = 1, nl0
if (l0_obs(i) < 0) cycle
chi2term = &
pow2((l0_freq_corr(i) - l0_obs(i))/l0_obs_sigma(i))
if (trace_okay .and. trace_chi2_seismo_frequencies_info) &
write(*,'(a30,2i6,99(1pe20.10))') 'chi2term l0', &
s% model_number, i, chi2term, &
l0_freq_corr(i), l0_obs(i), l0_obs_sigma(i)
chi2sum1 = chi2sum1 + chi2term
chi2N1 = chi2N1 + 1
end do
if (max_el >= 1) then
do i = 1, nl1
if (l1_obs(i) < 0) cycle
chi2term = &
pow2((l1_freq_corr(i) - l1_obs(i))/l1_obs_sigma(i))
if (trace_okay .and. trace_chi2_seismo_frequencies_info) &
write(*,'(a30,2i6,99(1pe20.10))') 'chi2term l1', &
s% model_number, i, chi2term, &
l1_freq_corr(i), l1_obs(i), l1_obs_sigma(i)
chi2sum1 = chi2sum1 + chi2term
chi2N1 = chi2N1 + 1
end do
end if
if (max_el >= 2) then
do i = 1, nl2
if (l2_obs(i) < 0) cycle
chi2term = &
pow2((l2_freq_corr(i) - l2_obs(i))/l2_obs_sigma(i))
if (trace_okay .and. trace_chi2_seismo_frequencies_info) &
write(*,'(a30,2i6,99(1pe20.10))') 'chi2term l2', &
s% model_number, i, chi2term, &
l2_freq_corr(i), l2_obs(i), l2_obs_sigma(i)
chi2sum1 = chi2sum1 + chi2term
chi2N1 = chi2N1 + 1
end do
end if
if (max_el >= 3) then
do i = 1, nl3
if (l3_obs(i) < 0) cycle
chi2term = &
pow2((l3_freq_corr(i) - l3_obs(i))/l3_obs_sigma(i))
if (trace_okay .and. trace_chi2_seismo_frequencies_info) &
write(*,'(a30,2i6,99(1pe20.10))') 'chi2term l3', &
s% model_number, i, chi2term, &
l3_freq_corr(i), l3_obs(i), l3_obs_sigma(i)
chi2sum1 = chi2sum1 + chi2term
chi2N1 = chi2N1 + 1
end do
end if
num_chi2_seismo_terms = chi2N1
chi2_frequencies = chi2sum1 ! /max(1,chi2N1)
end if
if (chi2_seismo_r_010_fraction > 0 .and. max_el >= 1) then
if (ratios_n == 0) then
write(*,*) 'ERROR: chi2_seismo_r_010_fraction > 0 but cannot evaluate r_010'
ierr = -1
return
end if
chi2sum1 = 0
do i=1,ratios_n
model_r01 = interpolate_ratio_r010( &
l0_obs(i + ratios_l0_first), ratios_l0_first, l0_freq, model_ratios_r01)
if (trace_okay .and. trace_chi2_seismo_ratios_info) &
write(*,2) 'r01 obs, model, interp, model - interp', &
i, ratios_r01(i), model_ratios_r01(i + ratios_l0_first), model_r01, &
model_ratios_r01(i + ratios_l0_first) - model_r01
model_r10 = interpolate_ratio_r010( &
l1_obs(i + ratios_l1_first), ratios_l1_first, l1_freq, model_ratios_r10)
if (trace_okay .and. trace_chi2_seismo_ratios_info) &
write(*,2) 'r10 obs, model, interp, model - interp', &
i, ratios_r10(i), model_ratios_r10(i + ratios_l1_first), model_r10, &
model_ratios_r10(i + ratios_l1_first) - model_r10
chi2term = &
pow2((model_r01 - ratios_r01(i))/sigmas_r01(i)) + &
pow2((model_r10 - ratios_r10(i))/sigmas_r10(i))
chi2sum1 = chi2sum1 + chi2term
if (trace_okay .and. trace_chi2_seismo_ratios_info) &
write(*,2) 'chi2 ratios terms r01 r10', i, chi2term, &
pow2((model_r01 - ratios_r01(i))/sigmas_r01(i)), &
pow2((model_r10 - ratios_r10(i))/sigmas_r10(i)), &
model_r01, model_r10
end do
n = 2*ratios_n
chi2_r_010_ratios = chi2sum1 ! /max(1,n)
end if
if (chi2_seismo_r_02_fraction > 0 .and. max_el >= 2) then
chi2sum1 = 0
n = 0
do i=2,nl0
if (sigmas_r02(i) == 0d0) cycle
model_r02 = interpolate_ratio_r02( &
l0_obs(i + ratios_l0_first), l0_freq, model_ratios_r02)
if (trace_okay .and. trace_chi2_seismo_ratios_info) &
write(*,2) 'r02 obs, model, interp, model - interp', &
i, ratios_r02(i), model_ratios_r02(i), model_r02, &
model_ratios_r02(i) - model_r10
chi2sum1 = chi2sum1 + &
pow2((model_r02 - ratios_r02(i))/sigmas_r02(i))
n = n+1
end do
if (n == 0) then
write(*,*) 'ERROR: chi2_seismo_r_02_fraction > 0 but cannot evaluate r_02'
ierr = -1
return
end if
chi2_r_02_ratios = chi2sum1 ! /max(1,n)
end if
chi2_seismo = &
chi2_seismo_r_010_fraction*chi2_r_010_ratios + &
chi2_seismo_r_02_fraction*chi2_r_02_ratios + &
chi2_seismo_freq_fraction*chi2_frequencies + &
chi2_seismo_delta_nu_fraction*chi2_delta_nu + &
chi2_seismo_nu_max_fraction*chi2_nu_max
chi2sum2 = 0
chi2N2 = 0
if (Teff_sigma > 0 .and. include_Teff_in_chi2_spectro) then
Teff = s% Teff
chi2term = pow2((Teff - Teff_target)/Teff_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term Teff', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (logL_sigma > 0 .and. include_logL_in_chi2_spectro) then
logL = s% log_surface_luminosity
chi2term = pow2((logL - logL_target)/logL_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term logL', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (logg_sigma > 0 .and. include_logg_in_chi2_spectro) then
chi2term = pow2((logg - logg_target)/logg_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term logg', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (FeH_sigma > 0 .and. include_FeH_in_chi2_spectro) then
chi2term = pow2((FeH - FeH_target)/FeH_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term FeH', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (logR_sigma > 0 .and. include_logR_in_chi2_spectro) then
chi2term = pow2((logR - logR_target)/logR_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term logR', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (age_sigma > 0 .and. include_age_in_chi2_spectro) then
chi2term = pow2((s% star_age - age_target)/age_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term age', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (surface_Z_div_X_sigma > 0 .and. include_surface_Z_div_X_in_chi2_spectro) then
chi2term = pow2((surface_Z_div_X - surface_Z_div_X_target)/surface_Z_div_X_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term surface_Z_div_X', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (surface_He_sigma > 0 .and. include_surface_He_in_chi2_spectro) then
chi2term = pow2((surface_He - surface_He_target)/surface_He_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term surface_He', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (Rcz_sigma > 0 .and. include_Rcz_in_chi2_spectro) then
chi2term = pow2((Rcz - Rcz_target)/Rcz_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term Rcz', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (csound_rms_sigma > 0 .and. include_csound_rms_in_chi2_spectro) then
chi2term = pow2((csound_rms - csound_rms_target)/csound_rms_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term csound_rms', s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (my_var1_sigma > 0 .and. include_my_var1_in_chi2_spectro) then
chi2term = pow2((my_var1 - my_var1_target)/my_var1_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term ' // trim(my_var1_name), s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (my_var2_sigma > 0 .and. include_my_var2_in_chi2_spectro) then
chi2term = pow2((my_var2 - my_var2_target)/my_var2_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term ' // trim(my_var2_name), s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
if (my_var3_sigma > 0 .and. include_my_var3_in_chi2_spectro) then
chi2term = pow2((my_var3 - my_var3_target)/my_var3_sigma)
if (trace_okay .and. trace_chi2_spectro_info) &
write(*,2) 'chi2_spectro_term ' // trim(my_var3_name), s% model_number, chi2term
chi2sum2 = chi2sum2 + chi2term
chi2N2 = chi2N2 + 1
end if
num_chi2_spectro_terms = chi2N2
chi2_spectro = chi2sum2 ! /max(1,chi2N2)
frac = chi2_seismo_fraction
chi2 = frac*chi2_seismo + (1-frac)*chi2_spectro
get_chi2 = chi2
if (chi2_seismo_fraction < 0 .or. chi2_seismo_fraction > 1) then
write(*,1) 'chi2_seismo_fraction', chi2_seismo_fraction
stop
end if
!if (is_bad_num(chi2)) stop 'get_chi2'
end function get_chi2
end module astero_support
-------------- next part --------------
# the following started as a copy of work_standard_makefile
# with changes to rule for $(STAR)
include $(MESA_DIR)/utils/makefile_header
#################################################################
ifndef STAR
STAR = star
endif
# define the PGSTAR files for astero
ifeq ($(USE_PGSTAR),YES)
PGSTAR_OBJS = pgstar_astero_plots.o
else
PGSTAR_OBJS = pgstar_astero_plots_stub.o
endif
ifeq ($(USE_GYRE),YES)
GYRE_OBJS = gyre_support.o
LOAD_GYRE = -lgyre $(LOAD_LAPACK) $(LOAD_BLAS)
else
GYRE_OBJS = gyre_support_stub.o
LOAD_GYRE =
endif
OBJS = \
astero_data.o adipls_support.o $(GYRE_OBJS) \
astero_support.o run_star_extras.o $(PGSTAR_OBJS) \
extras_support.o run_star_extras_astero.o run_star.o run.o
WORK_DIR = ..
WORK_SRC_DIR = $(WORK_DIR)/src
STAR_JOB_DIR = $(MESA_DIR)/star/job
ASTERO_SRC_DIR = $(MESA_DIR)/star/astero/src
ASTERO_DEFAULTS_DIR = $(MESA_DIR)/star/astero/defaults
$(STAR) : $(OBJS) adipls_support_procs.o
$(LOADER) $(FCopenmp) -o $(WORK_DIR)/$(STAR) \
$(OBJS) $(LOAD_MESA_STAR) $(LOAD_GYRE) \
-ladipls adipls_support_procs.o
#################################################################
# change this as necessary. see utils/makefile_header for definitions.
WORK_COMPILE = $(FC) $(FCbasic) $(FCopenmp) $(FCchecks) $(FCdebug) \
-I$(MESA_INCLUDE_DIR) -c $(FCfree)
run.o: $(WORK_SRC_DIR)/run.f
$(WORK_COMPILE) $<
run_star_extras.o: $(WORK_SRC_DIR)/run_star_extras.f
$(WORK_COMPILE) $<
gyre_support.o: $(ASTERO_SRC_DIR)/gyre_support.f
$(WORK_COMPILE) $<
gyre_support_stub.o: $(ASTERO_SRC_DIR)/gyre_support_stub.f
$(WORK_COMPILE) $<
pgstar_astero_plots.o: $(ASTERO_SRC_DIR)/pgstar_astero_plots.f
$(WORK_COMPILE) $<
pgstar_astero_plots_stub.o: $(ASTERO_SRC_DIR)/pgstar_astero_plots_stub.f
$(WORK_COMPILE) $<
extras_support.o: $(ASTERO_SRC_DIR)/extras_support.f
$(WORK_COMPILE) $<
run_star_extras_astero.o: $(ASTERO_SRC_DIR)/run_star_extras_astero.f
$(WORK_COMPILE) $<
adipls_support_procs.o: $(ASTERO_SRC_DIR)/adipls_support_procs.f
$(WORK_COMPILE) $<
adipls_support.o: $(ASTERO_SRC_DIR)/adipls_support.f
$(WORK_COMPILE) -I$(ASTERO_DEFAULTS_DIR) $<
#astero_support.o: $(ASTERO_SRC_DIR)/astero_support.f
astero_support.o: $(WORK_SRC_DIR)/astero_support.f
$(WORK_COMPILE) -I$(ASTERO_DEFAULTS_DIR) $<
astero_data.o: $(ASTERO_SRC_DIR)/astero_data.f
$(WORK_COMPILE) -I$(ASTERO_DEFAULTS_DIR) $<
%.o: $(STAR_JOB_DIR)/%.f
$(WORK_COMPILE) $<
clean:
- at rm -f *.o *.mod $(WORK_DIR)/$(STAR)
remk:
- at rm -f run.o $(WORK_DIR)/$(STAR)
More information about the Mesa-users
mailing list