> |
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:
|