function moffat2d, moff_dim, moffat_params, div=div
width=moff_dim/2

if NOT(keyword_set(div)) then div=1
div=fix(div)               ; div must be an integer

;m=2*width+1

xxn=(findgen(moff_dim) - moff_dim/2)/div  ; create x (first introduction
yyn=xxn                    ; create y
yyn(*)=1.0
x=xxn#yyn                  ; transform into matrices, x: left to right
y=transpose(x)             ; y: top to bottom


;xxn1=(findgen(m*div)-(m*div)/2)/div
;xxn1=m*(xxn-min(xxn))/(max(xxn)-min(xxn))-width - div/2

;stop

q=size(moffat_params)
a = moffat_params(0) & b = moffat_params(1) & c = moffat_params(2) & xc = moffat_params(3)
yc = moffat_params(4) & epsx = moffat_params(5) & epsy = moffat_params(6)        ; rename parameters
bg=moffat_params(7)

a2=a^2 & xa=x-xc & ya=y-yc
ra2=xa^2+ya^2 & eps=sqrt(epsx^2+epsy^2)         ; convenient substitutions

r2=xa^2*(1.+eps-epsx)+ya^2*(1.+eps+epsx)-2.*xa*ya*epsy  ; r^2 for rotated ellipses

;if (a2 eq 0.0) then begin                       ; safety checks
;     print, 'a^2 equal to zero!'
;     a2=a2+0.001
;endif
;if (min(r2) lt 0.0) then stop, 'Negative r^2!'

z = (1.+r2/a2)

moff=c*z^(-b)                   ; Evaluate the PSF core function

return, moff
end

;***********************************************