>	restart; with(plots): with(plottools):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

>	Rsys := 10.: m1 := 30.: m2 := 16.: Npart := 1000: Pimpact := 4: H3 := 32.: Rmin := Rsys/10: G:=1: Nt:=300: dt:=0.1: eps:=Rsys/50:

>	RAN := rand(-100..100):

Condiciones Iniciales

>

for i from 1 to Npart do   rtt:=0.;   for k while (rtt > Rsys or rtt <Rmin) do     xx:=(RAN()/100)*Rsys;     yy:=(RAN()/100)*Rsys;     zz:=0.;     rtt := sqrt(xx^2 + yy^2 + zz^2);   end do;   r3[i,1]:=xx;   r3[i,2]:=yy;   r3[i,3]:=0.;   v3[i,1]:= -sqrt(m1/rtt) * r3[i,2]/rtt;   v3[i,2]:= +sqrt(m1/rtt) * r3[i,1]/rtt;   v3[i,3]:=0. end do: r1[1,1] := 0.: v1[1,1]:=0: r1[1,2] := 0.: v1[1,2]:=0: r1[1,3] := 0.: v1[1,3]:=0: r2[1,1] := Pimpact: v2[1,1]:=0: r2[1,2] := 0.:      v2[1,2]:=0: r2[1,3] := H3:      v2[1,3]:=-4.5:

>

plt_1:=plot([seq([r3[i,1],r3[i,2]],i=1..Npart)],x=-Rsys-1..Rsys+1, y=-Rsys-1..Rsys+1, style=point, color=blue): cw_1:= disk ([r1[1,1],r1[1,2]], 0.25, color=red): display({plt_1, cw_1},axes=boxed,scaling=constrained,title=`t=0`); pltz_1:=plot([seq([r3[i,1],r3[i,3]],i=1..Npart)],x=-Rsys-1..Rsys+1, y=33..-5, style=point, color=blue): pert_1:= disk ([r2[1,1], r2[1,3]], 0.25, color=yellow): cwz_1:= disk ([r1[1,1],r1[1,3]],.45, color=red): display({pltz_1, pert_1, cwz_1},axes=boxed,scaling=constrained,title=`t=0`);

>

Aceleracion := proc (i, ac, pos)   local r13,r23,acel,r21;   global r1,r2,r3,G,m1,m2,eps;   if (ac=3) then     r13:=sqrt((pos[i,1]-r1[1,1])^2 + (pos[i,2]-r1[1,2])^2 + (pos[i,3]-r1[1,3])^2);     r23:=sqrt((pos[i,1]-r2[1,1])^2 + (pos[i,2]-r2[1,2])^2 + (pos[i,3]-r2[1,3])^2);     acel[1,1] := -G*m1*(pos[i,1]-r1[1,1])/(r13+eps)^3 - G*m2*(pos[i,1]-r2[1,1])/(r23+eps)^3;     acel[1,2] := -G*m1*(pos[i,2]-r1[1,2])/(r13+eps)^3 - G*m2*(pos[i,2]-r2[1,2])/(r23+eps)^3;     acel[1,3] := -G*m1*(pos[i,3]-r1[1,3])/(r13+eps)^3 - G*m2*(pos[i,3]-r2[1,3])/(r23+eps)^3;     return(acel);   end if;   if (ac=2) then     r21:=sqrt((pos[i,1]-r1[1,1])^2 + (pos[i,2]-r1[1,2])^2 + (pos[i,3]-r1[1,3])^2);     acel[1,1] := -G*m1*(pos[1,1]-r1[1,1])/(r21+eps)^3;     acel[1,2] := -G*m1*(pos[1,2]-r1[1,2])/(r21+eps)^3;     acel[1,3] := -G*m1*(pos[1,3]-r1[1,3])/(r21+eps)^3;     return(acel);   end if;   if (ac=1) then     r21:=sqrt((r2[i,1]-pos[1,1])^2 + (r2[i,2]-pos[1,2])^2 + (r2[i,3]-pos[1,3])^2);     acel[1,1] := G*m2*(r2[1,1]-pos[1,1])/(r21+eps)^3;     acel[1,2] := G*m2*(r2[1,2]-pos[1,2])/(r21+eps)^3;     acel[1,3] := G*m2*(r2[1,3]-pos[1,3])/(r21+eps)^3;     return(acel);   end if; end proc:

