Hi I am trying to find out this simple equation which gives the value of deposition in bends for different Stokes number

I am trying to plot I vs St. I is given by following equation and the graph that I wanna reproduce is attached here. For some reasons, I am not able to get exact same graph. Any help is appreciated. I think the problem might be from the value of root. myimagewithformula_and_desired_plot

my code is

from scipy.optimize import fsolve
from scipy import optimize

Ro = 8 
St = 0.1
dep = []
St_1 = []
    
while St < 1.7:    
        
    K = (4/Ro) * St  
        
    A = ((((1 + (K**2))**(1/2)) + 1)**(1/2))/1.41
     
    B = 0.5 * (K/A)
        
    def func(x):
        return((((A**2) + (B**2)*np.cos(A*x)*np.cos(B*x))  + ( 
               (A**3)*  np.sinh(A*x)*np.cos(B*x)) - ((B**3) *np.sin(B*x)*np.cosh(A*x)))) root = fsolve(func, 10)  
        
    print("This is root", root)  
        
    eta = (B**3) /((A**2) + (B**2)) *(np.cos(B*root)*np.sinh(A*root)) + np.sin(B*root)* 
        (np.sinh(A*root) + ((A**3) / ((A**2) + (B**2)) *np.cosh(A*root)) )

    BB = (eta - np.exp(root))**2

    CC = (eta + np.exp(root))**2

    ZC1 =(( 1 - (Ro**2 *( BB/CC))))

    ZC = (( abs(ZC1)))**(1/2)

    print("This is zc", ZC)

    m = (Ro**2 +1 )*ZC - 1/3 * ZC**3

    n = -1 + (1/(eta**2) *np.exp(2 *root))

    p = (1 + (1/(eta**2) *np.exp(2 *root))) *Ro

    o = ZC *((1- (ZC**2))**(1/2)) + np.arcsin(ZC1)

    Imp0 = 1/(3.14*Ro)*((m *n) + (Ro *p *o))

    imp = 1-Imp0

    dep.append(imp)

    St_1.append(St)

    St+= 0.1
fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1])
ax.plot(St_1, dep, ':.') #plot between stokes no vs deposition