function dydt=e_track(t,y,flag,L,E0,P,w)
%模拟电子在光场中的运动，ft是激光的波形函数，L是脉冲宽度
%激光朝z轴正方向传播。P=1线偏振0为圆偏振。w是圆频率
%电子运动方程d^2r/dt^2=-q/m/γ*(E+V×B),em为荷质比,lz为1/γ
%梁彬彬 编写于2011年5月21日
em=1.75882e11;C=2.99792458e8;
lz=(1-(y(2)^2+y(4)^2+y(6)^2)/C^2)^0.5;
ft=exp(-4*log(2)*(y(5)-C*t)^2/L^2);
if P==0
Ex=E0*ft*sqrt(1/2)*cos(w*(y(5)/C-t));
Ey=E0*ft*sqrt(1/2)*sin(w*(y(5)/C-t));
elseif P==1
Ex=E0*ft*cos(w*(y(5)/C-t));
Ey=0;
else
[];
end
Ez=0;
Bx=Ey/C;By=-Ex/C;Bz=0;
dydt=[y(2);
-em*lz*(Ex+y(4)*Bz-y(6)*By);
y(4);
-em*lz*(Ey+y(6)*Bx-y(2)*Bz);
y(6);
-em*lz*(Ez+y(2)*By-y(4)*Bx)];
再编写一个名为Thomson.m文件。存放计算操作指令和绘图命令。
function Thomson(sj,V,L,E0,P,w)
%sj为模拟时间一半。激光朝z轴正方向传播。V是电子初速
%L是脉冲宽度。E0是峰值强度。P=1线偏振0圆偏振。w是圆频率
%命令举例:Thomson(1.2e-13,[0,0,0],2e-5,5e13,0,5e14)
%梁彬彬 编写于2011年5月23日
close all;
tspan=[-sj,sj];y0=[0,V(1),0,V(2),0,V(3)];
options=odeset('RelTol',1e-9,'AbsTol',1e-12);
[t,y]=ode45('e_track',tspan,y0,options,L,E0,P,w);
figure(1);
plot3(y(:,1),y(:,3),y(:,5))
xlabel('x轴');ylabel('y轴');zlabel('z轴');
title('电子轨迹[單位:米]');
grid on;hold on;
%%
x=[-sj:sj/2e3:sj];C=2.99792458e8;
figure(2);
if P==0
z1=E0.*exp(-(C.*x).^2/L^2)*sqrt(1/2).*cos(-w.*x);
z2=E0.*exp(-(C.*x).^2/L^2)*sqrt(1/2).*sin(-w.*x);
plot3(z1,z2,x,'r');
title('光场电分量時間变化圖');
xlabel('Ex');ylabel('Ey');
zlabel('時間[單位:秒]');grid on;
elseif P==1
z1=E0.*exp(-(C.*x).^2/L^2).*cos(-w.*x);
z2=0;
plot(x,z1,'r');
title('Ex時間变化圖');
xlabel('時間[單位:秒]');grid on;
else
[];
end
%%
figure(3);
subplot(2,2,1);
plot(t,y(:,2),'b');
title('速度Vx的時間变化圖');xlabel('時間[單位:秒]');grid on;
subplot(2,2,2);
plot(t,y(:,4),'r');
title('速度Vy的時間变化圖');xlabel('時間[單位:秒]');grid on;
subplot(2,2,3);
plot(t,y(:,6),'g');
title('速度Vz的時間变化圖');xlabel('時間[單位:秒]');grid on;
subplot(2,2,4);
u=sqrt(y(:,2).^2+y(:,4).^2+y(:,6).^2);
plot(t,u,'k');
title('总速度V的時間变化圖');xlabel('時間[單位:秒]');grid on;
命令窗口输入：>> Thomson(1.2e-13,[0,0,0],2e-5,5e13,0,5e14)
得到的模拟图为：
※※※※※※
天地之道，以阴阳二气造化万物。是故易有太极，是生两仪。两仪生四象，四象生八卦。
