function [pimg_model, pimg_jac] = pimg_fit_modelfun(X, psf)
    % particle image model
    %% extract parameters of model
    centroid        = X(1:2);
    radius          = X(3);
    black           = X(4);
    white           = X(5);
    
    %% create synthetic image
    ksize               = size(psf.kernel);
    rmat		        = sqrt((psf.irel_aa(:) - centroid(1)).^2 + (psf.jrel_aa(:).' - centroid(2)).^2);
    [pimg_sharp, dPdr]	= pimg_step(rmat, radius, 1);
	pimg_sharp_pad      = padarray(pimg_sharp, floor(ksize/2), 'replicate', 'both');
	% convolve with kernel
    pimg_conv_up        = convn(pimg_sharp_pad, psf.kernel, 'valid');
	pimg_model   	    = psf_downsample(pimg_conv_up, psf.aa_factor);
    pimg_model          = (white-black) * pimg_model + black;
    %pimg_model=pimg_model(:);
    if nargout == 1
        return;
    end

    %% jacobian, if requested
    drdi                = (centroid(1) - psf.irel_aa(:)  ) .* rmat;
    drdj                = (centroid(2) - psf.jrel_aa(:).') .* rmat;
    dPdi                = (white-black) * dPdr .* drdi;
    dPdj                = (white-black) * dPdr .* drdj;
    dPdR                =-(white-black) * dPdr;
    % convolve with kernel
    dPdi_pad            = padarray(dPdi, floor(ksize/2), 0, 'both');
    dPdj_pad            = padarray(dPdj, floor(ksize/2), 0, 'both');
    dPdR_pad            = padarray(dPdR, floor(ksize/2), 0, 'both');
    dPdi_conv_up        = convn(dPdi_pad, psf.kernel, 'valid');
    dPdj_conv_up        = convn(dPdj_pad, psf.kernel, 'valid');
    dPdR_conv_up        = convn(dPdR_pad, psf.kernel, 'valid');
    dPdi_conv           = psf_downsample(dPdi_conv_up, psf.aa_factor);
    dPdj_conv           = psf_downsample(dPdj_conv_up, psf.aa_factor);
    dPdR_conv           = psf_downsample(dPdR_conv_up, psf.aa_factor);

    n_data              = numel(pimg_model);
    pimg_jac            = zeros(n_data, 5);
    pimg_jac(:,1)       = dPdi_conv(:); % dP/di0
    pimg_jac(:,2)       = dPdj_conv(:); % dP/dj0
    pimg_jac(:,3)       = dPdR_conv(:); % dP/dR
    pimg_jac(:,4)       = ones(n_data,1); % dP/dB
    pimg_jac(:,5)       = (pimg_model(:) - black) / (white-black); % dP/dW
end