function otf = sotf_diffraction_model(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]
   
   otf         = otf_opts;
   
   %% size 
   % precomputed function, has adequate resolution
   n_psi       = 200;
   n_phi       = 200;
   n_k         = 50;
   phi_lim     = [-20 20];
   psi_lim     = [-20 20];
   
   %% grid
   phi         = linspace(phi_lim(1), phi_lim(2), n_phi);
   psi         = linspace(psi_lim(1), psi_lim(2), n_psi);
   [phi_mat, psi_mat] = ndgrid(phi, psi);

   %% Lommel functions
   % equations 8.19
   Un       = zeros(n_phi, n_psi, 2);

   for n = 1 : 2
      for k = 0 : n_k-1
         Un(:,:,n)    = Un(:,:,n) ...
                        + (-1)^k * besselj(n+2*k, psi_mat) .* (phi_mat./psi_mat).^(n+2*k);
      end
   end

   %% point spread function
   % equation 3.18
   I                 = (Un(:,:,1).^2 + Un(:,:,2).^2) .* (2./phi_mat).^2;
   [dIdpsi, dIdphi]  = gradient(I, psi, phi);   % note switch due to silly behaviour of gradient
   gradI             = cat(3, dIdphi, dIdpsi);
   
   %% store
   otf.I             = I;
   otf.gradI         = gradI;
   otf.psi           = psi; % dimensionless radial coordinate
   otf.phi           = phi; % dimensionless axial coordinate
end