clear  %方形平凹鏡腔邊界元迭代法高斯分布反射率雙鏡傾斜計算（平面鏡光場分布、遠場分布、本征值）
for t1=0.01:0.01:0.13
    for t2=0.01:0.01:0.13
L=340;  %腔長
a=1;  %鏡寬
R1=500;    %凹鏡曲率半徑
W=10.6e-4;  %波長
N=a^2/(W*L)  %菲涅爾數
M=100;  %分元數
t=0.247   %傾斜角度
q1=t1*pi/180;  %傾斜角度的弧度數
q2=t2*pi/180;  %傾斜角度的弧度數
K=2*pi/W;   %K是波數。
G1X=zeros(M);  %X方向凹鏡向平鏡度越矩陣
G2X=zeros(M);    %X方向平鏡向凹鏡度越矩陣
G1Y=zeros(M);   %Y方向凹鏡向平鏡度越矩陣
G2Y=zeros(M);   %Y方向平鏡向凹鏡度越矩陣
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(q1)));
      G2X(m,n)=f1*exp(-(-a+2*a*n/M)^2/w1^2)*(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(q1)));
      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)+(-a+2*a*n/M)*sin(q2)));
      G2Y(m,n)=f1*exp(-(-a+2*a*n/M)^2/w1^2)*(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(q2)));
   end
end
GX=G1X*G2X;%算出為平鏡光場
GY=G1Y*G2Y;
%save 0000 G1X G2X G1Y G2Y
clear G1X G1Y G2X G2Y
%以上為總矩陣。
%save 0001 GX GY
E1X=ones(M,1);  %平鏡X方向光場分布
E1Y=ones(M,1);  %平鏡Y方向光場分布
for d=1:5000
    E1X=GX*E1X;
    E1X=E1X/max(abs(E1X));
    E1Y=GY*E1Y;
    E1Y=E1Y/max(abs(E1Y));
end
MX=max(abs(GX*E1X));  %X向本征值
MY=max(abs(GY*E1Y));  %Y向本征值
M1=MX*MY;  %總本征值
E=E1X*E1Y.';   %平鏡總光場分布
[maxcol,maxindexX]=max(E);  %極大值x坐標
[maxrow,maxindexY]=max(maxcol)  %極大值y坐標
for m2=1:M
indexX(m2)=-a+2*a*m2/M;
indexY(m2)=indexX(m2);
end
mesh(indexX,indexY,abs(E));  %光場幅度畫圖
pause   %程序暫停*******************************************
mesh(indexX,indexY,angle(E));  %畫圖
pause  %程序暫停*******************************************
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
QX=QXG*E1X;  %遠場X向分布
QY=QYG*E1Y;  %遠場Y向分布
Q=QX*QY.';   %遠場總分布
indexx=zeros(1,M);   %遠場分布坐標
 indexy=zeros(1,M);   %遠場分布坐標
for m1=1:M
    indexx(m1)=-15+15*m2/M;
    indexy(m1)=indexx;
end
mesh(indexx,indexy,abs(Q));  %遠場強度畫圖
end
end