Antenna Directivity Calculation by field data

Hi, everyone. I have a basic question about using MATLAB program to calculate the antenna directivity by using the field data.

The field data are stored in a N-by-M array of magnitude of the E field, N is the phi angle from 0 degree to 360 degree and M is the theta angle from 0 degree to 180 degree.

My calculation is basically converting the field magnitude to radiation intensity, then summing the radiation intensity to total radiated power, then dividing the radiation intensity to the average radiated power for directivity.

However, I always found some discrepancy between my calculations for isotropic source and other software results. Would anybody look at my program and give me some hints? Thanks in advance.

The MATLAB programs are shown below:

%%%%%%%%%%%%%%%%%%%%%%%%%%
function [RP_DirDB]=RP_Dir(RP_Mag)
% This function converts a bunch of field magitude to the power directivity

% Input:
% RP_Mag: 1-by-3 cells, each cell contains a radiation pattern
% (total, theta, phi-components) represented by field magitude.
% Output:
% RP_DirDB: 1-by-3 cells, each cell contains a radiation pattern
% (total, theta, phi-components) represented by power
% directivity.


[P,S]=TotalP(RP_Mag{1});
RP_DirDB{1}=EMag2DirDB(RP_Mag{1},P);

% The following commands are used for field sub-components
%RP_DirDB{2}=EMag2DirDB(RP_Mag{2},P);
%RP_DirDB{3}=EMag2DirDB(RP_Mag{3},P);

function [P,S]=TotalP(Mag)
% This function gets the normalized total radiated power from the total magitute of E
% fields

% Input:
% Mag: N_phi-by-N_theta matrix, contains the radiation pattern in
% absolute value.
% Output:
% P: Normalized total radiated power
% S: Integral surface area

% get the radiation intensity
Unor=abs(Mag).^2;

% Get the phi and theta data
[Nphi,Ntheta]=size(Mag);
Phi=linspace(0,2*pi,Nphi);
Dphi=2*pi/(Nphi-1);
Theta=linspace(0,pi,Ntheta);
Dtheta=pi/(Ntheta-1);

% Initialize the total radiated power P, integral area S, delta power dP
% and detal integral area dS.
P=0;
S=0;
dP=[];
dS=[];

% Calculate the total radiation power and radiation intensity
dS(1,1)=2*pi*(1-cos(Dtheta/2));
dP(1,1)=Unor(1,1)*dS(1,1);
P=P+dP(1,1);
S=S+dS(1,1);

dS(1,Ntheta)=2*pi*(1-cos(Dtheta/2));
dP(1,Ntheta)=Unor(1,Ntheta)*dS(1,Ntheta);
P=P+dP(1,Ntheta);
S=S+dS(1,Ntheta);

for i_theta=2:(Ntheta-1)
for i_phi=1:(Nphi-1)
dS(i_phi,i_theta)=Dphi*(cos(Theta(i_theta)-Dtheta/2)-cos(Theta(i_theta)+Dtheta/2));
dP(i_phi,i_theta)=Unor(i_phi,i_theta)*dS(i_phi,i_t heta);
P=P+dP(i_phi,i_theta);
S=S+dS(i_phi,i_theta);
end
end

function [DirDB]=EMag2DirDB(Mag,P)
% This function converts the absolute magnitude number to decibel number

% Input:
% Mag: N_phi-by-N_theta matrix, contains the radiation pattern in
% absolute value.
% Output:
% DirDB: N_phi-by-N_theta matrix, contains the radiation pattern in
% decibel value


% get the radiation intensity
Unor=abs(Mag).^2;

% Renormalize the radiation intensity to directivity
Dir=4*pi*Unor/P;

% Calculate the directivity in decibel
DirDB=10*log10(Dir);