python - 吹奏楽器のインピーダンス

Question

この論文https://hal.archives-ouvertes.fr/file/index/docid/683477/filename/clarinette-logique-8.pdfで説明されているモデルを再現しようとしています。

以下は、半径 a と長さ L の円柱の伝達行列を返すメソッドです。

# density of air
rho = 1.2250

# viscosity of air
eta = 18.5

# transfer matrix for a cylinder
def Tcyl(a, L):
    rv = a * math.sqrt(rho*omega/eta)
    Z0 = (rho * constants.c) / (constants.pi * a**2)
    Zc = Z0 * complex(1 + 0.369/rv, -(0.369/rv + 1.149/rv**2))
    gamma = k * complex(1.045/rv + 1.080/rv**2, 1 + 1.045/rv)
    return
       np.matrix([[cmath.cosh(gamma*L),Zc*cmath.sinh(gamma*L)],
       [(1.0/Zc)*cmath.sinh(gamma*L),cmath.cosh(gamma*L)]], dtype=complex)

入力インピーダンスは、次のように異なる周波数で計算されます。

for f in range(1,4000):
    
    # angular frequency
    omega = 2*constants.pi*f

    # wave number
    k = omega/constants.c

    # transfer matrix for a cylinder of radius 7.5mm and length 100mm
    T = Tcyl(7.5, 100.0)

    # radiation impedance
    Zr = complex(0,0)

    # [P,U] vector (pression and velocity)
    T = np.dot(T,np.matrix([[Zr],[complex(1,0)]], dtype=complex))

    # input impedance (Z = P/U)
    Z = T.item(0)/T.item(1)

再生周波数は、式Ｉｍ［Ｚ］＝０を満たす。ZI の虚数部をプロットすると、次の図が得られます:間違ったインピーダンス

期待される出力は次のようになるはずなので、これは明らかに間違っています:正しいインピーダンス

私は何を間違っていますか？ありがとうございました。