function otf = otf_diffraction_model_interpolant(otf_opts)
	% OTF structure contains
	%  resolution  resolution for integrating profile [pixels]
	%  Z_0         sensor plane position[mm]
	%  f           focal length         [mm]
	%  D           aperture diameter    [mm]
	%  pitch       pixel ptich          [um]
	%  lambda      wavelength           [um]
	%  phi         dimensionless axial distance
	%  psi         dimensionless radial coordinate
	%  I           diffraction pattern I(phi,psi)
	%  gradI       gradient of I with respect to [phi,psi]
	%  support     support function

	otf         = otf_opts;

	%% size 
	% precomputed function, has adequate resolution
	n_phi       = 100;
	n_k         = 50;
	phi_lim     = [-20 20];

	%% grid for evaluation of OTF
	% before integration over pixel
	sz			= size(otf_opts.support);
	irel_pad	= -sz(1)/2-0.5 : 0.1 : sz(1)/2+0.5;
	jrel_pad	= -sz(2)/2-0.5 : 0.1 : sz(2)/2+0.5;
	n_pad		= [length(irel_pad) length(jrel_pad)];
	% integration kernel
	kernel		= trapz(-0.5 : 0.1 : 0.5, eye(11), 2);
	irel		= conv(irel_pad, kernel, 'valid');
	jrel		= conv(jrel_pad, kernel, 'valid');
	
	%% grid
	phi			= linspace(phi_lim(1), phi_lim(2), n_phi);
	[imat, jmat, phimat] = ndgrid(irel_pad, jrel_pad, phi);
	rmat		= sqrt(imat.^2 + jmat.^2);
	psimat		= 2*pi*rmat*otf.pitch./otf_opts.lambda * otf_opts.D/(otf_opts.f*2);
	
	%% Lommel functions
	% equations 8.19
	Un       = zeros([n_pad n_phi 2]);

	for n = 1 : 2
		for k = 0 : n_k-1
			Un(:,:,:,n)    = Un(:,:,:,n) ...
				+ (-1)^k * besselj(n+2*k, psimat) .* (phimat./psimat).^(n+2*k);
		end
	end
	% fix bessel function at i=j=0
	% besselj(m, 0) / x^m = 2^-m / gamma(1+m)
	k		= 0:n_k-1;
	Un0		= zeros(2,1,n_phi);
	i0		= find(irel_pad==0);
	j0		= find(jrel_pad==0);
	for n = 1 : 2
		m			= n+2*k;
		Un0(n,1,:)	= sum(phimat(i0,j0,:).^m .* (-1).^k .* 2.^-m ./ gamma(1+m), 2);
	end
	Un(i0,j0,:,1) = Un0(1,1,:);
	Un(i0,j0,:,2) = Un0(2,1,:);

   %% point spread function
   % equation 3.18
   I_pad			= (Un(:,:,:,1).^2 + Un(:,:,:,2).^2) .* (2./phimat).^2;
   % convolve to integrate over the area of one pixel
   I				= convn(convn(I_pad, kernel, 'valid'), kernel.', 'valid');
   [dIdj, dIdi]		= gradient(I, jrel, irel, phi);   % note switch due to silly behaviour of gradient
   gradI            = cat(4, dIdi, dIdj);
   
   %% store
   otf.I            = I;
   otf.gradI        = gradI;
   otf.irel         = irel; % relative i coordinate coordinate
   otf.jrel			= jrel; % relative j coordinate
   otf.phi			= phi; % dimensionless axial coordinate
end