Polar plot in Matplotlib by mapping into Cartesian coordinate

I have a variable (P) which is a function of angle (theta):

In this equation the K is a constant, theta_p is equal to zero and I is the modified Bessel function of the first kind (order 0) which is defined as:

Now, I want to plot the P versus theta for different values of constant K. First I calculated the parameter I and then plug it into the first equation to calculate P for different angles theta. I mapped it into a Cartesian coordinate by putting :

x = P*cos(theta)

y = P*sin(theta)

Here is my python implementation using matplotlib and scipy when the constant k=2.0:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

def integrand(x, a, k):
   return a*np.exp(k*np.cos(x))

theta = (np.arange(0, 362, 2))
theta_p = 0.0

X = []
Y = []

for i in range(len(theta)):
    a = (1 / np.pi)
    k = 2.0
    Bessel = quad(integrand, 0, np.pi, args=(a, k))
    I = list(Bessel)[0]
    P = (1 / (np.pi * I)) * np.exp(k * np.cos(2 * (theta[i]*np.pi/180. - theta_p)))
    x = P*np.cos(theta[i]*np.pi/180.)
    y = P*np.sin(theta[i]*np.pi/180.)
    X.append(x)
    Y.append(y)

plt.plot(X,Y,  linestyle='-', linewidth=3, color='red')
axes = plt.gca()
plt.show()

I should get a set of graphs like the below figure for different K values:

(Note that the distributions were plotted on a circle of unit 1 to ease visualization)

However it seems like the graphs produced by the above code are not similar to the above figure. Any idea what is the issue with the above implementation? Thanks in advance for your help.

Here is how it looks like (for k=2):

The reference for these formulas are the equation 5 and 6 that you could find here

You had a mistake in your formula.

Your formula gives the delta of your function above a unit circle. So in your function to get the plot you want, simply add 1 to it.

Here is what you want, with some tidied up python. ...note you can do the whole calculation of the 'P' values as a numpy vector line, you don't need to loop over the indicies. ...also you can just do a polar plot directly in matplotlib - you don't need to transform it into cartesian.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

theta = np.arange(0, 2*np.pi+0.1, 2*np.pi/100)

def integrand(x, a, k):
   return a*np.exp(k*np.cos(x))

for k in np.arange(0, 5, 0.5):
    a = (1 / np.pi)
    Bessel = quad(integrand, 0, np.pi, args=(a, k))
    I = Bessel[0]    
    P = 1 + (1/(np.pi * I)) * np.exp(k * np.cos(2 * theta))
    plt.polar(theta, P)

plt.show()