tic
clear;
lambda=0.0004;%波長單位mm
theater0=0.00035;%高斯光束發散角單位rad
w0=lambda/pi/theater0;
A0=250;
A=-50;
B=-9802;
D=-0.02;
L=704;%ABCD光學系統光軸上光程
nx=50;%觀察平面上分100個點
ny=50;
nx1=50;%將束腰分成100個點
ny1=50;
r1=0.5;%單位mm
r2=100;%觀察圓面半徑單位mm
x0=linspace(-r1,r1,nx1);%束腰平面上的點
y0=linspace(-r1,r1,ny1);
x00=linspace(-r2,r2,nx);%觀察平面上的點
y00=linspace(-r2,r2,ny);
[X00,Y00]=meshgrid(x00,y00);
        for m1=1:nx
            for n1=1:ny
                a=0;
                b=0;
                for m11=1:nx1
                    for n11=1:ny1
                    x=-r1+(2.*m11-1)./nx1.*r1;
                    y=-r1+(2.*n11-1)./ny1.*r1;
                    e(m11,n11)=pi.*r1.^2./nx1./ny1.*exp((-(x.^2+y.^2))./w0./w0+i.*pi.*A./lambda./B.*(x.^2+y.^2)+i.*2.*pi./B./lambda.*(x.*x00(m1)+y.*y00(n1)));
                    end 
                end
                    a=real(e)+a;
                    b=imag(e)+b;
                    u=(a+i.*b).*A0./w0.*(-i)./lambda./B.*exp(i.*2.*pi./lambda.*L).*exp(i.*pi.*D./B.*(X00.^2+Y00.^2));
                    G=(real(u).^2+imag(u).^2)*10*30;
            end            
        end 
        surf(X00,Y00,G);
        xlabel('X(mm)','fontsize',12);
        ylabel('Y(mm)','fontsize',12);
        zlabel('Relative intensity','fontsize',12);
        toc