pro psf_1_gauss_func, arg, param, z, deriv

    x = arg.indx mod arg.nx
    y = arg.indx / arg.nx
    p = psf_1_gauss_param2struct(param)
    
    n = n_elements(x)

    z = dblarr(n)
    if n_params() ge 4 then begin
        deriv = dblarr(n,n_elements(param))
    end

    
    ; calculate intermediate variables and the function
    x2 = (x-p.xctr)^2.d
    y2 = (y-p.yctr)^2.d
    r2 = x2 + y2
    w2 = (p.fwhm)^2.d / ( 4. * alog(2.) )
    w = sqrt(w2)
    
    z_ctr1 = exp( - ( (r2/w2) < 100.d) )  ; unity at center
    z_tot1 = z_ctr1 / (w2 * !dpi)         ; for unit total
    z = p.total * z_tot1
    
    if  n_params() ge 4 then begin
        deriv[*,0] = p.total * z_tot1 * 2.d * (x-p.xctr)/w2
        deriv[*,1] = p.total * z_tot1 * 2.d * (y-p.yctr)/w2
        deriv[*,2] = z_tot1
        deriv[*,3] = 2.d * (p.total * z_tot1/p.fwhm) * ( (r2/w2) - 1.d )
    endif
end