该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

clear,clc

global steps L k a

lamda=input('波长lamda=');

L=input('腔长L=');

a=input('镜长a=');

eps_u=input('精度=');

k=2*pi/lamda;

steps=500;

x=linspace(-a,a,steps);

u_=ones(1,steps);

for m=1:30

for mm=1:steps

u0(mm)=QU(x(mm),u_);

end;

u_=u0/max(abs(u0));

u30(m,:)=u0;

end

flag=1; %自再现标志位

u30_flag=1; %判据矩阵维数标记为

N=30;

while(flag==1)

if max(abs(u30(:,steps/4)))-min(abs(u30(:,steps/4)))

flag=0;

else

N=N+1;

for mm=1:steps

u0(mm)=QU(x(mm),u_);

end;

u_=u0/max(abs(u0));

if u30_flag==31

u30_flag=1;

end

u30( u30_flag,:)=u0;

u30_flag=u30_flag+1;

end

end

N

subplot(2,1,1)

plot(x,abs(u0)/abs(u0(steps/2)))

xlabel('x');ylabel('相对振幅');

angle_u0=angle(u0)/pi*180;

angle_u0=angle_u0-angle_u0(steps/2);

subplot(2,1,2)

plot(x,angle_u0)

xlabel('x');ylabel('相对相位');

function y=QU(x,u)

global steps L k a

x_=linspace(-a,a,steps);

step_length=2*a/(steps-1);

y=sqrt(1i/L*exp(-1i*k*L))*sum(exp(-1i*k/2/L*(-x_+x).^2).*u)*step_length;