def fun(phi, theta):
    term_1 = (np.cos(phi) * np.cos(2.*theta))**2
    term_2 = 0.5 * np.sin(phi)**2
    return term_1 + term_2

import numpy as np
import matplotlib.pyplot as plt

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

theta = np.linspace(-halfpi, halfpi, 201)
phi   = np.linspace(-pi,     pi,     401)

Theta, Phi = np.meshgrid(theta, phi, indexing='ij')
print Theta.shape
print Phi.shape

F = fun(Phi, Theta)

plt.figure()
plt.imshow(F)
plt.colorbar()
plt.show()