s3 (Stochastic Spatial Simulator) is an X Windows software package for the interactive exploration of particle systems on one and two dimensional grids with at most eight states at each site. It offers a number of built-in models which can be simulated from a variety of initializations. In addition, s3 is designed so that the user can easily add new models by writing C code that describes the simulation loop and then adding the name of the model to a parameter structure. s3 must then be recompiled (using the accompanying make file) but (i) only the new model and one other file must be recompiled and (ii) no knowledge of the many X Windows library routines used in s3 is needed to add a model. This document explains:
There is a separate HTML document that describes the existing models and gives a brief introduction to interacting particle systems.
As our story begins, you have obtained a copy of a file called something like s3.tar.Z and you have it sitting in your home directory. The magic incantation
zcat s3.tar.Z | tar xf-
will uncompress the file and untar it all at once. If some reason this does not work you can fall back on the two step procedure
uncompress s3.tar.Z tar xf s3.tar
Untarring the file creates a directory s3 and several subdirectories. The program files come to you as ANSI C source code so the next step is to "make" the executable. If you have a Sparcstation with gcc and the X11 libraries in /usr/openwin/lib then you should be able to create the executable by
cd s3 make
If you have a different C compiler or the X11 libraries reside in a different place then you will need to edit the Makefile.
# Makefile for s3
#
SHELL=/bin/csh
# Set prefix directory for X11 libs and include files
# Uncomment one of the following. The idea is
# to have X11 libs in $(XPREFIX)/lib and
# the X11 header files in $(XPREFIX)/include.
# XPREFIX=/usr
# XPREFIX=/usr/X11
# XPREFIX=/usr/X11R5
# XPREFIX=/usr/local/X11
# XPREFIX=/usr/local/X11R5
# XPREFIX=/usr/openwin
XPREFIX=$(OPENWINHOME)
# These should now be okay.
INCLUDEDIRS= -I${XPREFIX}/include
LIBDIR=$(XPREFIX)/lib
# Keep the following if using gcc.
#
CC=gcc
DFLAGS= -g
#DLDFLAGS=-static
CMDFLAGS=-O2 ${INCLUDEDIRS}
CFLAGS=-O2 ${INCLUDEDIRS} ${DFLAGS}
LDFLAGS=-dynamic $(DLDFLAGS) -L$(LIBDIR) -lXaw -lXmu -lXt -lX11 -lXext -lm
# For SUN's ansi c compiler, comment out the previous three lines
# and uncomment the following three lines.
# CC=cc
# CFLAGS=-xO4 ${INCLUDEDIRS} ${DFLAGS}
# LDFLAGS=-Bdynamic -L$(LIBDIR) -lXaw -lXmu -lXt -lX11 -lXext -lm
OBJS = main.o graph.o parm.o states.o rand.o
# Here we go.
all: s3s $(OBJS) libs s3s
$(CC) -o s3 $(CMDFLAGS) -Iinclude $(OBJS) $(LIBS) $(LDFLAGS) $(LIBDIR)/libXmu.a
LIBS = inits/inits.a models/models.a sims/sims.a widgets/widgets.a
libs:
chmod 755 submake
./submake inits; cd inits; make "CFLAGS=${CFLAGS} -I../include" "CC=${CC}"
./submake models; cd models; make "CFLAGS=${CFLAGS} -I../include" "CC=${CC}"
./submake sims; cd sims; make "CFLAGS=${CFLAGS} -I../include" "CC=${CC}"
./submake widgets; cd widgets; make "CFLAGS=${CFLAGS} -I../include" "CC=${CC}"
s3s:
unset noclobber ;
echo '#!/bin/csh' > s3s ;
echo "" >> s3s ;
echo "setenv LD_LIBRARY_PATH ${XPREFIX}/lib" >> s3s ;
echo "" >> s3s ;
echo './s3' >> s3s ;
chmod 755 s3s
clean:
/bin/rm -f */*.a *.o */*.o *% *~{} */*% */*~{}; rm s3 s3s
The details of the Makefile are not important once you get the program to
run. However there is one feature we will use later:
make clean removes all the .a and .o files to force all the .c files to be recompiled. It also removes certain backup files that your editor might make.
Once you successfully start the program, a window will appear on the screen and be filled in with an array of colored rectangles. You are now looking at The Main Panel.
Pushing a button (here and in what follows this means positioning the cursor on the button and pushing the left mouse button) on the first row pops up an auxiliary panel.
The descriptions of the auxilary panels are given in correspondingly
labelled paragraphs below.
The second row of buttons more immediately affects the main panel.
"Quit" terminates s3.
"Reset" resets the time counter to 0 and reinitializes the cell
array.
"Start/Stop" starts the simulation (from the point where you left
off) when it is stopped and
stops the simulation when it is started. If you want to change parameters
it is not necessary to stop the simulation but it is easier to fill in the
menus if you do.
(i) Do not try to get rid of an auxillary panel by selecting Quit from the pull down menu on the menu bar, since that will Quit the entire application.
(ii) You can quit s3 by using Quit on the pull down menu on the main panel but you will get a message that an abnormal termination has occurred.
(iii) Any auxilary panel can be iconified by selecting Close from its pull down menu but in most cases you can simply use Close Panel instead.
(iv) The main panel can be iconified by selecting Close from it pull down menu and the simulation will keep running.
Size = The number of times steps that can be recorded.
TicksX = The number of tick marks to place on the horizontal axis of
the graph.
TicksY = The number of tick marks to place on the vertical axis of the
graph.
"Apply Parameters" must be pushed in order for the changes to the take effect.
"Reset" clears a graph panel of densities from past time steps.
Note: if, for instance, you change colors for a one dimensional simulation on the suggested 250 x 500 it will take several seconds to redraw the 125,000 rectangles in the array.
In addition to being changed by the states panel, the colors of a simulation can be set in the model_*.c files. See part 3 below. In fact our philosophy is that the main purpose of the states panel is to aid in finding good choices for the default colors. For this purpose it is useful to know the names of the colors you are looking at.
0 = Black 8 = Brown 16 = Tan 1 = DarkRed 9 = Red 17 = LightRed 2 = RedGray 10 = Orange 18 = LightOrange 3 = BlueGray 11 = GreenGray 19 = Yellow 4 = DarkGreen 12 = Green 20 = LightGreen 5 = DarkCyan 13 = Cyan 21 = LightCyan 6 = DarkBlue 14 = Blue 22 = LightBlue 7 = DarkMagenta 15 = Magneta 23 = LightMagenta
The shades of gray in 24--31 are simply called Gray0--Gray7. Some colors are called by two names: Pink = LightRed, White = Gray7, LightGray = Gray5, Gray = Gray2. In referring to these colors note that each word begins with a capital letter and recall C is "case sensitive" i.e., it does not think darkred and DarkRed are the same thing. If you are really fussy about how the colors look you can edit the RGB values in the function DefineColors at the end of main.c. For help in defining new colors look at the file /usr/X11/lib/rgb.txt which gives the RGB values for several hundred color names.
#include "../main.h" #include "model_GBT.h"".." means the parent of the current directory, i.e. s3, since the file is in s3/models. The header file main.h is needed to define the global variables GiNStates = the number of states in the model, cmap[8] = the colors assigned to the states, and the list which translates color names into numbers, e.g. #define White 31.
The next thing to do is to declare the model parameters, which are labelled static to make them private variables invisible to other modules.
static double Beta0 = 4.5, Beta0_default = 4.5; static double Delta0 = 1.0, Delta0_default = 1.0; static double Beta1 = 2.0, Beta1_default = 2.0; static double Delta1 = 1.0, Delta1_default = 1.0; static int Range = 4, Range_default = 4;Next we construct the menu which pops up when the model is chosen. GaParm = "global argument which is a PARM" and we suggest you retain this part of the name.
PARM GaParm_GBT [] = ${$
{"0 = grass, 1 = bushes, 2 = trees", NULL, param_unused, -1, -1, NULL, NULL, TRUE},
{"Higher state can give birth onto lower", NULL, param_unused, -1, -1, NULL, NULL, TRUE},
{"State 1 (Bushes):", NULL, param_unused, -1, -1, NULL, NULL, TRUE},
{" Linear birth rate: ", &Beta0, param_double, -1, 70, &Beta0_default, NULL, TRUE},
{" Constant death rate: ", &Delta0, param_double, -1, 70, &Delta0_default), NULL, TRUE},
{"State 2 (Trees):", NULL, param_unused, -1, -1, NULL, NULL, TRUE},
{" Linear birth rate: ", &Beta1, param_double, -1, 70, &Beta1_default[1]), NULL, TRUE},
{" Constant death rate: ", &Delta1, param_double, -1, 30,
&Delta1_default[1]), NULL, TRUE},
{" Range [0 = nearest nrb] : ", &Range, param_int, -1, 70, &Range_default[1]), NULL, TRUE},
{NULL, NULL, param_unused, -1, -1, NULL, NULL, TRUE},
}; }
To explain this, we first note that the structure PARM is defined in
parm.h as
typedef struct {
char *MpcName;
void *Value;
enum param_type paramType;
short stringSize;
short displaySize;
void *Default;
Widget MW;
int MbNewline;
} PARM; }
Here the naming conventions are M = member of a structure, p = pointer,
c = character, d = double, W = Widget, b = Boolean. So
MpcName is a string (pointer to a character) which is printed on
the model menu. Notice that extra blanks at the end of the string are
not ignored. This can be used to neaten the display.
Value is a null pointer, to allow for "arbritrary" data types.
param_type specifies the type of data that Value, the possible
types are:
indent param_double, param_int, param_string, parm_unused.
stringsize determines the number of characters allowed to define the
parameter value under consideration. If stringsize is set
equal to -1, s3 will handle this question automatically. The recommended
practice is to always set stringsize equal to 1.
displaysize determines whether or not a box will be created on the
menu, and how large this box will be. If displaysize is set equal to
-1, no box will be created. If displaysize is a positive integer, a
box of which holds that many characters will be created. A typical safe value
in this case is 30.
Default is a pointer to the default value of the parameter, which is
input by "Load Defualts" on the Model Panel. In most cases this should
agree with the initialization provided when the parameter is declared.
MW is a Widget, the basic object from which the windows and
menus are built. For the moment all you need to know that this is
a pointer to something and you should always fill in this blank
with NULL.
MbNewLine is a Boolean which is TRUE if this is the end of a line on
the model menu. By putting FALSE here you can have two
or more parameters on the same
line. In designing your menus you need to take into account the fact that
resulting widget inherits its width from the one to which it is
attached and extra text is clipped off.
Note that each line ends with a comma not a semi-colon. This punctuation
mistake or omitting one of the arguments will get the compiler quite
upset. In addition each PARM must end with the lines
{ NULL, NULL, param_unused, -1, -1, NULL, NULL, TRUE },
} ;
The first NULL is the important entry since it indicates that the end of the
menu definition has been reached. The }; on the last line is
needed to close the definition of Ga_Parm.
Next comes the function which is the simulation engine.
void FA2pcA4i_GBT (char *ApcArray, char *ApcAux,
int AiW, int AiH, int AiTime, int AiSteps)
The somewhat strange string "FA2pcA4i" means this is a Function that takes
as arguments 2 pointers to characters and 4 integers and returns nothing
(absence of letters between F and first A). In a simple model file this
precision in naming things is obnoxious but when you are dealing with
the source code for s3 which resides in a dozen .c files
such information is a useful debugging tool. You can rename
this function something like Rule_GBT if you desire but you have to remember
the name you choose and modify parts b and c below appropriately.
We will describe how to write this function in the next section. For the
moment we only want to metion that the code must contain
GiNStates = 3;
if (AiSteps==0){
cmap[0] = Black; cmap[1] = Green; cmap[2] = Tan;}
The first line sets a global variable GiNStates that contains the
number of states in the model. The last two lines set the color map
for the model but only when the function is called with AiSteps = 0.
The "if" (and some other code) is needed so that the default colors are
loaded only when you change the model and then push "Apply Parameters."
b. making the model_xyz.h file
In the case of GBT it is
#ifndef __model_GBT_h__
#define __model_GBT_h__
#include "model.h"
extern PARM GaParm_GBT [ ];
void FA2pcA4i_GBT (char *ApcSrc, char *ApcDst,
int AiW, int AiH, int AiTime, int AiSteps);
#endif
If your new model is xyz, you can simply replace each GBT by xyz. Of
course if you decided to rename the simulation engine Rule_GBT
or something else that name should go after the "void".
The workings of this file are not important for most users but if you
are curious:
The first two lines and the last one guarantee that the header file
model_GBT.h will only be included once when the program compiles.
model.h must be included to define PARM. The storage class
specifier extern tells the compiler to pass the variable name
GaParm_GBT to the linker. An extern is not
needed on the function since it is assumed to be extern unless
it is declared to be static.
c. Modifying model.c
Search for the line
# Models to include
and insert
#include "model_GBT.h"
somehere (it does not matter where) in the list of includes. Then add the
following line to the definition of GaModel
{"Grass Bushes Trees", GaParm_GBT, FA2pcA4i_GBT, CMType2D},}
The first entry is the name as you want it displayed on the
menu. The location of this line in the list determines the location of
this name in the menu.
The second entry is the name of the PARM structure.
The third entry is the name of the simulation engine. If you kept the
FA2pcA4i then you see one small reason for this notation. As you fill
in the line you see that your new function takes the same arguments as
the others. If you have changed the name of the simulation engine to
Rule_GBT, you should of course fill in your new name.
The fourth entry in the model type. The choices are
CMtype2D = two dimensions, continuous time
CMtype1D = one dimension, continuous time
CMtype2DSync = two dimensions, discrete time
CMtype1DSync = one dimension, discrete time
The meanings of these choices will be explained in the next section.
4. Writing New Models
We begin by explaining the arguments passed to the function which
is the simulation engine. As above,
A indicates that these are arguments to the function, pc = pointer
to a character, and i = integer. Two new abbreviations are Src = source,
Dst = destination.
ApcArray = the array of states
ApcAux = an auxillary arrray which may be useful someday
AiW = width of the array
AiH = height of the array
AiTime = the number of times we have updated the array. This is
needed in one dimensional simulations to write the new values in the
right place on the screen. It is ignored in continuous time.
AiSteps = the number of sites we update in continuous time
= AiW * AiH * "UpdateRate"/100
This is ignored in discrete time.
Note: The two pointers ApcArray and ApcAux currently point to
the same array. In later releases the Auxilary array will be used
in two dimensional discrete time simulations to aid in updating
the array or in nonhomogeneous models to describe the environment.
A. Two dimensional continuous time.
The simulation engine for GBT
begins by declaring variables
int rep, Diameter, i, iX, iY, iDir, iX2, iY2;
int nbrX[4] = { 1, 0, -1, 0 };
int nbrY[4] = { 0, -1, 0, 1 };
double r, max_rate, c10, c11, c20, c22;
Then continues with the two lines that all models must contain
GiNStates = 3;
if (AiSteps==0){
cmap[0] = Black; cmap[1] = Green; cmap[2] = Tan;}
For help with picking colors, see the list in the discussion of
the States Panel in Section 2.
To explain the rest of the code we have to describe the model.
The states 0 = vacant or grass, 1 = a bush, 2 = a tree.
(i) At rate Delta0 a 1 dies and becomes a 0.
(ii) At rate Beta0 a 1 gives birth to a new 1 which is sent to a
randomly chosen neighbor. If the neighbor is a 0 it becomes a 1.
(iii) At rate Delta1 a 2 dies and becomes a 0.
(iv) At rate Beta1 a 2 gives birth to a new 2 which is sent to a
randomly chosen neighbor. If the neighbor is a 1 or 0 it becomes a 2.
We will explain the phrase "randomly chosen neighbor" in a moment.
The first step is to note that 1's try to do something (give birth or
die) at rate Beta0 + Delta0, while 2's try to do something
(give birth or die) at rate Beta1 + Delta1, so we have to
compute the max_rate at which things happen:
max_rate = Beta0 + Delta0;
if (max_rate < Beta1 + Delta1)
max_rate = Beta1 + Delta1;
In writing a continuous time model we imagine each site is trying to do
something at rate max_rate so we need to compare the model rates with
max_rate to determine the probability a transition will occur:
c10 = Delta0/max_rate;
c11 = c10 + Beta0/max_rate;
c20 = Delta1/max_rate;
c22 = c20 + Beta1 /max_rate;
Again the reader will see the reasons for these definitions in a
minute. Next, the line
Diameter = 2*Range+1;
will help in determining the "randomly chosen neighbor."
Now we come to the main loop. AiSteps is the number of sites we
should update. Since we imagine each site is trying to
do something at rate max_rate they will all have the same chance of
being the next one chosen. RANDOM is an alias for a random long
integer (which is defined in rand.c and can be changed to the user's
favorite random number generator) so we call it twice:
once modulo AiW and once modulo AiH
to get a random point (iX,iY) in the grid:
for (i = 0; i < AiSteps; i++) {
iX = RANDOM % AiW;
iY = RANDOM % AiH;}
To accomodate one or two dimensional arrays our basic structure is always a
one dimensional array. In two dimensions
the element in row iX and column
iY is found in location iX + AiW * iY. The code for choosing the
random neighbor is:
if (Range==0) {
iDir = RANDOM & 0x3;
iX2 = (iX + nbrX [iDir] + AiW) % AiW;
iY2 = (iY + nbrY [iDir] + AiH) % AiH;
}
else
iDir = -Range + Diameter*UNIFORM;
iX2 = (iX + iDir +AiW) % AiW;
iDir = -Range + Diameter*UNIFORM;
iY2 = (iY + iDir +AiH) % AiH;
}
Here UNIFORM is an alias for a double uniformly distributed on
[0,1), again defined in rand.c.
If the value of Range is 0, then the neighbor is supposed to be one
of the four nearest neighbors of the first site. To pick it,
we "and" RANDOM with 3 to get a random number
iDir (integer Direction) that is 0, 1, 2, or 3 with equal probability.
nbrX[iDir] and nbrY[iDir] give the displacement of the new
point from the old one. To keep from running off the edge of the array
we use "periodic boundary conditions" that is we do our arithmetic
modulo AiW and AiH respectivley. Since many computers think
-1 mod 100 = -1 we instead add 100 first to make the number
positive and get the desired result (-1 + 100) mod 100 = 99. If the
value of Range is larger than 0, then the neighbor is supposed to be
selected at random from the sqaure of side Diameter centered at the
first site. Now iDir becomes an integer selected at random from
{-Range,...,+Range}, which is added to iX to obtain
the x-coordinate of the neighbor site. Another random iDir is added to
iY to obtain the y-coordinate, and thus the random neighbor
(iX2,iY2) is determined.
Returning to the code, what we do is based on what we find at the chosen sites:
switch (ApcArray [iX + AiW * iY]) {
case 0:
/* do nothing */
break;
case 1:
r = UNIFORM;
if (r < c10 ) /* die */
ApcArray [ iX + AiW * iY ] = 0 ;
else if (r < c11 ) /* give birth if possible */
{
if ( ApcArray [ iX2 + AiW * iY2 ] == 0 )
ApcArray [ iX2 + AiW * iY2 ] = 1 ;
}
break;
Deaths are supposed to occur to 1's at rate Delta0
while events occur at rate max_rate so we choose to have a death
event with probability c10, i.e., when r &#lt; c10.
Likewise, we choose to have a birth
event with probability Beta0/max_rate = c11 - c10, i.e., when
c10 &#lt;= r <lt; c11. Given this, we find the current value at the
neighbor (iX2, iY2), and set it = 1 if it is 0.
The case in which a 2 is chosen is similar but now we don't have to look
at neighbor but simply always set it = 2. If the neighbor is 0 or 1 it should
become 2; if it is already 2 we have done nothing.
case 2:
r = UNIFORM;
if (r < c20 ) /* die */
ApcArray [ iX + AiW * iY ] = 0 ;
else if (r < c22 ) /* give birth on anything */
{
ApcArray [ iX2 + AiW * iY2 ] = 2 ;
}
break;
Finally, three }'s are needed to close (i) the switch, (ii) the
for loop, and (iii) the function itself.
B. One dimensional continuous time.
We use the one dimensional
voter model, model_1DVM.c, as an example. In this model the
number of states, GdNStates, can take any value
between 2 and 8. The simulation engine begins by declaring variables
int i, j, k, iDir;
int y = AiTime % AiH;
int y2 = (AiTime + 1) % AiH;
int Range = GdRange;
int Nbrs = 2*Range;
Then continues with the
if (AiSteps==0){
for (i=0; i<GiNStates; i++) cmap[i] = 16+i;}
Here the colormap consists of the pastel colors 16--23. To explain the
rest of the code we have to describe the voter model. The states are
thought of as opinions. At rate 1 the voter at x decides to change her
opinion to that of a person chosen at random from x-r, ..., x-1,
x+1, ..., x+r where r is the range of the interaction.
The state of the sites 0 <= i < AiW at time AiTime are
recorded in the
array at i + y * AiW where y = AiTime % AiH, so the first
thing we do is copy the current contents of row y to row
y2 = (AiTime + 1) % AiH).
for (i = 0; i < AiW; i++) {
ApcArray [i + y2 * AiW] = ApcArray [i + y * AiW];
Now to update row y2:
for (i = 0; i < AiW; i++) {
j = RANDOM % AiW;
iDir = -Range + (RANDOM % Nbrs);
if (iDir==0) iDir = Range;
k = (j + AiW + iDir ) % AiW;
ApcArray [j + y2 * AiW] = ApcArray [k + y2 * AiW];
Inside the loop: (1) we pick a site j at random between 0 and AiW-1;
(2) makes iDir random between -Range and
Range-1 but we don't want to pick 0 and we do want to pick
Range so (3) corrects this. (4) picks the neighbor k using modulo
arithmetic to avoid running off the ends of the interval. (5) sets the
value at j to the one found at k.
C. Two dimensional discrete time.
We use the nonlinear voter model,
model_NLV.c, as an example. This time we need two new includes:
#include <stdlib.h>
#include <string.h>
The first is needed to allocate memory with malloc, the second is
needed for memcpy. The simulation engine
begins as usual by declaring variables
int i, x, y, x2, y2, n;
double r;
char *tmp;
int nbrX[4] = { 1, 0, -1, 0 };
int nbrY[4] = { 0, -1, 0, 1 };
double prob[6];
and continues with the two lines that all models must contain
GiNStates = 2
if (AiSteps==0){
cmap[0] = Black; cmap[1] = LightGreen;
In the nonlinear voter model, to determine the state of x at time t+1
we examine the state of x and its four neighbors at time t. If we find
k neighbors that are 1 then x will be 1 with probability prob[k].
To make the model symmetric under interchange of 0's and 1's we
require prob[5-i] = 1-prob[i].
prob[0]=0.0; prob[1]=p1; prob[2]=p2;
for (i=0; i<3; i++) prob[5-i] = 1-prob[i];
To calculate the state at time t+1 without disturbing the state at time t
we need to allocate a new array
tmp = (char *) malloc(AiW*AiH);
Then we go through the array one site at a time, count the number of 1's
we find and decide the new state:
for (x = 0; x &#lt; AiW; x++) {
for (y = 0; y &#lt; AiH; y++)
{
n= ApcArray [x + AiW*y];
for (i = 0; i &#lt; 4; i++)
{
x2 = (x + nbrX[i] + AiW) % AiW;
y2 = (y + nbrY[i] + AiH) % AiH;
n += ApcArray[x2 + AiW*y2];
}
/* default is set state to 0 */
tmp[x + AiW*y] = 0;
/* now check to see if should really have a 1 */
r = UNIFORM;
if (r &#lt; prob[n] ) tmp[x + AiW*y] = 1;
}
}
Finally we copy the new array in dst into ApcDst and free the
memory allocated. We have to free the memory every time to prepare for
the moment at which the simulation is terminated.
memcpy(ApcArray,tmp,AiW*AiH*sizeof(char));
free(tmp);
D. One dimensional discrete time.
Our example is model_1DDTCP.c,
the discrete time contact process. The simulation engine
begins as usual by declaring variables
int i, j, k;
int y = AiTime % AiH;
int y2 = (AiTime + 1) % AiH;
and continues with the two lines that all models must contain
GiNStates = 2
if (AiSteps==0){
cmap[0] = Black; cmap[1] = LightCyan;}
In the discrete time contact process 1 indicates an occupied site and
0 a vacant site. Occupied sites survive from time t to time t+1 with
probability GdDeath. If the 1 at x survives it gives birth
onto x+1 and onto x-1 with probability GdBirth each. We think of
occupied sites as the locations of plants so there can be at most one
individual per site.
/* determine survivors */
for (i = 0; i < AiW; i++) {
if ((ApcArray[i+y*AiW]==1) && (UNIFORM > GdDeath))
ApcArray[i+y2*AiW] = 1;
else ApcArray[i+y2*AiW]=0;
}
/* Now we determine births */
for (i = 0; i < AiW; i++)
{
j = (i -1 + AiW) % AiW;
if ((ApcArray[j+y2*AiW]==1) && (UNIFORM < GdBirth))
ApcArray[i+y2*AiW]=1;
j = (i +1 + AiW) % AiW;
if ((ApcArray[j+y2*AiW]==1) && (UNIFORM < GdBirth))
ApcArray[i+y2*AiW]=1;
}