(*     Module #21:  Motion of a Charged Particle in
       a Magnetic Field                            *)

(*   Program for constant motion of a positive charge in
     a constant magnetic field      *)

Clear[x,y,Vx,Vy];
x[0]=0.; Vx[0]=1.5*(10^7); y[0]=0.; Vy[0]=0.; 
Delta=0.5*(10^-10); Steps =1400;
q=1.6*(10^-19);m=1.7*(10^-27); B=1.;  
Do[Fx = q*Vy[n]*B; Fy = -q*Vx[n]*B;
   Vx[n+1] = Vx[n] + (Fx/m)*Delta;
   Vy[n+1] = Vy[n] + (Fy/m)*Delta; 
   x[n+1] = x[n] + Vx[n+1]*Delta;
   y[n+1] = y[n] + Vy[n+1]*Delta,
{n,0,Steps}]

(*   Make Graph of trajectory     *)
data = Table[{x[n],y[n]}, {n,0,Steps}];
plot1=ListPlot[data, AxesLabel->{"x","y"},AspectRatio->1]

(*    Program for motion of negatively charged particle
      in a constant magnetic field     *)
x[0]=0.; Vx[0]=1.5*(10^7); y[0]=0.; Vy[0]=0.; 
Delta=0.5*(10^-10); Steps =1400;
q=-1.6*(10^-19);m=1.7*(10^-27);B=1.;                          
Do[Fx=q*Vy[n]*B; Fy=-1*q*Vx[n]*B;
   Vx[n+1] = Vx[n] + (Fx/m)*Delta;
   Vy[n+1] = Vy[n] + (Fy/m)*Delta; 
   x[n+1] = x[n] + Vx[n+1]*Delta;
   y[n+1] = y[n] + Vy[n+1]*Delta,
{n,0,Steps}]

(*   Make graph of trajectory    *)
data = Table[{x[n],y[n]}, {n,0,Steps}];
plot2=ListPlot[data, AxesLabel->{"x","y"},AspectRatio->1]

(*     Motion of charge where magnetic field is
       twice as large as before     *)
x[0]=0.; Vx[0]=1.5*(10^7); y[0]=0.; Vy[0]=0.; 
Delta=0.5*(10^-10); Steps =650;
q=-1.6*(10^-19);m=1.7*(10^-27); B=2.;                          
Do[Fx = q*Vy[n]*B; Fy = -q*Vx[n]*B;
   Vx[n+1] = Vx[n] + (Fx/m)*Delta;
   Vy[n+1] = Vy[n] + (Fy/m)*Delta; 
   x[n+1] = x[n] + Vx[n+1]*Delta;
   y[n+1] = y[n] + Vy[n+1]*Delta,
   {n,0,Steps}]

(*     Make graph of tejectory     *)
data = Table[{x[n],y[n]}, {n,0,Steps}];
plot3=ListPlot[data, AxesLabel->{"x","y"},AspectRatio->1]

(*   Combined graph of cases of two different
     magnetic fields      *)
Show[plot2,plot3]

(*   A case where the magnetic field is switched off   *)

Clear[x,y,Vx,Vy];
x[0]=0.; Vx[0]=1.5*(10^7); y[0]=0.; Vy[0]=0.; 
Delta=0.5*(10^-10); Steps =750;
q=1.6*(10^-19);m=1.7*(10^-27); B=1.;  
Do[If[n<665,Fx = q*Vy[n]*B; Fy = -q*Vx[n]*B,Fx=0.;Fy=0.];
   Vx[n+1] = Vx[n] + (Fx/m)*Delta;
   Vy[n+1] = Vy[n] + (Fy/m)*Delta; 
   x[n+1] = x[n] + Vx[n+1]*Delta;
   y[n+1] = y[n] + Vy[n+1]*Delta,
{n,0,Steps}] 
data = Table[{x[n],y[n]}, {n,0,Steps}];
ListPlot[data, AxesLabel->{"x","y"},AspectRatio->1.]