slippage.m

来自「自由电子激光」· M 代码 · 共 51 行

betaz=0.990;%?????
Nb=400;
Np=100;
N=500;
Nw=100;
ru=0.02;
deltaz1=4*pi*ru*(1-betaz)/betaz;%??????%
deltaz2=4*pi*ru%*(1-betaz)/betaz;%??????%
G=4*pi*ru*Nw;
S=Nw/Nb;
K=S/G;
delta=0.09;

theta=zeros(Np,80000);
p=zeros(Np,80000);
A=zeros(N,Np);
A(1,:)=0.001+i*0.001;%???????????

Random=rand(1,100)*pi;
for j=1:1:800
    for k=1:1:100     
    l=(j-1)*100+k;
    theta(1,l)=Random(k)+(j-1)*pi;
    end
end

for j=1:1:Np-1
j
    theta(j+1,:)=theta(j,:)+deltaz2*p(j,:);
    for k=1:1:80000
        N1=ceil(k/200);
        p(j+1,k)=p(j,k)-deltaz2*(A(N1,j)*exp(i*theta(j,k))+(A(N1,j)*exp(i*theta(j,k)))');
    end
    for m=1:1:N-1
        if m<=400
            xie=1;
            A(m+1,j+1)=A(m,j+1)+deltaz1*xie*(mean(cos(theta(j+1,1+100*2*(m-1):100*2*m))-i*sin(theta(j+1,1+100*2*(m-1):100*2*m)))+i*delta*A(m,j+1));
        else
            xie=0;
            A(m+1,j+1)=A(m,j+1)+deltaz1*i*delta*A(m,j+1);
        end
    end
end
for k=1:1:N 
B(k)=A(k,Np)*(A(k,Np))';
end
plot(B)
for j=1:1:400
C(j)=mean(p(100,1+200*(j-1):200*j));
end
plot(C)