pythran compiled by g++ produces wrong results, OK with clang

_array.py

import numpy as np

#pythran export getArrayFactor (complex[:], float[:], float[:], float[:,:], float[:,:], float)
def getArrayFactor(W, px, py, theta, phi, L):
    N = px.shape[0]
    AF = np.zeros(theta.shape, dtype=complex)
    kx = np.sin(theta) * np.cos(phi) * 2.0 * np.pi / L
    ky = np.sin(theta) * np.sin(phi) * 2.0 * np.pi / L
    for n in range(N):
        AF += W[n] * np.exp(-1j * (kx * px[n] + ky * py[n]))
    magAF = np.abs(AF)
    return magAF

#pythran export latticeLayout (float, float, float, int, int, str)
def latticeLayout(L, DX, DY, Nx, Ny, CENTER_FLAG):
    dx = DY * L
    dy = DX * L

    # px = np.zeros((Nx * Ny, 1))
    # py = np.zeros((Nx * Ny, 1))
    px = np.zeros((Nx * Ny))
    py = np.zeros((Nx * Ny))
    tmpi = 0

    for ix in range(Nx):
        if ix % 2 == 0:
            offset = 0
        else:
            offset = 0.5 * dx

        for iy in range(Ny):
            px[tmpi] = ix * dx
            py[tmpi] = iy * dy + offset
            tmpi += 1

    if CENTER_FLAG == 'CENTER':
        mean_px = np.mean(px)
        mean_py = np.mean(py)
        px = px - mean_px
        py = py - mean_py

    return px, py

array.py

import numpy as np
from scipy import constants

F_des = 20.2 * 1e9
C = constants.speed_of_light
L = C / F_des  # wavelength in meters
patchR = 0.015 / 4
patchH = 0.000422
er = 3.66

#-------------------------- COORDINATE SPACE -----------------------------

# Uniform mesh in theta-phi coordinate space
theta_start = 0
theta_step = np.pi / 1000
theta_end = np.pi - theta_step
phi_start = -np.pi
phi_step = np.pi / 1000
phi_end = np.pi - phi_step
theta, phi = np.meshgrid(np.arange(theta_start, theta_end, theta_step),
                         np.arange(phi_start, phi_end, phi_step))

# equivalent non-uniform mesh in uv space
uGrid = np.sin(theta) * np.cos(phi)
vGrid = np.sin(theta) * np.sin(phi)

from patch import circular_patch
magEF =  circular_patch(patchR, patchH, L, er, theta, phi);

# array geometry
dx = 0.5  # 'dx (wavelengths)', '0.015';...
dy = 0.5  # 'dy (wavelengths)', '1.7321';...
Nx = 16  # 'Nx', '4';...
Ny = 16  # 'Ny', '4';...
theta_doa = np.pi / 3
phi_doa = np.pi / 6

# array geometry
def latticeLayout(L, DX, DY, Nx, Ny, CENTER_FLAG):
    dx = DY * L
    dy = DX * L

    # px = np.zeros((Nx * Ny, 1))
    # py = np.zeros((Nx * Ny, 1))
    px = np.zeros((Nx * Ny))
    py = np.zeros((Nx * Ny))
    tmpi = 0

    for ix in range(Nx):
        if ix % 2 == 0:
            offset = 0
        else:
            offset = 0.5 * dx

        for iy in range(Ny):
            px[tmpi] = ix * dx
            py[tmpi] = iy * dy + offset
            tmpi += 1

    if CENTER_FLAG == 'CENTER':
        mean_px = np.mean(px)
        mean_py = np.mean(py)
        px = px - mean_px
        py = py - mean_py

    return px, py

from _array import latticeLayout

def getArrayFactor(W, px, py, theta, phi, L):
    N = px.shape[0]
    AF = np.zeros(theta.shape, dtype=complex)
    kx = np.sin(theta) * np.cos(phi) * 2 * np.pi / L
    ky = np.sin(theta) * np.sin(phi) * 2 * np.pi / L
    for n in range(N):
        AF += W[n] * np.exp(-1j * (kx * px[n] + ky * py[n]))
    magAF = np.abs(AF)
    return magAF

#from _array import getArrayFactor
# -----------------------------ARRAY LAYOUT-------------------------------
px, py = latticeLayout(L, dx, dy, Nx, Ny, 'CENTER')

