clear  %方形平凹鏡腔邊界元法平鏡傾斜計(jì)算
L=140;
a=1;
R1=500; 
%以上為參數(shù)
W=10.6e-4;
N=a^2/(W*L)
M=200;
t=0.0134;
q=t*pi/180;
e=100;
%W是波長，M是分元總數(shù)。
K=2*pi/W;
%K是波數(shù)。
G1X=zeros(M);
G2X=zeros(M);
G1Y=zeros(M);
G2Y=zeros(M);
for m=1:M
   for n=1:M
      G1X(n,m)=(2*a/M)*sqrt((i/(W*L))*exp(-i*K*L))*exp(-i*K*(((-a+2*a*n/M)-(-a+2*a*m/M))^2/(2*L)-(-a+2*a*m/M)^2/(2*R1)+(-a+2*a*n/M)*sin(q)));
      G2X(m,n)=(2*a/M)*sqrt((i/(W*L))*exp(-i*K*L))*exp(-i*K*(((-a+2*a*m/M)-(-a+2*a*n/M))^2/(2*L)-(-a+2*a*m/M)^2/(2*R1)+(-a+2*a*n/M)*sin(q)));
      G1Y(n,m)=(2*a/M)*sqrt((i/(W*L))*exp(-i*K*L))*exp(-i*K*(((-a+2*a*n/M)-(-a+2*a*m/M))^2/(2*L)-(-a+2*a*m/M)^2/(2*R1)));
      G2Y(m,n)=(2*a/M)*sqrt((i/(W*L))*exp(-i*K*L))*exp(-i*K*(((-a+2*a*m/M)-(-a+2*a*n/M))^2/(2*L)-(-a+2*a*m/M)^2/(2*R1)));
   end
end
GX=G1X*G2X;%算出為平鏡光場
GY=G1Y*G2Y;
%save 0000 G1X G2X G1Y G2Y
clear G1X G1Y G2X G2Y
%以上為總矩陣。
%save 0001 GX GY
[eigenvectorGX,eigenvalueGX]=eig(GX);
[eigenvectorGY,eigenvalueGY]=eig(GY);
%以上是X方向，Y方向的本征值和本征矢。
clear GX GY
save 0006x eigenvectorGX eigenvectorGY
clear eigenvectorGX eigenvectorGY
evGX=zeros(e,1);
evGY=zeros(1,e);
for i=1:e
    evGX(i,1)=eigenvalueGX(i,i);
    evGY(1,i)=eigenvalueGY(i,i);
end
%save 0003 eigenvalueGX eigenvalueGY
clear eigenvalueGX eigenvalueGY
    sumEIG=evGX*evGY;
    msumEIG=abs(sumEIG);
%save 0002 msumEIG
%save 0005 sumEIG
%clear msumEIG sumEIG
%以上是部分總本征值的模。
%save 9999 L a R1 W M K
clear L a R1 W M K
%以上為數(shù)據(jù)存儲。