"""
Part of comet/pyhed/hed
Earth class for calculation of dipole (HED) fields for 1d layered earth.
The algorithms in method *calcFieldForLayer* of HED class is partly taken
from Kerry Key Dipole1D.f90 after the algorithms published in [Key2009G]_
(Appendix A).
Hankel factors of Hankelfc are based on the original values of Anderson (1990).
References:
-----------
.. [Key2009G] Key, K., 2009, 1D inversion of multicomponent,
multifrequency marine (CSEM) data: Methodology and synthetic
studies for resolving thin resistive layers: Geophysics.
"""
# COMET - COupled Magnetic resonance Electrical resistivity Tomography
# Copyright (C) 2019 Nico Skibbe
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# This program 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 General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
import numpy as np
import pygimli as pg
from comet import pyhed as ph
[docs]class World1D:
def __init__(self, rho=1000., thk=None, airspace_interface=0.0,
f=2000.):
""" This class contains and manages the 1d resistivity distribution.
z is defined positive upwards!
"""
# constants
self.e0 = 8.8541878176e-12
# initialising class attributes for overview
self.rho = None
self.thk = None
self.top_layer = None
self.air_inter = airspace_interface
self.freq = None
self.omega = None
self.setFrequency(f)
self.rho = None
self.sigma_complex = None
self.air_res = None
self.setRes(rho=rho, thk=thk, air_resistivty=1e12)
# dormant for now
self.tetm = 'tetm'
[docs] def setFrequency(self, freq):
""" Simple setter for frequency + implicit omega/sigma calculation """
old_omega = self.omega
self.freq = freq
self.omega = freq * 2. * np.pi
if old_omega is not None:
self._calcComplexSigma()
[docs] def setRes(self, rho=1000., thk=None, air_resistivty=1e13):
""" Sets the resisitvity model for the dipoles + calc sigma complex.
Parameters
----------
rho: float or array_like
Resistivity distribution in Ohm*m. Airspace is considered to have
0 Ohm*m. The first entry of rho correspond to the first layer
of the subsurface. The airspace interface is considered to be at
z = 0 m which simplyfies the calculations. For offsets in z,
a cordinate transformation has to be performed externally.
thk: float or array_like
Layer thicknesses of each subsurface layer except the last.
"""
self.air_res = air_resistivty
self.rho = np.append(air_resistivty, rho).astype(float)
self._calcComplexSigma()
if thk is None:
self.thk = np.empty(0)
else:
self.thk = np.atleast_1d(thk).astype(float)
self.__checkRhoThk()
# top level of layers
self.top_layer = np.cumsum(np.append([0.0], -self.thk)) +\
self.air_inter
def _calcComplexSigma(self):
self.sigma_complex = (1./self.rho) - 1j * self.omega * self.e0
[docs] def evalSrcIdx(self, src_depth):
""" Evaluates in which layer the source is considered. """
# src layer index (-self.toplayer due to z convention)
src_idx = np.searchsorted(-self.top_layer, -src_depth,
side='left')
# in case of depth == layer interface: source is considered to be in
# the top medium
return src_idx
def __checkRhoThk(self):
""" Check function: Assert num thicknesses = num res - 1.
Internal function, do not call from outside.
"""
len1 = len(self.rho)
len2 = len(self.thk)
if len1 != len2 + 2: # -1 for both air and ground halfspace
raise Exception('Found {} thicknesses for {} layer, expect {}.'
.format(len2, len1, len1 - 2))
[docs]class HED:
def __init__(self, src_z=-0.01, src_theta=0.0, src_ids=1.0, config=None,
timer=None, debug=False):
"""
This class is needed to calculate the electric and magnetic fields of
horizontal electric dipole sources.
It contains various internally used functions for the calcuation of
the recursive factors a, b, c, and d as stated in [Key2009G]_
(Appendix A).
References:
-----------
.. [Key2009G] Key, K., 2009, 1D inversion of multicomponent,
multifrequency marine (CSEM) data: Methodology and synthetic
studies for resolving thin resistive layers: Geophysics.
"""
self.debug = debug
if config:
self.config = config
else:
self.config = World1D()
self.src_z = src_z
ph.log.debug(f'src_z: {self.src_z}')
self.src_theta = src_theta
self.src_ids = src_ids
self.cartesian = None
self.cylindrical = None
self.e_field = None
self.b_field = None
self.hfc = None
self.lam = None
self.lam_c = None
self.lam_c2 = None
self.hfc = None
self.r_2d = None
self.lay_idx = None
self.src_idx = None
self.mu0 = 4e-7 * np.pi
# recursive parameters
# hed upward attenuation coefficients: defined at base of each layer
self.hed_a = None # eq. A-10
self.hed_c = None # eq. A-12
# hed downward attenuation coefficients: defined at top of each layer
self.hed_b = None # eq. A-11
self.hed_d = None # eq. A-13
if timer is not None:
self.timer = timer
else:
self.timer = ph.misc.NoneTimer()
[docs] def setTheta(self, theta):
self.src_theta = theta
[docs] def setCoords(self, cartesian, nodes=True, drop_tol=1e-2):
""" Sets coordinates of the receiver for calculation.
All calcualtions will be performed in cylindrically coordinates.
Parameters
----------
cartesian: np.ndarray
Cartesian coordinates (N points) of the receiver points with
shape (3, N). Z is defined positive upwards.
drop_tol: float
Tolerance in meter, where the horizontal src distance is capped
to ensure a safe division (singularity fix). Distances smaller
than drop_tol are distributed between droptol and 20% of the
first value outside the droptol. Raises Exception if all points
within drop_tol. Default value of 1cm.
"""
if isinstance(cartesian, str):
if cartesian.endswith('.npy'):
cartesian = np.load(cartesian)
elif cartesian.endswith(('.bms', '.h5')):
cartesian = pg.Mesh(cartesian)
else:
raise Exception('File format not implemented: {}'
.format(cartesian))
if isinstance(cartesian, pg.Mesh):
if nodes:
self.cartesian = cartesian.positions().array().T
else:
self.cartesian = cartesian.cellCenter().array().T
elif isinstance(cartesian, np.ndarray):
self.cartesian = np.array(cartesian, dtype=np.float64)
else:
raise Exception('Data type not understood: {}'
.format(type(cartesian)))
self.cylindrical = ph.misc.Ca2Cy(self.cartesian, drop_tol=drop_tol,
dipole_z=self.src_z)
# switch to y-directed dipole (np.pi/2.) + src direction (src_theta)
self.cylindrical[1] -= self.src_theta + np.pi/2.
if self.debug:
print('theta rx', np.rad2deg(self.cylindrical[1]))
[docs] def calculate(self):
"""Calculates the 1d layered earth recursive formula.
Calculates the recursive attenuation and reflection coefficients for
each layer on basis of the given set of cylindrical coordinates.
Fills the variables **R_p**, **R_m**, **r_p**, **r_s**, **S_p**,
**S_m**, **s_p**, **s_m**, **hem_a**, **hem_b**, **hem_c**, and
**hem_d**. The used formulas correspond to equations A-6 to A-13 in
[Key2009G]_ (Appendix A).
"""
if self.cylindrical is None:
raise Exception('No coords found. Use "setCoords" to define '
'reciever points.')
self.timer.update()
# %% initialize Hankel coefficients
# hold instance of Hankel class for calculation of lambda
# get Hankel factors and reshape for convolution
self.hfc = wer_201_2018()
self.hfc_j0 = self.hfc.j0[:, np.newaxis]
self.hfc_j1 = self.hfc.j1[:, np.newaxis]
# else:
# Christensen 101 filter -> turn around for convolution
# self.hfc = hankelfc()
# self.hfc_j0 = self.hfc.getFactors('j0')[::-1]
# self.hfc_j1 = self.hfc.getFactors('j1')[::-1]
# get indices and rho indices for each receiver and layer
self._getIndices()
# %% recursive part
self.e_field = np.zeros(self.cylindrical.shape, dtype=np.complex)
self.b_field = np.zeros(self.cylindrical.shape, dtype=np.complex)
for rx_lay, lay_indices in enumerate(self.lay_idx):
# overall loop: separate calculations for the points of each layer
# as there are no possibilities to cache, as the basic wavenumbers
# are different. The only time to gain here is due to more
# sophisticated broadcasting, but this comes with a lot of
# unnesessary empty spaces in some of the recursive arrays which
# cost a lot of memory and calulation time as well.
if self.debug:
print('#########################')
if rx_lay != self.src_idx:
print('Layer: {} ({} points)'
.format(rx_lay, self.lay_num[rx_lay]))
else:
print('Layer: {} (source) ({} points)'
.format(rx_lay, self.lay_num[rx_lay]))
print('#########################')
if self.lay_num[rx_lay] == 0:
# no receiver in layer num, no calculations.
continue
# what we need:
# Step 1/4: wavenumbers
# shape = (hankel, model)
# if ph.params.use_kk_hankel:
# Werthmüller 2018
lam = self.hfc.base[:, np.newaxis] /\
self.cylindrical[0, lay_indices]
# else:
# lam = self.hfc.getWavenumbes('j0')[:, np.newaxis] / \
# self.cylindrical[0, lay_indices]
# we need both lam and lam², so this is about memory vs time
lam_2 = lam ** 2
# complex wavenumber lam_c² = lam² - 1j * (2 * pi * mu0 * sigma)
# shape = (hankel, model, layer)
kk = - 1j * (self.config.omega * self.mu0) * \
self.config.sigma_complex
lam_c2 = lam_2[:, :, np.newaxis] + kk
lam_c = np.sqrt(lam_c2)
exp = np.zeros_like(lam_c)
exp[:, :, 1:-1] = np.exp(-lam_c[:, :, 1:-1] * self.thk[1:-1])
# Step 2/4: R+, R-, S+, S- for all layer
R_p, R_m, S_p, S_m = self.reflectionCoefficients(
rx_lay, lam, lam_2, lam_c, lam_c2, exp)
self.timer.tick('reflectionCoefficients: {:2.3f} sec')
# Step 3/4: Ay, Az (a, b, c, d) can be calulated recursively
e_field = self.calcFieldForLayer(rx_lay, lay_indices,
lam, lam_2, lam_c, lam_c2,
R_p, R_m, S_p, S_m, exp)
# Step 4/4: Append Field of layer rx_lay
self.e_field[:, lay_indices] = e_field
# Note: currently a bug occurs below the airspace where the z
# component is sign-flipped. This is fixed here explicitely,
# however, I need to find the bug in the real implementation.
# TODO!
ground = self.cylindrical[2] < self.config.air_inter
self.e_field[2, ground] *= -1.
# e-iwt conformity with custEM: added 19.11.2019
# switch real component of field
self.e_field = -self.e_field.real + 1j * self.e_field.imag
# scaling by dipole moment
self.e_field *= self.src_ids
# rotate back by theta
if not np.isclose(self.src_theta, 0):
ph.misc.rotate3(self.e_field, self.src_theta, axis='z', copy=False)
def _getIndices(self):
"""
Evaluates src-rec indices and receiver sorting for recursion
formulation.
"""
# how many layer above and below transmitter
self.src_idx = self.config.evalSrcIdx(self.src_z)
# z layer index for each point
self.lay_idx = ()
for ind in range(len(self.config.sigma_complex)):
if ind == 0:
# case 1/3: airspace (upper halfspace)
self.lay_idx += np.where(self.cylindrical[2] >=
self.config.air_inter)
elif ind < len(self.config.sigma_complex) - 1:
# case 2/3: all layer except last (N finite layers)
self.lay_idx += np.where(
np.logical_and(
self.cylindrical[2] < self.config.top_layer[ind - 1],
self.cylindrical[2] >= self.config.top_layer[ind]))
else:
# case 3/3: homogenous earth (lower halfspace)
self.lay_idx += np.where(
self.cylindrical[2] < self.config.top_layer[-1])
# number of receivers in each layer (can be zero)
self.lay_num = []
for arr in self.lay_idx:
self.lay_num.append(len(arr))
# number of receiver above source
self.above_num = int(np.sum(self.lay_num[:self.src_idx]))
# layer indices above source layer from top to bottom
self.above_idx = np.arange(0, self.src_idx, dtype=int)
# number of receiver below source
self.below_num = int(np.sum(self.lay_num[self.src_idx + 1:]))
# layer indices below source layer from bottom to top
self.below_idx = np.arange(self.src_idx + 1, len(self.lay_idx),
dtype=int)[::-1]
# number of receiver in source layer
self.src_num = int(self.lay_num[self.src_idx])
# layer thicknesses (1 value per layer for easy indexing)
self.thk = [np.inf]
self.thk.extend(self.config.thk)
self.thk.append(np.inf)
# layer top (1 value per layer for easy indexing)
self.z_top = [np.inf]
self.z_top.extend(self.config.top_layer)
[docs] def reflectionCoefficients(self, rx_lay, lam, lam_2, lam_c, lam_c2, exp):
""" Calculation of the general reflection coefficients R+, R-, S+,
and S- as stated in [Key2009G]_ (Appendix A, equations A-06 to A-09).
Computed from the air and halfspace, respectively, inward to the
source layer.
"""
# general rule: (a - b)/(a + b) = 2a / (a + b) - 1
# -->
# r- = (lamc_i - lamc_i-1)/(lamc_i + lamc_i-1)
# = 2*lamc_i/(lamc_i + lamc_i-1) - 1
# print_i = 99 # dipole1d - 1 number/index ->
# entspricht print - 1 beim plotten
R_p = [] # R+
S_p = [] # S+
R_m = [] # R-
S_m = [] # S-
# layer indices for downward recursion
top_2_src = np.append(self.above_idx, self.src_idx)
# layer indices for upward recursion
btm_2_src = np.append(self.below_idx, self.src_idx)
# common term used for calculation of r+ and r-
# (lamc[:-1] + lamc[1:])
common_r = lam_c[:, :, :-1] + lam_c[:, :, 1:]
# common term used for calculation of s+ and s-
# lamc[:-1] * sig[1:] + lamc[1:] * sig[:-1]
common_s = lam_c[:, :, :-1] * self.config.sigma_complex[1:] + \
lam_c[:, :, 1:] * self.config.sigma_complex[:-1]
# %% upward
for itr, lay in enumerate(btm_2_src):
# begin recursive calulation: from bottom to source (n -> src)
if itr == 0:
# R+ [N] = 0
R_p.append(np.zeros(lam.shape, dtype=complex))
# S+ [N] = 0
S_p.append(np.zeros(lam.shape, dtype=complex))
continue
# equation A-7: r+
# lam- = lamc(i) - lamc(i-1)
# lam+ = lamc(i) + lamc(i-1)
# r+: lam- / lam+
# == lam- * lam+ / (lam+)**2
# with 3rd binomial formula:
# lam- * lam+ = -i * omega * mu0 * (sig(i) - sig(i - 1))
# -> r+ = -i * omega * mu0 * (sig(i) - sig(i - 1)) / (lam+)**2
# this is after the code from Dipole1D.f90 and increases accuracy
# when it comes to very small wavenumbers
r_p = -1j * self.config.omega * self.mu0 * \
(self.config.sigma_complex[lay] -
self.config.sigma_complex[lay+1]) /\
common_r[:, :, lay]**2
# equation A-9: s+
s_p = (2 * lam_c[:, :, lay] *
self.config.sigma_complex[lay + 1] /
common_s[:, :, lay]) - 1
if itr == 1:
# R+[N] = 0 --> R+[N-1] = r+[n-1]
R_p.append(r_p)
# same for S+
S_p.append(s_p)
else:
# exp used in next recursive iteration
# exp term from last time and this time eq.
# exp(-lamc[i + 1] * h[i + 1])
exp_j1 = exp[:, :, lay + 1]
# equation A-6: R+
# (r+[i] + R+[i+1] * exp(-lamc[i+1] * h[i+1]) /
# 1 + r+[i] * R+[i+1] * exp(-lamc[i+1] * h[i+1]))
# * exp(-lamc[i] * h[i]) ! applied in next step
R_p[-1] = R_p[-1] * exp_j1
R_term = R_p[-1] * exp_j1
R_p.append((r_p + R_term) / (1 + r_p * R_term))
# equation A-8: S+
# (s+[i] + S+[i+1] * exp(-lamc[i+1] * h[i+1]) /
# 1 + s+[i] * S+[i+1] * exp(-lamc[i+1] * h[i+1]))
# * exp(-lamc[i] * h[i]) ! applied in next step
S_p[-1] = S_p[-1] * exp_j1
S_term = S_p[-1] * exp_j1
S_p.append((s_p + S_term) / (1 + s_p * S_term))
for _ in self.above_idx:
# R+, S+ only defined below source
# append empty array so that len(R/S) == number of layers
# and entries are at the correct indices
R_p.append(np.array([]))
S_p.append(np.array([]))
# N...0 -> 0...N, for indexing later on
S_p = S_p[::-1]
R_p = R_p[::-1]
# %% downward
for itr, lay in enumerate(top_2_src):
# begin recursive calulation: from top to source (0 -> src)
if itr == 0:
# R- [0] = 0
R_m.append(np.zeros(lam.shape, dtype=complex))
# S- [0] = 0
S_m.append(np.zeros(lam.shape, dtype=complex))
continue
# equation A-7: r-
# r-[i] see r+ for derivation (and account for index shift)
r_m = -1j * self.config.omega * self.mu0 * \
(self.config.sigma_complex[lay] -
self.config.sigma_complex[lay-1]) /\
common_r[:, :, lay - 1]**2
# equation A-9: s-
s_m = (2 * lam_c[:, :, lay] *
self.config.sigma_complex[lay - 1] /
common_s[:, :, lay - 1]) - 1
if itr == 1:
# R- [0] = 0 --> R- [1] = r- [1], same for S-
R_m.append(r_m)
S_m.append(s_m)
else:
# exp term from last time and this time eq.
# exp(-lamc[i - 1] * h[i - 1])
exp_j1 = exp[:, :, lay - 1]
# equation A-6: R-
# (r-[i] + R-[i-1] * exp(-lamc[i-1] * h[i-1]) /
# 1 + r-[i] * R-[i-1] * exp(-lamc[i-1] * h[i-1]))
# * exp(-lamc[i] * h[i]) ! applied in next step
R_m[-1] = R_m[-1] * exp_j1
R_term = R_m[-1] * exp_j1
R_m.append((r_m + R_term) / (1 + r_m * R_term))
# equation A-8: S-
# (s-[i] + S-[i-1] * exp(-lamc[i-1] * h[i-1]) /
# 1 + s-[i] * S-[i-1] * exp(-lamc[i-1] * h[i-1]))
# * exp(-lamc[i] * h[i]) ! applied in next step
S_m[-1] = S_m[-1] * exp_j1
S_term = S_m[-1] * exp_j1
S_m.append((s_m + S_term) / (1 + s_m * S_term))
for _ in self.below_idx:
# R-, S- only defined above source
# append empty array so that len(R/S) == number of layers
# and entries are at the correct indices
R_m.append(np.array([]))
S_m.append(np.array([]))
return R_p, R_m, S_p, S_m
[docs] def calcFieldForLayer(self, rx_lay, lay_indices, lam, lam_2, lam_c, lam_c2,
R_p, R_m, S_p, S_m, exp):
print_i = 99
# %% Evaluation of the attenuation coefficients a, b, c, d
self.timer.update()
if rx_lay == self.src_idx:
src_2_rx = [self.src_idx]
elif rx_lay > self.src_idx:
if self.debug:
print('target layer below source')
down = True
src_2_rx = np.arange(self.src_idx, rx_lay + 1)
else:
if self.debug:
print('target layer above source')
down = False
src_2_rx = np.arange(self.src_idx, rx_lay - 1, -1)
if self.debug:
print('calculate Potential recursively for layers: {}'
.format(src_2_rx))
for itr, lay in enumerate(src_2_rx):
# starting with the source layer i = j
if itr == 0:
# mu0 / (2 * lam_c[j])
mu_2_lamc = self.mu0 / (2. * lam_c[:, :, lay])
# mu0 / (2 * lam²)
mu_2_lam2 = self.mu0 / (2. * lam_2)
if lay == 0:
# case 1/3: src in airspace
if self.debug:
print('case 1/3: src in airspace')
# R- = S- = 0
# --> R+ * R- = 0
exp_btm = np.exp(-lam_c[:, :, lay] *
np.abs(self.z_top[lay + 1] - self.src_z))
coeff_a = R_p[lay] * mu_2_lamc * exp_btm
coeff_b = np.zeros(lam.shape, dtype=np.complex)
coeff_c = - S_p[lay] * mu_2_lam2 * exp_btm
coeff_d = np.zeros(lam.shape, dtype=np.complex)
elif lay == len(self.lay_num) - 1:
# case 2/3: src in lower halfspace
if self.debug:
print('case 2/3: src in lower halfspace')
# R+ = S+ = 0
# --> R+ * R- = 0
exp_top = np.exp(-lam_c[:, :, lay] *
np.abs(self.z_top[lay] - self.src_z))
coeff_a = np.zeros(lam.shape, dtype=np.complex)
coeff_b = R_m[lay] * mu_2_lamc * exp_top
coeff_c = np.zeros(lam.shape, dtype=np.complex)
coeff_d = S_m[lay] * mu_2_lam2 * exp_top
else:
# case 3/3: src in finite layer
if self.debug:
print('case 3/3: src in finite layer')
# lay == self.src_idx
exp_btm = np.exp(-lam_c[:, :, lay] *
np.abs(self.z_top[lay + 1] - self.src_z))
exp_top = np.exp(-lam_c[:, :, lay] *
np.abs(self.z_top[lay] - self.src_z))
# exponential term for R+, R-, S-, S- in source layer
exp_j = exp[:, :, lay]
Rp_Rm = 1. - R_p[lay] * R_m[lay] *\
exp_j * exp_j
Sp_Sm = 1. - S_p[lay] * S_m[lay] *\
exp_j * exp_j
# a (eq. A-10)
coeff_a = (exp_btm + R_m[lay] * exp_j * exp_top) * \
(R_p[lay] / Rp_Rm) * mu_2_lamc
# b (eq. A-11)
coeff_b = (R_p[lay] * exp_j * exp_btm + exp_top) * \
(R_m[lay] / Rp_Rm) * mu_2_lamc
# c (eq. A-12)
coeff_c = (-exp_btm + S_m[lay] * exp_j * exp_top) * \
(S_p[lay] / Sp_Sm) * mu_2_lam2
# d (eq. A-13)
coeff_d = (-S_p[lay] * exp_j * exp_btm + exp_top) * \
(S_m[lay] / Sp_Sm) * mu_2_lam2
else:
# source ready, now moving to rx layer
if down:
# moving downwards to rx layer
top_te = coeff_a + coeff_b * exp[:, :, lay - 1]
top_tm = coeff_c + coeff_d * exp[:, :, lay - 1]
if itr == 1:
# src term
src_ab = exp_btm * mu_2_lamc
src_cd = -exp_btm * mu_2_lam2
top_te += src_ab
top_tm += src_cd
coeff_b = top_te / (1. + R_p[lay] * exp[:, :, lay])
coeff_a = coeff_b * R_p[lay]
coeff_d = top_tm / (1. + S_p[lay] * exp[:, :, lay])
coeff_c = coeff_d * S_p[lay]
else:
# print('up to layer idx:', lay)
# moving upwards to rx layer
top_te = coeff_a * exp[:, :, lay + 1] + coeff_b
top_tm = coeff_c * exp[:, :, lay + 1] + coeff_d
if itr == 1:
src_ab = exp_top * mu_2_lamc
src_cd = exp_top * mu_2_lam2
top_te += src_ab
top_tm += src_cd
coeff_a = top_te / (1. + R_m[lay] * exp[:, :, lay])
coeff_b = coeff_a * R_m[lay]
coeff_c = top_tm / (1. + S_m[lay] * exp[:, :, lay])
coeff_d = coeff_c * S_m[lay]
self.timer.tick('coeff a, b, c, d: {:2.3f} sec')
if self.debug:
print('final coeff')
print('coeff a[{}]'.format(print_i + 1), coeff_a[print_i])
print('coeff b[{}]'.format(print_i + 1), coeff_b[print_i])
print('coeff c[{}]'.format(print_i + 1), coeff_c[print_i])
print('coeff d[{}]'.format(print_i + 1), coeff_d[print_i])
# %% field
# A = (0, A_y, d/dy(A_z)).T
# A_y: equation (A-2)
# d/dy(A_z): equation (A-3)
# E = iw * A + 1 / (mu * sigma) * grad(div(A))
# E_x = ( d2/dxdy(A_y) + d3/dxdydz(A_z) ) / (mu * sigma)
# E_y = iw * A_y + ( d2/dy²(A_y) + d3/dy²dz(A_z)) / (mu * sigma)
# E_z = iw * A_z + ( d2/dydz(A_y) + d3/dydz²(A_z) ) / (mu * sigma)
# B = rot(A)
# B_x = d2/dy²(A_z) - d/dz(A_y)
# B_y = - d2/dxdy(A_z)
# B_z = d/dx(A_y)
# cylindrical coords for rx layer: (r, phi, z)
local = self.cylindrical[:, lay_indices]
if self.debug:
print('--------len local {} ------------'.format(local.shape))
# calculation common terms
# exp(-lam (z - z_i) )
if rx_lay == 0:
# hack so that exp_b and exp_d are 0, which is correct
exp_i = np.zeros_like(lam_c[:, :, rx_lay])
else:
exp_i = np.exp(
-lam_c[:, :, rx_lay] * (-local[2] + self.z_top[rx_lay]))
# exp(lam (z - z_i+1) )
if rx_lay == len(self.lay_num) - 1:
# hack so that exp_a and exp_c are 0, which is correct
exp_i1 = np.zeros_like(lam_c[:, :, rx_lay])
else:
exp_i1 = np.exp(
lam_c[:, :, rx_lay] * (-local[2] + self.z_top[rx_lay + 1]))
exp_a = coeff_a * exp_i1
exp_b = coeff_b * exp_i
exp_c = coeff_c * exp_i1
exp_d = coeff_d * exp_i
if self.debug:
print('exp_a[{}]'.format(print_i + 1), exp_a[print_i])
print('exp_b[{}]'.format(print_i + 1), exp_b[print_i])
print('exp_c[{}]'.format(print_i + 1), exp_c[print_i])
print('exp_d[{}]'.format(print_i + 1), exp_d[print_i])
# building z derivatives before Hankel integration is applied
ab_plus = exp_a + exp_b
ab_minus = exp_a - exp_b
if self.debug:
print('exp_a + exp_b [{}]:'.format(print_i + 1), ab_plus[print_i])
print('exp_a - exp_b [{}]:'.format(print_i + 1), ab_minus[print_i])
if self.src_idx == rx_lay:
# only if source and receiver in the same layer
# explicit solves the Kronecker in (eq. A-2)
# and keeps calculation effort to minimum
src_exp = np.exp(-lam_c[:, :, self.src_idx] *
np.abs(local[2] - self.src_z))
src_term = self.mu0/(2.*lam_c[:, :, rx_lay]) * src_exp
if self.debug:
print('src_term(100)', src_term[print_i])
ay = ab_plus + src_term
# derivative of source_term is defined for z != src_z
# d/dz(exp(-lam_c * |z - src_z|)) = -lam_c * signum(z-src_z) *\
# exp(-lam_c * |z - src_z|)
# note the effect of the z convention inside the sign function
ddz_src_term = -self.mu0 / 2. * np.sign(-local[2] + self.src_z) *\
src_exp
if self.debug:
print('ddz_src_term(100)', ddz_src_term[print_i])
ddz_ay = ab_minus * lam_c[:, :, rx_lay] + ddz_src_term
else:
# A_y (inside Hankel) (eq. A-2)
ay = ab_plus
ddz_ay = ab_minus * lam_c[:, :, rx_lay]
if self.debug:
print('ay [{}]:'.format(print_i + 1), ay[print_i])
print('ddz_ay [{}]:'.format(print_i + 1), ddz_ay[print_i])
# A_z (inside Hankel) (eq. A-3)
az = (exp_c + exp_d) - lam_c[:, :, rx_lay] * ab_minus / lam_2
# d/dz(A_z)
ddz_az = lam_c[:, :, rx_lay] * (exp_c - exp_d) - \
ab_plus * lam_c2[:, :, rx_lay] / lam_2
if self.debug:
print('##########################################################')
print('gg', lam_c[print_i, :, rx_lay])
print('ce1', exp_c[print_i])
print('de2', exp_d[print_i])
print('ae1pbe2', ab_plus[print_i])
print('gg2', lam_c2[print_i, :, rx_lay])
print('lam2', lam_2[print_i])
print('ddz_az', ddz_az[print_i])
print('##########################################################')
print('az [{}]:'.format(print_i + 1), az[print_i])
print('ddz_az [{}]:'.format(print_i + 1), ddz_az[print_i])
# d²/dz²(A_z)
d2dz2_az = lam_c2[:, :, rx_lay] * az
if self.debug:
print('d2dz2_az [{}]:'.format(print_i + 1), d2dz2_az[print_i])
mu_sig_2pi = 1. /\
(self.mu0 * self.config.sigma_complex[rx_lay] * 2. * np.pi)
lam_3 = lam_2 * lam
if self.debug:
print('csigsite', self.config.sigma_complex[rx_lay])
# ay_ddz_az = (ay + ddz_az)
# %% Hankel Integration about here
# Integral(0...inf) ()
ex_j0 = np.sum(-lam_3 * (ay + ddz_az) * self.hfc_j0, 0, np.complex)
if self.debug:
print('ex_j0 [{}]:'.format(print_i + 1),
(-lam_3 * (ay + ddz_az))[print_i])
# Integral(0...inf) ()
ex_j1 = np.sum(2. * lam_2 / local[0] * (ay + ddz_az) * self.hfc_j1,
0, np.complex)
if self.debug:
print('ex_j1 [{}]:'.format(print_i + 1), (2. * lam_2 / local[0] *
(ay + ddz_az))[print_i])
# Integral(0...inf) ()
# 0.42993358039234764 (-np.pi off)
part1 = 1j * self.config.freq * ay * lam
# part2 = -mu_sig_2pi * (ay + ddz_az) * lam_3 * \
# np.sin(local[1] + np.pi/2.)**2
part2 = -mu_sig_2pi * (ay + ddz_az) * lam_3 * \
np.sin(local[1])**2
int_ey_j0 = part1 + part2
if self.debug:
print('ey_j0 [{}]:'.format(print_i + 1), int_ey_j0[print_i])
ey_j0 = np.sum(int_ey_j0 * self.hfc_j0, 0, np.complex)
# Integral(0...inf) ()
# int_ey_j1 = -mu_sig_2pi * (ay + ddz_az) * lam_2 *\
# np.cos(2.*(local[1] + np.pi/2.)) / local[0]
int_ey_j1 = -mu_sig_2pi * (ay + ddz_az) * lam_2 *\
np.cos(2.*(local[1])) / local[0]
ey_j1 = np.sum(int_ey_j1 * self.hfc_j1, 0, np.complex)
if self.debug:
print('ey_j1 [{}]:'.format(print_i + 1), int_ey_j1[print_i])
int_ez_j1 = (1j * self.config.omega * az *
self.config.sigma_complex[rx_lay] +
(ddz_ay + d2dz2_az) / self.mu0) * lam_2
if self.debug:
print('(ddz_ay + d2dz2_az) [{}]:'.format(print_i + 1),
(ddz_ay + d2dz2_az)[print_i])
print('int_ez_j1 [{}]:'.format(print_i + 1), int_ez_j1[print_i])
print('hfc_j0 [{}]:'.format(print_i + 1), '{:.12e}'
.format(self.hfc_j0[print_i, 0]))
print('hfc_j1 [{}]:'.format(print_i + 1), '{:.12e}'
.format(self.hfc_j1[print_i, 0]))
ez_j1 = np.sum(int_ez_j1 * self.hfc_j1,
0, np.complex)
if self.debug:
print('ez_j1:', ez_j1)
print('exh:', ex_j0 + ex_j1)
print('eyh:', ey_j0 + ey_j1)
# Hankel ready, except for a division by horizontal distance to tx
# return fields
e_field = np.zeros((3, self.lay_num[rx_lay]), dtype=np.complex)
# attention: e_field[0] = e_field[1] and e_field[1] = -e_field[0]
# due to convention of x-directed src
# this simply replaces the rotation of the field
# e_field[1] = np.sin(2. * (local[1] + np.pi/2.)) / (2. / mu_sig_2pi) *\
# (ex_j0 + ex_j1) / local[0]
e_field[1] = np.sin(2. * (local[1])) / (2. / mu_sig_2pi) *\
(ex_j0 + ex_j1) / local[0]
e_field[0] = -(ey_j0 + ey_j1) / local[0]
# e_field[2] = -np.sin(local[1] + np.pi/2.) /\
# (2. * np.pi) * ez_j1 / local[0]
e_field[2] = -np.sin(local[1]) /\
(2. * np.pi) * ez_j1 / local[0]
if self.debug:
print('jz:', e_field[2])
# jz1D(i) = jz1D(i)/ (sigsite(i) + II*eps*2*pi*ftx1D)
e_field[2] /= np.conjugate(self.config.sigma_complex[rx_lay])
if self.debug:
print('II*eps*2*pi*ftx1D', 1j * self.config.e0 * 2. *
np.pi * self.config.freq)
print('sigsite(i)',
np.conjugate(self.config.sigma_complex[rx_lay]))
print('ex:', e_field[0])
print('ey:', e_field[1])
print('ez:', e_field[2])
# TODO implement B fields as well
# B = rot(A)
# B_x = d2/dy²(A_z) - d/dz(A_y)
# B_y = - d2/dxdy(A_z)
# B_z = d/dx(A_y)
# b_field = np.zeros((3, self.lay_num[rx_lay]), dtype=np.complex)
self.timer.tick('fields: {:2.3f} sec')
return e_field
[docs]class hankelfc:
def __init__(self):
# filter coefficients for hankel transformation
# coefficients by Anderson (1980)
self.factors = {}
self.factors['sin'] = [np.array([
2.59526236e-7, 3.66544843e-7, 5.17830795e-7, 7.31340622e-7,
1.03322805e-6, 1.45918500e-6, 2.06161065e-6, 2.91137793e-6,
4.11357863e-6, 5.80876420e-6, 8.20798075e-6, 1.15895083e-5,
1.63778560e-5, 2.31228459e-5, 3.26800649e-5, 4.61329334e-5,
6.52101085e-5, 9.20390575e-5, 1.30122935e-4, 1.83620431e-4,
2.59656626e-4, 3.66311982e-4, 5.18141184e-4, 7.30717340e-4,
1.03392184e-3, 1.45742714e-3, 2.06292302e-3, 2.90599911e-3,
4.11471902e-3, 5.79042763e-3, 8.20004722e-3, 1.15192930e-2,
1.63039133e-2, 2.28257757e-2, 3.22249222e-2, 4.47864328e-2,
6.27329625e-2, 8.57059100e-2, 1.17418314e-1, 1.53632655e-1,
1.97717964e-1, 2.28849849e-1, 2.40311038e-1, 1.65409220e-1,
2.84701476e-3, -2.88016057e-1, -3.69097406e-1, -2.50107514e-2,
5.71811256e-1, -3.92261572e-1, 7.63280044e-2, 5.16233994e-2,
-6.48012082e-2, 4.89047141e-2, -3.26936331e-2, 2.10539842e-2,
-1.33862549e-2, 8.47124695e-3, -5.35123972e-3, 3.37796651e-3,
-2.13174466e-3, 1.34513833e-3, -8.48749612e-4, 5.35531006e-4,
-3.37898780e-4, 2.13200109e-4, -1.34520273e-4, 8.48765787e-5,
-5.35535069e-5, 3.37899801e-5, -2.13200365e-5, 1.34520337e-5,
-8.48765949e-6, 5.35535110e-6, -3.37899811e-6, 2.13200368e-6,
-1.34520338e-6, 8.48765951e-7, -5.35535110e-7, 3.37899811e-7],
np.float), 80, 40]
self.factors['cos'] = [np.array([
1.63740363e-7, 1.83719709e-7, 2.06136904e-7, 2.31289411e-7,
2.59510987e-7, 2.91176117e-7, 3.26704977e-7, 3.66569013e-7,
4.11297197e-7, 4.61483045e-7, 5.17792493e-7, 5.80972733e-7,
6.51862128e-7, 7.31401337e-7, 8.20645798e-7, 9.20779729e-7,
1.03313185e-6, 1.15919300e-6, 1.30063594e-6, 1.45933752e-6,
1.63740363e-6, 1.83719709e-6, 2.06136904e-6, 2.31289411e-6,
2.59510987e-6, 2.91176117e-6, 3.26704977e-6, 3.66569013e-6,
4.11297197e-6, 4.61483045e-6, 5.17792493e-6, 5.80972733e-6,
6.51862128e-6, 7.31401337e-6, 8.20645798e-6, 9.20779729e-6,
1.03313185e-5, 1.15919300e-5, 1.30063594e-5, 1.45933752e-5,
1.63740363e-5, 1.83719709e-5, 2.06136904e-5, 2.31289411e-5,
2.59510987e-5, 2.91176117e-5, 3.26704977e-5, 3.66569013e-5,
4.11297197e-5, 4.61483045e-5, 5.17792493e-5, 5.80972733e-5,
6.51862128e-5, 7.31401337e-5, 8.20645798e-5, 9.20779729e-5,
1.03313185e-4, 1.15919300e-4, 1.30063594e-4, 1.45933752e-4,
1.63740363e-4, 1.83719709e-4, 2.06136904e-4, 2.31289411e-4,
2.59510987e-4, 2.91176117e-4, 3.26704976e-4, 3.66569013e-4,
4.11297197e-4, 4.61483045e-4, 5.17792493e-4, 5.80972733e-4,
6.51862127e-4, 7.31401337e-4, 8.20645797e-4, 9.20779730e-4,
1.03313185e-3, 1.15919300e-3, 1.30063593e-3, 1.45933753e-3,
1.63740362e-3, 1.83719710e-3, 2.06136901e-3, 2.31289411e-3,
2.59510977e-3, 2.91176115e-3, 3.26704948e-3, 3.66569003e-3,
4.11297114e-3, 4.61483003e-3, 5.17792252e-3, 5.80972566e-3,
6.51861416e-3, 7.31400728e-3, 8.20643673e-3, 9.20777603e-3,
1.03312545e-2, 1.15918577e-2, 1.30061650e-2, 1.45931339e-2,
1.63734419e-2, 1.83711757e-2, 2.06118614e-2, 2.31263461e-2,
2.59454421e-2, 2.91092045e-2, 3.26529302e-2, 3.66298115e-2,
4.10749753e-2, 4.60613861e-2, 5.16081994e-2, 5.78193646e-2,
6.46507780e-2, 7.22544422e-2, 8.03873578e-2, 8.92661837e-2,
9.80670729e-2, 1.07049506e-1, 1.13757572e-1, 1.18327217e-1,
1.13965041e-1, 1.00497783e-1, 6.12958082e-2, -1.61234222e-4,
-1.11788551e-1, -2.27536948e-1, -3.39004453e-1, -2.25128800e-1,
8.98279919e-2, 5.12510388e-1, -1.31991937e-1, -3.35136479e-1,
3.64868100e-1, -2.34039961e-1, 1.32085237e-1, -7.56739672e-2,
4.52296662e-2, -2.78297002e-2, 1.73727753e-2, -1.09136894e-2,
6.87397283e-3, -4.33413470e-3, 2.73388730e-3, -1.72477355e-3,
1.08821012e-3, -6.86602007e-4, 4.33213523e-4, -2.73338487e-4,
1.72464733e-4, -1.08817842e-4, 6.86594042e-5, -4.33211523e-5,
2.73337984e-5, -1.72464607e-5, 1.08817810e-5, -6.86593962e-6,
4.33211503e-6, -2.73337979e-6, 1.72464606e-6, -1.08817810e-6,
6.86593961e-7, -4.33211503e-7, 2.73337979e-7, -1.72464606e-7],
np.float), 164, 122]
self.factors['j0'] = [np.array([
2.89878288e-7, 3.64935144e-7, 4.59426126e-7, 5.78383226e-7,
7.28141338e-7, 9.16675639e-7, 1.15402625e-6, 1.45283298e-6,
1.82900834e-6, 2.30258511e-6, 2.89878286e-6, 3.64935148e-6,
4.59426119e-6, 5.78383236e-6, 7.28141322e-6, 9.16675664e-6,
1.15402621e-5, 1.45283305e-5, 1.82900824e-5, 2.30258527e-5,
2.89878259e-5, 3.64935186e-5, 4.59426051e-5, 5.78383329e-5,
7.28141144e-5, 9.16675882e-5, 1.15402573e-4, 1.45283354e-4,
1.82900694e-4, 2.30258630e-4, 2.89877891e-4, 3.64935362e-4,
4.59424960e-4, 5.78383437e-4, 7.28137738e-4, 9.16674828e-4,
1.15401453e-3, 1.45282561e-3, 1.82896826e-3, 2.30254535e-3,
2.89863979e-3, 3.64916703e-3, 4.59373308e-3, 5.78303238e-3,
7.27941497e-3, 9.16340705e-3, 1.15325691e-2, 1.45145832e-2,
1.82601199e-2, 2.29701042e-2, 2.88702619e-2, 3.62691810e-2,
4.54794031e-2, 5.69408192e-2, 7.09873072e-2, 8.80995426e-2,
1.08223889e-1, 1.31250483e-1, 1.55055715e-1, 1.76371506e-1,
1.85627738e-1, 1.69778044e-1, 1.03405245e-1, -3.02583233e-2,
-2.27574393e-1, -3.62173217e-1, -2.05500446e-1, 3.37394873e-1,
3.17689897e-1, -5.13762160e-1, 3.09130264e-1, -1.26757592e-1,
4.61967890e-2, -1.80968674e-2, 8.35426050e-3, -4.47368304e-3,
2.61974783e-3, -1.60171357e-3, 9.97717882e-4, -6.26275815e-4,
3.94338818e-4, -2.48606354e-4, 1.56808604e-4, -9.89266288e-5,
6.24152398e-5, -3.93805393e-5, 2.48472358e-5, -1.56774945e-5,
9.89181741e-6, -6.24131160e-6, 3.93800058e-6, -2.48471018e-6,
1.56774609e-6, -9.89180896e-7, 6.24130948e-7, -3.93800005e-7,
2.48471005e-7, -1.56774605e-7, 9.89180888e-8, -6.24130946e-8],
np.float), 100, 60]
self.factors['j1'] = [np.array([
1.84909557e-13, 2.85321327e-13, 4.64471808e-13, 7.16694771e-13,
1.16670043e-12, 1.80025587e-12, 2.93061898e-12, 4.52203829e-12,
7.36138206e-12, 1.13588466e-11, 1.84909557e-11, 2.85321327e-11,
4.64471808e-11, 7.16694771e-11, 1.16670043e-10, 1.80025587e-10,
2.93061898e-10, 4.52203829e-10, 7.36138206e-10, 1.13588466e-9,
1.84909557e-9, 2.85321326e-9, 4.64471806e-9, 7.16694765e-9,
1.16670042e-8, 1.80025583e-8, 2.93061889e-8, 4.52203807e-8,
7.36138149e-8, 1.13588452e-7, 1.84909521e-7, 2.85321237e-7,
4.64471580e-7, 7.16694198e-7, 1.16669899e-6, 1.80025226e-6,
2.93060990e-6, 4.52201549e-6, 7.36132477e-6, 1.13587027e-5,
1.84905942e-5, 2.85312247e-5, 4.64449000e-5, 7.16637480e-5,
1.16655653e-4, 1.79989440e-4, 2.92971106e-4, 4.51975783e-4,
7.35565435e-4, 1.13444615e-3, 1.84548306e-3, 2.84414257e-3,
4.62194743e-3, 7.10980590e-3, 1.15236911e-2, 1.76434485e-2,
2.84076233e-2, 4.29770596e-2, 6.80332569e-2, 9.97845929e-2,
1.51070544e-1, 2.03540581e-1, 2.71235377e-1, 2.76073871e-1,
2.16691977e-1, -7.83723737e-2, -3.40675627e-1, -3.60693673e-1,
5.13024526e-1, -5.94724729e-2, -1.95117123e-1, 1.99235600e-1,
-1.38521553e-1, 8.79320859e-2, -5.50697146e-2, 3.45637848e-2,
-2.17527180e-2, 1.37100291e-2, -8.64656417e-3, 5.45462758e-3,
-3.44138864e-3, 2.17130686e-3, -1.36998628e-3, 8.64398952e-4,
-5.45397874e-4, 3.44122545e-4, -2.17126585e-4, 1.36997597e-4,
-8.64396364e-5, 5.45397224e-5, -3.44122382e-5, 2.17126544e-5,
-1.36997587e-5, 8.64396338e-6, -5.45397218e-6, 3.44122380e-6,
-2.17126543e-6, 1.36997587e-6, -8.64396337e-7, 5.45397218e-7],
np.float), 100, 60]
[docs] def getFactors(self, string):
""" Returns the requested set of Hankel factors.
Parameters
----------
string: [str]
Evaluates which Hankel factors the wavenumbers is calculated.
Possible choices are **sin**, **cos**, **j0**, or **j1**.
Returns
-------
factors: [np.ndarray]
Hankel factors.
"""
fac = self.factors[string][0]
return np.reshape(fac, (-1, 1))
[docs] def getWavenumbes(self, string):
""" Calcualtes the wavenumbers for the requested set of Hankel factors.
Parameters
----------
string: [str]
Evaluates for which Hankel factors the wavenumbers is calculated.
Possible choices are **sin**, **cos**, **j0**, or **j1**.
Returns
-------
wavenumbers: [np.ndarray]
Normed wavenumber factors for evaluation of the Hankel integral.
Divide by horizontal distance of the receiver to get
horizontal wavenumber lambda = sqrt(k_x² + k_y²).
"""
nc, nc0 = self.factors[string][1:]
fac_n = np.arange(nc0 - nc, nc0)
# 10 values per decade: 0.1 * np.log(10) = 0.23025850929940459
# predefined by the usage of these Hankel factors
lam_factor = np.exp(-fac_n * 0.1 * np.log(10))
return lam_factor
[docs]class wer_201_2018():
"""
Hankel factors after Werthmüller 2018 implemented from the empymod
package after consultation with Dieter Werthmüller.
Thank you very much!
The filter coefficient are published in:
Werthmüller, D., K. Key, and E. Slob, 2019, A tool for designing
digital filters for the Hankel and Fourier transforms in potential,
diffusive, and wavefield modeling: 84(2), F47-F56;
DOI: 10.1190/geo2018-0069.1
under the Apache 2.0 license.
"""
def __init__(self):
self.base = np.array([
8.653980893285999343e-04, 9.170399868578506730e-04,
9.717635708540675208e-04, 1.029752738345341345e-03,
1.091202360259058363e-03, 1.156318936280037154e-03,
1.225321288786770891e-03, 1.298441298197734088e-03,
1.375924682198844014e-03, 1.458031821470676028e-03,
1.545038634690233410e-03, 1.637237505747725754e-03,
1.734938266294200121e-03, 1.838469236921894991e-03,
1.948178330476121357e-03, 2.064434221206373505e-03,
2.187627583685520637e-03, 2.318172405660454734e-03,
2.456507379246002792e-03, 2.603097375137131825e-03,
2.758435004793547123e-03, 2.923042275846299189e-03,
3.097472346289414230e-03, 3.282311383351382318e-03,
3.478180533293271509e-03, 3.685738008752825322e-03,
3.905681300649099623e-03, 4.138749522080598965e-03,
4.385725892093550113e-03, 4.647440367666977670e-03,
4.924772432759196537e-03, 5.218654053788364042e-03,
5.530072811478751842e-03, 5.860075219597365645e-03,
6.209770241733286907e-03, 6.580333017937914017e-03,
6.973008813749223544e-03, 7.389117204870794542e-03,
7.830056511567909730e-03, 8.297308497682591086e-03,
8.792443350058264107e-03, 9.317124955107483966e-03,
9.873116490254248145e-03, 1.046228634904101083e-02,
1.108661441981135046e-02, 1.174819873906770597e-02,
1.244926254186268406e-02, 1.319216173291635694e-02,
1.397939280356631266e-02, 1.481360122115480543e-02,
1.569759031904567614e-02, 1.663433071714528685e-02,
1.762697030458526201e-02, 1.867884481811322647e-02,
1.979348905174003331e-02, 2.097464873531334692e-02,
2.222629312193481060e-02, 2.355262832652090660e-02,
2.495811146033087569e-02, 2.644746560896089893e-02,
2.802569570413705746e-02, 2.969810534264444649e-02,
3.147031460891126092e-02, 3.334827896114080786e-02,
3.533830924445728605e-02, 3.744709289831891358e-02,
3.968171642946591998e-02, 4.204968922592208086e-02,
4.455896879207719985e-02, 4.721798748965116282e-02,
5.003568087440296575e-02, 5.302151772380819111e-02,
5.618553185661335353e-02, 5.953835585119460899e-02,
6.309125677603129312e-02, 6.685617405236506106e-02,
7.084575957628110043e-02, 7.507342023504076645e-02,
7.955336296053971967e-02, 8.430064247129397115e-02,
8.933121186338736919e-02, 9.466197622039211612e-02,
1.003108494224146802e-01, 1.062968143451746283e-01,
1.126399866514113390e-01, 1.193616823889895873e-01,
1.264844896228644322e-01, 1.340323443416230609e-01,
1.420306108936851552e-01, 1.505061672234652703e-01,
1.594874951939309615e-01, 1.690047762990831426e-01,
1.790899930879975566e-01, 1.897770366412597776e-01,
2.011018204609656135e-01, 2.131024011570106513e-01,
2.258191063352318062e-01, 2.392946701171655144e-01,
2.535743767468323084e-01, 2.687062127671349665e-01,
2.847410282772538936e-01, 3.017327078129410922e-01,
3.197383514239506286e-01, 3.388184665571116749e-01,
3.590371713898611872e-01, 3.804624102975323052e-01,
4.031661821784713884e-01, 4.272247824042622599e-01,
4.527190592081252185e-01, 4.797346853730769523e-01,
5.083624461328500876e-01, 5.386985442530570767e-01,
5.708449233178135573e-01, 6.049096103082166609e-01,
6.410070786239048246e-01, 6.792586327676195523e-01,
7.197928159854947161e-01, 7.627458422329327359e-01,
8.082620539176789132e-01, 8.564944069583262376e-01,
9.076049847882753374e-01, 9.617655430324436594e-01,
1.019158086687097953e+00, 1.079975481742401877e+00,
1.144422103303022853e+00, 1.212714522384784388e+00,
1.285082233695329812e+00, 1.361768426844476521e+00,
1.443030803575896748e+00, 1.529142443766406734e+00,
1.620392723103028176e+00, 1.717088285521651159e+00,
1.819554073675149652e+00, 1.928134420893808709e+00,
2.043194208307565596e+00, 2.165120091018537973e+00,
2.294321797444361266e+00, 2.431233506198735572e+00,
2.576315305136151590e+00, 2.730054737463883274e+00,
2.892968440116886253e+00, 3.065603879901345419e+00,
3.248541193241110570e+00, 3.442395135709440002e+00,
3.647817147897402634e+00, 3.865497544561224075e+00,
4.096167834405143537e+00, 4.340603178295351583e+00,
4.599624994165753655e+00, 4.874103717369282940e+00,
5.164961725750846000e+00, 5.473176439271492555e+00,
5.799783604600056819e+00, 6.145880775710010013e+00,
6.512631002177978523e+00, 6.901266737578352739e+00,
7.313093981108034214e+00, 7.749496666359123154e+00,
8.211941311987910552e+00, 8.701981949908562441e+00,
9.221265347572652260e+00, 9.771536541883744320e+00,
1.035464470334362019e+01, 1.097254935013656407e+01,
1.162732693303372145e+01, 1.232117781324618733e+01,
1.305643365667540934e+01, 1.383556526940939513e+01,
1.466119090079532228e+01, 1.553608504199116247e+01,
1.646318774956322173e+01, 1.744561452546165370e+01,
1.848666678657500739e+01, 1.958984295904654616e+01,
2.075885023463468926e+01, 2.199761702862400270e+01,
2.331030618115176267e+01, 2.470132894631230158e+01,
2.617535981604930129e+01, 2.773735222865149197e+01,
2.939255521463919862e+01, 3.114653103598043415e+01,
3.300517387791185087e+01, 3.497472965617872376e+01,
3.706181700625468523e+01, 3.927344952507589682e+01,
4.161705934003131091e+01, 4.410052208441296528e+01,
4.673218336325475519e+01, 4.952088679849771324e+01,
5.247600374772711262e+01, 5.560746479634944706e+01,
5.892579312903905020e+01, 6.244213989259677078e+01,
6.616832166905808776e+01, 7.011686018497628936e+01,
7.430102439032442874e+01, 7.873487504841911289e+01,
8.343331198671111792e+01, 8.841212416722630962e+01,
9.368804274491759543e+01])
self.factor = np.array([1.0596741524693181])
self.j0 = np.array([
2.940900904253498815e+00, -1.601154970027019786e+01,
3.574488144287594338e+01, -3.775710631592443178e+01,
9.347313619702582344e+00, 1.903333998229986435e+01,
-1.266160545445113073e+01, -1.274822633210828471e+01,
1.499040669570122830e+01, 1.128114232815630835e+01,
-3.141814347069891511e+01, 2.185520874118305557e+01,
4.204688908817206361e+00, -2.039396651211594147e+01,
1.703214350560853418e+01, -3.858035665251962953e+00,
-5.664929771474184861e+00, 6.236013601339658763e+00,
-8.349084962847823643e-01, -5.076648434675173682e+00,
8.069437229646323928e+00, -7.693517618528387558e+00,
5.264540454351635645e+00, -2.378232295562936471e+00,
9.488986688091295696e-02, 1.211665432819921451e+00,
-1.641788465968915700e+00, 1.505783472315168181e+00,
-1.120376089133631847e+00, 7.202627618735454318e-01,
-4.293308092527763353e-01, 2.880519253847271255e-01,
-2.765909110132210857e-01, 3.540771325316399709e-01,
-4.702415607316920432e-01, 5.886079250132733032e-01,
-6.804407177384068639e-01, 7.365579045711867501e-01,
-7.516232294769148448e-01, 7.332518369461620278e-01,
-6.840940481215002089e-01, 6.145375558739786248e-01,
-5.263752170486778459e-01, 4.305607160730071659e-01,
-3.307504282674695872e-01, 2.429511351824290011e-01,
-1.749591466871897039e-01, 1.456366850874822316e-01,
-1.595894233786172844e-01, 2.280395161025021433e-01,
-3.413711397379035062e-01, 4.950015142033588056e-01,
-6.617826176259098414e-01, 8.243458750063681340e-01,
-9.468449720632085009e-01, 1.012979648696065826e+00,
-9.946084380621768029e-01, 8.932243609342348512e-01,
-7.019325270091921753e-01, 4.480078578272302381e-01,
-1.456381084707468188e-01, -1.604285443826692914e-01,
4.506359421881006022e-01, -6.820827325803067165e-01,
8.491728004462286705e-01, -9.252486609867097700e-01,
9.260533367474167443e-01, -8.402031897556576645e-01,
6.982570225067336045e-01, -4.946654104012083719e-01,
2.669962789856030194e-01, -1.042528939201309117e-02,
-2.319933928270803414e-01, 4.654858099944349514e-01,
-6.436180173401860882e-01, 7.780252822887177011e-01,
-8.272840865922600484e-01, 8.197748606729908794e-01,
-7.271116816505861502e-01, 5.992292042110317629e-01,
-4.191980233493379782e-01, 2.522483458861179417e-01,
-8.121597742330896597e-02, -2.652177515646544220e-02,
1.028351781199190462e-01, -8.958878867947789315e-02,
4.230973855067972356e-02, 8.503909842995018009e-02,
-2.142619566112575480e-01, 3.865561094837716705e-01,
-5.104727910435372662e-01, 6.353554205789004872e-01,
-6.678424785380696616e-01, 6.769052092946485910e-01,
-5.740789051340543514e-01, 4.533865101372042683e-01,
-2.311110970314779189e-01, 2.495118155831516429e-02,
2.508326503987150513e-01, -4.617568849058834579e-01,
7.064398492531979157e-01, -8.419497320020807862e-01,
9.894625411708843910e-01, -1.001995584015555441e+00,
1.027606495406535592e+00, -9.137586391106879979e-01,
8.347545540345707726e-01, -6.273557157264747497e-01,
4.895538883108097039e-01, -2.415916082478662963e-01,
1.027651673270286170e-01, 1.284590129111473078e-01,
-2.133464520148750654e-01, 3.796869867196227544e-01,
-3.706104772018420923e-01, 4.433984947033918766e-01,
-3.235750823461782666e-01, 2.994380356487793549e-01,
-7.888021235143724552e-02, -1.910472056257344828e-02,
3.043509756878240990e-01, -4.298489130677841108e-01,
7.233655581626995401e-01, -8.161936207950235556e-01,
1.054311731055406431e+00, -1.057859070622381603e+00,
1.189505276805579825e+00, -1.071819355511486993e+00,
1.077289848919564141e+00, -8.458743245226547636e-01,
7.443908003078970603e-01, -4.461507052923350813e-01,
2.868208498373531756e-01, 8.293917061440093941e-03,
-1.688393814622326516e-01, 3.930003549460732160e-01,
-5.148213607425790039e-01, 6.232961968619017412e-01,
-6.962255871378753014e-01, 6.749680427277989780e-01,
-7.157027353444623818e-01, 5.769543174180707945e-01,
-6.151602424632486299e-01, 3.892025319337197864e-01,
-4.506933666706669506e-01, 1.804687704243498336e-01,
-2.727244303574414830e-01, 1.174346932607069939e-02,
-1.137579568068424057e-01, -7.811421895993872488e-02,
1.779393134675293781e-02, -8.675295089134033022e-02,
1.335093001151206882e-01, -5.514694464669991220e-02,
2.441624548660136784e-01, -5.518933795483266236e-02,
3.328770195966135881e-01, -1.532145946078544430e-01,
3.579812932040986606e-01, -3.542418291888540516e-01,
3.138078445399646865e-01, -5.530414472417091165e-01,
2.996207005996078809e-01, -5.735406549188309944e-01,
4.429014477665748628e-01, -3.892397979507834505e-01,
6.319796721730875921e-01, -3.224362675105876264e-01,
5.280557016061210307e-01, -5.925152169592526885e-01,
3.097933933610434454e-01, -6.161398121500974989e-01,
5.824585099080735739e-01, -2.863603061506506120e-01,
5.579069742591284964e-01, -5.790177379949956737e-01,
1.778203303660013113e-01, -2.248159329174458654e-01,
4.570279405020428731e-01, -2.082694991220853942e-01,
-1.639463813834285966e-01, 1.219471051626438846e-01,
1.315907246465166380e-01, -1.561922951483753763e-01,
-9.000173966500184253e-02, 3.555651953426591794e-01,
-4.594671107387718334e-01, 4.051183609553802301e-01,
-2.826594312508529105e-01, 1.664628231342878961e-01,
-8.573386548433219179e-02, 3.940897144775467459e-02,
-1.632396319789397934e-02, 6.095399479649578345e-03,
-2.034197210251703809e-03, 5.959585007671103435e-04,
-1.489396688112877424e-04, 3.040429693056684047e-05,
-4.735729642331119566e-06, 4.982902654694162416e-07,
-2.645962918790746080e-08])
self.j1 = np.array([
-2.594301879688918743e-03, 2.339490069567079153e-02,
-9.478827695764932559e-02, 2.269100107268227362e-01,
-3.508900104631907935e-01, 3.515648778060092017e-01,
-1.994483426338835574e-01, 9.293484158449249674e-03,
8.048790331434139966e-02, -5.691508909482047296e-02,
2.075411236619714023e-02, -5.275228907131672418e-02,
1.348060780914397960e-01, -1.939017682982676627e-01,
1.854125006516780250e-01, -1.238162114378456441e-01,
5.435433909585837137e-02, -1.270327896980103649e-02,
7.560629247467166685e-03, -2.713648547163816441e-02,
5.383176198347366936e-02, -7.463014212375274070e-02,
8.413254388281748986e-02, -8.280914671407160754e-02,
7.393816085156304507e-02, -6.119942568117115594e-02,
4.746907894866293082e-02, -3.450474036399159977e-02,
2.313104032185092293e-02, -1.354138672902785445e-02,
5.604638952874172256e-03, 9.486899997092870969e-04,
-6.382314415042419746e-03, 1.091601370650388016e-02,
-1.467479730415811347e-02, 1.771829574440562591e-02,
-2.001990396213982823e-02, 2.152720536808908763e-02,
-2.214634467574278301e-02, 2.181839780171816040e-02,
-2.049006667874660528e-02, 1.819226618052163791e-02,
-1.498058384418661862e-02, 1.101097898088004151e-02,
-6.444381466006598655e-03, 1.527769500667217105e-03,
3.533045890696502513e-03, -8.470577184056853060e-03,
1.310891443998619434e-02, -1.722186749311952272e-02,
2.071070403809925978e-02, -2.340849145682623311e-02,
2.528895672008608583e-02, -2.620974028011485712e-02,
2.617463680378217042e-02, -2.501974545161460978e-02,
2.276381800709064221e-02, -1.923696679546721411e-02,
1.454384146700743972e-02, -8.605818651281786635e-03,
1.738844503682727928e-03, 5.947517515245405624e-03,
-1.385885800081448037e-02, 2.171703574368676753e-02,
-2.872814698074024550e-02, 3.461838232135941440e-02,
-3.858941536214401807e-02, 4.060565074712448042e-02,
-4.004001014668581715e-02, 3.725276070158690944e-02,
-3.184536604257427045e-02, 2.462644816656397312e-02,
-1.539955242663786951e-02, 5.437594699869644464e-03,
5.315602443775078317e-03, -1.514798704965643165e-02,
2.414409295479183482e-02, -3.032245442310742278e-02,
3.416403682006976389e-02, -3.375211066545635158e-02,
3.045425146826126819e-02, -2.272683240323554107e-02,
1.317397219169818418e-02, -6.119848802195663844e-04,
-1.115685958739305421e-02, 2.344848055710842955e-02,
-3.169397030448774938e-02, 3.817520279763731567e-02,
-3.809651169502087376e-02, 3.552394037620545952e-02,
-2.568176960339346032e-02, 1.487169880055576841e-02,
1.997965754252607837e-03, -1.644546122863936588e-02,
3.480415781283187349e-02, -4.677537576184215284e-02,
6.111820051573001178e-02, -6.590953690231306228e-02,
7.332112878498593667e-02, -6.941405787722150500e-02,
7.037058951810715168e-02, -5.921683716330301134e-02,
5.655416878814193554e-02, -4.097804417398832194e-02,
3.810136479256000241e-02, -2.058552251762190213e-02,
2.012009960638341116e-02, -1.904932224648239053e-03,
5.324133880430357950e-03, 1.357085636140308721e-02,
-5.592398994521355186e-03, 2.589415408390925710e-02,
-1.292224858788840886e-02, 3.561835458678683924e-02,
-1.700336052804993919e-02, 4.302272543873187499e-02,
-1.764501119714675936e-02, 4.776824187802861110e-02,
-1.404967795338326469e-02, 4.909105531804618811e-02,
-5.268363386320836297e-03, 4.641444115511134810e-02,
9.124170345995139680e-03, 3.991380510079723526e-02,
2.859680534219202416e-02, 3.073032173944831302e-02,
5.151734601415858261e-02, 2.097520864952298614e-02,
7.511832720419686638e-02, 1.354148686126329000e-02,
9.576070647813632319e-02, 1.140272504589405489e-02,
1.097360846564945647e-01, 1.641156037670641471e-02,
1.142607074214286172e-01, 2.822722328175052836e-02,
1.079741026261866604e-01, 4.397838549427537935e-02,
9.072532008032478668e-02, 5.865331390992817306e-02,
6.310580311949447185e-02, 6.592890500839264367e-02,
2.624768245556352575e-02, 5.954092122676180737e-02,
-1.789442088462365327e-02, 3.555984804371448149e-02,
-6.551973841104229146e-02, -4.489823587807304991e-03,
-1.089659607197719371e-01, -5.052377390460756346e-02,
-1.345089168583019634e-01, -8.286532978746213862e-02,
-1.229744228392625760e-01, -7.728334735138943368e-02,
-5.831402820224085293e-02, -2.088719549713474385e-02,
5.205792146829352207e-02, 6.258826423151815643e-02,
1.537593434905694390e-01, 9.979165672991972824e-02,
1.545936043232943868e-01, 1.442374184540378551e-02,
1.426608350698001064e-02, -1.462978689171700042e-01,
-1.285617304470880462e-01, -1.519560295563492092e-01,
-4.700754607069132507e-02, 1.106161196544740571e-01,
1.168039951766200457e-01, 2.366820009670577984e-01,
-1.283956636561735531e-01, 1.406024618127853579e-02,
-4.011626826391421763e-01, 2.656020926376266300e-01,
-1.291178852844470371e-01, 5.190017993719785450e-01,
-5.307617026829943851e-01, 2.388177316428383157e-01,
-4.213777919843004205e-01, 7.425023205388556757e-01,
-6.039491377476802203e-01, 2.966551251977549986e-01,
-3.157763414193298090e-01, 5.842678377857475347e-01,
-7.375386163194860289e-01, 6.321787190293964853e-01,
-3.871768116388242809e-01, 1.635153404860406612e-01,
-3.134030444480916111e-02, -2.007964312111037639e-02,
2.747060309825222202e-02, -1.969027573928272892e-02,
1.080160200622757964e-02, -4.917741573140644099e-03,
1.900872502605262596e-03, -6.228037791451519409e-04,
1.698280169898911708e-04, -3.718088395106743311e-05,
6.139140918052498045e-06, -6.796809441105533029e-07,
3.780807126974975489e-08])
if __name__ == '__main__':
pass
# The End