%% Monte Carlo Simulation of a Quenched Rod
% In this script we simulate the behavior of a triangular rod which has a
% dimensionless length of 2 on each side (equilateral triangle).  We shall 
% use a
% Brownian Dynamics simulation where we follow tracers as they execute a
% random walk through the domain.  
% 
% The simulation process is very simple: all we have to do is to give the
% tracers a normally distributed random kick with standard deviation
% (2dt)^.5 in each direction at every time step!
% Those who escape through the bottom or
% through the triangular sides are eliminated.  The average temperature of 
% the rod is just
% the total number of tracers at any instant divided by the number we
% started with!
%

n=10000;
xdots=2*rand(n,1)-1; %Create x,y coordinates in the range
ydots=2*rand(n,1);   % x=[-1,1] and y=[0,2]
dt=0.001; %Our discretization in time
t=[dt:dt:1.5];

x=[0:0.01:1];
height = sqrt(3).*(1-x); %We will need this to draw the triangle

% Look for those points less than the height. This will remove the points
% generated outside our geometry.

indexkeep = find(ydots<sqrt(3).*(1-abs(xdots))); %We find points inside the domain
num_points_inside = length(indexkeep); %number of valid points
n=num_points_inside;


xdots = xdots(indexkeep);
ydots = ydots(indexkeep);

% Plot this up to see the temperature distribution:
figure(1)
plot(x,height,'red',-x,height,'red')
hold on
plot(xdots,ydots,'xb')
hold off
axis([-1 1 0 sqrt(3)])
axis('equal')
title('Diffusion from a Triangular Rod')
drawnow

i=1;
temperature = zeros(size(t)); %initializing the temperature array

while i < length(t) & num_points_inside

%Propose a deltaX and deltaY random displacements
deltax = sqrt((2*dt))*randn(num_points_inside,1);
deltay = sqrt((2*dt))*randn(num_points_inside,1);

xdots = xdots+deltax;
ydots = ydots+deltay;

%Dissipation
index_inside = find(ydots<sqrt(3).*(1-abs(xdots)) & ydots>0);
xdots = xdots(index_inside);
ydots = ydots(index_inside);
num_points_inside = length(index_inside);

%Store average temperatures
temperature(i) = num_points_inside/n;

i=i+1;

figure(1) %Movie: comment this out to run faster!
plot(x,height,'red',-x,height,'red')
hold on
plot(xdots,ydots,'xb')
hold off
axis([-1 1 0 sqrt(3)])
axis('equal')
title('Diffusion from a Triangular Rod')
drawnow

end

figure(2)
plot(t, temperature);
grid on
xlabel('t')
ylabel('average temperature')
title('Average Temperature in a Triangular Rod as a Function of Time')
axis([0 .5 0 1])