huitu.zip_生物 光子_生物光子_生物光子学

clear all %%%%%%%%%%%%%%%%%%%% %A diffractive circle with more than 90% of energy in 700mm %of diffractive distance. The uint is millimeter in the program. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% lamd=0.532*10^-3;%% incident wavelength w=0.75;%% waist of incident beam R=4;%% limiting aperture of the element D=700;%% diffractive distance; r1=1.5;r2=3;%% inside and outside radius of circle L0=10;%%mm k=2*pi/lamd;%%%wave number; N=1024; %%%%%%%%%%%%% judging(analytic transformation) % Judging=(sqrt(N*lamd*D)<=L0); % if Judging==0 % disp('????ERROR'); % disp('......Fresnel Analytic Transformation is not satisfied'); % break; % end %%%%%%%%%%%%% x11=linspace(-L0/2,L0/2,N);y11=linspace(-L0/2,L0/2,N); [x1,y1]=meshgrid(x11,y11); J1=zeros(N); for m=1:N for n=1:N if x1(m,n)^2+y1(m,n)^2<=R^2 J1(m,n)=1; end end end A=J1; %%%%% fx=1/L0*(-N/2:N/2-1);fy=1/L0*(-N/2:N/2-1); [fx,fy]=meshgrid(fx,fy); %%%%% JJ=zeros(N); JJ(235:275,235:275)=1;JJ(491:531,235:275)=1;JJ(747:787,235:275)=1; JJ(235:275,491:531)=1;JJ(491:531,491:531)=1;JJ(747:787,491:531)=1; JJ(235:275,747:787)=1;JJ(491:531,747:787)=1;JJ(747:787,747:787)=1; a=sum(sum(A.^4))/sum(sum(JJ.^4)); J2=JJ*sqrt(sqrt(a)); %imagesc(J2);axis square;colormap(gray) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% G-S %pha0=load('C:\Documents and Settings\wjq\My Documents\MATLAB\pha90.txt'); pha0=2*pi*rand(N)-pi;M=0;CC=[];eta=0; while M<30 M=M+1; U1=A.*exp(i*pha0); temp=(fft2(U1));temp1=exp(i*angle(temp)); U2=ifft2(temp1.*J2);pha0=angle(U2);%A2=abs(U2);I2=abs(U2).^2; %Err(M)=sum(sum((I2-J2.^2 ).^2));eta(M)=sum(sum(I2.*JJ))/sum(sum(I2)); subplot(221);imagesc(x11,y11,A);axis square; title('入射光强度') % plotyy(1:M,Err,1:M,eta) % title('Err (left) and eta (right)');xlabel('M'); subplot(222);imagesc(x11,y11,J2);axis square;title('理想的像方平面光照') xlabel('x / mm');ylabel('y / mm');zlabel('Amplitude') subplot(223);imagesc(x11,y11,abs(temp));axis square;title('实际的像方光照情况') xlabel('x / mm');ylabel('y / mm');zlabel('Amplitude') %%%%%%%% %U2=J2.*exp(i*pha2); %tmep3=fft2(U2);temp4=exp(-i*k*D*(1-lamd^2/2*(fx.^2+fy.^2))); %U11=ifft2(tmep3.*temp4);pha0=angle(U11); subplot(224);imagesc(x11,y11,unwrap(unwrap(pha0),[],2));;axis square; title('入射光相位情况') colormap(gray(256)); %subplot(224);imagesc(x11,y11,pha0);axis square; title('phase of DOE') xlabel('x / mm');ylabel('y / mm') pause(0.01) end