|
|
[求助]
傅里叶光学MATLAB仿真的问题
来个大神吧!我搞了很多天实在搞不出才上这来问的!!!
要在MATLAB中模拟 光的衍射,用的 DavidGeorgeVoelz_Computational_Fourier_Optics_a_MATLAB_tutorial 这本书的代码和插图
![傅里叶光学MATLAB仿真的问题]()
![傅里叶光学MATLAB仿真的问题-1]()
实际实验中,observation plane处放置的是一个CCD相机
相机的像素尺寸为dx2 = 29.6 um
图片的大小为 512 x 512, M = 512
L2 = M*dx2 (1.52 cm)
dx2,M , L2 这些参数都是固定的,不能改变
传播距离 z = 5 cm,光的波长 lambda = 675 nm
运用夫琅禾费衍射公式
![傅里叶光学MATLAB仿真的问题-2]()
![傅里叶光学MATLAB仿真的问题-3]()
得到dx1 = lambda*z/L2 (dx1 为source plane 处的像素尺寸,为2.227um)
L1 = M*dx1 (L1为source plane的边长 1.14 mm,source plane处也是一个 M x M 大小的图片)
我的仿真是在source plane 处放置一个边长 D = 50 um 的图片(此处用矩孔代替),然后得到衍射图样,这个 D 也是固定不能变的
但是这样一来,D/dx1 = 22.4521 个像素,图片基本上看不出细节
w = D/2
仿真结果如下
![傅里叶光学MATLAB仿真的问题-4]()
![傅里叶光学MATLAB仿真的问题-5]()
我的问题是 dx2 ,M ,lambda ,z ,D 这些参数定死了,那么dx1 ,L1 这些参数根据公式也定死了,导致我插入的图片看不出细节(D/dx1 = 22.4521 个像素)
我想的是在source plane 处的 L1 能否缩小,然后M还是不变的,这样dx1就变小了,D/dx1就能更大些,插入的图片也能看到更多的细节?
在边长为D的矩形孔外都是0,所以L1的缩小在真实的实验中无关紧要
就像下面这样
![傅里叶光学MATLAB仿真的问题-6]()
%========================================
clear all;clc;
cm = 1e-2;
mm = 1e-3;
um = 1e-6;
nm = 1e-9;
%%
lambda = 675*nm;
k = 2*pi/lambda;
dx2 = 29.6*um;
M = 512;
z = 5*cm;
%%
L2 = M*dx2;
L1 = lambda*z/dx2;
dx1 = L1/M;
w = 25*um; % 这个 w = D/2
%%
x1 = -L1/2:dx1:L1/2-dx1;
y1 = x1;
[X1,Y1] = meshgrid(x1,y1);
u1 = rect(X1/(2*w)).*rect(Y1/(2*w));
figure(1) %irradiance image
imagesc(x1./mm,y1./mm,u1);
xlabel(\'x (mm)\'); ylabel(\'y (mm)\'); title(\'source plane z = 0\');
colormap(\'gray\');
axis square;
axis xy;
%%
[u2,L2]=propFF(u1,L1,lambda,z);
dx2 = L2/M;
x2 = -L2/2:dx2:L2/2-dx2; %obs ords
y2 = x2;
I2 = abs(u2.^2);
figure(2) %irradiance image
imagesc(x2./cm,y2./cm,nthroot(I2,3));
xlabel(\'x (cm)\'); ylabel(\'y (cm)\'); title(strcat(\'observation plane z = \',num2str(z/cm),\' cm\'));
colormap(\'gray\');
axis square;
axis xy;
%=================================================
function [out] = rect(x)
out=abs(x)<=1/2;
end
%=================================================
function[u2,L2]=propFF(u1,L1,lambda,z)
% propagation - Fraunhofer pattern
% assumes uniform sampling
% u1 - source plane field
% L1 - source plane side length
% lambda - wavelength
% z - propagation distance
% L2 - observation plane side length
% u2 - observation plane field
%
[M,N] = size(u1); %get input field array size
dx1 = L1/M; %source sample interval
k = 2*pi/lambda; %wavenumber
%
L2 = lambda*z/dx1; %obs sidelength
dx2 = lambda*z/L1; %obs sample interval
x2 = -L2/2:dx2:L2/2-dx2; %obs coords
[X2,Y2] = meshgrid(x2,x2);
%
c = 1/(j*lambda*z)*exp(j*k/(2*z)*(X2.^2+Y2.^2));
u2 = c.*ifftshift(fft2(fftshift(u1)))*dx1^2;
end
%====================================================
[ Last edited by shepherd2014 on 2015-3-19 at 10:31 ] |
|