Did I properly implement double integral equation in dblquad Python? Helicoidal Surface Theory

Question

I show the equations (simplified Helicoidal Surface Theory for implementation purposes) that I want to calculate numerically using Python.

and the code:

from scipy.integrate import quad
import numpy as np
from scipy import interpolate
from scipy.integrate import dblquad
import time
start_time = time.time()

input="-0.5 0.0 \
-0.3 0.9    \
0.0 0.8 \
0.3 0.4 \
0.5 0.02"

input_coordinates = np.genfromtxt(input.splitlines()).reshape(-1,2) # shape to 2 columns, any number of rows
x_coordinates = input_coordinates[:,0]
H_values = input_coordinates[:,1]
H_interpolation = interpolate.InterpolatedUnivariateSpline(x_coordinates, H_values)

def complex_dblquad(func, a, b, g, h, **kwargs):
    def real_func(z, x):
        return np.real(func(z, x))
    def imag_func(z, x):
        return np.imag(func(z, x))
    real_integral = dblquad(real_func, a, b, g, h, **kwargs)
    imag_integral = dblquad(imag_func, a, b, g, h, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

complex_integral = complex_dblquad(lambda z,x: np.sqrt(1+z*z)**2*(2/np.sqrt(1+z*z))**2*H_interpolation(x)*np.exp(1j*2/np.sqrt(1+z*z)*x), 0, 1, -0.5, 0.5)

print("Quad",complex_integral)

print("--- %s seconds ---" % (time.time() - start_time))

The question is - is it implemented correctly? If not - what is wrong?

Answer 1

It looks correct, but it is difficult to interpret it is not a good practice to expand all the variables in one single expression, you can use the use definitions that match your equations more easily.

def integrand_function(z, x):
  Mr = np.sqrt(1 + z**2)
  kx = 2/Mr;
  # merging the integrand of Psi in order to use the dblquad function.
  return 2 * Mr**2 * kx**2 H_interpolation(x)*np.exp(1j * kx * x)

Then your complex integral is evaluated with

complex_integral = complex_dblquad(integrand_function, 0, 1, -0.5, 0.5)

And hopefully you can prove more easily the correctness this way.