I'd like to add a cross-phase diagnostic, to plot the phase difference between two variables as a function of position and frequency. Should this be implemented as a standalone function or as a data array method?
i.e. if ds is a dataset:
result = crossphase(ds['Ne'], ds['phi'])
or
result = ds['Ne'].crossphase(ds['phi'])
?