Make PSF work with nd-input

[maxnoe]

Trying to make #2643 work with ND input.

Using this script to plot:

import astropy.units as u
import numpy as np
from ctapipe.instrument.optics import ComaModel
from ctapipe.instrument import SubarrayDescription
import matplotlib.pyplot as plt

subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst")

lst1 = subarray.select_subarray([1])
psf_model = ComaModel(
    subarray=lst1,
    asymmetry_params=[0.5, 10, 0.15],
    radial_scale_params=[0.015, -0.1, 0.06, 0.03],
    az_scale_params=[0.25, 7.5, 0.02],
)

edges_x = np.linspace(-0.1, 0.1, 250) * u.m
edges_y = np.linspace(-0.1, 0.1, 250) * u.m
centers_x = 0.5 * (edges_x[:-1] + edges_x[1:])
centers_y = 0.5 * (edges_y[:-1] + edges_y[1:])
x, y = np.meshgrid(centers_x, centers_y)

psf_center = psf_model.pdf(x, y, 0.0 * u.m, 0.0 * u.m)
psf_border = psf_model.pdf(x + 1 * u.m, y + 1 * u.m, 1 * u.m, 1 * u.m)

fig, (ax1, ax2) = plt.subplots(1, 2, layout="constrained")

ax1.pcolormesh(
    edges_x.to_value(u.m),
    edges_y.to_value(u.m),
    psf_center,
    cmap='inferno',
)

ax2.pcolormesh(
    edges_x.to_value(u.m),
    edges_y.to_value(u.m),
    psf_border,
    cmap='inferno',
)

ax1.set(
    aspect=1,
    title="PSF at (0 m, 0 m)",
)
ax2.set(
    aspect=1,
    title="PSF at (1 m, 1 m)",
)

plt.show()

This is the current result on this branch (The PSF at (1, 1) all evaluates to nan_:

[maxnoe]

I checked, the result is the same in #2643 when using loops to evaluate single points, so the vectorization here is correct and the results are weird.

[maxnoe]

@mexanick I pushed some changes to notation and fixed the order of paramters to polyval, with those changes, I now get this for the PSF plots:

Still not what I would expect... but no nans is an improvement I'd say