Impedance of a blowing instrument

I'm trying to reproduce the model described in this paper https://hal.archives-ouvertes.fr/file/index/docid/683477/filename/clarinette-logique-8.pdf.

Here is a method that returns the transfer matrix for a cylinder of radius a and length 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)

The input impedance is calculated for different frequencies as follow.

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)

The playing frequency satisfy the equation Im[Z]=0. When plotting the imaginary part of Z I get the following figure : wrong impedance

This is clearly wrong as the expected output should be something like this : correct impedance

What am I doing wrong? Thank you.

Solution

You have

    Z0 = (rho * constants.c) / (constants.pi * a**2)

The impedance of a clarinet depends on the speed of SOUND, not the speed of LIGHT. Replace constants.c with 343 and your results will be closer. I'm still not sure it's quite right, but closer.

As a clarinetist, I try to make people think my fingers move as fast as light, but it ain't so.