;  param = psf_1_gauss_initparam(d, b=b,x=x,y=y,f=f,w=w,rad=rad)

function psf_1_gauss_initparam, d, x, y, rad 

    b = 0.
    f = 1.
    w = 2.

    nxnyn,d,nx,ny
    if n_elements(rad) eq 0 then rad = 10.

    x = 0 > x < (nx-1)
    x0 = 0>(round(x)-rad)
    x1 = (round(x)+rad) < (nx-1)
    y = 0 > y < (ny-1)
    y0 = 0>(round(y)-rad)
    y1 = (round(y)+rad) < (ny-1)
    
    ; centroid on the subarray to get total
    centrod, d, x, y, rad, 0, 0, b, xbar, ybar, f

    w2 = f/ (!pi * d[round(x),round(y)] )
    if w2 gt 0 then begin
        w = sqrt( w2 )          ;  exp( - (x^2+y^2)/w^2 )
    endif else begin
        w = rad/2.
    endelse
    fwhm = 2. * sqrt(alog(2.)) * w

    return, [xbar,ybar,f,fwhm]

end