spherical aberration and Chromatic aberration

Asked by Alan AF on 20 Apr 2012

Latest activity Commented on by Eric on 19 Jun 2017

157 views (last 30 days)

How can we simulate spherical aberration in matlab? It might be by using the below code, but what are the parameters ??

PSF = fspecial(?); % create PSF
Blurred = imfilter(I,PSF,?,'conv');

Thank you

0 Comments

4 Answers

Link

Answer by Eric on 20 Apr 2012

Accepted Answer

The following code lets you simulate the PSF associated with spherical aberration (as well as diffraction). You need to specify the PSF sampling pitch, the wavelength, the aperture diameter, the system focal length, the amount of spherical aberration, and the PSF array size.

The idea is to calculate a wavefront with spherical aberration and then use the Fourier transform to calculate the PSF. You need to be careful with sampling, hence the assert() statement.

You can easily include other aberrations as well, if desired. Simply modify how the wavefront is calculated.

You can check that this code is working by setting RMS_SA to zero (i.e., simulate a diffraction-limited system, in which case the PSF is an Airy disk). Then the first null of the PSF should have a radius of 1.22*lambda*focal_length/aperture_diameter.

You'll have to do some work to simulate a spectrally broad system. In that case you could simulate a number of monochromatic PSFs and then sum them together. You would need to know how the spherical aberration varies as a function of wavelength (variation as the reciprocal of wavelength would be a good estimate - i.e., longer wavelengths have less aberration) as well as the spectral output of the source, the spectral transmittance of the optical system, and the spectral responsivity of the detector. The product of these spectral terms would be used as the weights in a weighted average to calculate the broadband PSF from the monochromatic PSFs.

Once you get the PSF you can convolve this with your scene to simulate imaging performance. imfilter() may work, but I like convnfft ( http://www.mathworks.com/matlabcentral/fileexchange/24504-fft-based-convolution) from the File Exchange. You could use conv2() as well.

Good luck,

Eric

%%Define parameters
clear all;close all;
psf_sampling = 0.5e-6;%focal plane sampling in meters
lambda = 0.6328e-6;%wavelength in meters
N = 256;
aperture_diameter = 0.0254;%meters; 1 inch
focal_length = 5*aperture_diameter;%meters
RMS_SA = 0.25;%RMS spherical aberration content in waves
%%Calculate pupil plane sampling
delta_fx = 1/(psf_sampling*N);
x_pupil = (-fix(N/2):fix((N-1)/2)) * delta_fx * lambda * focal_length;
[X_pupil,Y_pupil] = meshgrid(x_pupil);
R_pupil = sqrt(X_pupil.^2 + Y_pupil.^2);
R_norm = R_pupil/(aperture_diameter/2);%Pupil normalized to aperture diameter
assert(max(R_norm(:))>=sqrt(2),'Sampling is not sufficient to reconstruct the entire wavefront.');
%%Create wavefront
W = RMS_SA * sqrt(5) * (6*R_norm.^4 - 6*R_norm.^2 + 1);%Spherical Aberration wavefront
W(R_norm>1) = 0;
E = exp(1i*2*pi*W);%Complex amplitude
E(R_norm>1) = 0;%Impose aperture size
figure;imagesc(angle(E)/(2*pi));colorbar;title('Wavefront Phase (waves)');
%%Create point-spread function
psf = abs(fftshift(fft2(ifftshift(E)))).^2;
psf = psf/sum(psf(:));%Normalize to unity energy
x_psf = (-fix(N/2):fix((N-1)/2)) * psf_sampling;
figure;imagesc(x_psf*1e6,x_psf*1e6,psf);
title(sprintf('PSF with %.4f waves RMS of Spherical Aberration', RMS_SA));
xlabel(sprintf('Position (\\mum)'));
ylabel(sprintf('Position (\\mum)'));

11 Comments

Show 8 older comments

Eric on 23 Apr 2012

I sure don't know how to change things for a TEM. My training is in imaging systems that use photons rather than electrons.

I'm guessing that since it's called spherical aberration the cause is the same, though: a spherical surface is a better approximation to a parabolic surface on-axis than off. So the shape of the point-spread function would be the same. One could potentially measure a point source and tweak the RMS_SA term until the measured data agree reasonably with the simulated data.

For chromatic aberration I'll assume you mean longitudinal chromatic aberration. This is variation in focal length as a function of wavelength or electron energy (as opposed to lateral chromatic aberration, in which case there's a lateral shift of the focused spots as a function of wavelength/energy).

