import matplotlib.pyplot as plt
import numpy as np
import math

def pdf(b,g,beta,gamma,dbeta,dgamma):
    return (1/np.sqrt(2*np.pi*np.power(dgamma,2))) * (1/np.sqrt(2*np.pi*np.power(dbeta,2))) * np.exp(-np.power((g - gamma),2)/(2*np.power(dgamma,2))) * np.exp(-np.power((b - beta),2)/(2*np.power(dbeta,2))) 


def plot_isotope(ax,beta,gamma,dbeta,dgamma):

	g = np.arange(0, np.radians(61), np.radians(61)/120)
	b = np.arange(0,beta + 4*dbeta,0.01)
	b,g = np.meshgrid(b,g)
	Z = pdf(b,g,beta,gamma,dbeta,dgamma)

	ctf = ax.contourf(g,b,Z,40,cmap='jet')
	ax.scatter(gamma,beta,c='k',s=5)

	ax.set_thetamax(60)
	ax.set_rlabel_position(-22.5)  # get radial labels away from plotted line
	ax.grid(True)

beta = 0.4
gamma = 30 # degrees
dbeta = 0.1
dgamma = 5 # degrees

fig = plt.figure()
ax = fig.add_subplot(111,projection='polar')
fig.subplots_adjust(right=1,top=1,left=0.0,bottom=0.05)

plot_isotope(ax,beta,np.radians(gamma),dbeta,np.radians(dgamma))

ax.set_xlabel(r'$\beta$',fontsize=12)
ax.text(np.radians(30),ax.get_rmax()*1.2,r"$\gamma$",rotation=-60,ha='center',va='center',fontsize=12)
ax.tick_params('both',labelsize=13)

plt.show()