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?