How to integrate over a surface with non-uniform spherical coordinate?

Dear community,

I need to integrate the P over the entire cone surface composed of theta and phi. The numerical value of P(theta,phi) is known, and the ux and uy in the direction cosine are known too. Since the theta and phi are non-uniform in spherical coordinate, I tried the below two integration methods, why the results are so different (2.649e-12 Vs. 5.7498e-12)? I am not sure which one is correct. If both are inaccurate, can you please give me a best solution for this issue?

clear all;
load data.mat;

% Method 1:
% ux is Mx1, uy is Nx1
[Ux, Uy] = meshgrid(ux, uy);
Ux = Ux';  %  M×N
Uy = Uy';  

% compute theta and phi
Uz = real(sqrt(1 - Ux.^2 - Uy.^2));
theta = acos(Uz);  % M×N
phi = atan2(Uy, Ux);  % M×N

% compute integration ∬ P(θ,ϕ) sinθ dθ dϕ
integrand = P.* sin(theta) ;

% compute dθ and dϕ
% 1. compute dθ 
dtheta = zeros(size(theta));
dtheta(2:end, :) = diff(theta, 1, 1);  % (M-1)×N
dtheta(1, :) = dtheta(2, :);  % boundary

% 2. compute dphi 
dphi = zeros(size(phi));
dphi(:, 2:end) = diff(phi, 1, 2);  % M×(N-1)
dphi(:, 1) = dphi(:, 2);  

% 3. sum
weight = sin(theta) .* dtheta .* dphi;
result = sum(integrand .* weight, 'all');


% Mehtod 2:
% compute integration ∬ P(θ,ϕ) sinθ dθ dϕ
integrand = P.* sin(theta) ;

% generate the uniform mesh of theta_uniform and phi_uniform
theta_uniform = linspace(min(theta(:)), max(theta(:)), 200);
phi_uniform = linspace(min(phi(:)), max(phi(:)), 200);
[THETA, PHI] = meshgrid(theta_uniform, phi_uniform);

% 2D interpolation
F = scatteredInterpolant(theta(:), phi(:), integrand(:));
P_interp = F(THETA, PHI);

% 3. trapz integration of uniform mesh
result2 = trapz(phi_uniform, trapz(theta_uniform, P_interp, 1), 2);

Answer 1

[Mathieu NOE]

hello

seems to me in the first method you are applying sin(theta) twice (one only is needed)

integrand = P.* sin(theta)

weight = sin(theta) .* dtheta .* dphi;

so if I remove one of the two , both results are now closer :

result = 5.2374e-12

result2 = 5.7498e-12

load data.mat;
warning off

% Method 1:
% ux is Mx1, uy is Nx1
[Ux, Uy] = meshgrid(ux, uy);
Ux = Ux';  %  M×N
Uy = Uy';  
% compute theta and phi
Uz = real(sqrt(1 - Ux.^2 - Uy.^2));
theta = acos(Uz);  % M×N
phi = atan2(Uy, Ux);  % M×N
% compute integration ∬ P(θ,ϕ) sinθ dθ dϕ
% integrand = P.* sin(theta) ;
% compute dθ and dϕ
% 1. compute dθ 
dtheta = zeros(size(theta));
dtheta(2:end, :) = diff(theta, 1, 1);  % (M-1)×N
dtheta(1, :) = dtheta(2, :);  % boundary
% 2. compute dphi 
dphi = zeros(size(phi));
dphi(:, 2:end) = diff(phi, 1, 2);  % M×(N-1)
dphi(:, 1) = dphi(:, 2);  
% 3. sum
weight = sin(theta) .* dtheta .* dphi;
% result = sum(integrand .* weight, 'all')
result = sum(P .* weight, 'all')
% Mehtod 2:
% compute integration ∬ P(θ,ϕ) sinθ dθ dϕ
integrand = P.* sin(theta) ;
% generate the uniform mesh of theta_uniform and phi_uniform
theta_uniform = linspace(min(theta(:)), max(theta(:)), 200);
phi_uniform = linspace(min(phi(:)), max(phi(:)), 200);
[THETA, PHI] = meshgrid(theta_uniform, phi_uniform);
% 2D interpolation
F = scatteredInterpolant(theta(:), phi(:), integrand(:));
P_interp = F(THETA, PHI);
% 3. trapz integration of uniform mesh
result2 = trapz(phi_uniform, trapz(theta_uniform, P_interp, 1), 2)

[Mathieu NOE]

hello again

problem solved ?

[Mathieu NOE]

hello again

problem solved ?