# Bibliotecas
import numpy as np
import scipy.special as sp
import matplotlib.pyplot as plt

# Função de Theodorsen
def Ck(k):
    y = (sp.hankel2(1,k))/((sp.hankel2(1,k))+(1j*(sp.hankel2(0,k))))
    return y

# Exibição
rf = np.linspace(1e-5,1.0)
Tf = Ck(rf)
fig, axs = plt.subplots(2)
axs[0].set_title('Magnitude')
axs[1].set_title('Fase')
axs[0].plot(rf,np.abs(Tf))
axs[1].plot(rf,((np.arctan((np.imag(Tf))/(np.real(Tf))))*(180/np.pi)))
plt.show()