こちらの文献を参考に作成してみました。
http://library.jsce.or.jp/jsce/open/00550/2010/16-01-0030.pdf
元のWolf氏の文献を見つけらず「z0＝近似円錐地盤の頂点高さ」にどのような数値を設定すべきか不明です。
import math
import cmath
import scipy.special
G = 1.8E+07
VS = 100.0
poisson = 0.33
R = 6.0
Z0 = 15.0
for i in range(0,80):
freq = float(i) * 0.1
omega = 2.0 * math.pi * freq
a0 = omega * R / VS
K = 8.0 * G * R / ( 2.0 - poisson )
ca0 = math.pi * ( 2.0 -poisson ) / 8.0
a0cao = a0 * ca0
Zsway = K * complex(1.0, a0cao) / ( G * R )
mu = 0.3 * math.pi * ( poisson - 1.0 / 3.0 )
Kt = 8.0 * G * R**3.0 / ( 3.0 * ( 1.0 - poisson ) )
a02 = a0**2.0
kta0 = 1.0 - 1.0 / 3.0 * mu / math.pi * Z0 / R * a02 - 1.0 / 3.0 * a02 / ( ( 2.0 * R / Z0 )**2.0 + a02 )
cta0 = Z0 / ( 6.0 * R) * a02 / ( ( 2.0 * R / Z0 )**2.0 + a02 )
Zrock = Kt * complex( kta0, a0 * cta0 ) / ( G * R**3.0 )
print(a0, ',', Zsway.real, ',', Zsway.imag, ',', Zrock.real, ',', Zrock.imag)