MATLAB Answers

Trouble indexing matrices to simulate a particle splitting

2 views (last 30 days)
Joseph Schmidt
Joseph Schmidt on 3 Aug 2019
Answered: Pujitha Narra on 6 Aug 2019
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
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
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
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
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
%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
t = t + dt; %progress time for ALL species not EACH species
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.


Walter Roberson
Walter Roberson on 4 Aug 2019
"a charged particle (charge=+20 called dimers), after a certain time, will split into a neutral particle and a smaller charged particle (charge=-10 called monomer)."
In order for that to work, then because of conservation of charge, you would have to produce two monomers and they would have to travel backwards in time. You should be expecting self-interference from future particle production, and you need to model it all as a quantum system.
Joseph Schmidt
Joseph Schmidt on 4 Aug 2019
The quantity charge is actually suppose to be charge/mass. And even then the values (20 and -10) were simply placeholder values to test the code. No quantum system needed for this case, just a big particle splitting into two lighter particles but maintaining charge.
Any idea on what's going on with the index overlapping problem?

Sign in to comment.

Answers (1)

Pujitha Narra
Pujitha Narra on 6 Aug 2019
Can you provide the values for npmax,npeaks,Epx,Epy,Epz,dt to help reproduce the problem?


Sign in to comment.

Community Treasure Hunt

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

Start Hunting!