Inconsistent behaviour in AnalyticBeam.interp()

[ronniyjoseph]

Hi folks, I was playing around with AnalyticBeam and noticed a difference between how interp() deals with the az_array and za arraysfor Gaussian beams vs the Airy and uniform. The 'uniform' and 'airy beam' are pretty much agnostic to the shapes of the input az_array and input za_array, but when using the 'gaussian' beam these need to be flattened before passing to .interp(), caused by

https://github.com/RadioAstronomySoftwareGroup/pyuvsim/blob/1ef8cdbfd98fc6f49b598f75d759852407601ebf/pyuvsim/analyticbeam.py#L151

Minimum working example:

import pyuvsim
import numpy as np

#define some interp range
phi = np.linspace(0,np.pi,5)
theta = np.linspace(0,np.pi/2,6)
freqs = np.linspace(135,165,7)*1e6

#create 2D grid for 
pp,tt = np.meshgrid(phi,theta)

#uniform beam test
uniform = pyuvsim.AnalyticBeam(type='uniform')
uniform.efield_to_power()
uniform_beam = airy_disk.interp(pp, tt, freqs)

#airy beam test
airy_disk = pyuvsim.AnalyticBeam(type='airy', diameter=4.5)
airy_disk.efield_to_power()
airy_beam = airy_disk.interp(pp, tt, freqs)

Works fine, but with

# gaussian beam test
gaussian = pyuvsim.AnalyticBeam(type='gaussian', diameter=6)
gaussian.efield_to_power()
gauss_beam = gaussian.interp(pp,tt,freqs)

You get a classic shape mismatch error between the za_array and frequency_array in line 151.

gaussian = pyuvsim.AnalyticBeam(type='gaussian', diameter=6)
gaussian.efield_to_power()
gauss_beam = gaussian.interp(pp.flatten(),tt.flatten(),freqs)

Flattening like the above makes my problem go away (which also seems more in-line with the pyuvsim documentation). It seems that numpy.meshgrid actually flattens nd-arrays (not really in the numpy documentation) hence why the non-flattened case works for the Airy beam.

I'll gladly accept that I simply shouldn't be passing multi-dimensional arrays to .interp(). An alternative would be a np.meshgrid implementation in the gaussian interpolation, and some flattening afterwards?

[bhazelton]

@ronniyjoseph Thanks for pointing this out. I don't think we want a meshgrid implementation because most of the time our sources do not occur on a mesh 😆. I think we always want to have a flat array coming in, so we'll add some checks to catch when that's not happening and give a better error message.