www.pudn.com > ch9_ex.rar > ex925.m

function ex925 
clear all,clf,hold off 
 
e=1.6e-19;        %Charge of electron, coulomb 
m=9.1e-31;        %Electron mass, kg 
B=[0;0;0.1]*e/m;        %Maganetic field strength, tesla 
E=[0;2e4;0]*e/m;        %Electic field strength, vplt/meter 
h=0.5e-12;      %Time step, second 
v(:,1)=1e5*[-10;2;0.1];     %Initial velocity of electron, meter/sec 
t(1)=0;     %Time initialization, sec 
xyz(:,1)=[0;0;0]; 
 
for i=2:4000 
    t(i)=h*(i-1); 
    k1=vxv_(v(:,i-1),B)+E; 
    k2=vxv_(v(:,i-1)+0.5*h*k1,B)+E; 
    k3=vxv_(v(:,i-1)+0.5*h*k2,B)+E; 
    k4=vxv_(v(:,i-1)+h*k3,B)+E; 
    v(:,i)=v(:,i-1)+h*(k1+2*k2+2*k3+k4)/6; 
    xyz(:,i)=xyz(:,i-1)+0.5*(v(:,i-1)+v(:,i))*h; 
end 
v' 
 
figure(1) 
comet3(xyz(1,:),xyz(2,:),xyz(3,:)) 
 
figure(2) 
plot3(xyz(1,:),xyz(2,:),xyz(3,:)) 
xlabel('x'); ylabel('y'); zlabel('z'); 
view([0,0,1]) 
 
figure(3) 
plot3(v(1,:),v(2,:),v(3,:)) 
xlabel('Vx'); ylabel('Vy'); zlabel('Vz'); 
 
figure(4) 
plot(t,v(1,:),t,v(2,:),t,v(3,:)*100) 
xlabel('t'); ylabel('Velocities'); 
text(t(90),v(1,90),'Vx') 
text(t(1200),v(2,1200),'Vy') 
text(t(500),v(3,500)*100,'Vz*100') 
 
function c=vxv_(a,b) 
  c=[a(2)*b(3)-a(3)*b(2); 
      -a(1)*b(3)+a(3)*b(1); 
      a(1)*b(2)-a(2)*b(1)];