基於MATLAB菲涅爾衍射模擬

format long 
clear 
n1=1.45 
n2=1.447 
a=4.5  
k=(a^2)*(n1^2-n2^2); 
u1=(a^2)*(n1^2); 
w1=(a^2)*(n2^2);  
delta1=(n1^2+n2^2)/(2*n1^2); 
delta2=(n1^2-n2^2)/(2*n1^2); 
i=1; 
n=n2;  
for V=0:0.01:6     
    k02=V^2/k;     
    k01=sqrt(k02);      
    for BeiTa=n:0.00001:n1         
        %n          
        nn=n+0.00001          
        U2=u1*k02-(a^2)*k02*(BeiTa^2);         
        U=sqrt(U2);          
        W2=(a^2)*k02*(BeiTa^2)-w1*k02;         
        W=sqrt(W2);         
        if(U==0 || W==0)             
            break;              
            %disp('sss')         
        else
            z1=n1^2*besselj(1,U)/(U*besselj(0,U));             
            z2=n2^2*besselk(1,W)/(W*besselk(0,W));              
            %z3=1/(U^2)+delta1*(1/W^2-z2)-sqrt(delta2^2*(1/W^2-z2)^2+(BeiTa/n1)^2*(V/(U*W))^4);              
            if(abs(z1+z2)<0.01)                 
                x(i)=V;                 
                y(i)=BeiTa;                 
                i=i+1;                 
                z1;                 
                z2;                   
                %disp('nnn')                 
                n=BeiTa;                 
                break;             
            end
        end
    end
end
plot(x,y);  
axis([0 6.5 1.447 1.45]);