%% 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