Time-Delay beamforming delay computation

Question

I'm trying to implement a beamforming algorithm in python and I'm stuck with the time-delay computation for a rectangular array.

I'm using the following equation to compute the delays based on the desired steering angle : $$d_{i,j} = \frac{1}{c}\sin(\theta)(\cos(\phi)x_{i,j} + \sin(\phi)y_{i,j})$$

Where $d_{i,j}$ is the delay for the element $(i,j)$, $(\theta,\phi)$ is the azimuth/elevation steering angle, and $x_{i,j},y_{i,j}$ are the coordinates of the element $(i,j)$.

I am using Matlab's PhasedArray toolbox as a ground truth to verify my work, but I'm getting different results.

Here is the Matlab code I use to get the time delays :

array = phased.URA('Size', [4,4],...
                   'ElementSpacing', 0.042);
delay = phased.ElementDelay('SensorArray',array, 'PropagationSpeed', 340.0);
tau = delay([-30;-20])*fs

And here the Python equivalent :

ang_dft = np.array([-30, -20], float)

myArray = pa.URA(4, 4, 0.042, 0.042)
TDBeamformer = bf.TimeDelayBeamformer(myArray, ang_dft, 340.0, 8000.0)

myArray.plot_array()
TDBeamformer.compute_delays()

print(TDBeamformer.sec2spl(TDBeamformer.delay))

plt.show()

Where the compute_delays method is the following :

def compute_delays(self):
    for x in range(0, self.array.sizeX):
        for y in range(0, self.array.sizeY):
            self.delay[y, x] = (1/self.c *
                                np.sin(np.pi / 180. * self.steeringAngle[0]) *
                                (np.cos(np.pi / 180. * self.steeringAngle[1]) * self.array.coordinatesX[y, x] +
                                 np.sin(np.pi / 180. * self.steeringAngle[1]) * self.array.coordinatesY[y, x]))

Could you help me pinpoint where the problem lies? Is my equation wrong (I found it from this website)? Is it a geometry problem in the code?