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.
Good luck,
Eric
clear all;close all;
psf_sampling = 0.5e-6;
lambda = 0.6328e-6;
N = 256;
aperture_diameter = 0.0254;
focal_length = 5*aperture_diameter;
RMS_SA = 0.25;
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);
assert(max(R_norm(:))>=sqrt(2),'Sampling is not sufficient to reconstruct the entire wavefront.');
W = RMS_SA * sqrt(5) * (6*R_norm.^4 - 6*R_norm.^2 + 1);
W(R_norm>1) = 0;
E = exp(1i*2*pi*W);
E(R_norm>1) = 0;
figure;imagesc(angle(E)/(2*pi));colorbar;title('Wavefront Phase (waves)');
psf = abs(fftshift(fft2(ifftshift(E)))).^2;
psf = psf/sum(psf(:));
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)'));
Best Answer