a11111111.m

多模光纤模式分布 采用了一种新的方法 欢迎指正

clc;
clear;
xd=1;     %空间高斯脉宽x向:um
yd=1;       %空间高斯脉宽y向
width_x=20;   %窗口宽度:um
width_y=20;   %窗口宽度:um
E0=1;
xa=-width_x;      %%um-x窗口宽度,可改变模式数量
xb=width_x;
ya=-width_y;      %%um-y窗口宽度,可改变模式数量
yb=width_y;
n=100;        %空间采样点数
xh=(xb-xa)/n; %空间步长
yh=(yb-ya)/n;
x=xa:xh:xb;        %%um--x向
y=ya:yh:yb;        %%um--y向
lamida=1.55;     %%um
C=3e8;     %%%光速

[x,y]=meshgrid(xa:xh:xb,ya:yh:yb);    %将x,y向坐标转变成数组，以便于计算二维高斯和绘图
E=E0.*exp(-(x.^2+y.^2)/(xd^2+yd^2));  %假设透镜前的空间场分布

% figure(1)
% mesh(x,y,E);
% xlabel('x向空间距离（um）');
% ylabel('y向空间距离（um）');
% zlabel('归一化光源场强度');
% title('光源强度分布');
% axis([xa,xb, ya,yb, 0,E0]);
% 
f=E0*xh*fftshift(fft2(E));
F=abs(f).^2;

fx_max=1/xh;
fy_max=1/yh;
fx=-fx_max/2:fx_max/n:fx_max/2;   %直接傅氏变换横坐标fx,单位1/um;
fy=-fy_max/2:fy_max/n:fy_max/2;   %直接傅氏变换横坐标fy,单位1/um;

cosine_x=lamida.*fx;   %cos(theta_x)=lamida.*fx，故现横坐标cos(theta_x)
cosine_y=lamida.*fy;   %cos(theta_y)=lamida.*fy，故现横坐标cos(theta_y)
[fx,fy]=meshgrid(-fx_max/2:fx_max/n:fx_max/2);
[cosine_x,cosine_y]=meshgrid(-lamida.*fx_max/2:lamida.*fx_max/n:lamida.*fx_max/2);

figure(2)
mesh(fx,fy,F)
xlabel('频率(1/um)');
ylabel('频率（1/um）');
zlabel('幅值');
title('角谱')
% axis([-fx_max/2,fx_max/2, -fy_max/2,fy_max/2]);%未定义z轴值
view(-37.5,30);   %三维图视角default

figure(3)
mesh(cosine_x,cosine_y,F)
xlabel('cosine（x）');
ylabel('cosine（y）');
zlabel('幅值');
title('角谱');

for i=1:n+1
    Nx1=length(find(cosine_x(i,:)<=-1))+1;%%sin(theta)=-1点的后一点处的数组位置    关心的区间[-1,1]--out of memory,故[-0.6,0.6]!!!!!!!!!!!
    Nx2=length(find(cosine_x(i,:)<=1));%%sin(theta)=1点的前一点处的数组位置
end
for j=1:n+1
    Ny1=length(find(cosine_y(:,j)<=-1))+1;%%sin(theta)=-1点的后一点处的数组位置    关心的区间[-1,1]--out of memory,故[-0.6,0.6]!!!!!!!!!!!
    Ny2=length(find(cosine_y(:,j)<=1));%%sin(theta)=1点的前一点处的数组位置
end
Nx=Nx2-Nx1+1;  %%cosine_x在[-1,1]区间的点个数
Ny=Ny2-Ny1+1;  %%cosine_y在[-1,1]区间的点个数
for i=1:Nx
    for j=1:Ny
        p(i,j)=F(i+Nx1-1,j+Ny1-1);   %%% 将F在[-1,1]区间的能量分布赋给p
    end
end

total_colume=sum(p);    %可得数组的列和
total=sum(total_colume);%%%%%总能量

for i=1:Nx
    for j=1:Ny
        pp(i,j)=p(i,j)./total; %%% 将cosine_y在[-1,1]区间的能量分布权重赋给p赋给cosine_y1
    end
end
for i=1:Nx
    for j=1:Ny
        cosine_x1(i,j)=cosine_x(i+Nx1-1,j+Ny1-1);%%% 将cosine_x在[-1,1]区间的坐标赋给cosine_x1
        cosine_y1(i,j)=cosine_y(i+Nx1-1,j+Ny1-1);%%% 将cosine_y在[-1,1]区间的坐标赋给cosine_y1
    end
end

theta_x=acos(cosine_x1).*180/pi;        %%%%角度谱theta-intensity
theta_y=acos(cosine_y1).*180/pi;        %%%%角度谱theta-intensity

figure(5)
mesh(theta_x,theta_y,pp);
xlabel('theta-x （度）');
ylabel('theta-y （度）');
zlabel('归一化权值');
title('角谱');
% % axis([-lamida.*fx_max/2,+lamida.*fx_max/2, -lamida.*fy_max/2,+lamida.*fy_max/2]);%未定义z轴值
% view(-37.5,30);   %三维图视角default
% view(90,0);