i = i-esima particula m3
ac = indicador para la particula que se esta moviendo
dt = paso de tiempo
vel = velocidad
pos = posicion

>

Sapito := proc (pp, ac, dtt, vel, pos)   local acc, v12, rr2, vv2;   acc := Aceleracion (pp, ac, pos);   v12[1,1]:= vel[pp,1] + dtt*acc[1,1]/2;   v12[1,2]:= vel[pp,2] + dtt*acc[1,2]/2;   v12[1,3]:= vel[pp,3] + dtt*acc[1,3]/2;   rr2[1,1]:= pos[pp,1] + dtt*v12[1,1];   rr2[1,2]:= pos[pp,2] + dtt*v12[1,2];   rr2[1,3]:= pos[pp,3] + dtt*v12[1,3];   acc:= Aceleracion (1, ac, rr2);   vv2[1,1]:= v12[1,1] + dtt*acc[1,1];   vv2[1,2]:= v12[1,2] + dtt*acc[1,2];   vv2[1,3]:= v12[1,3] + dtt*acc[1,3];   return (rr2,vv2); end proc:

>

for i from 1 to Nt do   for k from 1 to Npart do     (tmp1,tmp2):=Sapito(k, 3, dt, v3, r3);     r3[k,1]:=tmp1[1,1]; v3[k,1]:=tmp2[1,1];     r3[k,2]:=tmp1[1,2]; v3[k,2]:=tmp2[1,2];     r3[k,3]:=tmp1[1,3]; v3[k,3]:=tmp2[1,3];   end do;   (tmp1,tmp2):=Sapito(1, 2, dt, v2, r2);   r2[1,1]:=tmp1[1,1]; v2[1,1]:=tmp2[1,1];   r2[1,2]:=tmp1[1,2]; v2[1,2]:=tmp2[1,2];   r2[1,3]:=tmp1[1,3]; v2[1,3]:=tmp2[1,3];   (tmp1,tmp2):=Sapito(1, 1, dt, v1, r1);   r1[1,1]:=tmp1[1,1]; v1[1,1]:=tmp2[1,1];   r1[1,2]:=tmp1[1,2]; v1[1,2]:=tmp2[1,2];   r1[1,3]:=tmp1[1,3]; v1[1,3]:=tmp2[1,3];   if ((i mod 5) = 0) then     print('Pasos',i);     cw[i]:= disk ([r1[1,1],r1[1,2]], 0.25, color=red):     plt[i]:=plot([seq([r3[ii,1],r3[ii,2]], ii=1..Npart)],x=-15..25, y=-20..20, style=point,title=convert(i,string), color=blue):     anima[i]:=display(cw[i],plt[i]);   end if;   if ((i mod 5) = 0) then      pltz[i]:=plot([seq([r3[i,1],r3[i,3]],i=1..Npart)],x=-15..25, y=32..-25, style=point, color=blue, title=convert(i,string)):      pert[i]:= disk ([r2[1,1], r2[1,3]], 0.25, color=yellow):      cwz[i]:= disk ([r1[1,1],r1[1,3]],.45, color=red):      animaz[i]:=display({pltz[i], pert[i], cwz[i]},axes=boxed,scaling=constrained);    end if; end do:

>	display({seq(animaz[i*5],i=1..30)}, axes=boxed, insequence=true, scaling=constrained); display({seq(anima[i*5],i=1..30)}, axes=boxed, insequence=true, scaling=constrained);

>

>