Commit 4105594e authored by Janosch Preuß's avatar Janosch Preuß
Browse files

add repro for dtn-freq analysis for helio

parent e0889e1f
Pipeline #233785 passed with stage
in 14 minutes and 16 seconds
......@@ -257,7 +257,12 @@ runs postprocessing steps for computing the cross-covariances.
When the computations have finished the following output data will be contained in the folder.
* Fig 6.2 (__a__): The power spectrum based on the Atmo model *power-computed-Atmo.png*.
* Fig. 6.2 (__b__): The power spectrum based on the VAL-C model *power-computed-VALC.png*.
* Fig. 6.2 (__b__): The power spectrum based on the VAL-C model *power-computed-VALC.png*.
* Fig. 6.2 (__c__): Samples of absolute value of dtn for for the Atmo model as function of frequency and harmonic degree *abs-dtn-Atmo-model.png*.
* Fig. 6.2 (__d__): Samples of absolute value of dtn for for the VAL-C model as function of frequency and harmonic degree *abs-dtn-VALC-model.png*.
* Fig. 6.2 (__e__): The values of the dtn function for the Atmo model at l=200 are available in the file *dtn-freq-Atmo-model-l200.dat*. The first column *mHz* is the frequency in mHz.
The second column *dtnReal* contains the real part of dtn as a function of frequency while the last column *dtnImag* contains the imaginary part.
* Fig. 6.2 (__f__): The values of the dtn function for the VAL-C model at l=200 are available in the file *dtn-freq-VALC-model-l200.dat*. The columns are labelled as in Fig. 6.2 (__e__).
* Fig. 6.5 (__b__): Time-distance diagram based on MDI observations *Time-Distance-observed.png*.
* Fig. 6.5 (__c__): Time-distance diagram based on Atmo model *Time-Distance-Atmo.png*.
* Fig. 6.5 (__d__): Time-distance diagram based on VAL-C model *Time-Distance-VALC.png*.
......@@ -277,6 +282,12 @@ If you would like to reproduce the latex plots of Fig. 6.4 and Fig. 6.6 with the
copy this data from the folder `scripts/6_helio_tbc/obs_helio` into the folder `data/6_helio_tbc/obs_helio` and follow the instructions for Fig. 6.4 and Fig. 6.6 below
for compiling the latex figures.
### Fig. 6.2
Change to directory `plots/6_helio_tbc/obs_helio`.
* (__e__): Run `lualatex dtn-freq-Atmo.tex`.
* (__f__): Run `lualatex dtn-freq-VALC.tex`.
### Fig. 6.3
Change to directory `plots/6_helio_tbc/obs_helio`.
......@@ -815,7 +826,11 @@ the script `power-spectrum.py` consists of two steps.
python3 power-spectrum.py "VALC"
python3 power-spectrum.py "Atmo"
Afterwards please continue with the instructions given for the case of already provided data, see e.g.[here](#IntroHelioDownload). The "download" flag
First of all, these script actually compute the poles of the dtn functions with respect to the frequency which are displayed in the lower
panels of Fig. 6.2 (__e__) and Fig. 6.2 (__f__). They are stored in the files *dtn-freq-Atmo-model-l200.dat* and *dtn-freq-VALC-model-l200.dat*.
The columns *poles_real* and *poles_imag* contain the real and imaginary parts of the poles, respectively.
After finishing these computations please continue with the instructions given for the case of already provided data, see e.g.[here](#IntroHelioDownload). The "download" flag
in these instructions can now be omitted.
### Fig. 6.3 (__b__)
......
poles_real poles_imag
5.178914558324499318e+00 -2.412779778345396694e-03
5.179646709982960395e+00 -1.531267507390972988e-03
5.179263447391017294e+00 -2.011903501403160960e-03
5.180096309913728270e+00 -9.046586145244077306e-04
5.178587969158392923e+00 -2.759382765158180633e-03
poles_real poles_imag
5.514461047255593051e+00 -1.134824217875595630e-07
6.627343821549520619e+00 -3.147537491842500908e-08
4.248551745685618108e+00 -2.761798749118265158e-07
8.039810098007452055e+00 -1.945924054384935379e-07
\documentclass[tikz]{standalone}
\usetikzlibrary{spy,shapes,shadows,calc}
\usepackage{amsmath}
\usepackage{physics}
\usepackage{pgfplots}
\pgfplotsset{compat=1.3}
\usepackage{amsmath}
\DeclareFontFamily{OT1}{pzc}{}
\DeclareFontShape{OT1}{pzc}{m}{it}{<-> s * [1.10] pzcmi7t}{}
\DeclareMathAlphabet{\mathpzc}{OT1}{pzc}{m}{it}
\pgfplotsset{
legend style = {font=\small}
}
\begin{document}
\begin{tikzpicture}[scale = 0.8]
\begin{axis}[
name = dtnplot,
height = 5.5cm,
width = 8.5cm,
every axis plot/.append style={thick},
axis y line*=left,
xlabel= { mHz },
scaled y ticks={base 10:-3},
%legend pos = north east,
legend style={at={(0.05,0.7)},anchor=west},
xmin = 0,
xmax = 8.3,
ymin = -3500,
ymax = 3500,
restrict y to domain=-1000000:1000000,
label style={at={(axis description cs:0.5,-0.08)},anchor=north},
]
\addplot[blue,line width=2pt]
table[mark=none,x=mHz,y=dtnReal] {../../../data/6_helio_tbc/obs_helio/dtn-freq-Atmo-model-l200.dat};
\addplot[red,line width=2pt]
table[mark=none,x=mHz,y=dtnImag] {../../../data/6_helio_tbc/obs_helio/dtn-freq-Atmo-model-l200.dat};
\legend{$\Re $ ,$\Im$}
\end{axis}
\begin{axis}[
name=polesplot,
at={($(dtnplot.south)-(0,1cm)$)},
anchor = north,
height = 3.5cm,
width = 8.5cm,
every axis plot/.append style={thick},
axis y line*=left,
xticklabels = {,,},
ylabel = {$\Im$},
y label style={at={(axis description cs:0.05,0.8)},anchor=east},
xlabel = {$\Re$},
x label style={at={(axis description cs:0.4,0.15)},anchor=east},
legend pos = north east,
xmin = 0,
xmax = 8.3,
ymin = -3e-3,
ymax = 1e-7,
]
\addplot[lightgray,very thick,mark=x, only marks,mark options={scale=2}]
table[mark=none,x=poles_real,y=poles_imag] {../../../data/6_helio_tbc/obs_helio/poles_freq_dtn_Atmo_l_200.dat};
\addplot[mark=none,black,samples=2,very thin,domain=0:8.3] {0.0};
\legend{poles}
\end{axis}
\end{tikzpicture}
\end{document}
\documentclass[tikz]{standalone}
\usetikzlibrary{spy,shapes,shadows,calc}
\usepackage{amsmath}
\usepackage{physics}
\usepackage{pgfplots}
\pgfplotsset{compat=1.3}
\usepackage{amsmath}
\DeclareFontFamily{OT1}{pzc}{}
\DeclareFontShape{OT1}{pzc}{m}{it}{<-> s * [1.10] pzcmi7t}{}
\DeclareMathAlphabet{\mathpzc}{OT1}{pzc}{m}{it}
\pgfplotsset{
legend style = {font=\small}
}
\begin{document}
\begin{tikzpicture}[scale = 0.8]
\begin{axis}[
name = dtnplot,
height = 5.5cm,
width = 8.5cm,
every axis plot/.append style={thick},
axis y line*=left,
xlabel= { mHz },
%legend pos = north east,
%y tick label style={xshift=2em},
scaled ticks=true,
scaled y ticks={base 10:-3},
legend style={at={(0.05,0.675)},anchor=west},
xmin = 0,
xmax = 8.3,
ymin = -3500,
ymax = 3500,
restrict y to domain=-1000000:1000000,
label style={at={(axis description cs:0.5,-0.08)},anchor=north},
]
\addplot[blue,line width=2pt]
table[mark=none,x=mHz,y=dtnReal] {../../../data/6_helio_tbc/obs_helio/dtn-freq-VALC-model-l200.dat};
\addplot[red,line width=2pt]
table[mark=none,x=mHz,y=dtnImag] {../../../data/6_helio_tbc/obs_helio/dtn-freq-VALC-model-l200.dat};
\legend{$\Re $ ,$\Im$}
\end{axis}
\begin{axis}[
name=polesplot,
at={($(dtnplot.south)-(0,1cm)$)},
anchor = north,
height = 3.5cm,
width = 8.5cm,
every axis plot/.append style={thick},
axis y line*=left,
xticklabels = {,,},
ylabel = {$\Im$},
y label style={at={(axis description cs:0.05,0.8)},anchor=east},
xlabel = {$\Re$},
x label style={at={(axis description cs:0.4,0.15)},anchor=east},
legend pos = north east,
xmin = 0,
xmax = 8.3,
ymin = -5e-7,
ymax = 5e-7,
]
\addplot[lightgray,very thick,mark=x, only marks,mark options={scale=2}]
table[mark=none,x=poles_real,y=poles_imag] {../../../data/6_helio_tbc/obs_helio/poles_freq_dtn_VALC_l_200.dat};
\addplot[mark=none,black,samples=2,very thin,domain=0:8.3] {0.0};
\legend{poles}
\end{axis}
\end{tikzpicture}
\end{document}
......@@ -29,6 +29,7 @@ figures:
cd ${current_dir}/6_helio_tbc/obs_helio; lualatex dtn-approx-VALC-3p5mHz-damping-real.tex; lualatex dtn-approx-VALC-3p5mHz-damping-imag.tex; lualatex dtn-atmo-vs-VALC-8mHz-real.tex; lualatex dtn-atmo-vs-VALC-8mHz-imag.tex; latexmk -c
cd ${current_dir}/6_helio_tbc/obs_helio; lualatex FWHM_expl.tex; lualatex power-law-damping-f.tex; lualatex power-cut.tex; lualatex cross-tt.tex; lualatex rays-modelS.tex; lualatex double-ridge-filter.tex; latexmk -c
cd ${current_dir}/6_helio_tbc/tbc_Atmo; lualatex dtn-Atmo-3mHz.tex; lualatex dtn-Atmo-5p3mHz.tex; lualatex dtn-Atmo-6p5mHz.tex; latexmk -c
cd ${current_dir}/6_helio_tbc/obs_helio; lualatex dtn-freq-Atmo.tex; lualatex dtn-freq-VALC.tex; latexmk -c
cd ${current_dir}/6_helio_tbc/tbc_Atmo; lualatex cross-covariance-low-theta14.tex; lualatex cross-covariance-low-theta14.tex; lualatex cross-covariance-low-theta28.tex; lualatex cross-covariance-low-theta28.tex;
cd ${current_dir}/6_helio_tbc/tbc_Atmo; lualatex cross-covariance-high-theta36.tex; lualatex cross-covariance-high-theta36.tex; lualatex cross-covariance-high-theta72.tex; lualatex cross-covariance-high-theta72.tex; latexmk -c
cd ${current_dir}/6_helio_tbc/tbc_VALC; lualatex cross-covariance-low-theta.tex; lualatex cross-covariance-high-theta.tex; latexmk -c
......
......@@ -5,11 +5,23 @@ import matplotlib.colors as colors
import sys
import urllib.request
plt.rc('legend', fontsize=18)
plt.rc('axes', titlesize=18)
plt.rc('axes', labelsize=18)
plt.rc('xtick', labelsize=14)
plt.rc('ytick', labelsize=14)
plt.rcParams['savefig.dpi'] = 300
plt.rcParams["figure.dpi"] = 100
if ( len(sys.argv) > 1 and sys.argv[1] == "download"):
print("Downloading data from Gro.data catalogue")
file_atmo, header_atmo = urllib.request.urlretrieve("https://test.data.gro.uni-goettingen.de/api/access/datafile/:persistentId?persistentId=hdl:21.T11987/UCVKVX/OWKZL8","power-spectrum-Atmo-whitw.out")
file_VALC, header_VALC = urllib.request.urlretrieve("https://test.data.gro.uni-goettingen.de/api/access/datafile/:persistentId?persistentId=hdl:21.T11987/UCVKVX/LBCRB5","power-spectrum-VALC-meshed.out")
file_MDI, header_MDI = urllib.request.urlretrieve("https://test.data.gro.uni-goettingen.de/api/access/datafile/:persistentId?persistentId=hdl:21.T11987/UCVKVX/YV8MGG","mdi_medium_l.np")
file_atmo_dtn, header_atmo_dtn = urllib.request.urlretrieve("https://test.data.gro.uni-goettingen.de/api/access/datafile/:persistentId?persistentId=hdl:21.T11987/UCVKVX/ERLOZT","dtn_nr_Whitw_intro.out")
file_VALC_dtn, header_VALC_dtn = urllib.request.urlretrieve("https://test.data.gro.uni-goettingen.de/api/access/datafile/:persistentId?persistentId=hdl:21.T11987/UCVKVX/OJAAEF","dtn_nr_VALC_intro.out")
print("done")
......@@ -39,14 +51,53 @@ index_1mhz = np.where( f_hzs > 1.0)[0][0]
#f_hzs = np.linspace(0.001,8.3,300)
f_hzs = (10**-3)*f_hzs
plt.rc('legend', fontsize=18)
plt.rc('axes', titlesize=18)
plt.rc('axes', labelsize=18)
plt.rc('xtick', labelsize=14)
plt.rc('ytick', labelsize=14)
dtn_nr_VALC = np.loadtxt('dtn_nr_VALC_intro.out',dtype="complex",delimiter=',')
dtn_nr_Atmo = np.loadtxt("dtn_nr_Whitw_intro.out",dtype="complex",delimiter=',')
idx_l = 200
for dtn_nr,dtn_str in zip([dtn_nr_VALC, dtn_nr_Atmo],["VALC","Atmo"]):
#plt.figure(figsize=(6,5),constrained_layout=True)
plt.figure(figsize=(6,3),constrained_layout=True)
#im = plt.pcolormesh(L_range,10**3*f_hzs,np.abs(dtn_nr)/np.max( np.abs(dtn_nr) ),
#im = plt.pcolormesh(L_range,10**3*f_hzs,np.imag(dtn_nr) ,
im = plt.pcolormesh(L_range,10**3*f_hzs,np.abs(dtn_nr) ,
cmap=plt.get_cmap('viridis'),
norm=colors.LogNorm( vmin= 302 ,vmax= 14792 )
#norm=colors.LogNorm( vmin= np.abs(dtn_nr).min() ,vmax= np.abs(dtn_nr).max() )
)
plt.axvline(x=idx_l ,color='r',linestyle='--',linewidth=3,ymin= 10**3*f_hzs[10] , ymax= 10**3*f_hzs[-10] )
plt.rc('legend', fontsize=16)
plt.rc('axes', titlesize=16)
plt.rc('axes', labelsize=16)
plt.xlabel("Harmonic degree $\ell$")
plt.ylabel("Frequency (mHz)")
#plt.title("Model S + PML")
#plt.title("Model S + smoothed VAL-C + PML")
cbar = plt.colorbar(im,)
#cbar.set_ticks([1e3,5e3,1e4])
#cbar.set_ticklabels(['$10^3$','$10^4$'])
plt.savefig("abs-dtn-{0}-model.png".format(dtn_str),transparent=True,dpi=300)
plt.show()
plt.plot(10**3*f_hzs,dtn_nr[:,idx_l].real,label="real")
plt.plot(10**3*f_hzs,dtn_nr[:,idx_l].imag,label="imag")
#plt.title("Real")
plt.show()
#plt.plot(10**3*f_hzs,dtn_nr[:,idx_l].imag)
#plt.title("Imag")
#plt.show()
fname_dtn = "dtn-freq-{0}-model-l{1}.dat".format(dtn_str,idx_l)
np.savetxt(fname=fname_dtn,
X=np.transpose([ 10**3*f_hzs, dtn_nr[:,idx_l].real , dtn_nr[:,idx_l].imag ]),
header = 'mHz dtnReal dtnImag',
comments=''
)
plt.rcParams['savefig.dpi'] = 300
plt.rcParams["figure.dpi"] = 100
####################################################
# load reference power spectrum
......
......@@ -69,6 +69,8 @@ if atmospheric_model == "VALC":
else:
a = 1.0007126
l_idx = 200 #output dtn at freq for l = l_idx
Nmax = 5 # maximal dofs for learned IE in radial direction
alpha_decay = 1/60 # exponential decay factor of weights
Nrs = list(range(Nmax))
......@@ -552,6 +554,74 @@ def get_ODE_DtN(dtn_nr,use_PML):
np.savetxt('ODE_dtn_nr.out',dtn_nr,delimiter=',')
def get_ODE_DtN_omega_squared(omega_squared,lam_continuous,use_PML=False):
dtn_nr = []
mesh_1D = Make1DMesh_givenpts(mesh_above_surface)
fes_DtN = H1(mesh_1D, complex=True, order=order_ODE, dirichlet=[1])
u_DtN,v_DtN = fes_DtN.TnT()
gfu_DtN = GridFunction (fes_DtN)
f_DtN = LinearForm(fes_DtN)
f_DtN.Assemble()
r_DtN = f_DtN.vec.CreateVector()
frees_DtN = [i for i in range(1,fes_DtN.ndof)]
if use_PML:
pml_phys = get_PML(f_hz)
mesh_1D.SetPML(pml_phys,definedon=1)
for idx_lami,lami in enumerate(lam_continuous):
pot_1D = pot_rho_B - omega_squared/pot_c_B**2
a_DtN = BilinearForm (fes_DtN, symmetric=False)
a_DtN += SymbolicBFI((x**2)*grad(u_DtN)*grad(v_DtN))
a_DtN += SymbolicBFI(( lami.item()*a**2 + pot_1D*x**2 )*u_DtN*v_DtN)
a_DtN.Assemble()
gfu_DtN.vec[:] = 0.0
gfu_DtN.Set(1.0, BND)
r_DtN.data = f_DtN.vec - a_DtN.mat * gfu_DtN.vec
gfu_DtN.vec.data += a_DtN.mat.Inverse(freedofs=fes_DtN.FreeDofs(),inverse=solver) * r_DtN
rows_DtN,cols_DtN,vals_DtN = a_DtN.mat.COO()
A_DtN = csr_matrix((vals_DtN,(rows_DtN,cols_DtN)))
val = -P_DoFs(A_DtN,[0],frees_DtN).dot(gfu_DtN.vec.FV().NumPy()[frees_DtN])-P_DoFs(A_DtN,[0],[0]).dot(gfu_DtN.vec.FV().NumPy()[0])
dtn_nr.append(-val.item()/a**2)
if use_PML:
mesh_1D.UnSetPML(1)
return dtn_nr
if atmospheric_model == "VALC":
from GRPF import Poles_and_Roots
from cmath import sqrt as csqrt
def fun(z):
DtN_list = get_ODE_DtN_omega_squared(z,[lam_continuous[l_idx]],use_PML=False)
return DtN_list[0]
param = {} # for searching poles
param["Tol"] = 1e-9
param["visual"] = 1
param["ItMax"] = 25
param["NodesMax"] = 500000
param["xrange"] = [305983634275718.6, 1317450785328390.8]
param["yrange"] = [ -5e12 , 5e12]
param["h"] = 3e11
result = Poles_and_Roots(fun,param)
poles = 10**3*np.sqrt(result["poles"])/ (2*pi*RSun)
print("poles = ", poles)
fname_analytic_poles = "poles_freq_dtn_VALC_l_{0}.dat".format(l_idx)
header_str_analytic_poles = "poles_real poles_imag"
plot_collect_analytic_poles = [poles.real,poles.imag]
np.savetxt(fname= fname_analytic_poles ,
X= np.transpose(plot_collect_analytic_poles ) ,
header= header_str_analytic_poles ,
comments='')
def get_atmo_dtn(dtn_nr):
print("Getting Atmo dtn numbers from whittaker function")
......@@ -585,6 +655,59 @@ def get_atmo_dtn(dtn_nr):
np.savetxt('dtn_nr_Whitw.out',dtn_nr,delimiter=',')
def get_atmo_dtn_sigma_squared(sigma_squared):
from ngs_arb import Whitw_deriv_over_Whitw,Whitw_deriv_over_Whitw_par,Whitw_deriv_over_Whitw_spar
relevant_ls = [l_idx]
k_squared = sigma_squared/c_infty**2 - 0.25*alpha_infty**2
k_omega = sqrt_2branch(k_squared)
eta = alpha_infty/(2*k_omega)
chi = -1j*eta
z_arg = -2*1j*k_omega*a
kappas_l = [chi for l in relevant_ls ]
mus_l = [l+0.5 for l in relevant_ls ]
zs_l = [z_arg for l in relevant_ls ]
#wh_qout = np.array(Whitw_deriv_over_Whitw( kappas_l,mus_l,zs_l,prec))
wh_qout = np.array(Whitw_deriv_over_Whitw_par( kappas_l,mus_l,zs_l,prec))
DtN_atmo = 1/a + 2*1j*k_omega*wh_qout
dtn_nr = [DtN_atmo[0]]
return dtn_nr
if atmospheric_model == "Atmo":
l_idx = 200
from GRPF import Poles_and_Roots
from cmath import sqrt as csqrt
def fun(z):
DtN_list = get_atmo_dtn_sigma_squared(z)
return DtN_list[0]
param = {} # for searching poles
param["Tol"] = 1e-9
param["visual"] = 1
#param["ItMax"] = 25
param["ItMax"] = 15
param["NodesMax"] = 500000
#param["xrange"] = [305983634275718.6, 1317450785328390.8]
param["xrange"] = [5.10e14, 5.15e14]
#param["yrange"] = [ -5e14 , 5e14]
#param["yrange"] = [ -1e14 , 1e1]
param["yrange"] = [ -6.0e11 , 1e3]
param["h"] = 9e10
#param["h"] = 5e12
result = Poles_and_Roots(fun,param)
poles = 10**3*np.sqrt(result["poles"])/ (2*pi*RSun)
print("poles = ", poles)
fname_analytic_poles = "poles_freq_dtn_Atmo_l_{0}.dat".format(l_idx)
header_str_analytic_poles = "poles_real poles_imag"
plot_collect_analytic_poles = [poles.real,poles.imag]
np.savetxt(fname= fname_analytic_poles ,
X= np.transpose(plot_collect_analytic_poles ) ,
header= header_str_analytic_poles ,
comments='')
def get_atmo_nonlocal(dtn_nr):
for idx_f,f_hz in enumerate(f_hzs):
print("idx = {0}, f_hz = {1}".format(idx_f,f_hz))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment