'Errors in the output of scipy's residue function

I have a pade approximation of an admittance function (Y) with coefficients a and b. I then tried to do a partial fractions expansion (with scipy's residue function) using coefficients a and b in order to obtain the residues (r), the poles (k) and the direct term (k). However, when I try to recompute Y_calc from the obtained values of r, p, and k evaluated at every frequency, the results were off as shown in the attached plot. Also the coefficients obtained using scipy's invres were quite different from the original. I have tried to troubleshoot but cannot pinpoint what might be the source of the error. Kindly assist. My code is provided below.

# import libraries
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sps

# frequency array
F = np.array([2.00e+05, 1.35e+05, 9.15e+04, 6.19e+04, 4.19e+04, 2.83e+04,
       1.92e+04, 1.30e+04, 8.77e+03, 5.93e+03, 4.01e+03, 2.72e+03,
       1.84e+03, 1.24e+03, 8.40e+02, 5.69e+02, 3.85e+02, 2.60e+02,
       1.76e+02, 1.19e+02, 8.05e+01, 5.45e+01, 3.69e+01, 2.49e+01,
       1.69e+01, 1.14e+01, 7.72e+00, 5.22e+00, 3.53e+00, 2.39e+00,
       1.61e+00, 1.09e+00, 7.40e-01, 5.00e-01])

# admittance (Y) array
Y = np.array([0.03635013+0.00070057j, 0.03622358+0.00225245j,0.03602907+0.00347189j, 0.03546447+0.00521688j,
        0.03422052+0.00714964j, 0.03223397+0.00959102j,0.02929454+0.01175767j, 0.02495787+0.01340181j,
        0.02015355+0.01391555j, 0.0154913 +0.01303236j,0.01178335+0.01132911j, 0.00902089+0.00930647j,
       0.00709459+0.00743243j, 0.00575937+0.00586753j,0.0048024 +0.0045821j , 0.00406744+0.00359585j,
       0.0034929 +0.00282663j, 0.0030694 +0.00218837j,0.00275143+0.00166534j, 0.00251572+0.00125786j,
       0.0023522 +0.00092682j, 0.00224678+0.00067734j,0.00218716+0.00049425j, 0.002148  +0.00035753j,
       0.00212246+0.00026348j, 0.00209845+0.00020584j,0.00208301+0.00016725j, 0.00206836+0.00014835j,
       0.00204372+0.00014058j, 0.00201088+0.0001369j ,0.00198245+0.00013782j, 0.00194749+0.00013568j,
       0.00191001+0.00013381j, 0.00187448+0.00012814j])

# numerator
a =  np.array([3.01279728e-19, 1.70866222e-18,
       0.00000000e+00, 0.00000000e+00, 1.09422259e-08, 5.39497245e-07,
       2.66646345e-05, 1.92254916e-04, 2.21808439e-03, 1.62562196e-03])
# denominator
b = np.array([8.40390817e-18, 0.00000000e+00, 0.00000000e+00, 7.79149876e-11,
       2.69403425e-07, 2.05058225e-05, 3.08320721e-03, 1.01323926e-01,
       9.83433619e-01, 1.00000000e+00])

def pade(f,a,b):
    """
    function to compute Y from the pade coefficients
    """
    Yf = np.polyval(a,np.sqrt(1j*f))/np.polyval(b,np.sqrt(1j*f))
    return Yf

# compute the partial fraction expansion
rpk = sps.residue(a, b) 

recalculated_Y = np.zeros(len(F), dtype=complex)
for idx, freq in enumerate(F):
    # adding the terms of the partial fraction expansion
    recalculated_Y[idx] = np.sum((rpk[0]/(np.sqrt(1j*freq) - rpk[2]))) + rpk[2]

# compute Yf using the pade function
Yf = pade(F, a, b)

# get back the original coefficients (which are different from the original a and b)
coeffs = sps.invres(rpk[0], rpk[2], rpk[2])
coeffs

# (array([3.58499548e-02+0.00000000e+00j, 0.00000000e+00-4.44089210e-16j,
#         1.98951966e-12-7.10542736e-15j, 3.32374977e+05+3.12866177e-10j,
#         1.14923919e+09-2.90572643e-07j, 6.41959948e+10+0.00000000e+00j,
#         3.17288504e+12-2.44140625e-04j, 2.28768464e+13-1.95312500e-02j,
#         2.63934868e+14+5.93750000e-01j, 1.93436426e+14-7.26562500e-01j]),
#  array([1.00000000e+00+0.00000000e+00j, 0.00000000e+00+0.00000000e+00j,
#         5.82076609e-11+0.00000000e+00j, 9.27128022e+06+8.84756446e-09j,
#         3.20569216e+10-7.62939453e-06j, 2.44003410e+12+2.44140625e-04j,
#         3.66877784e+14+3.12500000e-02j, 1.20567626e+16+1.00000000e+00j,
#         1.17020986e+17+8.00000000e+00j, 1.18992257e+17+0.00000000e+00j]))

# plt.plot(Y.real, Y.imag, 'b-', color = "red", linewidth = 2.0, label = "expt")
# plt.plot(Yf.real, Yf.imag, 'b-', linewidth = 2.0, label = "pade_fit")
# plt.scatter(recalculated_Y.real, recalculated_Y.imag, color = "orange", label = "calculated") # computed transfer function
# plt.legend(loc = "best")