Trouble indexing matrices to simulate a particle splitting

Question

Joseph Schmidt 2019-8-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/474751-trouble-indexing-matrices-to-simulate-a-particle-splitting

回答： Pujitha Narra 2019-8-6

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.

2 个评论
显示无隐藏无

Walter Roberson 2019-8-4

"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 2019-8-4

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?

请先登录，再进行评论。

请先登录，再回答此问题。