Scipy: Acelerando o cálculo de uma integral complexa 2D

Eu quero calcular repetidamente uma integral complexa bidimensional usando dblquad de scipy.integrate. Como o número de avaliações será bastante alto, gostaria de aumentar a velocidade de avaliação do meu código.

Dblquad não parece ser capaz de lidar com integrandos complexos. Assim, eu dividi o integrando complexo em uma parte real e uma parte imaginária:

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, k e lam são variáveis ​​defind antecipadamente. Para avaliar a integral sobre a área de um círculo com raio ra utilizo os seguintes 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]

De acordo com o profiler, os dois integrantes são avaliados cerca de 35.000 vezes. O cálculo total demora cerca de um segundo, o que é muito longo para o aplicativo que tenho em mente.

Eu sou um iniciante em computação científica com Python e Scipy e ficaria feliz com comentários que apontam maneiras de melhorar a velocidade de avaliação. Existem maneiras de reescrever os comandos nas funções integrand_real e integrand_complex que poderiam levar a melhorias significativas na velocidade?

Faz sentido compilar essas funções usando ferramentas como o Cython? Se sim: Qual ferramenta melhor se encaixaria nessa aplicação?