建屋側面地盤のインピーダンスを簡易計算する手法としてNOVAKのばねがよく知られています。
こちらのPDFにその理論式が掲載されています。
これをもとにpythonで作成したプログラムが以下です。
import math
import cmath
import scipy.special
G = 1.0E+09
VS = 700.0
poisson = 0.3
R = 10.0
HJ = 1.0
h = 0.00
VP = VS * math.sqrt( ( 2.0 * ( 1.0 - poisson ) ) / ( 1.0 - 2.0 * poisson ) )
for i in range(1,2000):
freq = float(i) * 0.01
omega = 2.0 * math.pi * freq
a0 = omega * R / VS
b0 = omega * R / VP
Z = cmath.sqrt( 1.0 + 2.0 * h * complex( 0.0, 1.0 ))
a0Z = a0 * complex( 0.0, 1.0 ) / Z
b0Z = b0 * complex( 0.0, 1.0 ) / Z
# B = math.sqrt( ( 2.0 * ( 1.0 - poisson ) ) / ( 1.0 - 2.0 * poisson ) )
# b0Z = a0Z / B
# 第２種変形ベッセル関数
k0a0 = scipy.special.kv(0, a0Z)
k0b0 = scipy.special.kv(0, b0Z)
k1a0 = scipy.special.kv(1, a0Z)
k1b0 = scipy.special.kv(1, b0Z)
p1 = - math.pi * HJ * G * a0**2.0
p2 = 4.0 * k1b0 * k1a0 + a0Z * k1b0 * k0a0 + b0Z * k0b0 * k1a0
p3 = b0Z * k0b0 * k1a0 + a0Z * k1b0 * k0a0 + b0Z * a0Z * k0b0 * k0a0
kuj = p1 * p2 / p3
print(freq, ',', kuj.real, ',', kuj.imag)