function []=Fox()
clear;
clc;
lam=input('波长(nm)=');
turn=input('渡越次数=');
f=input('原始分布(0 方波/1 正弦)=');
lamda=lam*10^(-9);
L=100*lamda;
a=25*lamda;
k=2*pi/lamda;
N=500;%将积分区间分成 N 小段
x=linspace(-a,a,N+1);%x 轴上等间距插入 N+1 个点
if f==0
u=ones(1,N+1);%初始条件
else
u=sin(pi*x/a);
end
for m=0:turn
A=abs(u);%相对振幅
phase=angle(u)/pi*180;%转换成角度
phase=phase-phase((N+2)/2);%相对于中心位置处的相位差
%---------------画图部分---------------%
subplot(2,1,1),plot(x,A,'b');
xlabel('x/m'),ylabel('相对振幅'),title(['渡越次数：' num2str(m)]);
box on,grid on;
subplot(2,1,2),plot(x,phase,'r');
xlabel('x/m'),ylabel('相对相位/度'),title(['相对相位分布']);
box on,grid on;
drawnow;
u=calc(u);
end
%------------函数部分------------%
function [u2] = calc(u1)
u2=zeros(1,N+1);
for s=1:N+1
for t=0:N
u2(s)=u2(s)+sqrt(1i/lamda/L*exp(-1i*k*L))*(2*a/N)*exp(-1i*k/2/L*(2*a/N*(s-t))^2)*u1(t+1);
end
end
u2=u2/max(abs(u2));%归一化
end
end
