code: def function(self, x, y, amp, alpha, beta, center_x, center_y):
          
          x_shift = x - center_x
          y_shift = y - center_y
          return amp * (1. + (x_shift**2+y_shift**2)/alpha**2)**(-beta)
docs: returns Moffat profile