clear
clc
npoint=2000;
w0=0.5e-3;
n1=60;
lam=1550e-9;
l=1;
NA=1.32;
nt=1.518;
R=1;
a=linspace(-n1*w0,n1*w0,npoint);
b=linspace(-n1*w0,n1*w0,npoint);
[A,B]=meshgrid(a,b);
rho=sqrt(A.^2+B.^2);
apt=(rho>=0)&(rho<=0.03);
theta=2*mod((atan2(B,A)+pi),2*pi);
phi=zeros(npoint,npoint);
T1=pi/100; %周期
%%
detax1=5e-3;
detay1=0e-3;
%%
detax2=(2.5*sqrt(2))*1e-3;
detay2=(2.5*sqrt(2))*1e-3;
%%
detax3=0e-3;
detay3=5e-3;
%%
detax4=-(2.5*sqrt(2))*1e-3;
detay4=(2.5*sqrt(2))*1e-3;
%%
detax5=-5e-3;
detay5=0e-3;
%%
detax6=-(2.5*sqrt(2))*1e-3;
detay6=-(2.5*sqrt(2))*1e-3;
%%
detax7=0e-3;
detay7=-5e-3;
%%
detax8=(2.5*sqrt(2))*1e-3;
detay8=-(2.5*sqrt(2))*1e-3;
%%
% detax1=0;detay1=0;detax2=0;detay2=0;
% detax3=0;detay3=0;detax4=0;detay4=0;
% detax5=0;detay5=0;detax6=0;detay6=0;
% detax7=0;detay7=0;detax8=0;detay8=0;
%%
l_1=0;l_2=1;l_3=2;l_4=3;l_5=4;l_6=5;l_7=6;l_8=7;
%%
for m=1:npoint
for n=1:npoint
if (theta(m,n)-floor(theta(m,n)/T1)*T1)>0 && (theta(m,n)-floor(theta(m,n)/T1)*T1)<=T1/8
phi(m,n)=2*pi/lam*NA/R/nt*(A(m,n).*detax1+B(m,n).*detay1)+l_1*theta(m,n);
elseif (theta(m,n)-floor(theta(m,n)/T1)*T1)>T1/8 && (theta(m,n)-floor(theta(m,n)/T1)*T1)<=2*T1/8
phi(m,n)=2*pi/lam*NA/R/nt*(A(m,n).*detax2+B(m,n).*detay2)+l_2*theta(m,n);
elseif (theta(m,n)-floor(theta(m,n)/T1)*T1)>2*T1/8 && (theta(m,n)-floor(theta(m,n)/T1)*T1)<=3*T1/8
phi(m,n)=2*pi/lam*NA/R/nt*(A(m,n).*detax3+B(m,n).*detay3)+l_3*theta(m,n);
elseif (theta(m,n)-floor(theta(m,n)/T1)*T1)>3*T1/8 && (theta(m,n)-floor(theta(m,n)/T1)*T1)<=4*T1/8
phi(m,n)=2*pi/lam*NA/R/nt*(A(m,n).*detax4+B(m,n).*detay4)+l_4*theta(m,n);
elseif (theta(m,n)-floor(theta(m,n)/T1)*T1)>4*T1/8 && (theta(m,n)-floor(theta(m,n)/T1)*T1)<=5*T1/8
phi(m,n)=2*pi/lam*NA/R/nt*(A(m,n).*detax5+B(m,n).*detay5)+l_5*theta(m,n);
elseif (theta(m,n)-floor(theta(m,n)/T1)*T1)>5*T1/8 && (theta(m,n)-floor(theta(m,n)/T1)*T1)<=6*T1/8
phi(m,n)=2*pi/lam*NA/R/nt*(A(m,n).*detax6+B(m,n).*detay6)+l_6*theta(m,n);
elseif (theta(m,n)-floor(theta(m,n)/T1)*T1)>6*T1/8 && (theta(m,n)-floor(theta(m,n)/T1)*T1)<=7*T1/8
phi(m,n)=2*pi/lam*NA/R/nt*(A(m,n).*detax7+B(m,n).*detay7)+l_7*theta(m,n);
elseif (theta(m,n)-floor(theta(m,n)/T1)*T1)>7*T1/8 && (theta(m,n)-floor(theta(m,n)/T1)*T1)<=8*T1/8
phi(m,n)=2*pi/lam*NA/R/nt*(A(m,n).*detax8+B(m,n).*detay8)+l_8*theta(m,n);
end
end
end
phi=phi.*apt;
%%
w0=5e-3;
z0=1;
z=10;
k=2*pi/lam;
w=w0*sqrt(1+(lam*z0/pi/w0^2)^2);
%%
dl=2*n1*w0;
dt=dl/m;
fx=linspace(-0.5/dt,0.5/dt,m); %x方向频率
fy=linspace(-0.5/dt,0.5/dt,m); %y方向频率
fx(1,101)=1/dt/m; %消去0点
fy(1,101)=1/dt/m;
cn=1e-14; %湍流强度
[Fx,Fy]=meshgrid(fx,fy);
l0=2e-4; kl=3.3/l0; L0=50;%内外湍流尺度
hh=((2*pi.*Fx).^2+(2*pi.*Fy).^2)/(kl.^2);
hs=(1+1.802.*sqrt(hh)-0.254.*hh.^(-7/12)).*exp(hh);
S=((2*pi.*Fx).^2+(2*pi.*Fy).^2+1/L0.^2).^(-11/6);
ss=(2*pi/2/dl).^2.*k.^2.*2.*pi.*0.033.*cn.*z.*hs.*S;
ee=exp(-(rho/w).^2).*exp(-1i*l_4*theta);% ee=lg_1+1i*lg_1;
%%
R=randn(m,m)+1i*randn(m,m); %随机复矩阵
%功率谱反演法
ph = fftshift(fft2(R.*sqrt(ss))); %等效相位屏
ee1=ee.*exp(1i*abs(ph));
fx=-1/2/dt:1/dl:1/2/dt-1/dl; %x方向频率
[Fx,Fy]=meshgrid(fx,fx);
H=exp(-1i*pi*lam*z*(Fx.^2+Fy.^2));%传递函数
H=fftshift(H);
u11=fft2(fftshift(ee1));
u2=H.*u11; %频域处理
ee=500*ifftshift(ifft2(u2));
%%
% gauss=exp(-(rho/w).^2).*exp(-1i*l_4*theta);
gauss=fftshift(fft2(ee.*exp(1i*phi)));
%%
figure
mesh(abs(gauss))
view(0,90)
colormap hot
20181022.rar_Photonic crystal_matlab_一维光子晶体_传输矩阵_光子晶体
版权申诉
5星 · 超过95%的资源 65 浏览量
2022-07-15
14:51:46
上传
评论 4
收藏 9KB RAR 举报
weixin_42653672
- 粉丝: 93
- 资源: 1万+
评论8