2D Peak

[Mark Dean]

We fairly commonly need to fit 2D peaks, such a Lorentzian or Lorentzian squared or Gaussian peaks.

This can be done in lmfit, but my experience is that people often find it pretty hard. 

Would you be interested in me contributing an example for this here:

https://lmfit.github.io/lmfit-py/examples/index.html

If so, I'd further prefer to add a function and model with form similar to that copied below.  Would a PR on this topic be welcomed? 

def twodpeak(x, y, height=1., centerx=0., centery=0., sigmax=1., sigmay=1.,

             rotation=0., power=1., form='gaussian'):

    """Return a two dimensional gaussian or lorentzian raised to a power.

    ``form`` sets the functional form

    form='gaussian' (default)

        =height*exp(-R**2)**power

    form='lorentzian' (default)

        =height*(1/(1+R**2))**power

    Where radius, ``R`` is defined in terms of ``centerx`` and ``centery``

    with widths ``sigmax`` and ``sigmay``. The peak can be rotated by choosing

    the value of ``rotation`` in radians.

    xp = (x - centerx)*cos(rotation) - (y - centery)*sin(rotation)

    yp = (x - centerx)*sin(rotation) + (y - centery)*cos(rotation)

    R = (xp/sigmax)**2 + (yp/sigmay)**2

    Intensity is defined by ``height`` i.e. the value at ``centerx``=``centery``=0

    """

    xp = (x - centerx)*cos(rotation) - (y - centery)*sin(rotation)

    yp = (x - centerx)*sin(rotation) + (y - centery)*cos(rotation)

    R = (xp/sigmax)**2 + (yp/sigmay)**2

    if form.lower() == 'gaussian':

        return height*exp(-R**2)**power

    elif form.lower() == 'lorentzian':

        return height*(1/(1+R**2))**power

    else:

        msg = ("Invalid value '{}' for argument 'form'".format(form)

               + "; should be 'gaussian' or 'lorentzian'")

        raise ValueError(msg)

[Matt Newville]

Hi Mark, 

I think it's a great idea to add a fit to a 2D peak as an example. Please do!

I have to admit a bias against labeling randomly-defined peak shapes as 'Gaussian', 'Lorentzian', etc and a bias toward using unit-normalized peak-shapes with consistent naming conventions.  

For examples working with 2D arrays, I would also suggest that this would be a fine opportunity to showcase the use of `numpy.meshgrid` and/or related methods.    That is, a 2D Gaussian for a 2D array with 'x' and 'y' units other than "pixel" can be done nicely as

    import numpy as np

    from lmfit.lineshapes import gaussian 

    xx, yy = meshgrid(np.arange(512)*0.075, np.arange(1024)*0.075)

    gauss2d = gaussian(xx, center=15.7, sigma=4.5)*gaussian(yy, center=52.1, sigma=6)

Once you can do that, model data as a 2d gaussian seems much easier, well at least to me ;).

Does that seem reasonable to you?

--Matt

 

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/5e5d2c8e-cd24-4405-a50e-28c4d4b7df85%40googlegroups.com.

--

--Matt Newville <newville at cars.uchicago.edu> 630-252-0431