function [pimg, jac] = pimg_model(param_vec, PSF, irel_aa, jrel_aa, aa_mtx)
	% [pimg, jac] = pimg_model(param_vec, irel, jrel, PSF, aa_mtx)
	% 
	% evaluate particle image model for least squares fitting
	%
	% Input:
	%	param_vec		fit parameters (image space coordinates)
	%	param_vec(1)	radius
	%	param_vec(2)	i coordinate
	%	param_vec(3)	j coordinate
	%	param_vec(4)	black level
	%	param_vec(5)	white level
	%
	%	PSF				point spread function kernel
	%	irel_aa			i coordinates of pixels in upsampled image
	%	jrel_aa			j coordinates of pixels in upsampled image
	%	aa_mtx			matrix representing downsampling operation 

	%% extract parameters 
	R			= param_vec(1);
	im_ctr		= param_vec(2:3);
	black		= param_vec(4);
	white		= param_vec(5);
	ksize		= size(PSF,1);
	assert(ismatrix(PSF) && size(PSF,2) == ksize, 'PSF kernel must be square');
	assert(mod(ksize,2) == 0, 'PSF kernel size must be odd');
	
	%% evaluate sharp image
	r2mat		= (irel_aa(:) - im_ctr(1)).^2 + (jrel_aa(:).' - im_ctr(2)).^2;
	pimg_sharp	= r2mat >= R^2;
	pimg_sharp  = (white-black)*pimg_sharp + black;

	%% convolve with PSF
	pimg_pad	= padarray(pimg_sharp, floor(ksize/2*[1 1]), 'replicate', 'both');
	pimg_conv	= conv2(pimg_pad, PSF, 'valid');
	pimg		= aa_mtx * pimg_conv(:);
end