Quiero calcular repetidamente una integral compleja bidimensional usando dblquad de scipy.integrate. Como la cantidad de evaluaciones será bastante alta, me gustaría aumentar la velocidad de evaluación de mi código.

Dblquad no parece ser capaz de manejar integrands complejos. Por lo tanto, he dividido el complejo integrand en una parte real y otra imaginaria:

def integrand_real(x, y):
    R1=sqrt(x**2 + (y-y0)**2 + z**2)
    R2=sqrt(x**2 + y**2 + zxp**2)
    return real(exp(1j*k*(R1-R2)) * (-1j*z/lam/R2/R1**2) * (1+1j/k/R1))

def integrand_imag(x,y):
    R1=sqrt(x**2 + (y-y0)**2 + z**2)
    R2=sqrt(x**2 + y**2 + zxp**2)
    return imag(exp(1j*k*(R1-R2)) * (-1j*z/lam/R2/R1**2) * (1+1j/k/R1))

y0, z, zxp, ky lam son variables definidas de antemano. Para evaluar la integral sobre el área de un círculo con radio ra, uso los siguientes comandos:

from __future__ import division
from scipy.integrate import dblquad
from pylab import *

def ymax(x):
    return sqrt(ra**2-x**2)

lam = 0.000532
zxp = 5.
z = 4.94
k = 2*pi/lam
ra = 1.0

res_real = dblquad(integrand_real, -ra, ra, lambda x: -ymax(x), lambda x: ymax(x))
res_imag = dblquad(integrand_imag, -ra, ra, lambda x: -ymax(x), lambda x: ymax(x))
res = res_real[0]+ 1j*res_imag[0]

Según el perfilador, los dos integrandes se evalúan unas 35000 veces. El cálculo total toma aproximadamente un segundo, lo que es demasiado largo para la aplicación que tengo en mente.

Soy un principiante de la computación científica con Python y Scipy y me encantaría recibir comentarios que indiquen formas de mejorar la velocidad de evaluación. ¿Hay formas de volver a escribir los comandos en las funciones integrand_real y integrand_complex que podrían llevar a mejoras de velocidad significativas?

¿Tendría sentido compilar esas funciones usando herramientas como Cython? En caso afirmativo: ¿Qué herramienta se ajustaría mejor a esta aplicación?