How to produce lorentzian 2D sources in python

I am trying to produce 2D sources and add them to an image. Currently I'm generating gaussian sources through Gaussian2D (astropy):

min_x = int(rapix[i]) - 300
max_x = int(rapix[i]) + 300
min_y = int(decpix[i]) - 300
max_y = int(decpix[i]) + 300
y, x = np.mgrid[min_y:max_y, min_x:max_x] #HERE I CREATE THE GRID TO PUT THE SOURCES ON

fakesource = Gaussian2D(intensity, rapix[i], decpix[i], dimension, dimension)(x, y)

where the intensity and dimensions have been defined previously. I'd like to produce Lorentzian sources instead of Gaussians, but I did not find anything similar to Gaussian2D. Which is the best way to do it?

Solution

Lorentz1D exists, though no Lorentz2D, but you can define one following the guide for implementing custom models. Here's a basic example:

from astropy.modeling import Fittable2DModel, Parameter


class Lorentz2D(Fittable2DModel):
    amplitude = Parameter()
    x_0 = Parameter()
    y_0 = Parameter()
    fwhm = Parameter

    @staticmethod
    def evaluate(x, y, amplitude, x_0, y_0, fwhm):
        hwhm = fwhm / 2.0
        return (amplitude * hwhm /
                ((x - x_0)**2 + (y - y_0)**2 + hwhm**2)**1.5)

I don't know if this is exactly how your 2D Lorentzian model is defined; I just adapated this definition from Wikipedia. But you can modify this example as-needed.

Then, if you think this would be valuable to others, you might consider submitting it as a contribution.