%% Simulation of an explosion % In this script we simulate the autocatalytic instability leading to a % nuclear explosion. Simplistically, we will take the neutrons to be point % particles which are diffusing through a fissile material. At each time % step they take a random step (diffusion) in the x, y, and z directions. % The length of this step is a normally distributed random variable with % standard deviation (2*D*dt)^.5 where D is the diffusion coefficient. % Because we've rendered everything dimensionless, the diffusion % coefficient is just 1. They also have a probability of duplicating % themselves if a random number is less than dt*alpha (which is a parameter % much less than one). Particles which leave the sphere of dimensionless % radius of one are eliminated from the list of particles. % % We shall use a value of alpha which is a bit greater than the critical % value of pi^2 so that it blows up (eventually). We start with our % particles distributed uniformly over the sphere. Because most of them % are close to the edges, the total number initially decreases a bit (a % uniform distribution is not the most unstable mode). Eventually, % however, we wind up with the distribution corresponding to the most % unstable mode and the number of particles then increases exponentially in % time. % % Note that while this is written for a specific application (a bomb), % similar approaches can be used for many other autocatalytic problems: % viral spreading is a good example! % We kick off the simulation with a hundred neutrons. n=100; x=2*rand(n,1)-1; y=2*rand(n,1)-1; z=2*rand(n,1)-1; r=(x.^2+y.^2+z.^2).^.5; % We need to only keep the particles which are inside the sphere! ikeep=find(r<1); x=x(ikeep); y=y(ikeep); z=z(ikeep); ntot=length(x); %The total number of particles inside the sphere % We select a value of dt and alpha: alpha = 15; dt = 0.0002; t = 0.0001; %We intialize time so that the graphics (title) looks better. while t < 0.6 % We duplicate particles where a random number between 0 and 1 exceeds the % threshold. Note that this requires the threshold to be much less than 1. dup = rand(size(x)); idup = find(dup