Field-Focus Plot

field_focus_plot.m

function field_focus_plot(W, Qt, myaxes, lambda, ua, uc)

h = linspace(myaxes(3),myaxes(4),200)/uc;
hsq = h.*h;
Qsq = (4/9)*W(1)^2 + hsq.*((2/3)*(W(2))^2 + hsq.*(1/2)*W(3)^2);
A = (4/3)*W(1)+hsq*(W(4)+0.5*W(3));
eps = lambda/ua;
Qtsq = (Qt/(eps))^2;
radical = (Qtsq-Qsq)/2;
idx = find(radical>=0);
zeta = 2*lambda/ua^2;
x = sqrt(radical(idx));
y = h(idx);
m = length(x);
Am = -A(idx);
xp = Am-x;
xp(m+(1:m)) = fliplr(Am + x);
yp = y;
yp(m+(1:m)) = fliplr(y);
plot(xp*zeta,yp*uc,'--k');
hold on;
plot(Am*zeta,y*uc,'k');
hold off;
axis(myaxes);
xlabel('focal shift (mm)');
ylabel('image height (deg)');

field_focus_image.m

function rgb = field_focus(W, Qt, myaxes, lambda, ua, uc)

N=256;
h = linspace(myaxes(3),myaxes(4),N)/uc;
zeta = 2*lambda/ua^2;
A = linspace(myaxes(1),myaxes(2),N)/zeta;
k = ones(N,1);
hsq = h.*h;
Qsq = (4/9)*W(1)^2 + hsq.*((2/3)*(W(2))^2 + hsq.*(1/2)*W(3)^2);
Am = (4/3)*W(1)+hsq*(W(4)+0.5*W(3));
out = (k*Qsq)'+2.0*(k*A+(k*Am)').^2;
eps = lambda/ua;
Qtsq = (Qt/(eps))^2;
out = flipud(out/Qtsq);
m = out<1;
b = zeros(size(out));
out = m-m.*out;
rgb = cat(3,out,out,1-m);
imshow(rgb);

Example

>> W = [ 1 2 5 4 ];
>> ua = 0.05;
>> uc = 10;
>> lambda = 0.55e-3;
>> et = 25e-3;
>> myaxes = [ -2.5 0.5 0 8 ];
>> field_focus_plot(W, et, myaxes, lambda, ua, uc);
>> out = field_focus(W, et, myaxes, lambda, ua, uc);
>> imwrite(out,'field_focus.tif');

Field-Focus plot

Field-Focus image

Maintained by John Loomis, last updated 28 Nov 2001