# beam weight vector
N = px.shape[0]
kx = np.sin(theta_doa) * np.cos(phi_doa) * 2 * np.pi / L
ky = np.sin(theta_doa) * np.sin(phi_doa) * 2 * np.pi / L
W = np.conj(np.exp(-1j * (kx * px + ky * py)))
W /= np.linalg.norm(W)

# -----COMPUTE ARRAY FACTOR
magAF = getArrayFactor(W, px, py, theta, phi, L)
from _array import getArrayFactor as _getArrayFactor
magAF2 = _getArrayFactor (W, px, py, theta, phi, L)
assert np.allclose (magAF, magAF2)

# # Include effect of element factor into array factor
# magAF = magAF * magEF

# import matplotlib.pyplot as plt
# plt.figure(3)
# plt.grid()
# #plt.hold()
# plt.contour(theta, phi, magAF)
# plt.show()

patch.py

# Constants
C = 299792458  # speed of light in m/s
F_des = 20.2  # design frequency in GHz
patchR = 0.015  # circular patch radius in meters
patchH = 0.000422  # circular patch height in meters
er = 3.66  # dielectric constant
PATCH = 'CIRC'

# Calculate wavelength in meters
L = C / (F_des * 1e9)

# Declare rectangular patch width and length in wavelengths
# (these variables are not defined in the original code snippet)
pW = None
pL = None

import numpy as np

# Constants
theta_start = 0
theta_step = np.pi/1000
theta_end = np.pi-theta_step
phi_start = -np.pi
phi_step = np.pi/1000
phi_end = np.pi-phi_step

# Uniform mesh in theta-phi coordinate space
theta, phi = np.meshgrid(np.arange(theta_start, theta_end+theta_step, theta_step),
                         np.arange(phi_start, phi_end+phi_step, phi_step))

# equivalent non-uniform mesh in uv space
uGrid = np.sin(theta) * np.cos(phi)
vGrid = np.sin(theta) * np.sin(phi)

import numpy as np
import scipy.special as sp

def rect_patch(pWf, pLf, L_des, L_op, er, theta, phi):
    # Patch dimensions are based on the design freq/wavelength
    pW = pWf * L_des / np.sqrt(er)  # patch width
    pL = pLf * L_des / np.sqrt(er)  # patch Length

    # Array is exposed to radiation of operational freq/wavelength
    K = 2 * np.pi / L_op

    # Response calculations
    tmp1 = np.sinc(0.5 * K * pW * np.sin(theta) * np.sin(phi) / np.pi)
    tmp2 = np.cos(0.5 * K * pL * np.sin(theta) * np.cos(phi))

    E_theta = tmp1 * tmp2 * np.cos(phi)
    E_phi = -tmp1 * tmp2 * np.cos(theta) * np.sin(phi)
    magEF = np.sqrt(E_theta * E_theta + E_phi * E_phi)

    # Assuming a infinite XY groundplane
    win = theta <= np.pi / 2
    magEF = magEF * win

    return magEF

def circular_patch(a, h, L_op, er, theta, phi):
    # effective radius
    ae = a * np.sqrt(1 + (2 * h / (np.pi * a * er)) * (np.log(0.5 * np.pi * a / h) + 1.7726))

    k0 = 2 * np.pi / L_op

    sin_theta = np.sin(theta)

    bfn0 = sp.jn(0, k0 * ae * sin_theta)
    bfn2 = sp.jn(2, k0 * ae * sin_theta)
    Jp_02 = bfn0 - bfn2
    J_02 = bfn0 + bfn2

    # ((k0 * ae * exp(-1i*k0*r))/(2*r)) is a constant wrt theta & phi , so is
    # not included below.
    E_theta = np.abs(-1j * (np.cos(phi) * Jp_02))

    E_phi = np.abs(1j * (np.cos(theta) * np.sin(phi) * J_02))

    magEF = np.sqrt(E_theta * E_theta + E_phi * E_phi)

    return magEF

# Include element factor if flag is set
if PATCH == 'CIRC':
    magEF = circular_patch(patchR, patchH, L, er, theta, phi)
elif PATCH == 'SQUARE':
    magEF = rect_patch(pW, pL, L, L, er, theta, phi)