Figure 10.14:

Tracking ten particles undergoing different growth-rate dispersion mechanisms; intrinsic growth-rate dispersion (left), and growth-dependent dispersion (right).

Code for Figure 10.14

Text of the GNU GPL.

main.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
% comparing intrinsic growth rate dispersion to growth dependent dispersion
% jbr 10/5/04

nhalf = 5;
npart = 2*nhalf;

k = 10;

rng(0)
rng(0)


% particles have a distribution of growth rates
g0 = 1;
std = 0.1;
size0 = 10;
ntime = 15;
time = (0:ntime-1)';
x(:,1) = size0*ones(npart,1);
g = std*randn(npart,1) + g0;

for i = 1:ntime-1
  x(:,i+1) = x(:,i) + g;
end

% particles have same growth rate; choose which one grows at random
% particle increases size by increment delta

delta = 0.1;
nsim = 1500;
ipart = fix(rand(nsim,1)*npart)+1;
xran = size0*ones(npart,1);
simtim = (0:nsim-1)';
tau(1) = 0;

for i = 1:nsim-1
  tau(i+1) = tau(i) -log(rand)/(k*npart);
  part = ipart(i);
  xran(:,i+1) = xran(:,i);
  xran(part,i+1) = xran(part,i) + delta;
end

%sort the particles by size; split into two groups, large and small,
%and grow them again

[size, label] = sort(x(:,ntime));

y1(:,1) = x(label(1:nhalf),ntime);
g1 = g(label(1:nhalf));
y2(:,1) = x(label(nhalf+1:npart),ntime);
g2 = g(label(nhalf+1:npart));

for i = 1:ntime-1
  y1(:,i+1) = y1(:,i) + g1;
  y2(:,i+1) = y2(:,i) + g2;
end

[size, label] = sort(xran(:,nsim));
y1ran(:,1) = xran(label(1:nhalf),nsim);
y2ran(:,1) = xran(label(nhalf+1:npart),nsim);
simtim = (0:(nsim-1)/2)';
tau1(1) = 0;
tau2(1) = 0;

ipart1 = fix(rand(nsim,1)*nhalf)+1;
ipart2 = fix(rand(nsim,1)*nhalf)+1;
for i = 1:(nsim-1)/2
  tau1(i+1) = tau1(i) -log(rand)/(k*npart)*2;
  tau2(i+1) = tau2(i) -log(rand)/(k*npart)*2;
  part1 = ipart1(i);
  part2 = ipart2(i);
  y1ran(:,i+1) = y1ran(:,i);
  y2ran(:,i+1) = y2ran(:,i);
  y1ran(part1,i+1) = y1ran(part1,i) + delta;
  y2ran(part2,i+1) = y2ran(part2,i) + delta;
end

table1 = [time,x'];
table2 = [tau',xran'];

save dispmodel.dat table1 table2;

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
figure(1);
plot (table1(:,1), table1(:,2:10));
% TITLE dispmodel1

figure(2);
plot (table2(:,1), table2(:,2:10));
% TITLE dispmodel2
end % PLOTTING