You'll have to take into account the polychromatic considerations I list above. One way to handle this would be to model each monochromatic PSF with a different amount of defocus. You could include a defocus term by adding

W_defocus = RMS_defocus * sqrt(3)*(2*R_norm.^2 - 1);

W_defocus(R_norm>1) = 0;

and then sum W and W_defocus. In this code RMS_defocus is the RMS magnitude of the amount of defocus to add in. You could vary this as a function of wavelength to induce a wavelength-/energy-dependent blur. In an optical system I believe RMS_defocus generally varies parabolically with wavelength and you need to have an accurate model of the system to know this behavior.

Note that you would not want to vary the focal length as a function of wavelength/energy in my code. That would simulate making an observation at the different focal planes for the different wavelength/energy values. That does not accurately model the system, since in reality all measurements are made at a single focal plane.

Good luck,

Eric

Eric on 23 May 2012

The reference for the equation for W is:

R. J. Noll, "Zernike Polynomials and Atmospheric-Turbulence," J Opt Soc Am 66, 207-211 (1976).

online at:

http://apps.webofknowledge.com/InboundService.do?SID=4AKJi8LAJLJLFkHocME&product=WOS&UT=A1976BK23700002&SrcApp=EndNote&DestFail=http%3A%2F%2Fwww.webofknowledge.com&Init=Yes&action=retrieve&Func=Frame&customersID=ResearchSoft&SrcAuth=ResearchSoft&IsProductCode=Yes&mode=FullRecord

If you look in Table 1 you can see that Z11 (which has n = 4 and m = 0) is the equation I’ve given.

Note that there are other definitions of Zernike polynomials. Many forms omit the sqrt(5) term I've included, in which case the coefficient is the peak-to-valley wavefront error. Other definitions also change the indexing of the terms. Whenever you're working with Zernike polynomials it is absolutely critical to specify which scaling/indexing convention you are using. Since I specify this as spherical aberration (and not Z11 or some other indexing notation) and since I've specified the equation explicitly, I've satisfied that requirement.

The rest of the code is pretty straightforward Fourier Optics, for which I recommend Goodman’s Introduction to Fourier Optics.

RMS_SA is a scaling factor that specifies the root-mean-square of the wavefront. See

http://mathworld.wolfram.com/Root-Mean-Square.html

for what is meant by RMS. If you calculate the RMS of the wavefront over the aperture, you will get RMS_SA. Making RMS_SA bigger increases the amount of spherical aberration. Setting it to zero allows you to model a diffraction-limited optical system.

Let obj.Map be the phase of the surface in waves. Assume that obj.Map is NaN outside of the aperture. I calculate the RMS via

idx = find(~isnan(obj.Map));

output = sqrt(sum(sum((obj.Map(idx)).^2))/numel(idx));

This implements Equation 2 of the aforementioned Mathworld article.

By expressing phase in waves, I mean that you need to multiply this number by 2*pi to express the phase in radians.

Hopefully this helps.

-Eric

Eric on 25 Nov 2013

I received from questions from "Hadi" via personal email but am unable to reply via email. Hopefully Hadi will check here for answers. To wit...

Question: "What do you mean by "Sampling Plan", on line 2; psf_sampling = 0.5e-6;%focal plane sampling in meters"

Answer: This is the center-to-center spacing of the pixels at the focal plane in microns. The pixels are assumed to be square.

Question: On line 15, why do u say: "assert(max(R_norm(:))>=sqrt(2),'Sampling is not sufficient to reconstruct the entire wavefront.');"

Answer: This is a rather inelegant solution, but this code was rather quick-and-dirty. Traditional Fourier sampling properties apply. Variables delta_fx, x_pupil, etc. calculate the spatial sampling of the pupil plane based on the defined sampling of the focal plane. R_norm is the radial coordinate of the pupil plane normalized to the specified aperture radius. If the normalized pupil sampling doesn't go all the way up to 1 at the edge of the sampling array, you have not sampled the entire pupil. This means that at the corner (where R_norm(:) is a max, the value for R_norm should be sqrt(2).

A more elegant solution would be to calculate the optical sampling ratio (frequently denoted Q). Set

Q = lambda * focal_length/(aperture_diameter*psf_sampling);

You should also be able to test that Q>=1 for proper sampling. Q>=2 is generally preferred to avoid aliasing effects.

Question: "-- Is there any other ways to impose Aperture? : E(R_norm>1) = 0;%Impose aperture size"

Answer: You are free to define whatever aperture shape and amplitude you like. I have made the pupil an unobscured circular aperture whose amplitude is equal to 1 on the inside, 0 on the outside. The phase is defined by 2*pi*W. You could apodize the aperture by multiplying the E I calculate by an appropriate function. You could add an obscuration by setting E equal to 0 on the interior as well.

Question: "- Could you please define, the second command parameters : figure;imagesc(angle(E)/(2*pi));colorbar;title('Wavefront Phase (waves)'); i mean "angle(E)/(2*pi)"

