Векторы индукции неоднородного магнитного и электрического полей

>       restart:
 
#=====================================================
move_particle := proc(M,S,T)
#=====================================================
  local num_frames,num_moves_per_frame,n,dt,
        movements,xmin,xmax,ymin,ymax,qmax,
        plot_background,plot_legend,plot_stationary,plot_movements;
 
  num_frames := 20;
  num_moves_per_frame := 5;
  n := num_frames*num_moves_per_frame;
  dt := T / n;
 
  movements := compute_movements(M,S,dt,n);
  compute_bounds(movements,M[2],S,
                 'radius','xmin','xmax','ymin','ymax','qmax');
 
  plot_background := draw_background(xmin,xmax,ymin,ymax);
  plot_legend     := draw_legend(xmin,xmax,ymin,ymax,qmax);
  plot_stationary := draw_stationary(S,radius,qmax);
  plot_movements  := draw_movements(movements,radius,qmax,M[2],num_moves_per_frame);
 
  plots[display]([plot_movements,plot_stationary,plot_legend,plot_background],
                 scaling=constrained,axes=none,
                 title=`Movement of a Charged Particle Under an Electric Field`,
                 titlefont=[TIMES,BOLD,12]);
end:
 
#=====================================================
compute_movements := proc(M,S,dt,n)
#=====================================================
  local m,q,config,movements,i;
 
  m  := M[1];
  q  := M[2];
  config     := M[3..6];
  movements  := config[1..2];
 
  for i from 1 by 1 to n do
    config    := compute_a_movement(q,m,config,S,dt);
    movements := movements,config[1..2];
  od;
 
  return [movements];
end:
 
#=====================================================
compute_a_movement := proc(q,m,C,S,dt)
#=====================================================
  local px,py,vx,vy,Enet,Fnetx,Fnety,
        ax,ay,dpx,dpy,new_px,new_py,new_vx,new_vy;
 
  px := C[1];
  py := C[2];
  vx := C[3];
  vy := C[4];
 
  Enet := compute_electric_field(px,py,S);
 
  Fnetx := Enet[1]*q;
  Fnety := Enet[2]*q;
 
  ax := Fnetx/m;
  ay := Fnety/m;
 
  dpx := vx*dt + 1/2*ax*dt*dt;
  dpy := vy*dt + 1/2*ay*dt*dt;
 
  new_px := px + dpx;
  new_py := py + dpy;
 
  new_vx := vx + ax*dt;
  new_vy := vy + ay*dt;
 
  return [new_px,new_py,new_vx,new_vy];
end:
 
#=====================================================
compute_electric_field := proc(px,py,S)
#=====================================================
  local Enetx,Enety,s,dx,dy,r,c;
 
  Enetx := 0.0;
  Enety := 0.0;
 
  for s in S do
    dx := px-s[2];
    dy := py-s[3];
 
    r  := sqrt(dx*dx + dy*dy);
    c  := 8.99E9*s[1] / (r*r*r);
 
    Enetx := Enetx + c * dx;
    Enety := Enety + c * dy;
  od;
 
  return [Enetx,Enety];
end:
 
#=====================================================
compute_bounds := proc(movements,q,S,radiusp,xminp,xmaxp,yminp,ymaxp,qmaxp)
#=====================================================
  local movement,s,side_length,xcenter,ycenter,
        radius,xmin,xmax,ymin,ymax,qmax;
 
  xmin := movements[1][1];
  xmax := movements[1][1];
  ymin := movements[1][2];
  ymax := movements[1][2];
 
  for movement in movements do
    xmin := min(xmin,movement[1]);
    xmax := max(xmax,movement[1]);
    ymin := min(ymin,movement[2]);
    ymax := max(ymax,movement[2]);
  od;
 
  for s in S do
    xmin := min(xmin,s[2]);
    xmax := max(xmax,s[2]);
    ymin := min(ymin,s[3]);
    ymax := max(ymax,s[3]);
  od;
 
  side_length  := max(abs(xmax-xmin),abs(ymax-ymin))/(1-4/30);
  radius       := side_length/30;
  xcenter      := (xmin+xmax)/2;
  ycenter      := (ymin+ymax)/2;
  xmin := xcenter-side_length/2;
  xmax := xcenter+side_length/2;
  ymin := ycenter-side_length/2;
  ymax := ycenter+side_length/2;
 
  qmax := abs(q);
  for s in S do
    qmax := max(qmax,abs(s[1]));
  od;
 
  radiusp := radius;
  xminp   := xmin;
  xmaxp   := xmax;
  yminp   := ymin;
  ymaxp   := ymax;
  qmaxp   := qmax;
 
  return;
end:
 
#=====================================================
draw_background := proc(xmin,xmax,ymin,ymax)
#=====================================================
  return plottools[rectangle]([xmin,ymax],[xmax,ymin],color=black);
end:
 
#=====================================================
draw_legend := proc(xmin,xmax,ymin,ymax,qmax)
#=====================================================
  local k,width,plot_legend,i,c;
 
  k := 18;
  width := (xmax-xmin)/k;
 
  plot_legend := TEXT([xmin-width,ymin-width],`-`,
                      FONT(COURIER,BOLD,11)),
                 TEXT([xmax+width,ymin-width],`+`,
                      FONT(COURIER,BOLD,11));
 
  for i from 1 to k do
    c := get_color(qmax,-qmax+2*qmax/(k-1)*(i-1));
    plot_legend := plot_legend,plottools[rectangle]([xmin+width*(i-1),ymin-width/2],
                                                    [xmin+width*i,    ymin-width/2*3],
                                                    color=c);
  od;
 
  return plot_legend;
end:
 
#=====================================================
draw_stationary := proc(S,radius,qmax)
#=====================================================
  local plot_stationary,s;
 
  plot_stationary := [];
  for s in S do
    plot_stationary := [op(plot_stationary),
                        plottools[disk]([s[2],s[3]],radius,color=get_color(qmax,s[1]))];
  od;
 
  return PLOT(op(plot_stationary));
end:
 
#=====================================================
draw_movements := proc(movements,radius,qmax,q,num_moves_per_frame)
#=====================================================
  local D,plot_movements,i;
 
  D := plottools[disk]([0,0],radius,color=get_color(qmax,q));
  plot_movements := [];
 
  for i from 1 by num_moves_per_frame to nops(movements) do
    plot_movements := [op(plot_movements),
                       [shift_disk(D,movements[i][1],movements[i][2])]];
  od;
 
  return PLOT(ANIMATE(op(plot_movements)));
end:
 
#=====================================================
shift_disk := proc(D,x,y)
#=====================================================
  return POLYGONS(map(xy -> [xy[1]+x,xy[2]+y],op(1,D)),op(2,D));
end:
 
#=====================================================
get_color := proc(qmax,q)
#=====================================================
  local color;
 
  if q > 0 then
    color := COLOR(RGB,1,1-q/qmax,0)
  else
    color := COLOR(RGB,0,1-abs(q)/qmax,1)
  fi;
 
  return color;
end:

>       move_particle(M,S,T);