Commit c6691b82 authored by Matthias Steffen's avatar Matthias Steffen
Browse files

test_ff123.py revised

test_ff123.py produces the same results as IDL script test_ff123.pro
parent 2c87ba41
;pro test_ff1
;
; use ff1 to see how a perfect LTE Spite plateau look after non-LTE corrections
;
; --- stars with log g=4.0, Teff = 5500 ... 6500 K, [Fe/H]= -1.5 ... -3.0
;
; --- sample of metallixities ---
; --- Spite plateau ---
z_test= -3.0 + 0.01*findgen(151)
; --- Full range ---
;z_test= -3.0 + 0.02*findgen(151)
nz = n_elements(z_test)
;
; --- sample of gravities ---
; --- Spite plateau ---
g_test = [4.0]
; --- Full range ---
;g_test = [3.5,4.0,4.5]
ng = n_elements(g_test)
;
; --- sample of temperatures ---
; --- Spite plateau ---
t_test= 5500.0 + 5.0*findgen(201)
; --- Full range ---
;t_test= 5000.0 + 10.0*findgen(151)
nt= n_elements(t_test)
; --- sample of Li abundances ---
a_test = [2.2]
na= n_elements(a_test)
;
a1corr=fltarr(na,nt,ng,nz)
a3corr=fltarr(na,nt,ng,nz)
;
a1inp=fltarr(na,nt,ng,nz)
t1inp=fltarr(na,nt,ng,nz)
g1inp=fltarr(na,nt,ng,nz)
z1inp=fltarr(na,nt,ng,nz)
;
a3inp=fltarr(na,nt,ng,nz)
t3inp=fltarr(na,nt,ng,nz)
g3inp=fltarr(na,nt,ng,nz)
z3inp=fltarr(na,nt,ng,nz)
;
W1FF2=fltarr(na,nt,ng,nz)
W3FF2=fltarr(na,nt,ng,nz)
a1ff3=fltarr(na,nt,ng,nz)
a3ff3=fltarr(na,nt,ng,nz)
for iz=0,nz-1 do begin
for ig=0,ng-1 do begin
for it=0,nt-1 do begin
for ia=0,na-1 do begin
t1inp[ia,it,ig,iz]=t_test[it]
g1inp[ia,it,ig,iz]=g_test[ig]
z1inp[ia,it,ig,iz]=z_test[iz]
a1inp[ia,it,ig,iz]=a_test[ia]
corr1 = ff1z_1d(t_test[it], g_test[ig], z_test[iz], a_test[ia])
a1corr[ia,it,ig,iz] = corr1
t3inp[ia,it,ig,iz]=t_test[it]
g3inp[ia,it,ig,iz]=g_test[ig]
z3inp[ia,it,ig,iz]=z_test[iz]
a3inp[ia,it,ig,iz]=a_test[ia]
corr3 = ff1z_3d(t_test[it], g_test[ig], z_test[iz], a_test[ia])
a3corr[ia,it,ig,iz] = corr3
endfor
endfor
endfor
endfor
;
a4ps, filename='test_ff1.ps', /half
!p.font=0
!x.thick=3
!y.thick=3
;
title1='Spite plateau corrections'
xrng1=[6600,5400]
xrng2=[-3.2,-1.3]
;
;title1='Full range corrections'
;xrng1=[6600,4900]
;xrng2=[-3.2,0.2]
;
plot, [5500,6500], [1.8,2.6], xtitle='Teff [K]', ytitle='A(Li)', $
title=title1, xrange=xrng1, xstyle=1, ystyle=1, /nodata
oplot, xrng1, [2.2,2.2], line=2
oplot, t1inp, a1inp+a1corr, psym=psymuser(102,/fill), symsize=0.25,$
col=!ctv.blue
oplot, t1inp, a3inp+a3corr, psym=psymuser(102,/fill), symsize=0.25,$
col=!ctv.red
legend, ['FFI (1D_NLTE)', 'FFI (3D_NLTE)'], line=[0,0], thick=[3,1], $
psym=[8,8], color=[!ctv.blue,!ctv.red], xali=1.0, yali=0.0
;
plot, [-3.2,-1.3], [1.8,2.6], xtitle='[Fe/H]', ytitle='A(Li)', $
title=title1, xrange=xrng2, xstyle=1, ystyle=1, /nodata
oplot, xrng2, [2.2,2.2], line=2
oplot, z1inp, a1inp+a1corr, psym=psymuser(102,/fill), symsize=0.35,$
col=!ctv.blue
oplot, z1inp, a3inp+a3corr, psym=psymuser(102,/fill), symsize=0.35,$
col=!ctv.red
legend, ['FFI (1D_NLTE)', 'FFI (3D_NLTE)'], line=[0,0], thick=[3,1], $
psym=[8,8], color=[!ctv.blue,!ctv.red], xali=1.0, yali=0.0
;
!p.font=-1
!x.thick=1
!y.thick=1
;
a4ps, /close
;
;
end ; test_ff1
; pro test_ff123
;
qiso = 0.000D0
z_tab = [-3.0D0, -2.0D0, -1.0D0, -0.5D0, 0.0D0]
nz = n_elements(z_tab)
t_tab = [5000.0D0, 5500.0D0, 5900.0D0, 6300.0D0 ,6500.0D0]
nt = n_elements(t_tab)
g_tab = [3.5D0, 4.0D0, 4.5D0]
ng = n_elements(g_tab)
a_tab = 1.0D0 + 0.2D0*dindgen(11)
na = n_elements(a_tab)
;
get_lun, unit
openw, unit, 'test_ff123_idl.out'
;
t1=' A*(Li)'
t2=' del_1'
t3=' A1(Li)'
t4=' EW1'
t5=' A1_inv'
t6=' del_3'
t7=' A3(Li)'
t8=' EW3'
t9=' A3_inv'
tt= t1+t2+t3+t4+t5+t6+t7+t8+t9
l1='--------'
ll=strjoin(replicate(l1,9))
for iz=0,nz-1 do begin
zinp=z_tab[iz]
for ig=0,ng-1 do begin
ginp=g_tab[ig]
for it=0,nt-1 do begin
tinp=t_tab[it]
printf, unit, 'Teff =', tinp, ' K, logg =', ginp, ', [Fe/H] =', zinp, $
', 6Li/7Li =', qiso, format='(A,F8.2,A,F5.2,A,F6.2,A,F7.4)'
printf, unit, ll
printf, unit, tt
printf, unit, ll
for ia=0,na-1 do begin
ainp=a_tab[ia]
ff1, tinp, ginp, zinp, ainp, qiso, ff1_out
aa1=ainp+ff1_out[0]
ff2, tinp, ginp, zinp, aa1, qiso, ff2_out
ew1=ff2_out[0]
ff3, tinp, ginp, zinp, ew1, qiso, ff3_out
aa1_inv=ff3_out[0]
aa3=ainp+ff1_out[1]
ff2, tinp, ginp, zinp, aa3, qiso, ff2_out
ew3=ff2_out[1]
ff3, tinp, ginp, zinp, ew3, qiso, ff3_out
aa3_inv=ff3_out[1]
printf, unit, ainp, ff1_out[0], aa1, ew1, aa1_inv, $
ff1_out[1], aa3, ew3, aa3_inv, format='(9F8.4)'
endfor
printf, unit, ll
printf, unit
endfor
endfor
endfor
printf, unit
printf, unit
;
ntst=21
z_tst = -3.0D0 + 0.15D0*dindgen(ntst)
t_tst = 6500.0D0 - 75.0D0*dindgen(ntst)
g_tst = 3.5D0 + 0.05*dindgen(ntst)
na=11
a_tst = 1.0D0 + 0.2D0*dindgen(na)
;
for ix=0,ntst-1 do begin
tinp=t_tst[ix]
ginp=g_tst[ix]
zinp=z_tst[ix]
printf, unit, 'Teff =', tinp, ' K, logg =', ginp, ', [Fe/H] =', zinp, $
', 6Li/7Li =', qiso, format='(A,F8.2,A,F5.2,A,F6.2,A,F7.4)'
printf, unit, ll
printf, unit, tt
printf, unit, ll
for ia=0,na-1 do begin
ainp=a_tst[ia]
ff1, tinp, ginp, zinp, ainp, qiso, ff1_out
aa1=ainp+ff1_out[0]
ff2, tinp, ginp, zinp, aa1, qiso, ff2_out
ew1=ff2_out[0]
ff3, tinp, ginp, zinp, ew1, qiso, ff3_out
aa1_inv=ff3_out[0]
aa3=ainp+ff1_out[1]
ff2, tinp, ginp, zinp, aa3, qiso, ff2_out
ew3=ff2_out[1]
ff3, tinp, ginp, zinp, ew3, qiso, ff3_out
aa3_inv=ff3_out[1]
printf, unit, ainp, ff1_out[0], aa1, ew1, aa1_inv, $
ff1_out[1], aa3, ew3, aa3_inv, format='(9F8.4)'
endfor
printf, unit, ll
printf, unit
endfor
;
close, unit
free_lun, unit
;
end ; test_ff123
import numpy as np
from ff_all import ff1, ff2, ff3
from ff_all import ff1, ff2, ff3
#
t_test=np.array([5000.0, 5500.0, 5900.0, 6300.0, 6500.0])
nt=t_test.size
g_test=np.array([3.5, 4.0, 4.5])
ng=g_test.size
z_test=np.array([-3.0, -2.0, -1.0, -0.5, 0.0])
nz=z_test.size
a_test=np.array([1.0, 2.0, 3.0])
na=a_test.size
qiso = 0.000
z_tab = np.array([-3.0, -2.0, -1.0, -0.5, 0.0])
nz = z_tab.size
t_tab = np.array([5000.0, 5500.0, 5900.0, 6300.0 ,6500.0])
nt = t_tab.size
g_tab = np.array([3.5, 4.0, 4.5])
ng = g_tab.size
a_tab = 1.0 + 0.2*np.arange(11)
na = a_tab.size
#
qiso=0.042
print("qiso : {:8.3f}".format(qiso))
for it in range(nt):
Teff=t_test[it]
t1=' A*(Li)'
t2=' del_1'
t3=' A1(Li)'
t4=' EW1'
t5=' A1_inv'
t6=' del_3'
t7=' A3(Li)'
t8=' EW3'
t9=' A3_inv'
tt= t1+t2+t3+t4+t5+t6+t7+t8+t9
l1='--------'
ll=l1*9
fmt1='Teff ={0:8.2f} K, logg ={1:5.2f}, [Fe/H] ={2:6.2f}, 6Li/7Li ={3:7.4f}'
fmt2='{:8.4f}'*9
f=open('test_ff123_py.out', 'w')
for iz in range(nz):
zinp=z_tab[iz]
for ig in range(ng):
logg=g_test[ig]
for iz in range(nz):
z=z_test[iz]
ginp=g_tab[ig]
for it in range(nt):
tinp=t_tab[it]
print(fmt1.format(tinp, ginp, zinp, qiso), file=f)
print(ll, file=f)
print(tt, file=f)
print(ll, file=f)
for ia in range(na):
ALi=a_test[ia]
res1 =ff1( Teff, logg, z, ALi, qiso)
res2 =ff2( Teff, logg, z, ALi, qiso)
EW=np.mean(res2)
print("Teff : {:8.1f}".format(Teff), "K,",
"log g: {:8.1f},".format(logg),
"[Fe/H]: {:8.1f},".format(z),
"log EW: {:12.8f}".format(EW))
res3 =ff3( Teff, logg, z, EW, qiso)
print("Results of FFI, FFII, FFIII:")
txt1="delta_1DNLTE delta_3DNLTE"
txt2=" ALi_1DNLTE ALi_3DNLTE"
txt3="logEW_1DNLTE logEW_3DNLTE"
print(txt1,txt2,txt3)
print("{:12.8f} {:12.8f}".format(res2[0],res2[1]),
"{:12.8f} {:12.8f}".format(res1[0],res1[1]),
"{:12.8f} {:12.8f}".format(res3[0],res3[1]))
ainp=a_tab[ia]
ff1_out=ff1(tinp, ginp, zinp, ainp, qiso)
aa1=ainp+ff1_out[0]
ff2_out=ff2(tinp, ginp, zinp, aa1, qiso)
ew1=ff2_out[0]
ff3_out=ff3(tinp, ginp, zinp, ew1, qiso)
aa1_inv=ff3_out[0]
aa3=ainp+ff1_out[1]
ff2_out=ff2(tinp, ginp, zinp, aa3, qiso)
ew3=ff2_out[1]
ff3_out=ff3(tinp, ginp, zinp, ew3, qiso)
aa3_inv=ff3_out[1]
print(fmt2.format(ainp, ff1_out[0], aa1, ew1, aa1_inv,
ff1_out[1], aa3, ew3, aa3_inv), file=f)
print(ll, file=f)
print( file=f)
#
print( file=f)
print( file=f)
#
ntst=21
z_tst = -3.0 + 0.15*np.arange(ntst)
t_tst = 6500.0 - 75.0*np.arange(ntst)
g_tst = 3.5 + 0.05*np.arange(ntst)
na=11
a_tst = 1.0 + 0.2*np.arange(na)
#
for ix in range(ntst):
tinp=t_tst[ix]
ginp=g_tst[ix]
zinp=z_tst[ix]
print(fmt1.format(tinp, ginp, zinp, qiso), file=f)
print(ll, file=f)
print(tt, file=f)
print(ll, file=f)
for ia in range(na):
ainp=a_tst[ia]
ff1_out=ff1(tinp, ginp, zinp, ainp, qiso)
aa1=ainp+ff1_out[0]
ff2_out=ff2(tinp, ginp, zinp, aa1, qiso)
ew1=ff2_out[0]
ff3_out=ff3(tinp, ginp, zinp, ew1, qiso)
aa1_inv=ff3_out[0]
aa3=ainp+ff1_out[1]
ff2_out=ff2(tinp, ginp, zinp, aa3, qiso)
ew3=ff2_out[1]
ff3_out=ff3(tinp, ginp, zinp, ew3, qiso)
aa3_inv=ff3_out[1]
print(fmt2.format(ainp, ff1_out[0], aa1, ew1, aa1_inv,
ff1_out[1], aa3, ew3, aa3_inv), file=f)
print(ll, file=f)
print( file=f)
#
f.close()
;pro test_ff2
;
; compare input EWs to ff2 output at all grid points
;
; --- grid of metallicities ---
z_test = [0.0, -0.5, -1.0, -2.0, -3.0]
nz = n_elements(z_test)
; --- grid of gravities ---
g_test = [3.5, 4.0, 4.5]
ng = n_elements(g_test)
; --- grid of temperatures ---
t_test = [5000.0, 5500.0, 5900.0, 6300.0, 6500.0]
nt= n_elements(t_test)
; --- grid of Li abundances ---
a_test = [1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0]
na= n_elements(a_test)
;
a1inp=fltarr(na,nt,ng,nz)
t1inp=fltarr(na,nt,ng,nz)
g1inp=fltarr(na,nt,ng,nz)
z1inp=fltarr(na,nt,ng,nz)
W1inp=fltarr(na,nt,ng,nz)
;
a3inp=fltarr(na,nt,ng,nz)
t3inp=fltarr(na,nt,ng,nz)
g3inp=fltarr(na,nt,ng,nz)
z3inp=fltarr(na,nt,ng,nz)
W3inp=fltarr(na,nt,ng,nz)
;
W1ff2=fltarr(na,nt,ng,nz)
W3ff2=fltarr(na,nt,ng,nz)
;
; --- read input EWs ---
get_lun, unit
openr, unit, 'LogEW_1DNLTE.txt'
for it=0,nt-1 do begin
for ig=0,ng-1 do begin
for iz=0,nz-1 do begin
for ia=0,na-1 do begin
readf, unit, ttab, gtab, ztab, atab, wtab
t1inp[ia,it,ig,iz]=ttab
g1inp[ia,it,ig,iz]=gtab
z1inp[ia,it,ig,iz]=ztab
a1inp[ia,it,ig,iz]=atab
W1inp[ia,it,ig,iz]=10^wtab
w1x=ff2_1d(ttab, gtab, ztab, atab)
W1ff2[ia,it,ig,iz]=10^w1x
endfor
endfor
endfor
endfor
close, unit
openr, unit, 'LogEW_3DNLTE.txt'
for it=0,nt-1 do begin
for ig=0,ng-1 do begin
for iz=0,nz-1 do begin
for ia=0,na-1 do begin
readf, unit, ttab, gtab, ztab, atab, wtab
t3inp[ia,it,ig,iz]=ttab
g3inp[ia,it,ig,iz]=gtab
z3inp[ia,it,ig,iz]=ztab
a3inp[ia,it,ig,iz]=atab
W3inp[ia,it,ig,iz]=10^wtab
w3x=ff2_3d(ttab, gtab, ztab, atab)
W3ff2[ia,it,ig,iz]=10^w3x
endfor
endfor
endfor
endfor
close, unit
free_lun, unit
;
err1 = (W1ff2-W1inp)/W1inp
err3 = (W3ff2-W3inp)/W3inp
;
a4ps, filename='test_ff2.ps', /half
!p.font=0
!x.thick=3
!y.thick=3
;
plot, 100.0*err1, xtitle='index', ytitle='EW_FF2-EW_inp [%]', thick=3
oplot, 100.0*err3, col=!ctv.red
legend, ['1D_NLTE', '3D_NLTE'], line=[0,0], thick=[3,1], $
color=[!ctv.black,!ctv.red]
;
!p.font=-1
!x.thick=1
!y.thick=1
;
a4ps, /close
;
err1rms = sqrt(avg(err1^2))
err3rms = sqrt(avg(err3^2))
; --- Maximum error ---
ix1 = where(ABS(err1) eq MAX(ABS(err1)))
ix3 = where(ABS(err3) eq MAX(ABS(err3)))
;
get_lun, unit
openw, unit, 'test_ff2.out'
;
printf, unit, 'RMS mean error of 1D EWs = ' +string(100.0*err1rms) + ' %'
printf, unit, 'Maximum error of 1D EWs = ' +string(100.0*err1[ix1]) + ' %'
printf, unit, 'at parameters'
printf, unit, 'Teff ='+ string(t1inp[ix1])
printf, unit, 'log g ='+ string(g1inp[ix1])
printf, unit, '[Fe/H]='+ string(z1inp[ix1])
printf, unit, 'A(Li) ='+ string(a1inp[ix1])
;
printf, unit, 'RMS mean error of 3D EWs = ' +string(100.0*err3rms) + ' %'
printf, unit, 'Maximum error of 3D EWs = ' +string(100.0*err3[ix3]) + ' %'
printf, unit, 'at parameters'
printf, unit, 'Teff ='+ string(t3inp[ix3])
printf, unit, 'log g ='+ string(g3inp[ix3])
printf, unit, '[Fe/H]='+ string(z3inp[ix3])
printf, unit, 'A(Li) ='+ string(a3inp[ix3])
;
close, unit
free_lun, unit
;
end ; test_ff2
;pro test_ff3
;
; compare ff3 output with known A(Li) of input EWs at all grid points
;
; --- grid of metallicities ---
z_test = [0.0, -0.5, -1.0, -2.0, -3.0]
nz = n_elements(z_test)
; --- grid of gravities ---
g_test = [3.5, 4.0, 4.5]
ng = n_elements(g_test)
; --- grid of temperatures ---
t_test = [5000.0, 5500.0, 5900.0, 6300.0, 6500.0]
nt= n_elements(t_test)
; --- grid of Li abundances ---
a_test = [1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0]
na= n_elements(a_test)
;
a1inp=fltarr(na,nt,ng,nz)
t1inp=fltarr(na,nt,ng,nz)
g1inp=fltarr(na,nt,ng,nz)
z1inp=fltarr(na,nt,ng,nz)
W1inp=fltarr(na,nt,ng,nz)
;
a3inp=fltarr(na,nt,ng,nz)
t3inp=fltarr(na,nt,ng,nz)
g3inp=fltarr(na,nt,ng,nz)
z3inp=fltarr(na,nt,ng,nz)
W3inp=fltarr(na,nt,ng,nz)
;
a1ff3=fltarr(na,nt,ng,nz)
a3ff3=fltarr(na,nt,ng,nz)
;
; --- read inout EWs ---
get_lun, unit
openr, unit, 'LogEW_1DNLTE.txt'
for it=0,nt-1 do begin
for ig=0,ng-1 do begin
for iz=0,nz-1 do begin
for ia=0,na-1 do begin
readf, unit, ttab, gtab, ztab, atab, wtab
t1inp[ia,it,ig,iz]=ttab
g1inp[ia,it,ig,iz]=gtab
z1inp[ia,it,ig,iz]=ztab
a1inp[ia,it,ig,iz]=atab
W1inp[ia,it,ig,iz]=10^wtab
a1x=ff3_1d(ttab, gtab, ztab, 10^wtab)
a1ff3[ia,it,ig,iz]=a1x
endfor
endfor
endfor
endfor
close, unit
openr, unit, 'LogEW_3DNLTE.txt'
for it=0,nt-1 do begin
for ig=0,ng-1 do begin
for iz=0,nz-1 do begin
for ia=0,na-1 do begin
readf, unit, ttab, gtab, ztab, atab, wtab
t3inp[ia,it,ig,iz]=ttab
g3inp[ia,it,ig,iz]=gtab
z3inp[ia,it,ig,iz]=ztab
a3inp[ia,it,ig,iz]=atab
W3inp[ia,it,ig,iz]=10^wtab
a3x=ff3_3d(ttab, gtab, ztab, 10^wtab)
a3ff3[ia,it,ig,iz]=a3x
endfor
endfor
endfor
endfor
close, unit
free_lun, unit
;
err1 = a1ff3-a1inp
err3 = a3ff3-a3inp
;
a4ps, filename='test_ff3.ps', /half
!p.font=0
!x.thick=3
!y.thick=3
;
plot, err1, xtitle='index', ytitle='A(Li)_FF3 - A(Li)_inp [dex]', thick=3
oplot, err3, col=!ctv.red
legend, ['1D_NLTE', '3D_NLTE'], line=[0,0], thick=[3,1], $
color=[!ctv.black,!ctv.red], xali=1.0, yali=0.0
;
!p.font=-1
!x.thick=1
!y.thick=1
;
a4ps, /close
;
;
err1rms = sqrt(avg(err1^2))
err3rms = sqrt(avg(err3^2))
; --- Maximum error ---
ix1 = where(ABS(err1) eq MAX(ABS(err1)))
ix3 = where(ABS(err3) eq MAX(ABS(err3)))
;
get_lun, unit
openw, unit, 'test_ff3.out'
;
printf, unit, 'RMS mean error of 1D A(Li) = ' +string(err1rms) + ' dex'
printf, unit, 'Maximum error of 1D A(Li) = ' +string(err1[ix1]) + ' dex'
printf, unit, 'at parameters'
printf, unit, 'Teff ='+ string(t1inp[ix1])
printf, unit, 'log g ='+ string(g1inp[ix1])
printf, unit, '[Fe/H]='+ string(z1inp[ix1])
printf, unit, 'A(Li) ='+ string(a1inp[ix1])
printf, unit, 'Ew = '+ string(w1inp[ix1])
;
printf, unit, 'RMS mean error of 3D A(Li) = ' +string(err3rms) + ' dex'
printf, unit, 'Maximum error of 3D A(Li) = ' +string(err3[ix3]) + ' dex'
printf, unit, 'at parameters'
printf, unit, 'Teff ='+ string(t3inp[ix3])
printf, unit, 'log g ='+ string(g3inp[ix3])
printf, unit, '[Fe/H]='+ string(z3inp[ix3])
printf, unit, 'A(Li) ='+ string(a1inp[ix3])
printf, unit, 'Ew = '+ string(w3inp[ix3])
;
close, unit
free_lun, unit
;
end ; test_ff3