Answer: angle(E) is using the built-in Matlab angle() function to calculate the phase of complex array E. In this case E is the complex pupil function (see Goodman's text for a discussion of this).

Question: "On Create point-spread function section, what do u mean by; figure;imagesc(x_psf*1e6,x_psf*1e6,psf);"

Answer: This is just opening a figure and then plotting the point-spread function. The first two arguments passed to imagesc() are the sampling vectors in the point-spread function plane. The sampling is the same in both the X and Y directions, hence the repeated use of x_psf. I multiply x_psf by 1e6 to convert from meters to microns. The third argument, psf, is the calculated point-spread function.

Hopefully this helps.

Best regards, Eric

Hadi on 25 Nov 2013

DEAR ERIC,

HI,

I received your valuable answer,thanks again. it is so helpful. As you mentioned, to simulate other aberrations, i should add the aberrations, and calculate the whole phase, other steps are the same as you define on the code, Am I right?

tanx again Hadi

Eric on 26 Nov 2013

Hadi,

That's right. Add other aberrations to the variable W and you should be set.

-Eric

Link

Answer by Hadi on 26 Nov 2013

Eric, hi again; The all problem is, I want to simulate aberrations by Zernike polynomials to 15th sentence, for a 16" telescope, but i think, somewhere there are some defeats on code, could you please help me to solve the problem? this is the code:

   clear all;
   close all;
   psf_sampling = 0.5e-6;%focal plane sampling in meters
   lambda = 0.6328e-6;%wavelength in meters
   N=256;
teta=45;% the angle to compute
aperture_diameter = 0.407;%meters; 16 inch
focal_length = 5*aperture_diameter;%meters

%RMS_SA = 0.15;%RMS spherical aberration content in waves

% Calculate Pupil Plane Sampling

   delta_fx = 1/(psf_sampling*N);
   x_pupil = (-fix(N/2):fix((N-1)/2)) * delta_fx * lambda * focal_length;
   [X_pupil,Y_pupil] = meshgrid(x_pupil);
   R_pupil = sqrt(X_pupil.^2 + Y_pupil.^2);
   R_norm = R_pupil/(aperture_diameter/2);%Pupil normalized to aperture diameter
   assert(max(R_norm(:))>=sqrt(2),'Sampling is not sufficient to reconstruct the entire 
wavefront.');

%% Create Wavefront & Aberations

%z0 is equal to 1

z1 =  (2 * R_norm) *cos(teta);%Tip Aberration
z1(R_norm>1)=0;
z2 =  (2*R_norm) *sin(teta);%Tilt Aberration
z2(R_norm>1)=0;
z3 = sqrt(3) * (2* R_norm.^2 - 1);% De-focus Aberration
z3(R_norm>1)=0;
z4 =  sqrt(6) * (R_norm.^2) * sin(2*teta);% Astigmatism Aberration
z4(R_norm>1)=0;
z5 = sqrt(6) * (R_norm.^2) * cos(2*teta);% Astigmatism Aberration
z5(R_norm>1)=0;
z6 = sqrt(8) * (3 * R_norm.^3 - 2 * R_norm) * sin(teta);% Coma Aberration
z6(R_norm>1)=0;
z7 =  sqrt(8) * (3 * R_norm.^3 - 2 * R_norm) * cos(teta);% Coma Aberration
z7(R_norm>1)=0;
z8 =  sqrt(8) * (R_norm.^3) * sin(3*teta);% Trefoil aberration
z8(R_norm>1)=0;
z9 =  sqrt(8) * (R_norm.^3) * cos(3*teta);% Trefoil aberration
z9(R_norm>1)=0;
z10 =  sqrt(5) * (6 * R_norm.^4 - 6 * R_norm.^2 + 1);%Spherical Aberration wavefront
z10(R_norm>1) = 0;
z11 =  sqrt(10) * (4 * R_norm.^4 - 3 * R_norm.^2) * cos(2 * teta);% ----- Aberration
z11(R_norm>1) = 0;
z12 =  sqrt(10) * (4 * R_norm.^4 - 3 * R_norm.^2) * sin(2 * teta);% ----- Aberration
z12(R_norm>1) = 0;
z13 =  sqrt(10) * (R_norm.^4) * cos(4 * teta)% ----- Aberration
z13(R_norm>1) = 0;
z14 = sqrt(10) * (R_norm.^4) * sin(4 * teta)% ----- Aberration
z14(R_norm>1) = 0;
    E= exp(1i*2*pi*z1) + exp(1i*2*pi*z2) + exp(1i*2*pi*z3) + exp(1i*2*pi*z4) + exp(1i*2*pi*z5) + exp(1i*2*pi*z6) + exp(1i*2*pi*z7) + exp(1i*2*pi*z8) + exp(1i*2*pi*z9) + exp(1i*2*pi*z10) + exp(1i*2*pi*z11) + exp(1i*2*pi*z12) + exp(1i*2*pi*z13) + exp(1i*2*pi*z14);%Complex amplitude for first 15 zernike polynomials
    E(R_norm>1) = 0;%Impose aperture size
   figure;imagesc(angle(E)/(2*pi));colorbar;title('Wavefront Phase (waves)');

%% Create Point-Spread Function

   psf = abs(fftshift(fft2(ifftshift(E)))).^2;
   psf = psf/sum(psf(:));%Normalize to unity energy
x_psf = (-fix(N/2):fix((N-1)/2)) * psf_sampling;
figure;imagesc(x_psf*1e6,x_psf*1e6,psf);
title(sprintf('PSF Aberration'));
xlabel(sprintf('Position (\\mum)'));
ylabel(sprintf('Position (\\mum)'));

I write only for specific angle(45), and the weights for all aberrations are equal to 1. Could u tell me this code is OK or not?

Hadi

2 Comments

Eric on 27 Nov 2013

Hadi,

I emailed you back but see you also posted here. I'll put a short answer here for posterity's sake.

I see two problems with this code:

1. teta should be set equal to atan2(Y_pupil,X_pupil) and is hence an array. This is the angular coordinate of the pupil plane. This means that all of the subsequent multiplications of teta terms by R_norm terms need to use .* rather than *. Once you do this you can plot the z1, z2, z3, etc. terms and they should look like the appropriate wavefronts.

2. Setting all of the aberrations equal to 1 creates a very aberrated wavefront. The spatial variations are severe and you need to be very careful to avoid aliasing. Even when increasing the array size to 2048 I still think the wavefront is aberrated in this case.

Good luck,

Eric

Hadi on 29 Nov 2013

Hi,

DEAR ERIC

First, Excuse me. I only check the mathworks page.

Second, YOUR answers are amazing, Thanks again for the helpful replies.

You are right the coefficients extremely effect on aberrations, but for the first step i used 1 for all of them. Yesterday i found a function, Created by Dick Brunson, on 12/1/97. This code generates all Zernike polynomial, and size of row or col of the zern output array.

The code is here:

function [zern, mask] = genzern(nz,arraysize,mask)
% function [zern, mask] = genzern(nz,arraysize,mask)
%
% function to generate a matrix of phase values from a single
% Zernike Polynomial coefficient.  Enter genzern with no parameters
% to view the definitions of first 15 Zernike Polynomial terms.
%
% Created by:
% Dick Brunson
% 12/1/97
%
% Inputs:
%
%  nz               number of Zernike Coefficient desired(if <0 print out
%                   definition of Zernike Polynomial
%  arraysize        size of row or col of the zern output array.
%
%  mask[optional]   array of 1's & 0's to define boundaries of Zernike 
%                   output array. If undefined, output array is max  circle 
%                   that can be inscribed inside the array.   
%
% Outputs:
%  zern             output array of phase values
%  mask             pupil function or definition of good phase point area
%                   (defined by mask input or default)
%
% Zernike Polynomials definiton from "Optical Shop Testing", Daniel Malacara, 1992
% pg 465
%
if nargin < 1 
   disp('Zernike #    n    m    Zernike                    Definition')
   disp('   1         0    0       1                       Piston')
   disp('   2         1    0    r*sin(phi)                 Tilt Y')
   disp('   3         1    1    r*cos(phi)                 Tilt X')
   disp('   4         2    0    r^2*sin(2*phi)             Astig 1st ord. 45 deg')
   disp('   5         2    1    2*r^2 - 1                  Defocus')
   disp('   6         2    2    r^2*cos(2*phi)             Astig 1st ord. 0 deg')
   disp('   7         3    0    r^3*sin(3*phi)             Trifoil 30 deg')
   disp('   8         3    1    (3*r^3 - 2*r)sin(phi)      Coma Y')
   disp('   9         3    2    (3*r^3 - 2*r)cos(phi)      Coma X')
   disp('  10         3    3    r^3*cos(3*phi)             Trifoil 0 deg')
   disp('  11         4    0    r^4*sin(4*phi)             Tetrafoil 22.5 deg')
   disp('  12         4    1    (4*r^4 - 3*r^2)sin(2*phi)  Astig 2nd ord 45 deg')
   disp('  13         4    2    6*r^4 - 2*r^2 - 1          Spherical Aberration')
   disp('  14         4    3    (4*r^4 - 3*r^2)cos(2*phi)  Astig 2nd ord 0 deg')
   disp('  15         4    4    r^4*cos(4*phi)             Tetrafoil 0 deg')
   return
end
xv = [1:arraysize]';
yv = [1:arraysize];
x = [];
y = [];
for i = 1:arraysize
   x = [x, xv];
   y = [y; yv];
end
center = 0.5 + arraysize/2;
r = sqrt((x - center).^2 + (y - center).^2);
r = r/(center - 1);
if nargin < 3
   % define inscribed circle for mask
   router = r(arraysize/2,1);
   mask = r <= router;
   rmax_new = max(max(r))/router;
   r = r * rmax_new/max(max(r));
else
   % rescale r array for outer radius of input mask
   maskmax = max(max(mask));
   maskmin = min(min(mask));
   thresh = maskmin + 0.25;
   mask = mask > thresh;
   % outer radius of input mask
   router = max(max(r .* mask));
   % rescale r array so that mask outer radius = 1.
   r = r/router;
end
   phi = atan2((x - center),(y - center));
   zern =zeros(arraysize);
%  compute indices m & n for Zernikes per definition
   csum = cumsum(0:nz);
   nindx = sum(csum < nz);
   mx = csum(nindx);
   n = nindx - 1;
   m = nz - mx - 1;
   %
   %  Zernike generating function from "Optical Shop Testing", Daniel Malacara, 1992
   %
   %              ___m
   %    (n-2m)    \          s                                (n-2s)
   %   R     =    /__    (-1)  [(n-s)!/(s!(m-s)!(n-m-s)!)] * r        
   %    n         s=0 
   %              ___m
   %              \              q(s)
   %         =    /__    p(s) * r        
   %              s=0
   % 
   %  Where the Zernike Polynomial is defined by:
   %
   %            ___k   ___n
   %            \      \         n-2m   sin(n-2m)phi
   % W(r,phi) = /__    /__  Anm Rn      or   
   %             n=0    m=0             cos(n-2m)phi
   %
   %  where sine function is used for n-2m>0 & cosine for n-2m<=0
   %
   for s = 0:m
      if (n - s) < 0 | (m - s) < 0 | (n - m - s) < 0 
         break
      end
      p = (-1)^s * prod(1:(n - s))/(prod(1:s) * prod(1:(m - s)) *   prod(1:(n - m - s)));
      q = n - 2*s;
      zern = zern + p * r .^ q;
%      disp(num2str([nz,n,m,n-2*m,p,q]))
   end
   if n - 2*m > 0
      zern = zern .* sin((n - 2*m)*phi);
   else
      zern = zern .* cos((n - 2*m)*phi);
   end
zern = zern .* mask;
return

I tried to use this function on simulation code, so we can select the number of Zernike Coefficient (and the aberrations) desired in code, and increase the precision of Zernike polynomials, on wavefront. But there are some issues, the first is how to relate aperture diameter to arraysize, the second is how to apply the coefficient weights? For the second issue, i think we can insert the weights by a "For Loop" in the function.

What is your idea?

thanks again

Hadi

Link

Answer by Eric on 3 Dec 2013

Hadi,

I'll have to take a look at this code later. I am reasonably sure that Malacara's Zernikes use a different ordering and scaling coefficient than the "Noll" Zernikes I use. That's one thing to be careful of. I'll have to look through it carefully to see how to apply the aperture size.

You might also be interested in the article at:

http://www.mathworks.com/company/newsletters/articles/analyzing-lasik-optical-data-using-zernike-functions.html

-Eric

2 Comments

dong on 7 Apr 2017

Hi Eric,

I have a question about simulation the PSF of an anisotropic point with half-plane pupil function. I know a vector point-spread function “Ipsf1” which is consisted of complex numbers, half-plane pupil function “pupil.” The calculating process is: 1. inverse Fourier transform “psf1” to get “z”; 2. multiple z with half-plane pupil function “pupil” to get “zp”; 3. Fourier transforms “zp” to get a new PSF “zz”.

The important point that the Ipsf1 is a vector point-spread function. The Ipsf0_1, zp, and zz should be consist of complex numbers. I think there is some problem in my Matlab code for Fourier transform. What should I do for the vector point-spread function Fourier transform? Is it necessary for me to use vectorial pupil function to calculate the new PSF? Could you give me some suggestion?

The Matlab code is:

[num,txt,raw]=xlsread('ipsf0_1.xlsx'); R=cellfun(@(x) ~isnumeric(x) isnan(x),raw); % Find non-numeric cells [m n]=size(raw); data=str2double(raw®); Ipsf0_1=reshape(data,m,n); %read the Ipsf0_1

N = 93; [x y]=meshgrid(-(N-1)/2:(N-1)/2); [theta rho]=cart2pol(x,y); Q=find(rho<30); pupil =zeros(N,N); % construct the pupil function pupil(Q)=1; [a,b]=size(pupil); for q=1:a/2 for p=1:b if pupil(q,p)>0 pupil(q,p)=0; % construct half-plane pupil function end end end

z=fftshift(fft2(Ipsf0_1)); z=abs(z); zp=z.*pupil; zz=ifftshift(ifft2(zp));

Ipsf1=abs(Ipsf0_1); Ipupil=pupil; Izz=abs(zz);

subplot(1,3,1);imshow(Ipsf1./max(max(Ipsf1))); subplot(1,3,2);imshow(Ipupil); subplot(1,3,3);imshow(Izz./max(max(Izz)));

The key code is: z=fftshift(fft2(Ipsf0_1)); z=abs(z); zp=z.*pupil; zz=ifftshift(ifft2(zp));

Walter Roberson on 7 Apr 2017

dong,

Eric will not be notified about new comments, so I recommend that you contact him by way of his profile https://www.mathworks.com/matlabcentral/profile/authors/210783-eric

Link

Answer by Eric on 3 Dec 2013

I looked over this function quickly and I think you want to set the "mask" input value of genzern() to R_norm<=1 from my code.

To implement this code you'll need to implement a for loop. The input argument nz is a scalar for the Zernike index and so this function returns only a single wavefront. You'll need to loop over the indices for the polynomials of interest and multiply the zern output of the function by the appropriate scalar. Sum these scaled wavefronts together to get the total wavefront.

Again, be careful about scaling. The normalization of these Zernikes is different than the ones you have in your earlier post. See the caveat in my original answer about being very specific about how Zernike coefficients are defined. There are lots of definitions. Malacara's are as valid as others, but just always make sure you know which definition you're working with. For what it's worth, I prefer Noll's definition from R. J. Noll, "Zernike Polynomials and Atmospheric-Turbulence," J Opt Soc Am 66, 207-211 (1976). The nice thing about this definition is that the root-sum-square of the Zernike coefficients is equal to the RMS wavefront error for an unobscured system. This makes comparing Zernike magnitudes of different terms considerably easier.

-Eric

1 Comment

Eric on 19 Jun 2017

The assert statement is checking that the sampling allows reconstruction of the entire wavefront. If this assertion fails, then the pupil sampling is such that you cannot reconstruct the entirety of the wavefront area. This can happen if the pixel pitch is relatively large. Assuming you don’t want to adjust the aperture diameter, you need to increase the values in R_norm.

For instance, if I set the psf_sampling to 10e-6 meters, the code from the page results in max(R_norm(:)) of 0.4475.

A simple fix is to calculate the PSF at a finer resolution and then bin that PSF to the desired low-resolution. So for the case where you want a PSF sampling of 10 microns, you could calculate the PSF at a sampling of 2.0 microns. Calculate the psf as shown in the script. Then do the following:

kernel = ones(5,5)/5^2;
psf_lowres_all = conv2(psf, kernel, 'same');
psf_lowres = psf_lowres_all(3:5:end,3:5:end);

In this case you’re doing convolution to spatially sum 5x5 neighborhoods of the high-resolution PSF (5 is the correct value because it’s the ratio of the desired low-resolution pixel spacing to the calculated spacing). This convolution produces 25 possible low-resolution PSFs, each at a different registration of the optical PSF to the pixel grid. In the last line I’ve selected the one in the middle. In the indexing you skip by 5 because you want to get sums of adjacent blocks (i.e., not sliding blocks). You can change the indexing to select another possible low-resolution PSF. For instance, you could also use

psf_lowres = psf_lowres_all(1:5:end,4:5:end);

to select another PSF.

You’ll notice when this is done that the wavefront calculation area is zero-padded to be larger than the specified system aperture. You can check this by plotting abs(E).

You are now following this question