clear  %方形平凹鏡腔邊界元法（平面鏡光場分布、遠場分布、本征值）
L=140;  %腔長
a=1;  %鏡寬
R1=500;   %凹鏡曲率半徑
W=10.6e-4;  %波長
N=a^2/(W*L)  %菲涅爾數(shù)
e=100;  %本征值個數(shù)
M=200;  %分元數(shù)
K=2*pi/W;  %波數(shù)
G1X=zeros(M);   %X方向凹鏡向平鏡度越矩陣
G2X=zeros(M);    %X方向平鏡向凹鏡度越矩陣
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)));
      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)));
   end
end
G1Y=G1X;   %Y方向凹鏡向平鏡度越矩陣
G2Y=G2X;   %Y方向平鏡向凹鏡度越矩陣
GX=G1X*G2X;   %X方向平鏡-凹鏡-平鏡度越矩陣
GY=G1Y*G2Y;   %Y方向平鏡-凹鏡-平鏡度越矩陣
%save 0000 G1X G2X G1Y G2Y
clear G1X G1Y G2X G2Y
%save 0001 GX GY
[eigenvectorGX,eigenvalueGX]=eig(GX);   %X方向的本征值和本征矢。
[eigenvectorGY,eigenvalueGY]=eig(GY);   %Y方向的本征值和本征矢。
clear GX GY
indexx=zeros(M);   %光場分布坐標
 indexy=zeros(M);   %光場分布坐標
for m1=1:M
    indexx(m1)=-a+2*a*m1/M;
    indexy(m1)=indexx(m1);
end
save 0006 eigenvectorGX eigenvectorGY indexx indexy
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);
    for m2=1:M
    for n2=1:M
        QXG(m2,n2)=(2*a/M)*exp(-i*K*((-15*W/(2*a)+15*2*W*m2/(2*a*M))*(-a+2*a*n2/M)))*sqrt(exp(-i*K*D));
        QYG(m2,n2)=(2*a/M)*exp(-i*K*((-15*W/(2*a)+15*2*W*m2/(2*a*M))*(-a+2*a*n2/M)))*sqrt(exp(-i*K*D));
    end
end
indexxf=zeros(1,M);   %遠場分布坐標
 indexyf=zeros(1,M);   %遠場分布坐標
for m1=1:M
    indexxf(m1)=-15+15*m2/M;
    indexyf(m1)=indexx;
end
save 0007  indexxf indexyf QXG QYG
save 0002 msumEIG
%save 0005 sumEIG
%clear msumEIG sumEIG
%以上是部分總本征值的模。
%save 9999 L a R1 W M K
clear L a R1 W M K