2 views (last 30 days)

I'm trying to write a code where a charged particle (charge/mass=+20 called dimers), after a certain time, will split into a neutral particle and a smaller charged particle (charge/mass=-10 called monomer). I split position and velocity into x,y,z component matrices along with q/m and fragmentation time matrices for each species (monomers, dimers, and neutrals).

The way the code works is 1 dimer is supposed to be injected per time step, but whenever the fragmentation time of that dimer is reached it will turn into a monomer and neutral particle. I got the injection and pos/vel update portions to work, but the fragmentation portion is causing a lot of trouble. Every time a fragmentation time is reached I move the dimer data (pos,vel) to the end of the monomer and neutral matrices and move the last dimer data to the destroyed dimer to overwrite it.

However, when I run the code the simulation appears to work correctly but then after about 17 cycles the monomers and neutrals will jump back to the original point where they split and rewrite over the correct pos/vel from before.

The error seems to be in indexing in passing information from the dimer to the neutrals and monomers but I'm not sure why it's happening.

The code is shown below with the main while loop stopping the simulation once a certain number of particles are injected. Inside the while loop are three for loops that inject particles, update vel/pos, and check for fragmentation (this is the loop causing trouble).

f_ind = [0,0,0];

qms0 = [-10,20]; %assign charge/mass values

t=0;

np = 0; %number of particles total

npnew = [1,1]; %!injected particles for each species

%_______________________Main loop__________________________________________

while (np+sum(npnew) <= npmax) %if particles injected is less than max particle limit

%INJECT PARTICLES

for k=1:npeaks %cycle through monomers (1) and dimers (2)

for i=1:npnew(k)

posx(i+f_ind(k),k)=0.0; %x-pos

posy(i+f_ind(k),k)=i*3; %y-pos

posz(i+f_ind(k),k)=k; %z-pos

velx(i+f_ind(k),k) = 0.01;

vely(i+f_ind(k),k) = 0.0;

velz(i+f_ind(k),k) = 0.0;

q_m(i+f_ind(k),k)=qms0(k); %assign q/m

if (k==2) %for dimers

tfrag(i+f_ind(k),k)=t + 0.05; %only set a fragmentation time for dimers

end

end

f_ind(k) = f_ind(k) + npnew(k); %add to index to keep species particle count

np = np + npnew(k); %add particle count of each injected species to total particle count

end

%UPDATE VELOCITY AND POSITION

for k=1:(npeaks+1) %cycle through monomer, dimers, and neutrals

for i=1:f_ind(k) %only cycle up to the last index of each species(saves time not cycling empty spaces)

velx(i,k) = velx(i,k) + Epx*(q_m(i,k))*dt;

vely(i,k) = vely(i,k) + Epy*(q_m(i,k))*dt;

velz(i,k) = velz(i,k) + Epz*(q_m(i,k))*dt;

posx(i,k) = posx(i,k) + velx(i,k)*dt; %x-pos

posy(i,k) = posy(i,k) + vely(i,k)*dt; %y-pos

posz(i,k) = posz(i,k) + velz(i,k)*dt; %z-pos

end

end

%Fragment dimer if it reaches fragmentation time (ERROR IN THIS PORTION)

for i=1:f_ind(2)

if((t >= tfrag(i,2))&&(tfrag(i,2)>0))

%pass dimer information to monomer

posx(f_ind(1)+1,1) = posx(i,2);

posy(f_ind(1)+1,1) = posy(i,2);

posz(f_ind(1)+1,1) = posz(i,2);

velx(f_ind(1)+1,1) = velx(i,2);

vely(f_ind(1)+1,1) = vely(i,2);

velz(f_ind(1)+1,1) = velz(i,2);

%assign monomer charge and update counter

q_m(f_ind(1)+1,1) = -10;

f_ind(1) = f_ind(1) + 1; %update monomer index

%pass dimer information to neutral

posx(f_ind(3)+1,3) = posx(i,2);

posy(f_ind(3)+1,3) = posy(i,2);

posz(f_ind(3)+1,3) = posz(i,2);

velx(f_ind(3)+1,3) = velx(i,2);

vely(f_ind(3)+1,3) = vely(i,2);

velz(f_ind(3)+1,3) = velz(i,2);

%assign neutral charge and update counter

q_m(f_ind(3)+1,3) = 0;

f_ind(3) = f_ind(3) + 1; %update neutral index

%destroy dimer and replace with the last dimer

posx(i,2) = posx(f_ind(2),2);

posy(i,2) = posy(f_ind(2),2);

posz(i,2) = posz(f_ind(2),2);

velx(i,2) = velx(f_ind(2),2);

vely(i,2) = vely(f_ind(2),2);

velz(i,2) = velz(f_ind(2),2);

q_m(i,2) = q_m(f_ind(2),2);

tfrag(i,2) = tfrag(f_ind(2),2);

f_ind(2) = f_ind(2) - 1; %update dimer index

end

end

t = t + dt; %progress time for ALL species not EACH species

end

The picture below shows the dimers being injected then breaking up into neutrals and monomers once their time is reached. The simulation works up to a point, but then it appears to double back and put particles at the dimer fragmentation point.

Pujitha Narra
on 6 Aug 2019

Can you provide the values for npmax,npeaks,Epx,Epy,Epz,dt to help reproduce the problem?

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 2 Comments

## Direct link to this comment

https://nl.mathworks.com/matlabcentral/answers/474751-trouble-indexing-matrices-to-simulate-a-particle-splitting#comment_731627

⋮## Direct link to this comment

https://nl.mathworks.com/matlabcentral/answers/474751-trouble-indexing-matrices-to-simulate-a-particle-splitting#comment_731627

## Direct link to this comment

https://nl.mathworks.com/matlabcentral/answers/474751-trouble-indexing-matrices-to-simulate-a-particle-splitting#comment_731733

⋮## Direct link to this comment

https://nl.mathworks.com/matlabcentral/answers/474751-trouble-indexing-matrices-to-simulate-a-particle-splitting#comment_731733

Sign in to comment.