Problem with the number of phase fronts when attempting to generate OAM

Problem with the number of phase fronts when attempting to generate OAM:

In brief, an OAM (Orbital Angular Momentum) beam with a specific mode l should present exactly l number of helical phase fronts in its radiation pattern:

(a) Propagating plane waves, (b) Propagating OAM waves of mode 3, (c) Wave front of plane waves and (d) Wave front of OAM waves of mode 3

To generate the OAM, we can use a uniform circular array. Corresponding to a number of equal finite azimuthal intervals, each antenna element is placed on a circle with equiangular spacing. All the antenna elements are fed by the same signal, but with an incremental phase shift of  , such that the  antenna element will have a phase of where N is the total number of the antenna elements and l is the OAM mode index.

So I try to implement the exact same logic in EMPro: 

My geometry:

(The code used to generate the geometry is based on the phased array antenna example in EMPro. )

When plotting the result, the excitation of the ports are set as follows:

# Ntx: the number of port

# l: OAM mode, the number of expected phase fronts

 exc = { p.name: 1/math.sqrt(Ntx)*complex(math.cos(2*math.pi*i*l/Ntx),math.sin(2*math.pi*i*l/Ntx)) for i, p in enumerate(empro.activeProject.circuitComponents()) if p.port }

 # Load the far field result

farFieldResult = results.farField(excitation=exc, sensorName=SENSOR_NAME)

This way, for mode 1 and Ntx = 10, the excitation profile looks like this: (which is correct)

However, when I export the far-field data and plot it in MATLAB, the number of phase fronts is always doubled, as shown in the three following figures (for mode 1-3 respectively), while radiation magnitude pattern is correct (the angle of the main lobe matches the theoretical value)

Mode 1

Mode 2

Mode 3

Here is the code I used to plot these figures:

emproData  = zeros(height(FarField),4);

emproData(:,1) = FarField.RealEtheta + i*FarField.ImagEtheta;

emproData(:,2) = FarField.RealEphi + i*FarField.ImagEphi;

emproData(:,3) = rad2deg(FarField.theta);

emproData(:,4) = rad2deg(FarField.phi);

pat = sqrt(power(emproData(:,1),2) + power(emproData(:,2),2));

theta = emproData(:,3);

phi = emproData(:,4);

subplot(1,2,1)

patternCustom(abs(pat),theta,phi);

title(strcat('Efield Magnitude Pattern',{' '},str));

view(90,0);

subplot(1,2,2)

patternCustom(unwrap(angle(pat)),theta,phi) ; % in radians

title(strcat('Efield Phase Pattern in Radians', {' '},str, {' '}, 'bottom view'));

It would be rather helpful if you could help me find what is wrong in my setup!