matlab仿真基模高斯光束经聚焦透镜聚焦到衍射光栅后，再经衍射光栅衍射后，衍射光由一块同焦距的透镜再次汇聚，画出焦面后光场在x方向上的分布曲线
具体内容如下：
相对于一块焦距为f（举例情况下不作具体约束，任意，后面的符号不作说明，可自定义验证，给出程序及图形结果）的聚焦透镜前一倍焦距处让一基模高斯光束正对该透镜的中心垂直入射，即E0=exp（-x^2/w0^2），入射波长为λ0，让光栅放置在该透镜的后一倍焦距处，光栅刻线d，该聚焦后的光束在光栅上的入射角度为α0（与光栅法线夹角），（这里会有光栅的相位因子对光场的影响，同时还有原入射面与光栅面之间不是同一个坐标系，所以还会有坐标变换的影响，两者需要同时考虑），这之后衍射光再次被一个焦距为f的透镜聚焦，这时候可以观察这个焦面后的光场分布，z=0m、0.1m等不同位置处的光场在x方向上的分别曲线。
[这里需要说明一下，这里坐标系为x0z，光轴方向定位z轴（即传输距离），且光前进的方向为z的正半轴，垂直于z轴定为x轴（也可以理解为空间上的水平方向），不考虑y轴方向光场分布。]
因为我自己在仿真过程中，我仿真的是多个源场光源入射在焦面后不同位置处的光场分布，但是随着传输距离的增加出现了光斑不重叠的现象，想弄明白是哪个地方出问题了，是fft变换错误了还是坐标变换有问题？
以下这是光路图及本人的代码，感谢各位的解答。
（光路图）
（仿真结果问题，各子光束的光强分布不重叠，即中心峰不重叠，这里图片的xy没细节标记，在这图片里x是样本点序号，y是光强计算值）
（这是理想的光强分布情况，中心峰位置一致，只是光斑大小在不同传输距离后不一样罢了
！提示一下，已经被我注释的地方大家就不用管了，还有就是，大家可能看不懂我的代码，所以我还是比较喜欢大家能自己编程一下，按照我一开始的内容设置进行一个光束的光强分布仿真，我想看一下大家是怎么仿真的，实在感谢！
close all
clear
clc
%% about self-Adapative spectral combining beam
%------set light number-----%
light_number=[5 10 15 20 25 30 35 40 45 50];
shunxu = 10;
F_number=light_number(shunxu); %sample number
center_number=(F_number+1)/2; %center sample number
%% ------Initial parameter-----%
w0=3.00e-4; %waist radius
w02=w0*w0; %square of waist
dp=5.0e-3; %adjacent array interval
f=0.3; %focal length
a0=36.25/180*pi; %incidence angle of central beam
r0=1.064e-6; %incidence wavelength of central beam
d=1/960/1000; %diffraction grating
beta=asin(r0/d-sin(a0)); %diffraction angle
circle=1; %circle of caculate M^2
dz=[0.0 1.0 3.0]; %the z of different waist radius
length_dz=length(dz); %the length of z matrix
%-------------Other deductive parameters and settings------------%
for m=1:F_number
r(m)=r0+(m-center_number)*dp*d*cos(a0)/f; %Incident Wavelength of subbeam
k(m)=2*pi/r(m); %Wave Vector
p(m)=(m-center_number)*dp; %Offset interval
p2(m)=p(m)*p(m); %Offset interval squared
a(m)=asin(r(m)/d-sin(beta)); %Incident angle of subbeam
delta(m)=sin(a(m))-sin(beta);%delta0=(sina-sinb)/cosb
end
%% ----------- source field----------%
B=f;C=-1/f;D=0;
beishu=19;beishu1=1;
dxx=d/(2^(beishu1));
E0=1;
M=2^(beishu+beishu1);
u=-M/2:1:M/2-1;
%-------x''o''z'' coords-----%
xx=u*dxx*cos(a0);L=M*dxx*cos(a0);
%-------坐标变换------%
zz1=2*f+xx*p(1)/f; zz2=2*f+xx*p(2)/f; zz3=2*f+xx*p(3)/f; zz4=2*f+xx*p(4)/f;
zz5=2*f+xx*p(5)/f; zz6=2*f+xx*p(6)/f; zz7=2*f+xx*p(7)/f; zz8=2*f+xx*p(8)/f;
zz9=2*f+xx*p(9)/f; zz10=2*f+xx*p(10)/f; zz11=2*f+xx*p(11)/f; zz12=2*f+xx*p(12)/f;
zz13=2*f+xx*p(13)/f; zz14=2*f+xx*p(14)/f; zz15=2*f+xx*p(15)/f; zz16=2*f+xx*p(16)/f;
zz17=2*f+xx*p(17)/f; zz18=2*f+xx*p(18)/f; zz19=2*f+xx*p(19)/f; zz20=2*f+xx*p(20)/f;
zz21=2*f+xx*p(21)/f; zz22=2*f+xx*p(22)/f; zz23=2*f+xx*p(23)/f; zz24=2*f+xx*p(24)/f;
zz25=2*f+xx*p(25)/f; zz26=2*f+xx*p(26)/f; zz27=2*f+xx*p(27)/f; zz28=2*f+xx*p(28)/f;
zz29=2*f+xx*p(29)/f; zz30=2*f+xx*p(30)/f; zz31=2*f+xx*p(31)/f; zz32=2*f+xx*p(32)/f;
zz33=2*f+xx*p(33)/f; zz34=2*f+xx*p(34)/f; zz35=2*f+xx*p(35)/f; zz36=2*f+xx*p(36)/f;
zz37=2*f+xx*p(37)/f; zz38=2*f+xx*p(38)/f; zz39=2*f+xx*p(39)/f; zz40=2*f+xx*p(40)/f;
zz41=2*f+xx*p(41)/f; zz42=2*f+xx*p(42)/f; zz43=2*f+xx*p(43)/f; zz44=2*f+xx*p(44)/f;
zz45=2*f+xx*p(45)/f; zz46=2*f+xx*p(46)/f; zz47=2*f+xx*p(47)/f; zz48=2*f+xx*p(48)/f;
zz49=2*f+xx*p(49)/f; zz50=2*f+xx*p(50)/f; %zz51=2*f+xx*p(51)/f;
%ABCD矩阵当中的A值
A1=2-zz1/f; A2=2-zz2/f; A3=2-zz3/f; A4=2-zz4/f; A5=2-zz5/f; A6=2-zz6/f;
A7=2-zz7/f; A8=2-zz8/f; A9=2-zz9/f; A10=2-zz10/f; A11=2-zz11/f; A12=2-zz12/f;
A13=2-zz13/f; A14=2-zz14/f; A15=2-zz15/f; A16=2-zz16/f; A17=2-zz17/f; A18=2-zz18/f;
A19=2-zz19/f; A20=2-zz20/f; A21=2-zz21/f; A22=2-zz22/f; A23=2-zz23/f; A24=2-zz24/f;
A25=2-zz25/f; A26=2-zz26/f; A27=2-zz27/f; A28=2-zz28/f; A29=2-zz29/f; A30=2-zz30/f;
A31=2-zz31/f; A32=2-zz32/f; A33=2-zz33/f; A34=2-zz34/f; A35=2-zz35/f; A36=2-zz36/f;
A37=2-zz37/f; A38=2-zz38/f; A39=2-zz39/f; A40=2-zz40/f; A41=2-zz41/f; A42=2-zz42/f;
A43=2-zz43/f; A44=2-zz44/f; A45=2-zz45/f; A46=2-zz46/f; A47=2-zz47/f; A48=2-zz48/f;
A49=2-zz49/f; A50=2-zz50/f; %A51=2-zz51/f;
%---------xoz???---------%
L1=(L/cos(a0))*cos(beta);
dx1=dxx*cos(beta);
x=-L1/2:dx1:L1/2-dx1;
fx=-1/(2*dx1):1/L1:1/(2*dx1)-1/L1;
%-----------------circle of M2--------------%
%part 1
ez1p1 =E0.*sqrt(k(1)*w02./(1i.*(2*B-1i*k(1).*A1*w02)));
ez2p1 =E0.*sqrt(k(2)*w02./(1i.*(2*B-1i*k(2).*A2*w02)));
ez3p1 =E0.*sqrt(k(3)*w02./(1i.*(2*B-1i*k(3).*A3*w02)));
ez4p1 =E0.*sqrt(k(4)*w02./(1i.*(2*B-1i*k(4).*A4*w02)));
ez5p1 =E0.*sqrt(k(5)*w02./(1i.*(2*B-1i*k(5).*A5*w02)));
ez6p1 =E0.*sqrt(k(6)*w02./(1i.*(2*B-1i*k(6).*A6*w02)));
ez7p1 =E0.*sqrt(k(7)*w02./(1i.*(2*B-1i*k(7).*A7*w02)));
ez8p1 =E0.*sqrt(k(8)*w02./(1i.*(2*B-1i*k(8).*A8*w02)));
ez9p1 =E0.*sqrt(k(9)*w02./(1i.*(2*B-1i*k(9).*A9*w02)));
ez10p1=E0.*sqrt(k(10)*w02./(1i.*(2*B-1i*k(10).*A10*w02)));
ez11p1=E0.*sqrt(k(11)*w02./(1i.*(2*B-1i*k(11).*A11*w02)));
ez12p1=E0.*sqrt(k(12)*w02./(1i.*(2*B-1i*k(12).*A12*w02)));
ez13p1=E0.*sqrt(k(13)*w02./(1i.*(2*B-1i*k(13).*A13*w02)));
ez14p1=E0.*sqrt(k(14)*w02./(1i.*(2*B-1i*k(14).*A14*w02)));
ez15p1=E0.*sqrt(k(15)*w02./(1i.*(2*B-1i*k(15).*A15*w02)));
ez16p1=E0.*sqrt(k(16)*w02./(1i.*(2*B-1i*k(16).*A16*w02)));
ez17p1=E0.*sqrt(k(17)*w02./(1i.*(2*B-1i*k(17).*A17*w02)));
ez18p1=E0.*sqrt(k(18)*w02./(1i.*(2*B-1i*k(18).*A18*w02)));
ez19p1=E0.*sqrt(k(19)*w02./(1i.*(2*B-1i*k(19).*A19*w02)));
ez20p1=E0.*sqrt(k(20)*w02./(1i.*(2*B-1i*k(20).*A20*w02)));
ez21p1=E0.*sqrt(k(21)*w02./(1i.*(2*B-1i*k(21).*A21*w02)));
ez22p1=E0.*sqrt(k(22)*w02./(1i.*(2*B-1i*k(22).*A22*w02)));
ez23p1=E0.*sqrt(k(23)*w02./(1i.*(2*B-1i*k(23).*A23*w02)));
ez24p1=E0.*sqrt(k(24)*w02./(1i.*(2*B-1i*k(24).*A24*w02)));
ez25p1=E0.*sqrt(k(25)*w02./(1i.*(2*B-1i*k(25).*A25*w02)));
ez26p1=E0.*sqrt(k(26)*w02./(1i.*(2*B-1i*k(26).*A26*w02)));
ez27p1=E0.*sqrt(k(27)*w02./(1i.*(2*B-1i*k(27).*A27*w02)));
ez28p1=E0.*sqrt(k(28)*w02./(1i.*(2*B-1i*k(28).*A28*w02)));
ez29p1=E0.*sqrt(k(29)*w02./(1i.*(2*B-1i*k(29).*A29*w02)));
ez30p1=E0.*sqrt(k(30)*w02./(1i.*(2*B-1i*k(30).*A30*w02)));
ez31p1=E0.*sqrt(k(31)*w02./(1i.*(2*B-1i*k(31).*A31*w02)));
ez32p1=E0.*sqrt(k(32)*w02./(1i.*(2*B-1i*k(32).*A32*w02)));
ez33p1=E0.*sqrt(k(33)*w02./(1i.*(2*B-1i*k(33).*A33*w02)));
ez34p1=E0.*sqrt(k(34)*w02./(1i.*(2*B-1i*k(34).*A34*w02)));
ez35p1=E0.*sqrt(k(35)*w02./(1i.*(2*B-1i*k(35).*A35*w02)));
ez36p1=E0.*sqrt(k(36)*w02./(1i.*(2*B-1i*k(36).*A36*w02)));
ez37p1=E0.*sqrt(k(37)*w02./(1i.*(2*B-1i*k(37).*A37*w02)));
ez38p1=E0.*sqrt(k(38)*w02./(1i.*(2*B-1i*k(38).*A38*w02)));
ez39p1=E0.*sqrt(k(39)*w02./(1i.*(2*B-1i*k(39).*A39*w02)));
ez40p1=E0.*sqrt(k(40)*w02./(1i.*(2*B-1i*k(40).*A40*w02)));
ez41p1=E0.*sqrt(k(41)*w02./(1i.*(2*B-1i*k(41).*A41*w02)));
ez42p1=E0.*sqrt(k(42)*w02./(1i.*(2*B-1i*k(42).*A42*w02)));
ez43p1=E0.*sqrt(k(43)*w02./(1i.*(2*B-1i*k(43).*A43*w02)));
ez44p1=E0.*sqrt(k(44)*w02./(1i.*(2*B-1i*k(44).*A44*w02)));
ez45p1=E0.*sqrt(k(45)*w02./(1i.*(2*B-1i*k(45).*A45*w02)));
ez46p1=E0.*sqrt(k(46)*w02./(1i.*(2*B-1i*k(46).*A46*w02)));
ez47p1=E0.*sqrt(k(47)*w02./(1i.*(2*B-1i*k(47).*A47*w02)));
ez48p1=E0.*sqrt(k(48)*w02./(1i.*(2*B-1i*k(48).*A48*w02)));
ez49p1=E0.*sqrt(k(49)*w02./(1i.*(2*B-1i*k(49).*A49*w02)));
ez50p1=E0.*sqrt(k(50)*w02./(1i.*(2*B-1i*k(50).*A50*w02)));
%ez51p1=E0.*sqrt(k(51)*w02./(1i.*(2*B-1i*k(51).*A51*w02)));
%part 2
ez1p2=exp(1i*k(1).*zz1); ez2p2=exp(1i*k(2).*zz2); ez3p2=exp(1i*k(3).*zz3);
ez4p2=exp(1i*k(4).*zz4); ez5p2=exp(1i*k(5).*zz5); ez6p2=exp(1i*k(6).*zz6);
ez7p2=exp(1i*k(7).*zz7); ez8p2=exp(1i*k(8).*zz8); ez9p2=exp(1i*k(9).*zz9);
ez10p2=exp(1i*k(10).*zz10);ez11p2=exp(1i*k(11).*zz11);ez12p2=exp(1i*k(12).*zz12);
ez13p2=exp(1i*k(13).*zz13);ez14p2=exp(1i*k(14).*zz14);ez15p2=exp(1i*k(15).*zz15);
ez16p2=exp(1i*k(16).*zz16);ez17p2=exp(1i*k(17).*zz17);ez18p2=exp(1i*k(18).*zz18);
ez19p2=exp(1i*k(19).*zz19);ez20p2=exp(1i*k(20).*zz20);ez21p2=exp(1i*k(21).*zz21);
ez22p2=exp(1i*k(22).*zz22);ez23p2=exp(1i*k(23).*zz23);ez24p2=exp(1i*k(24).*zz24);
ez25p2=exp(1i*k(25).*zz25);ez26p2=exp(1i*k(26).*zz26);ez27p2=exp(1i*k(27).*zz27);
ez28p2=exp(1i*k(28).*zz28);ez29p2=exp(1i*k(29).*zz29);ez30p2=exp(1i*k(30).*zz30);
ez31p2=exp(1i*k(31).*zz31);ez32p2=exp(1i*k(32).*zz32);ez33p2=exp(1i*k(33).*zz33);
ez34p2=exp(1i*k(34).*zz34);ez35p2=exp(1i*k(35).*zz35);ez36p2=exp(1i*k(36).*zz36);
ez37p2=exp(1i*k(37).*zz37);ez38p2=exp(1i*k(38).*zz38);ez39p2=exp(1i*k(39).*zz39);
ez40p2=exp(1i*k(40).*zz40);ez41p2=exp(1i*k(41).*zz41);ez42p2=exp(1i*k(42).*zz42);
ez43p2=exp(1i*k(43).*zz43);ez44p2=exp(1i*k(44).*zz44);ez45p2=exp(1i*k(45).*zz45);
ez46p2=exp(1i*k(46).*zz46);ez47p2=exp(1i*k(47).*zz47);ez48p2=exp(1i*k(48).*zz48);
ez49p2=exp(1i*k(49).*zz49);ez50p2=exp(1i*k(50).*zz50);%ez51p2=exp(1i*k(51).*zz51);
%part 3
ez1p3=exp(((-2*B*p(1)+1i*k(1)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(1).*A1*w02))).^2....
-p2(1)/w02);
ez2p3=exp(((-2*B*p(2)+1i*k(2)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(2).*A2*w02))).^2....
-p2(2)/w02);
ez3p3=exp(((-2*B*p(3)+1i*k(3)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(3).*A3*w02))).^2....
-p2(3)/w02);
ez4p3=exp(((-2*B*p(4)+1i*k(4)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(4).*A4*w02))).^2....
-p2(4)/w02);
ez5p3=exp(((-2*B*p(5)+1i*k(5)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(5).*A5*w02))).^2....
-p2(5)/w02);
ez6p3=exp(((-2*B*p(6)+1i*k(6)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(6).*A6*w02))).^2....
-p2(6)/w02);
ez7p3=exp(((-2*B*p(7)+1i*k(7)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(7).*A7*w02))).^2....
-p2(7)/w02);
ez8p3=exp(((-2*B*p(8)+1i*k(8)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(8).*A8*w02))).^2....
-p2(8)/w02);
ez9p3=exp(((-2*B*p(9)+1i*k(9)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(9).*A9*w02))).^2....
-p2(9)/w02);
ez10p3=exp(((-2*B*p(10)+1i*k(10)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(10).*A10*w02))).^2....
-p2(10)/w02);
ez11p3=exp(((-2*B*p(11)+1i*k(11)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(11).*A11*w02))).^2....
-p2(11)/w02);
ez12p3=exp(((-2*B*p(12)+1i*k(12)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(12).*A12*w02))).^2....
-p2(12)/w02);
ez13p3=exp(((-2*B*p(13)+1i*k(13)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(13).*A13*w02))).^2....
-p2(13)/w02);
ez14p3=exp(((-2*B*p(14)+1i*k(14)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(14).*A14*w02))).^2....
-p2(14)/w02);
ez15p3=exp(((-2*B*p(15)+1i*k(15)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(15).*A15*w02))).^2....
-p2(15)/w02);
ez16p3=exp(((-2*B*p(16)+1i*k(16)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(16).*A16*w02))).^2....
-p2(16)/w02);
ez17p3=exp(((-2*B*p(17)+1i*k(17)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(17).*A17*w02))).^2....
-p2(17)/w02);
ez18p3=exp(((-2*B*p(18)+1i*k(18)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(18).*A18*w02))).^2....
-p2(18)/w02);
ez19p3=exp(((-2*B*p(19)+1i*k(19)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(19).*A19*w02))).^2....
-p2(19)/w02);
ez20p3=exp(((-2*B*p(20)+1i*k(20)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(20).*A20*w02))).^2....
-p2(20)/w02);
ez21p3=exp(((-2*B*p(21)+1i*k(21)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(21).*A21*w02))).^2....
-p2(21)/w02);
ez22p3=exp(((-2*B*p(22)+1i*k(22)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(22).*A22*w02))).^2....
-p2(22)/w02);
ez23p3=exp(((-2*B*p(23)+1i*k(23)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(23).*A23*w02))).^2....
-p2(23)/w02);
ez24p3=exp(((-2*B*p(24)+1i*k(24)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(24).*A24*w02))).^2....
-p2(24)/w02);
ez25p3=exp(((-2*B*p(25)+1i*k(25)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(25).*A25*w02))).^2....
-p2(25)/w02);
ez26p3=exp(((-2*B*p(26)+1i*k(26)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(26).*A26*w02))).^2....
-p2(26)/w02);
ez27p3=exp(((-2*B*p(27)+1i*k(27)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(27).*A27*w02))).^2....
-p2(27)/w02);
ez28p3=exp(((-2*B*p(28)+1i*k(28)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(28).*A28*w02))).^2....
-p2(28)/w02);
ez29p3=exp(((-2*B*p(29)+1i*k(29)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(29).*A29*w02))).^2....
-p2(29)/w02);
ez30p3=exp(((-2*B*p(30)+1i*k(30)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(30).*A30*w02))).^2....
-p2(30)/w02);
ez31p3=exp(((-2*B*p(31)+1i*k(31)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(31).*A31*w02))).^2....
-p2(31)/w02);
ez32p3=exp(((-2*B*p(32)+1i*k(32)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(32).*A32*w02))).^2....
-p2(32)/w02);
ez33p3=exp(((-2*B*p(33)+1i*k(33)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(33).*A33*w02))).^2....
-p2(33)/w02);
ez34p3=exp(((-2*B*p(34)+1i*k(34)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(34).*A34*w02))).^2....
-p2(34)/w02);
ez35p3=exp(((-2*B*p(35)+1i*k(35)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(35).*A35*w02))).^2....
-p2(35)/w02);
ez36p3=exp(((-2*B*p(36)+1i*k(36)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(36).*A36*w02))).^2....
-p2(36)/w02);
ez37p3=exp(((-2*B*p(37)+1i*k(37)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(37).*A37*w02))).^2....
-p2(37)/w02);
ez38p3=exp(((-2*B*p(38)+1i*k(38)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(38).*A38*w02))).^2....
-p2(38)/w02);
ez39p3=exp(((-2*B*p(39)+1i*k(39)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(39).*A39*w02))).^2....
-p2(39)/w02);
ez40p3=exp(((-2*B*p(40)+1i*k(40)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(40).*A40*w02))).^2....
-p2(40)/w02);
ez41p3=exp(((-2*B*p(41)+1i*k(41)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(41).*A41*w02))).^2....
-p2(41)/w02);
ez42p3=exp(((-2*B*p(42)+1i*k(42)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(42).*A42*w02))).^2....
-p2(42)/w02);
ez43p3=exp(((-2*B*p(43)+1i*k(43)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(43).*A43*w02))).^2....
-p2(43)/w02);
ez44p3=exp(((-2*B*p(44)+1i*k(44)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(44).*A44*w02))).^2....
-p2(44)/w02);
ez45p3=exp(((-2*B*p(45)+1i*k(45)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(45).*A45*w02))).^2....
-p2(45)/w02);
ez46p3=exp(((-2*B*p(46)+1i*k(46)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(46).*A46*w02))).^2....
-p2(46)/w02);
ez47p3=exp(((-2*B*p(47)+1i*k(47)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(47).*A47*w02))).^2....
-p2(47)/w02);
ez48p3=exp(((-2*B*p(48)+1i*k(48)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(48).*A48*w02))).^2....
-p2(48)/w02);
ez49p3=exp(((-2*B*p(49)+1i*k(49)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(49).*A49*w02))).^2....
-p2(49)/w02);
ez50p3=exp(((-2*B*p(50)+1i*k(50)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(50).*A50*w02))).^2....
-p2(50)/w02);
%ez51p3=exp(((-2*B*p(51)+1i*k(51)*w02.*xx)./sqrt(2*B*w02.*(2*B-1i*k(51).*A51*w02))).^2....
%-p2(51)/w02);
%----------------------------------------------------------------------------------%
Eu1=ez1p1.*ez1p2.*ez1p3; Eu2=ez2p1.*ez2p2.*ez2p3; Eu3=ez3p1.*ez3p2.*ez3p3;
Eu4=ez4p1.*ez4p2.*ez4p3; Eu5=ez5p1.*ez5p2.*ez5p3; Eu6=ez6p1.*ez6p2.*ez6p3;
Eu7=ez7p1.*ez7p2.*ez7p3; Eu8=ez8p1.*ez8p2.*ez8p3; Eu9=ez9p1.*ez9p2.*ez9p3;
Eu10=ez10p1.*ez10p2.*ez10p3;Eu11=ez11p1.*ez11p2.*ez11p3;Eu12=ez12p1.*ez12p2.*ez12p3;
Eu13=ez13p1.*ez13p2.*ez13p3;Eu14=ez14p1.*ez14p2.*ez14p3;Eu15=ez15p1.*ez15p2.*ez15p3;
Eu16=ez16p1.*ez16p2.*ez16p3;Eu17=ez17p1.*ez17p2.*ez17p3;Eu18=ez18p1.*ez18p2.*ez18p3;
Eu19=ez19p1.*ez19p2.*ez19p3;Eu20=ez20p1.*ez20p2.*ez20p3;Eu21=ez21p1.*ez21p2.*ez21p3;
Eu22=ez22p1.*ez22p2.*ez22p3;Eu23=ez23p1.*ez23p2.*ez23p3;Eu24=ez24p1.*ez24p2.*ez24p3;
Eu25=ez25p1.*ez25p2.*ez25p3;Eu26=ez26p1.*ez26p2.*ez26p3;Eu27=ez27p1.*ez27p2.*ez27p3;
Eu28=ez28p1.*ez28p2.*ez28p3;Eu29=ez29p1.*ez29p2.*ez29p3;Eu30=ez30p1.*ez30p2.*ez30p3;
Eu31=ez31p1.*ez31p2.*ez31p3;Eu32=ez32p1.*ez32p2.*ez32p3;Eu33=ez33p1.*ez33p2.*ez33p3;
Eu34=ez34p1.*ez34p2.*ez34p3;Eu35=ez35p1.*ez35p2.*ez35p3;Eu36=ez36p1.*ez36p2.*ez36p3;
Eu37=ez37p1.*ez37p2.*ez37p3;Eu38=ez38p1.*ez38p2.*ez38p3;Eu39=ez39p1.*ez39p2.*ez39p3;
Eu40=ez40p1.*ez40p2.*ez40p3;Eu41=ez41p1.*ez41p2.*ez41p3;Eu42=ez42p1.*ez42p2.*ez42p3;
Eu43=ez43p1.*ez43p2.*ez43p3;Eu44=ez44p1.*ez44p2.*ez44p3;Eu45=ez45p1.*ez45p2.*ez45p3;
Eu46=ez46p1.*ez46p2.*ez46p3;Eu47=ez47p1.*ez47p2.*ez47p3;Eu48=ez48p1.*ez48p2.*ez48p3;
Eu49=ez49p1.*ez49p2.*ez49p3;Eu50=ez50p1.*ez50p2.*ez50p3;%Eu51=ez51p1.*ez51p2.*ez51p3;
%-------fftshift---------%
FF1=ifftshift(fft(fftshift(Eu1.*exp(1i*delta(1).*u*dxx))));
FF2=ifftshift(fft(fftshift(Eu2.*exp(1i*delta(2).*u*dxx))));
FF3=ifftshift(fft(fftshift(Eu3.*exp(1i*delta(3).*u*dxx))));
FF4=ifftshift(fft(fftshift(Eu4.*exp(1i*delta(4).*u*dxx))));
FF5=ifftshift(fft(fftshift(Eu5.*exp(1i*delta(5).*u*dxx))));
FF6=ifftshift(fft(fftshift(Eu6.*exp(1i*delta(6).*u*dxx))));
FF7=ifftshift(fft(fftshift(Eu7.*exp(1i*delta(7).*u*dxx))));
FF8=ifftshift(fft(fftshift(Eu8.*exp(1i*delta(8).*u*dxx))));
FF9=ifftshift(fft(fftshift(Eu9.*exp(1i*delta(9).*u*dxx))));
FF10=ifftshift(fft(fftshift(Eu10.*exp(1i*delta(10).*u*dxx))));
FF11=ifftshift(fft(fftshift(Eu11.*exp(1i*delta(11).*u*dxx))));
FF12=ifftshift(fft(fftshift(Eu12.*exp(1i*delta(12).*u*dxx))));
FF13=ifftshift(fft(fftshift(Eu13.*exp(1i*delta(13).*u*dxx))));
FF14=ifftshift(fft(fftshift(Eu14.*exp(1i*delta(14).*u*dxx))));
FF15=ifftshift(fft(fftshift(Eu15.*exp(1i*delta(15).*u*dxx))));
FF16=ifftshift(fft(fftshift(Eu16.*exp(1i*delta(16).*u*dxx))));
FF17=ifftshift(fft(fftshift(Eu17.*exp(1i*delta(17).*u*dxx))));
FF18=ifftshift(fft(fftshift(Eu18.*exp(1i*delta(18).*u*dxx))));
FF19=ifftshift(fft(fftshift(Eu19.*exp(1i*delta(19).*u*dxx))));
FF20=ifftshift(fft(fftshift(Eu20.*exp(1i*delta(20).*u*dxx))));
FF21=ifftshift(fft(fftshift(Eu21.*exp(1i*delta(21).*u*dxx))));
FF22=ifftshift(fft(fftshift(Eu22.*exp(1i*delta(22).*u*dxx))));
FF23=ifftshift(fft(fftshift(Eu23.*exp(1i*delta(23).*u*dxx))));
FF24=ifftshift(fft(fftshift(Eu24.*exp(1i*delta(24).*u*dxx))));
FF25=ifftshift(fft(fftshift(Eu25.*exp(1i*delta(25).*u*dxx))));
FF26=ifftshift(fft(fftshift(Eu26.*exp(1i*delta(26).*u*dxx))));
FF27=ifftshift(fft(fftshift(Eu27.*exp(1i*delta(27).*u*dxx))));
FF28=ifftshift(fft(fftshift(Eu28.*exp(1i*delta(28).*u*dxx))));
FF29=ifftshift(fft(fftshift(Eu29.*exp(1i*delta(29).*u*dxx))));
FF30=ifftshift(fft(fftshift(Eu30.*exp(1i*delta(30).*u*dxx))));
FF31=ifftshift(fft(fftshift(Eu31.*exp(1i*delta(31).*u*dxx))));
FF32=ifftshift(fft(fftshift(Eu32.*exp(1i*delta(32).*u*dxx))));
FF33=ifftshift(fft(fftshift(Eu33.*exp(1i*delta(33).*u*dxx))));
FF34=ifftshift(fft(fftshift(Eu34.*exp(1i*delta(34).*u*dxx))));
FF35=ifftshift(fft(fftshift(Eu35.*exp(1i*delta(35).*u*dxx))));
FF36=ifftshift(fft(fftshift(Eu36.*exp(1i*delta(36).*u*dxx))));
FF37=ifftshift(fft(fftshift(Eu37.*exp(1i*delta(37).*u*dxx))));
FF38=ifftshift(fft(fftshift(Eu38.*exp(1i*delta(38).*u*dxx))));
FF39=ifftshift(fft(fftshift(Eu39.*exp(1i*delta(39).*u*dxx))));
FF40=ifftshift(fft(fftshift(Eu40.*exp(1i*delta(40).*u*dxx))));
FF41=ifftshift(fft(fftshift(Eu41.*exp(1i*delta(41).*u*dxx))));
FF42=ifftshift(fft(fftshift(Eu42.*exp(1i*delta(42).*u*dxx))));
FF43=ifftshift(fft(fftshift(Eu43.*exp(1i*delta(43).*u*dxx))));
FF44=ifftshift(fft(fftshift(Eu44.*exp(1i*delta(44).*u*dxx))));
FF45=ifftshift(fft(fftshift(Eu45.*exp(1i*delta(45).*u*dxx))));
FF46=ifftshift(fft(fftshift(Eu46.*exp(1i*delta(46).*u*dxx))));
FF47=ifftshift(fft(fftshift(Eu47.*exp(1i*delta(47).*u*dxx))));
FF48=ifftshift(fft(fftshift(Eu48.*exp(1i*delta(48).*u*dxx))));
FF49=ifftshift(fft(fftshift(Eu49.*exp(1i*delta(49).*u*dxx))));
FF50=ifftshift(fft(fftshift(Eu50.*exp(1i*delta(50).*u*dxx))));
%FF51=fft(fftshift(Eu51.*exp(1i*delta(51).*u*dxx)));
for c2=2:2 %c2取不同，可以得到不同观察面处的光场分布
z=dz(c2);
%--------H(fx)-------%
Hfx1=exp(1i*k(1)*z).*exp(-1i*pi*r(1)*z.*fx.*fx);
Hfx2=exp(1i*k(2)*z).*exp(-1i*pi*r(2)*z.*fx.*fx);
Hfx3=exp(1i*k(3)*z).*exp(-1i*pi*r(3)*z.*fx.*fx);
Hfx4=exp(1i*k(4)*z).*exp(-1i*pi*r(4)*z.*fx.*fx);
Hfx5=exp(1i*k(5)*z).*exp(-1i*pi*r(5)*z.*fx.*fx);
Hfx6=exp(1i*k(6)*z).*exp(-1i*pi*r(6)*z.*fx.*fx);
Hfx7=exp(1i*k(7)*z).*exp(-1i*pi*r(7)*z.*fx.*fx);
Hfx8=exp(1i*k(8)*z).*exp(-1i*pi*r(8)*z.*fx.*fx);
Hfx9=exp(1i*k(9)*z).*exp(-1i*pi*r(9)*z.*fx.*fx);
Hfx10=exp(1i*k(10)*z).*exp(-1i*pi*r(10)*z.*fx.*fx);
Hfx11=exp(1i*k(11)*z).*exp(-1i*pi*r(11)*z.*fx.*fx);
Hfx12=exp(1i*k(12)*z).*exp(-1i*pi*r(12)*z.*fx.*fx);
Hfx13=exp(1i*k(13)*z).*exp(-1i*pi*r(13)*z.*fx.*fx);
Hfx14=exp(1i*k(14)*z).*exp(-1i*pi*r(14)*z.*fx.*fx);
Hfx15=exp(1i*k(15)*z).*exp(-1i*pi*r(15)*z.*fx.*fx);
Hfx16=exp(1i*k(16)*z).*exp(-1i*pi*r(16)*z.*fx.*fx);
Hfx17=exp(1i*k(17)*z).*exp(-1i*pi*r(17)*z.*fx.*fx);
Hfx18=exp(1i*k(18)*z).*exp(-1i*pi*r(18)*z.*fx.*fx);
Hfx19=exp(1i*k(19)*z).*exp(-1i*pi*r(19)*z.*fx.*fx);
Hfx20=exp(1i*k(20)*z).*exp(-1i*pi*r(20)*z.*fx.*fx);
Hfx21=exp(1i*k(21)*z).*exp(-1i*pi*r(21)*z.*fx.*fx);
Hfx22=exp(1i*k(22)*z).*exp(-1i*pi*r(22)*z.*fx.*fx);
Hfx23=exp(1i*k(23)*z).*exp(-1i*pi*r(23)*z.*fx.*fx);
Hfx24=exp(1i*k(24)*z).*exp(-1i*pi*r(24)*z.*fx.*fx);
Hfx25=exp(1i*k(25)*z).*exp(-1i*pi*r(25)*z.*fx.*fx);
Hfx26=exp(1i*k(26)*z).*exp(-1i*pi*r(26)*z.*fx.*fx);
Hfx27=exp(1i*k(27)*z).*exp(-1i*pi*r(27)*z.*fx.*fx);
Hfx28=exp(1i*k(28)*z).*exp(-1i*pi*r(28)*z.*fx.*fx);
Hfx29=exp(1i*k(29)*z).*exp(-1i*pi*r(29)*z.*fx.*fx);
Hfx30=exp(1i*k(30)*z).*exp(-1i*pi*r(30)*z.*fx.*fx);
Hfx31=exp(1i*k(31)*z).*exp(-1i*pi*r(31)*z.*fx.*fx);
Hfx32=exp(1i*k(32)*z).*exp(-1i*pi*r(32)*z.*fx.*fx);
Hfx33=exp(1i*k(33)*z).*exp(-1i*pi*r(33)*z.*fx.*fx);
Hfx34=exp(1i*k(34)*z).*exp(-1i*pi*r(34)*z.*fx.*fx);
Hfx35=exp(1i*k(35)*z).*exp(-1i*pi*r(35)*z.*fx.*fx);
Hfx36=exp(1i*k(36)*z).*exp(-1i*pi*r(36)*z.*fx.*fx);
Hfx37=exp(1i*k(37)*z).*exp(-1i*pi*r(37)*z.*fx.*fx);
Hfx38=exp(1i*k(38)*z).*exp(-1i*pi*r(38)*z.*fx.*fx);
Hfx39=exp(1i*k(39)*z).*exp(-1i*pi*r(39)*z.*fx.*fx);
Hfx40=exp(1i*k(40)*z).*exp(-1i*pi*r(40)*z.*fx.*fx);
Hfx41=exp(1i*k(41)*z).*exp(-1i*pi*r(41)*z.*fx.*fx);
Hfx42=exp(1i*k(42)*z).*exp(-1i*pi*r(42)*z.*fx.*fx);
Hfx43=exp(1i*k(43)*z).*exp(-1i*pi*r(43)*z.*fx.*fx);
Hfx44=exp(1i*k(44)*z).*exp(-1i*pi*r(44)*z.*fx.*fx);
Hfx45=exp(1i*k(45)*z).*exp(-1i*pi*r(45)*z.*fx.*fx);
Hfx46=exp(1i*k(46)*z).*exp(-1i*pi*r(46)*z.*fx.*fx);
Hfx47=exp(1i*k(47)*z).*exp(-1i*pi*r(47)*z.*fx.*fx);
Hfx48=exp(1i*k(48)*z).*exp(-1i*pi*r(48)*z.*fx.*fx);
Hfx49=exp(1i*k(49)*z).*exp(-1i*pi*r(49)*z.*fx.*fx);
Hfx50=exp(1i*k(50)*z).*exp(-1i*pi*r(50)*z.*fx.*fx);
%Hfx51=fftshift(exp(1i*k(51)*z).*exp(-1i*pi*r(51)*z.*fx.*fx));
%------ifft----%
E1=ifftshift(ifft(ifftshift(Hfx1.*FF1))); E2=ifftshift(ifft(ifftshift(Hfx2.*FF2)));
E3=ifftshift(ifft(ifftshift(Hfx3.*FF3))); E4=ifftshift(ifft(ifftshift(Hfx4.*FF4)));
E5=ifftshift(ifft(ifftshift(Hfx5.*FF5))); E6=ifftshift(ifft(ifftshift(Hfx6.*FF6)));
E7=ifftshift(ifft(ifftshift(Hfx7.*FF7))); E8=ifftshift(ifft(ifftshift(Hfx8.*FF8)));
E9=ifftshift(ifft(ifftshift(Hfx9.*FF9))); E10=ifftshift(ifft(ifftshift(Hfx10.*FF10)));
E11=ifftshift(ifft(ifftshift(Hfx11.*FF11)));E12=ifftshift(ifft(ifftshift(Hfx12.*FF12)));
E13=ifftshift(ifft(ifftshift(Hfx13.*FF13)));E14=ifftshift(ifft(ifftshift(Hfx14.*FF14)));
E15=ifftshift(ifft(ifftshift(Hfx15.*FF15)));E16=ifftshift(ifft(ifftshift(Hfx16.*FF16)));
E17=ifftshift(ifft(ifftshift(Hfx17.*FF17)));E18=ifftshift(ifft(ifftshift(Hfx18.*FF18)));
E19=ifftshift(ifft(ifftshift(Hfx19.*FF19)));E20=ifftshift(ifft(ifftshift(Hfx20.*FF20)));
E21=ifftshift(ifft(ifftshift(Hfx21.*FF21)));E22=ifftshift(ifft(ifftshift(Hfx22.*FF22)));
E23=ifftshift(ifft(ifftshift(Hfx23.*FF23)));E24=ifftshift(ifft(ifftshift(Hfx24.*FF24)));
E25=ifftshift(ifft(ifftshift(Hfx25.*FF25)));E26=ifftshift(ifft(ifftshift(Hfx26.*FF26)));
E27=ifftshift(ifft(ifftshift(Hfx27.*FF27)));E28=ifftshift(ifft(ifftshift(Hfx28.*FF28)));
E29=ifftshift(ifft(ifftshift(Hfx29.*FF29)));E30=ifftshift(ifft(ifftshift(Hfx30.*FF30)));
E31=ifftshift(ifft(ifftshift(Hfx31.*FF31)));E32=ifftshift(ifft(ifftshift(Hfx32.*FF32)));
E33=ifftshift(ifft(ifftshift(Hfx33.*FF33)));E34=ifftshift(ifft(ifftshift(Hfx34.*FF34)));
E35=ifftshift(ifft(ifftshift(Hfx35.*FF35)));E36=ifftshift(ifft(ifftshift(Hfx36.*FF36)));
E37=ifftshift(ifft(ifftshift(Hfx37.*FF37)));E38=ifftshift(ifft(ifftshift(Hfx38.*FF38)));
E39=ifftshift(ifft(ifftshift(Hfx39.*FF39)));E40=ifftshift(ifft(ifftshift(Hfx40.*FF40)));
E41=ifftshift(ifft(ifftshift(Hfx41.*FF41)));E42=ifftshift(ifft(ifftshift(Hfx42.*FF42)));
E43=ifftshift(ifft(ifftshift(Hfx43.*FF43)));E44=ifftshift(ifft(ifftshift(Hfx44.*FF44)));
E45=ifftshift(ifft(ifftshift(Hfx45.*FF45)));E46=ifftshift(ifft(ifftshift(Hfx46.*FF46)));
E47=ifftshift(ifft(ifftshift(Hfx47.*FF47)));E48=ifftshift(ifft(ifftshift(Hfx48.*FF48)));
E49=ifftshift(ifft(ifftshift(Hfx49.*FF49)));E50=ifftshift(ifft(ifftshift(Hfx50.*FF50)));
%E51=ifftshift(ifft(Hfx51.*FF51));
%------I(x,z)------%
I1=abs(E1).^2+abs(E2).^2+abs(E3).^2+abs(E4).^2+abs(E5).^2....
+abs(E6).^2+abs(E7).^2+abs(E8).^2+abs(E9).^2+abs(E10).^2....
+abs(E11).^2+abs(E12).^2+abs(E13).^2+abs(E14).^2+abs(E15).^2....
+abs(E16).^2+abs(E17).^2+abs(E18).^2+abs(E19).^2+abs(E20).^2....
+abs(E21).^2+abs(E22).^2+abs(E23).^2+abs(E24).^2+abs(E25).^2....
+abs(E26).^2+abs(E27).^2+abs(E28).^2+abs(E29).^2+abs(E30).^2....
+abs(E31).^2+abs(E32).^2+abs(E33).^2+abs(E34).^2+abs(E35).^2....
+abs(E36).^2+abs(E37).^2+abs(E38).^2+abs(E39).^2+abs(E40).^2....
+abs(E41).^2+abs(E42).^2+abs(E43).^2+abs(E44).^2+abs(E45).^2....
+abs(E46).^2+abs(E47).^2+abs(E48).^2+abs(E49).^2+abs(E50).^2;%....
%+abs(E51).^2;
end
%以下是五十个子光束单一的光强分布曲线和合束光束的分布曲线
figure
plot(abs(E1).^2);hold on;plot(abs(E2).^2);hold on;plot(abs(E3).^2);hold on
plot(abs(E4).^2);hold on;plot(abs(E5).^2);hold on;plot(abs(E6).^2);hold on
plot(abs(E7).^2);hold on;plot(abs(E8).^2);hold on;plot(abs(E9).^2);hold on
plot(abs(E10).^2);hold on;plot(abs(E11).^2);hold on;plot(abs(E12).^2);hold on
plot(abs(E13).^2);hold on;plot(abs(E14).^2);hold on;plot(abs(E15).^2);hold on
plot(abs(E16).^2);hold on;plot(abs(E17).^2);hold on;plot(abs(E18).^2);hold on;
plot(abs(E19).^2);hold on;plot(abs(E20).^2);hold on;plot(abs(E21).^2);hold on;
plot(abs(E22).^2);hold on;plot(abs(E23).^2);hold on;plot(abs(E24).^2);hold on;
plot(abs(E25).^2);hold on;plot(abs(E26).^2);hold on;plot(abs(E27).^2);hold on;
plot(abs(E28).^2);hold on;plot(abs(E29).^2);hold on;plot(abs(E30).^2);hold on;
plot(abs(E31).^2);hold on;plot(abs(E32).^2);hold on;plot(abs(E33).^2);hold on;
plot(abs(E34).^2);hold on;plot(abs(E35).^2);hold on;plot(abs(E36).^2);hold on;
plot(abs(E37).^2);hold on;plot(abs(E38).^2);hold on;plot(abs(E39).^2);hold on;
plot(abs(E40).^2);hold on;plot(abs(E41).^2);hold on;plot(abs(E42).^2);hold on;
plot(abs(E43).^2);hold on;plot(abs(E44).^2);hold on;plot(abs(E45).^2);hold on;
plot(abs(E46).^2);hold on;plot(abs(E47).^2);hold on;plot(abs(E48).^2);hold on;
plot(abs(E50).^2);hold on;plot(abs(E49).^2);hold on;%plot(abs(E42).^2);hold on;
%合束光束光强分布曲线
plot(I1);hold on
%axis([500000,550000,0,25.0])
不好意思，各位大佬，第一次用这个，不知道怎么上传matlab的代码，传了个c语言的，希望不会影响你们的阅读