Better way to calculate double-scattering diffraction using cartesian product of arrays?
You can replace:
d12 = np.array([a+b for a, b in list(itertools.product(p1.T, p2.T))])

with something like:
p1 = p1.T
p2 = p2.T
p3 = p3.T
d12 = p1[:,np.newaxis,:] + p2[np.newaxis,:,:]
d12 = my_d12.reshape((len(p1)*len(p2),2))

I find it most of the times easier to use the first index of an array for point_index and the second index for the dimensions, hence the .T's
With the use of the magic index np.newaxis at the right places you can create numpy array's of shape (M,N) with normal operators acting on arrays of shape (M) and (N).
With the reshape method d12 changed again to the shape in your original solution.