Numerical integration methods for surface integral

Hello. I have to perform this integral: and this other one that is basically the same: where:

• a is a known complex vector that is the FFT of the derivative of a real laboratory mesuerement
• are known
• S is a circumference surface with radius R (known)

Since a is a complex numerical vector, I tried to use cumtrapz and integral functions to compute this integral. The two methods give me worng results and the results are different depending on the used method.

Moreover, I have seen that the final result of the method 2 with the integral function is the same for the first integral (pa_a) and for the second one (p).

Could you please help me to find the right way to perform this calculation?

%% DATA
% Input
y = data;

% Compute acceleration from measured velocity
a = diff(y);
A = fft(a);

% Measurements data
Fs = 44100;     % Measure sampling frquency
f_start = 2;    % Initial measure frequency [Hz]
f_end = 5000;   % Final measure frequency [Hz]

% Parameters
rho_0 = 1.225; % kg/m3
ra = 1; % m 
rc = 0; % m 
p_0 = 1; % Standard pressure
c = 340; % Sound speed m/s

R = 0.254;   % Radius
Sc = pi*R^2; % Membrane surface [m2]

f = linspace(f_start, f_end, length(diff(fft(y)))); % Frequency vector
k = 2*pi*f/c; % Wavenumber

%% Integral computation
% METHOD 1: trapezoidal integration
domain1 = linspace(0, 2*pi, length(A));
domain2 = linspace(0, R, length(A));

pa_a_cumtrapz = rho_0/(2*pi) * cumtrapz(domain1, cumtrapz(domain2, A/abs(ra-rc)));
p_cumtrapz = rho_0/(2*pi) * cumtrapz(domain1, cumtrapz(domain2, A/abs(ra-rc) .* exp(-1i*k*abs(ra-rc))));

%------------------------

% METHOD 2: Double integral (I cannot use integral2 because I have a vector and not a matrix)
% First internal integral computation
integrand = @(om) A/abs(ra-rc);
temp = integral(integrand, 0, R, "ArrayValued",true);

% Second integral computation
integrand2 = @(om) temp;

pa_a_integral = rho_0/(2*pi).* integral(integrand2, 0, 2*pi, "ArrayValued",true);
p_integral = rho_0/(2*pi) * integral(integrand2, 0, 2*pi, "ArrayValued",true) .* exp(-1i*k*abs(ra-rc));

%% INTEGRATION RESULT PLOTS
figure;
loglog(abs(pa_a_cumtrapz));
hold on
loglog(abs(p_cumtrapz)*0.3) % *0.3 is a shift to check differences between two spectra
legend ('pa_a cumtrapz', 'p cumtrapz low shifted')
xlim([f_start f_end]);
grid minor

%------------------------

figure;
loglog(abs(pa_a_integral));
hold on
loglog(abs(p_integral)*0.3) % *0.3 is a shift to check differences between two spectra
legend ('pa_a integral', 'p integral low shifted')
xlim([f_start f_end]);
grid minor