Next: Listing for disreal3d.c Up: Computer programs for Previous: Code for random

## Listing for hydreal3d.c

```#define AGRATE 0.25        /* Probability of gypsum absorption by CSH */

/* routine to select a new neighboring location to (xloc, yloc, zloc) */
/* for a diffusing species */
/* Returns a prime number flag indicating direction chosen */
/* Calls ran1 */
/* Called by movecsh, extettr, extfh3, movegyp, extafm, moveettr, */
/* extpozz, movefh3, movech, extc3ah6, movec3a */
int moveone(xloc,yloc,zloc,act,sumold)
int *xloc,*yloc,*zloc,*act,sumold;
{
int plok,sumnew,xl1,yl1,zl1,act1;

sumnew=1;
/* store the input values for location */
xl1=(*xloc);
yl1=(*yloc);
zl1=(*zloc);
act1=(*act);

/* Choose one of six directions (at random) for the new */
/* location */
plok=6.*ran1(seed);
if((plok>5)||(plok<0)){plok=5;}

switch (plok){
case 0:
xl1-=1;
act1=1;
if(xl1<0){xl1=(SYSIZE-1);}
if(sumold%2!=0){sumnew=2;}
break;
case 1:
xl1+=1;
act1=2;
if(xl1>=SYSIZE){xl1=0;}
if(sumold%3!=0){sumnew=3;}
break;
case 2:
yl1-=1;
act1=3;
if(yl1<0){yl1=(SYSIZE-1);}
if(sumold%5!=0){sumnew=5;}
break;
case 3:
yl1+=1;
act1=4;
if(yl1>=SYSIZE){yl1=0;}
if(sumold%7!=0){sumnew=7;}
break;
case 4:
zl1-=1;
act1=5;
if(zl1<0){zl1=(SYSIZE-1);}
if(sumold%11!=0){sumnew=11;}
break;
case 5:
zl1+=1;
act1=6;
if(zl1>=SYSIZE){zl1=0;}
if(sumold%13!=0){sumnew=13;}
break;
default:
break;
}

/* Return the new location */
*xloc=xl1;
*yloc=yl1;
*zloc=zl1;
*act=act1;
/* sumnew returns a prime number indicating that a specific direction */
/* has been chosen */
return(sumnew);
}

/* routine to move a diffusing CSH species */
/* Inputs: current location (xcur,ycur,zcur) and flag indicating if final */
/* step in diffusion process */
/* Returns flag indicating action taken (reaction or diffusion/no movement) */
/* Calls moveone */
/* Called by hydrate */
int movecsh(xcur,ycur,zcur,finalstep)
int xcur,ycur,zcur,finalstep;
{
int xnew,ynew,znew,plok,action,sumback,sumin,check;

action=0;
/* Store current location of species */
xnew=xcur;
ynew=ycur;
znew=zcur;
sumin=1;
sumback=moveone(&xnew,&ynew,&znew,&action,sumin);

if(action==0){printf("Error in value of action \n");}
check=mic[xnew][ynew][znew];

/* if new location is solid C3S, C2S, or CSH, then convert */
/* diffusing CSH species to solid CSH */
if((check==C3S)||(check==C2S)||(check==CSH)||(finalstep==1)){
mic[xcur][ycur][zcur]=CSH;
/* decrement count of diffusing CSH species */
/* and increment count of solid CSH */
count[DIFFCSH]-=1;
count[CSH]+=1;
action=0;
}

if(action!=0){
/* if diffusion step is possible, perform it */
if(check==POROSITY){
mic[xcur][ycur][zcur]=POROSITY;
mic[xnew][ynew][znew]=DIFFCSH;
}
else{
/* indicate that diffusing CSH species remained */
/* at original location */
action=7;
}
}
return(action);
}

/* routine to return count of number of neighboring pixels for pixel */
/* (xck,yck,zck) which are not phase ph1, ph2, or ph3 which are input as */
/* parameters */
/* Calls no other routines */
/* Called by extettr, extfh3, extch, extafm, extpozz, extc3ah6 */
int edgecnt(xck,yck,zck,ph1,ph2,ph3)
int xck,yck,zck,ph1,ph2,ph3;
{
int ixe,iye,ize,edgeback,x2,y2,z2,check;

/* counter for number of neighboring pixels which are not ph1, ph2, or ph3 */
edgeback=0;

/* Examine all pixels in a 3*3*3 box centered at (xck,yck,zck) */
/* except for the central pixel */
for(ixe=(-1);ixe<=1;ixe++){
x2=xck+ixe;
for(iye=(-1);iye<=1;iye++){
y2=yck+iye;
for(ize=(-1);ize<=1;ize++){

if((ixe!=0)||(iye!=0)||(ize!=0)){
z2=zck+ize;
/* adjust to maintain periodic boundaries */
if(x2<0){x2=(SYSIZE-1);}
else if(x2>=SYSIZE){x2=0;}
if(y2<0){y2=(SYSIZE-1);}
else if(y2>=SYSIZE){y2=0;}
if(z2<0){z2=(SYSIZE-1);}
else if(z2>=SYSIZE){z2=0;}
check=mic[x2][y2][z2];
if((check!=ph1)&&(check!=ph2)&&(check!=ph3)){
edgeback+=1;
}
}
}
}
}
/* return number of neighboring pixels which are not ph1, ph2, or ph3 */
return(edgeback);
}
```

```/* routine to add extra ettringite when gypsum reacts with */
/* in a fashion to preserve needle growth */
/* Returns flag indicating action taken */
/* Calls moveone and edgecnt */
/* Called by movegyp and movec3a */
int extettr(xpres,ypres,zpres)
int xpres,ypres,zpres;
{
int check,newact,multf,numnear,sump,xchr,ychr,zchr,fchr,i1,plok;
float pneigh,ptest;
long int tries;

/* first try 6 neighboring locations until      */
/*	a) successful                           */
/*	b) all 6 sites are tried and full or    */
/*	c) 500 tries are made                   */
fchr=0;
sump=1;
/* Note that 30030 = 2*3*5*7*11*13 */
/* indicating that all six sites have been tried */
for(i1=1;((i1<=500)&&(fchr==0)&&(sump!=30030));i1++){

/* determine location of neighbor (using periodic boundaries) */
xchr=xpres;
ychr=ypres;
zchr=zpres;
newact=0;
multf=moveone(&xchr,&ychr,&zchr,&newact,sump);
if(newact==0){printf("Error in value of action \n");}

check=mic[xchr][ychr][zchr];

/* if neighbor is porosity, and conditions are favorable */
/* based on number of neighboring ettringite, C3A, or C4AF */
/* pixels then locate the ettringite there */
if(check==POROSITY){
numnear=edgecnt(xchr,ychr,zchr,ETTR,C3A,C4AF);
pneigh=(float)(numnear+1)/26.0;
pneigh*=pneigh;
ptest=ran1(seed);
if(pneigh>=ptest){
mic[xchr][ychr][zchr]=ETTR;
fchr=1;
}
}
else{
sump*=multf;
}
}

/* if no neighbor available, locate ettringite at random location */
/* in pore space in contact with at least another ettringite */
/* or aluminate surface  */
tries=0;
while(fchr==0){
tries+=1;
newact=7;
/* generate a random location in the 3-D system */
xchr=(int)((float)SYSIZE*ran1(seed));
ychr=(int)((float)SYSIZE*ran1(seed));
zchr=(int)((float)SYSIZE*ran1(seed));
if(xchr>=SYSIZE){xchr=0;}
if(ychr>=SYSIZE){ychr=0;}
if(zchr>=SYSIZE){zchr=0;}

check=mic[xchr][ychr][zchr];
/* if location is porosity, locate the ettringite there */
if(check==POROSITY){
numnear=edgecnt(xchr,ychr,zchr,ETTR,C3A,C4AF);
/* be sure that at least one neighboring pixel */
/* is ettringite, or aluminate clinker */
if((tries>5000)||(numnear<26)){
mic[xchr][ychr][zchr]=ETTR;
fchr=1;
}
}
}
return(newact);
}

/* routine to add extra FH3 when gypsum reacts with */
/* C4AF at location (xpres,ypres,zpres) */
/* Called by movegyp and moveettr */
/* Calls moveone and edgecnt */
void extfh3(xpres,ypres,zpres)
int xpres,ypres,zpres;
{
int multf,numnear,sump,xchr,ychr,zchr,check,fchr,i1,plok,newact;
long int tries;

/* first try 6 neighboring locations until      */
/*	a) successful                           */
/*	b) all 6 sites are tried and full or    */
/*	c) 500 tries are made                   */
fchr=0;
sump=1;
for(i1=1;((i1<=500)&&(fchr==0)&&(sump!=30030));i1++){

/* choose a neighbor at random */
xchr=xpres;
ychr=ypres;
zchr=zpres;
newact=0;
multf=moveone(&xchr,&ychr,&zchr,&newact,sump);
if(newact==0){printf("Error in value of newact in extfh3 \n");}
check=mic[xchr][ychr][zchr];

/* if neighbor is porosity   */
/* then locate the FH3 there */
if(check==POROSITY){
mic[xchr][ychr][zchr]=FH3;
fchr=1;
}
else{
sump*=multf;
}
}

/* if no neighbor available, locate FH3 at random location */
/* in pore space in contact with at least one FH3 */
tries=0;
while(fchr==0){
tries+=1;
/* generate a random location in the 3-D system */
xchr=(int)((float)SYSIZE*ran1(seed));
ychr=(int)((float)SYSIZE*ran1(seed));
zchr=(int)((float)SYSIZE*ran1(seed));
if(xchr>=SYSIZE){xchr=0;}
if(ychr>=SYSIZE){ychr=0;}
if(zchr>=SYSIZE){zchr=0;}
check=mic[xchr][ychr][zchr];
/* if location is porosity, locate the FH3 there */
if(check==POROSITY){
numnear=edgecnt(xchr,ychr,zchr,FH3,FH3,DIFFFH3);
/* be sure that at least one neighboring pixel */
/* is FH3 or diffusing FH3 */
if((numnear<26)||(tries>5000)){
mic[xchr][ychr][zchr]=FH3;
fchr=1;
}
}
}
}

/* routine to add extra CH when gypsum reacts with */
/* C4AF */
/* Called by movegyp and moveettr */
/* Calls edgecnt */
void extch()
{
int numnear,sump,xchr,ychr,zchr,fchr,i1,plok,check;
long int tries;

fchr=0;
tries=0;
/* locate CH at random location */
/* in pore space in contact with at least another CH */
while(fchr==0){
tries+=1;
/* generate a random location in the 3-D system */
xchr=(int)((float)SYSIZE*ran1(seed));
ychr=(int)((float)SYSIZE*ran1(seed));
zchr=(int)((float)SYSIZE*ran1(seed));
if(xchr>=SYSIZE){xchr=0;}
if(ychr>=SYSIZE){ychr=0;}
if(zchr>=SYSIZE){zchr=0;}
check=mic[xchr][ychr][zchr];

/* if location is porosity, locate the CH there */
if(check==POROSITY){
numnear=edgecnt(xchr,ychr,zchr,CH,DIFFCH,CH);
/* be sure that at least one neighboring pixel */
/* is CH or diffusing CH */
if((numnear<26)||(tries>5000)){
mic[xchr][ychr][zchr]=CH;
fchr=1;
}
}
}
}

/* routine to move a diffusing gypsum species */
/* from current location (xcur,ycur,zcur) */
/* Returns action flag indicating response taken */
/* Called by hydrate */
/* Calls moveone, extettr, extch, and extfh3 */
int movegyp(xcur,ycur,zcur,finalstep)
int xcur,ycur,zcur,finalstep;
{
int check,xnew,ynew,znew,plok,action,nexp,iexp;
int xexp,yexp,zexp,newact,sumold,sumgarb;
float pexp,pext;

sumold=1;

/* First be sure that a diffusing gypsum species is located at xcur,ycur,zcur */
if(mic[xcur][ycur][zcur]!=DIFFGYP){
action=0;
return(action);
}

/* Determine new coordinates (periodic boundaries are used) */
xnew=xcur;
ynew=ycur;
znew=zcur;
action=0;
sumgarb=moveone(&xnew,&ynew,&znew,&action,sumold);
if(action==0){printf("Error in value of action in movegyp \n");}
check=mic[xnew][ynew][znew];

/* if new location is CSH, check for absorption of gypsum */
if((check==CSH)&&((float)count[ABSGYP]<(gypabsprob*(float)count[CSH]))){
pexp=ran1(seed);
if(pexp<AGRATE){
/* update counts for absorbed and diffusing gypsum */
count[ABSGYP]+=1;
count[DIFFGYP]-=1;
mic[xcur][ycur][zcur]=ABSGYP;
action=0;
}
}

/* if new location is C3A or diffusing C3A, execute conversion */
/* to ettringite (including necessary volumetric expansion) */
if((check==C3A)||(check==DIFFC3A)){
/* Convert diffusing gypsum to an ettringite pixel */
action=0;
mic[xcur][ycur][zcur]=ETTR;
count[DIFFGYP]-=1;
if(check==DIFFC3A){
/* decrement count of diffusing C3A species */
count[DIFFC3A]-=1;
}

/* determine if C3A should be converted to ettringite */
/* 1 unit of gypsum requires 0.40 units of C3A */
/* and should form 3.30 units of ettringite */
pexp=ran1(seed);
nexp=2;
if(pexp<=0.40){
mic[xnew][ynew][znew]=ETTR;
nexp=1;
}
else{
/* maybe someday, use a new FIXEDC3A here */
/* so it won't dissolve later */
if(check==C3A){
mic[xnew][ynew][znew]=C3A;
}
else{
count[DIFFC3A]+=1;
mic[xnew][ynew][znew]=DIFFC3A;
}
nexp=2;
}

/* create extra ettringite pixels to maintain volume stoichiometry */
/* xexp, yexp, and zexp hold coordinates of most recently added ettringite */
/* species as we attempt to grow a needle like structure */
xexp=xcur;
yexp=ycur;
zexp=zcur;
for(iexp=1;iexp<=nexp;iexp++){
newact=extettr(xexp,yexp,zexp);
/* update xexp, yexp and zexp as needed */
switch (newact){
case 1:
xexp-=1;
if(xexp<0){xexp=(SYSIZE-1);}
break;
case 2:
xexp+=1;
if(xexp>=SYSIZE){xexp=0;}
break;
case 3:
yexp-=1;
if(yexp<0){yexp=(SYSIZE-1);}
break;
case 4:
yexp+=1;
if(yexp>=SYSIZE){yexp=0;}
break;
case 5:
zexp-=1;
if(zexp<0){zexp=(SYSIZE-1);}
break;
case 6:
zexp+=1;
if(zexp>=SYSIZE){zexp=0;}
break;
default:
break;
}
}

/* probabilistic-based expansion for last ettringite pixel */
pexp=ran1(seed);
if(pexp<=0.30){
newact=extettr(xexp,yexp,zexp);
}
}

/* if new location is C4AF execute conversion */
/* to ettringite (including necessary volumetric expansion) */
if(check==C4AF){
mic[xcur][ycur][zcur]=ETTR;
count[DIFFGYP]-=1;

/* determine if C4AF should be converted to ettringite */
/* 1 unit of gypsum requires 0.575 units of C4AF */
/* and should form 3.30 units of ettringite */
pexp=ran1(seed);
nexp=2;
if(pexp<=0.575){
mic[xnew][ynew][znew]=ETTR;
nexp=1;
pext=ran1(seed);
/* Addition of extra CH */
if(pext<0.2584){
extch();
}
pext=ran1(seed);
/* Addition of extra FH3 */
if(pext<0.5453){
extfh3(xnew,ynew,znew);
}

}
else{
/* maybe someday, use a new FIXEDC4AF here */
/* so it won't dissolve later */
mic[xnew][ynew][znew]=C4AF;
nexp=2;
}

/* create extra ettringite pixels to maintain volume stoichiometry */
/* xexp, yexp and zexp hold coordinates of most recently added ettringite */
/* species as we attempt to grow a needle like structure */
xexp=xcur;
yexp=ycur;
zexp=zcur;
for(iexp=1;iexp<=nexp;iexp++){
newact=extettr(xexp,yexp,zexp);
/* update xexp, yexp and zexp as needed */
switch (newact){
case 1:
xexp-=1;
if(xexp<0){xexp=(SYSIZE-1);}
break;
case 2:
xexp+=1;
if(xexp>=SYSIZE){xexp=0;}
break;
case 3:
yexp-=1;
if(yexp<0){yexp=(SYSIZE-1);}
break;
case 4:
yexp+=1;
if(yexp>=SYSIZE){yexp=0;}
break;
case 5:
zexp-=1;
if(zexp<0){zexp=(SYSIZE-1);}
break;
case 6:
zexp+=1;
if(zexp>=SYSIZE){zexp=0;}
break;
default:
break;
}
}

/* probabilistic-based expansion for last ettringite pixel */
pexp=ran1(seed);
if(pexp<=0.30){
newact=extettr(xexp,yexp,zexp);
}
action=0;
}

/* if last diffusion step and no reaction, convert back to */
/* solid gypsum */
if((action!=0)&&(finalstep==1)){
action=0;
count[DIFFGYP]-=1;
mic[xcur][ycur][zcur]=GYPSUM;
}

if(action!=0){
/* if diffusion is possible, execute it */
if(check==POROSITY){
mic[xcur][ycur][zcur]=POROSITY;
mic[xnew][ynew][znew]=DIFFGYP;
}
else{
/* indicate that diffusing gypsum remained at */
/* original location */
action=7;
}
}
return(action);
}

/* routine to add extra AFm phase when diffusing ettringite reacts */
/* with C3A (diffusing or solid) at location (xpres,ypres,zpres) */
/* Called by moveettr and movec3a */
/* Calls moveone and edgecnt */
void extafm(xpres,ypres,zpres)
int xpres,ypres,zpres;
{
int check,sump,xchr,ychr,zchr,fchr,i1,plok,newact,numnear;
long int tries;

/* first try 6 neighboring locations until      */
/*	a) successful                           */
/*	b) all 6 sites are tried or             */
/*	c) 100 tries are made                   */
fchr=0;
sump=1;
for(i1=1;((i1<=100)&&(fchr==0)&&(sump!=30030));i1++){

/* determine location of neighbor (using periodic boundaries) */
xchr=xpres;
ychr=ypres;
zchr=zpres;
newact=0;
sump*=moveone(&xchr,&ychr,&zchr,&newact,sump);
if(newact==0){printf("Error in value of newact in extafm \n");}
check=mic[xchr][ychr][zchr];

/* if neighbor is porosity, locate the AFm phase there */
if(check==POROSITY){
mic[xchr][ychr][zchr]=AFM;
fchr=1;
}
}

/* if no neighbor available, locate AFm phase at random location */
/* in pore space */
tries=0;
while(fchr==0){
tries+=1;
/* generate a random location in the 3-D system */
xchr=(int)((float)SYSIZE*ran1(seed));
ychr=(int)((float)SYSIZE*ran1(seed));
zchr=(int)((float)SYSIZE*ran1(seed));
if(xchr>=SYSIZE){xchr=0;}
if(ychr>=SYSIZE){ychr=0;}
if(zchr>=SYSIZE){zchr=0;}
check=mic[xchr][ychr][zchr];

/* if location is porosity, locate the extra AFm there */
if(check==POROSITY){
numnear=edgecnt(xchr,ychr,zchr,AFM,C3A,C4AF);
/* Be sure that at least one neighboring pixel is */
/* Afm phase, C3A, or C4AF */
if((tries>5000)||(numnear<26)){
mic[xchr][ychr][zchr]=AFM;
fchr=1;
}
}
}
}

/* routine to move a diffusing ettringite species */
/* currently located at (xcur,ycur,zcur) */
/* Called by hydrate */
/* Calls moveone, extch, extfh3, and extafm */
int moveettr(xcur,ycur,zcur,finalstep)
int xcur,ycur,zcur,finalstep;
{
int check,xnew,ynew,znew,plok,action,nexp,iexp;
int xexp,yexp,zexp,newact,sumold,sumgarb;
float pexp,pafm,pgrow;

/* First be sure a diffusing ettringite species is located at xcur,ycur,zcur */
if(mic[xcur][ycur][zcur]!=DIFFETTR){
action=0;
return(action);
}

/* Determine new coordinates (periodic boundaries are used) */
xnew=xcur;
ynew=ycur;
znew=zcur;
action=0;
sumold=1;
sumgarb=moveone(&xnew,&ynew,&znew,&action,sumold);
if(action==0){printf("Error in value of action in moveettr \n");}
check=mic[xnew][ynew][znew];

/* if new location is C4AF, execute conversion */
/* to AFM phase (including necessary volumetric expansion) */
if(check==C4AF){
/* Convert diffusing ettringite to AFM phase */
mic[xcur][ycur][zcur]=AFM;
count[DIFFETTR]-=1;

/* determine if C4AF should be converted to Afm */
/* or FH3- 1 unit of ettringite requires 0.348 units */
/* of C4AF to form 1.278 units of Afm, */
/* 0.0901 units of CH and 0.1899 units of FH3 */
pexp=ran1(seed);

if(pexp<=0.278){
mic[xnew][ynew][znew]=AFM;
pafm=ran1(seed);
/* 0.3241= 0.0901/0.278 */
if(pafm<=0.3241){
extch();
}
pafm=ran1(seed);
/* 0.4313= ((.1899-(.348-.278))/.278)   */
if(pafm<=0.4313){
extfh3(xnew,ynew,znew);
}
}
else if (pexp<=0.348){
mic[xnew][ynew][znew]=FH3;
}
action=0;
}

/* if new location is C3A or diffusing C3A, execute conversion */
/* to AFM phase (including necessary volumetric expansion) */
if((check==C3A)||(check==DIFFC3A)){
/* Convert diffusing ettringite to AFM phase */
action=0;
mic[xcur][ycur][zcur]=AFM;
count[DIFFETTR]-=1;

if(check==DIFFC3A){
/* decrement count of diffusing C3A species */
count[DIFFC3A]-=1;
}

/* determine if C3A should be converted to AFm */
/* 1 unit of ettringite requires 0.2424 units of C3A */
/* and should form 1.278 units of AFm phase */
pexp=ran1(seed);
if(pexp<=0.2424){
mic[xnew][ynew][znew]=AFM;
pafm=(-0.1);
}
else{
/* maybe someday, use a new FIXEDC3A here */
/* so it won't dissolve later */
if(check==C3A){
mic[xnew][ynew][znew]=C3A;
}
else{
count[DIFFC3A]+=1;
mic[xnew][ynew][znew]=DIFFC3A;
}
pafm=(0.278-0.2424)/(1.0-0.2424);
}

/* probabilistic-based expansion for new AFm phase pixel */
pexp=ran1(seed);
if(pexp<=pafm){
extafm(xcur,ycur,zcur);
}
}

/* Check for conversion back to solid ettringite */
if(check==ETTR){
pgrow=ran1(seed);
if(pgrow<=ETTRGROW){
mic[xcur] [ycur] [zcur]=ETTR;
action=0;
count[DIFFETTR]-=1;
}
}

/* if last diffusion step and no reaction, convert back to */
/* solid ettringite */
if((action!=0)&&(finalstep==1)){
action=0;
count[DIFFETTR]-=1;
mic[xcur][ycur][zcur]=ETTR;
}

if(action!=0){
/* if diffusion is possible, execute it */
if(check==POROSITY){
mic[xcur][ycur][zcur]=POROSITY;
mic[xnew][ynew][znew]=DIFFETTR;
}
else{
/* indicate that diffusing ettringite remained at */
/* original location */
action=7;
}
}
return(action);
}

/* routine to add extra pozzolanic CSH when CH reacts at */
/* pozzolanic surface (e.g. silica fume) located at (xpres,ypres,zpres) */
/* Called by movech */
/* Calls moveone and edgecnt */
void extpozz(xpres,ypres,zpres)
int xpres,ypres,zpres;
{
int check,sump,xchr,ychr,zchr,fchr,i1,plok,action,numnear;
long int tries;

/* first try 6 neighboring locations until      */
/*	a) successful                           */
/*	b) all 6 sites are tried or             */
/*	c) 100 tries are made                   */
fchr=0;
sump=1;
for(i1=1;((i1<=100)&&(fchr==0)&&(sump!=30030));i1++){

/* determine location of neighbor (using periodic boundaries) */
xchr=xpres;
ychr=ypres;
zchr=zpres;
action=0;
sump*=moveone(&xchr,&ychr,&zchr,&action,sump);
if(action==0){printf("Error in value of action in extpozz \n");}
check=mic[xchr][ychr][zchr];

/* if neighbor is porosity, locate the pozzolanic CSH there */
if(check==POROSITY){
mic[xchr][ychr][zchr]=POZZ;
fchr=1;
}
}

/* if no neighbor available, locate pozzolanic CSH at random location */
/* in pore space */
tries=0;
while(fchr==0){
tries+=1;
/* generate a random location in the 3-D system */
xchr=(int)((float)SYSIZE*ran1(seed));
ychr=(int)((float)SYSIZE*ran1(seed));
zchr=(int)((float)SYSIZE*ran1(seed));
if(xchr>=SYSIZE){xchr=0;}
if(ychr>=SYSIZE){ychr=0;}
if(zchr>=SYSIZE){zchr=0;}
check=mic[xchr][ychr][zchr];
/* if location is porosity, locate the extra pozzolanic CSH there */
if(check==POROSITY){
numnear=edgecnt(xchr,ychr,zchr,POZZ,CSH,POZZ);
/* Be sure that one neighboring species is CSH or */
/* pozzolanic material */
if((tries>5000)||(numnear<26)){
mic[xchr][ychr][zchr]=POZZ;
fchr=1;
}
}
}
}
```

```/* routine to move a diffusing FH3 species */
/* from location (xcur,ycur,zcur) with nucleation probability nucprob */
/* Called by hydrate */
/* Calls moveone */
int movefh3(xcur,ycur,zcur,finalstep,nucprob)
int xcur,ycur,zcur,finalstep;
float nucprob;
{
int check,xnew,ynew,znew,plok,action,sumold,sumgarb;
float pgen;

/* first check for nucleation */
pgen=ran1(seed);

if((nucprob>=pgen)||(finalstep==1)){
action=0;
mic[xcur][ycur][zcur]=FH3;
count[DIFFFH3]-=1;
}
else{

/* determine new location (using periodic boundaries) */
xnew=xcur;
ynew=ycur;
znew=zcur;
action=0;
sumold=1;
sumgarb=moveone(&xnew,&ynew,&znew,&action,sumold);
if(action==0){printf("Error in value of action in movefh3 \n");}
check=mic[xnew][ynew][znew];

/* check for growth of FH3 crystal */
if(check==FH3){
mic[xcur][ycur][zcur]=FH3;
count[DIFFFH3]-=1;
action=0;
}

if(action!=0){
/* if diffusion is possible, execute it */
if(check==POROSITY){
mic[xcur][ycur][zcur]=POROSITY;
mic[xnew][ynew][znew]=DIFFFH3;
}
else{
/* indicate that diffusing FH3 species */
/* remained at original location */
action=7;
}
}
}
return(action);
}

/* routine to move a diffusing CH species */
/* from location (xcur,ycur,zcur) with nucleation probability nucprob */
/* Called by hydrate */
/* Calls moveone and extpozz */
int movech(xcur,ycur,zcur,finalstep,nucprob)
int xcur,ycur,zcur,finalstep;
float nucprob;
{
int check,xnew,ynew,znew,plok,action,sumgarb,sumold;
float pexp,pgen;

/* first check for nucleation */
pgen=ran1(seed);
if((nucprob>=pgen)||(finalstep==1)){
action=0;
mic[xcur][ycur][zcur]=CH;
count[DIFFCH]-=1;
}
else{

/* determine new location (using periodic boundaries) */
xnew=xcur;
ynew=ycur;
znew=zcur;
action=0;
sumold=1;
sumgarb=moveone(&xnew,&ynew,&znew,&action,sumold);
if(action==0){printf("Error in value of action in movech \n");}
check=mic[xnew][ynew][znew];

/* check for growth of CH crystal */
if(check==CH){
mic[xcur][ycur][zcur]=CH;
count[DIFFCH]-=1;
action=0;
}

/* check for pozzolanic reaction */
/* 36.41 units of CH can react with 27 units of S */
if((check==POZZ)&&(npr<=(int)((float)nfill*1.35))){
mic[xcur][ycur][zcur]=POZZ;
/* update counter of number of diffusing CH */
/* which have reacted pozzolanically */
npr+=1;
count[DIFFCH]-=1;
action=0;
/* allow for extra pozzolanic CSH as needed */
pexp=ran1(seed);
/* 0.4795=(80.87-36.41-27)/36.41 */
if(pexp<=0.4795){
extpozz(xcur,ycur,zcur);
}
}
if(action!=0){
/* if diffusion is possible, execute it */
if(check==POROSITY){
mic[xcur][ycur][zcur]=POROSITY;
mic[xnew][ynew][znew]=DIFFCH;
}
else{
/* indicate that diffusing CH species */
/* remained at original location */
action=7;
}
}
}
return(action);
}

/* routine to add extra C3AH6 when diffusing C3A nucleates or reacts at */
/* C3AH6 surface at location (xpres,ypres,zpres) */
/* Called by movec3a */
/* Calls moveone and edgecnt */
void extc3ah6(xpres,ypres,zpres)
int xpres,ypres,zpres;
{
int check,sump,xchr,ychr,zchr,fchr,i1,plok,action,numnear;
long int tries;

/* First try 6 neighboring locations until      */
/* 	a) successful                           */
/*	b) all 6 sites are tried or             */
/*	c) 100 random attempts are made         */
fchr=0;
sump=1;
for(i1=1;((i1<=100)&&(fchr==0)&&(sump!=30030));i1++){

/* determine new coordinates (using periodic boundaries) */
xchr=xpres;
ychr=ypres;
zchr=zpres;
action=0;
sump*=moveone(&xchr,&ychr,&zchr,&action,sump);
if(action==0){printf("Error in action value in extc3ah6 \n");}
check=mic[xchr][ychr][zchr];

/* if neighbor is pore space, convert it to C3AH6 */
if(check==POROSITY){
mic[xchr][ychr][zchr]=C3AH6;
fchr=1;
}
}

/* if unsuccessful, add C3AH6 at random location in pore space */
tries=0;
while(fchr==0){
tries+=1;
xchr=(int)((float)SYSIZE*ran1(seed));
ychr=(int)((float)SYSIZE*ran1(seed));
zchr=(int)((float)SYSIZE*ran1(seed));
if(xchr>=SYSIZE){xchr=0;}
if(ychr>=SYSIZE){ychr=0;}
if(zchr>=SYSIZE){zchr=0;}
check=mic[xchr][ychr][zchr];

if(check==POROSITY){
numnear=edgecnt(xchr,ychr,zchr,C3AH6,C3A,C3AH6);
/* Be sure that new C3AH6 is in contact with */
/* at least one C3AH6 or C3A */
if((tries>5000)||(numnear<26)){
mic[xchr][ychr][zchr]=C3AH6;
fchr=1;
}
}
}
}

/* routine to move a diffusing C3A species */
/* from location (xcur,ycur,zcur) with nucleation probability of nucprob */
/* Called by hydrate */
/* Calls extc3ah6, moveone, extettr, and extafm */
int movec3a(xcur,ycur,zcur,finalstep,nucprob)
int xcur,ycur,zcur,finalstep;
float nucprob;
{
int check,xnew,ynew,znew,plok,action,sumgarb,sumold;
int xexp,yexp,zexp,nexp,iexp,newact;
float pgen,pexp,pafm,pgrow;

/* First be sure that a diffusing C3A species is at (xcur,ycur,zcur) */
if(mic[xcur][ycur][zcur]!=DIFFC3A){
action=0;
return(action);
}

/* Check for nucleation into solid C3AH6 */
pgen=ran1(seed);

if((nucprob>=pgen)||(finalstep==1)){
action=0;
mic[xcur][ycur][zcur]=C3AH6;
/* decrement count of diffusing C3A species */
count[DIFFC3A]-=1;
/* allow for probabilistic-based expansion of C3AH6 */
/* crystal to account for volume stoichiometry */
pexp=ran1(seed);
if(pexp<=0.69){
extc3ah6(xcur,ycur,zcur);
}
}
else{
/* determine new coordinates (using periodic boundaries) */
xnew=xcur;
ynew=ycur;
znew=zcur;
action=0;
sumold=1;
sumgarb=moveone(&xnew,&ynew,&znew,&action,sumold);
if(action==0){printf("Error in value of action in movec3a \n");}
check=mic[xnew][ynew][znew];

/* check for growth of C3AH6 crystal */
if(check==C3AH6){
pgrow=ran1(seed);
/* Try to slow down growth of C3AH6 crystals to */
/* promote ettringite and Afm formation */
if(pgrow<=C3AH6GROW){
mic[xcur][ycur][zcur]=C3AH6;
count[DIFFC3A]-=1;
action=0;
/* allow for probabilistic-based expansion of C3AH6 */
/* crystal to account for volume stoichiometry */
pexp=ran1(seed);
if(pexp<=0.69){
extc3ah6(xcur,ycur,zcur);
}
}
}

/* examine reaction with diffusing gypsum to form ettringite */
/* Only allow reaction with diffusing gypsum */
if(check==DIFFGYP){
/* convert diffusing gypsum to ettringite */
mic[xnew][ynew][znew]=ETTR;
/* decrement counts of diffusing C3A and gypsum */
count[DIFFC3A]-=1;
count[DIFFGYP]-=1;
action=0;

/* convert diffusing C3A to solid ettringite or else leave as a diffusing C3A */
pexp=ran1(seed);
nexp=2;
if(pexp<=0.40){
mic[xcur][ycur][zcur]=ETTR;
nexp=1;
}
else{
mic[xcur][ycur][zcur]=DIFFC3A;
count[DIFFC3A]+=1;
/* indicate that diffusing species remains in current location */
action=7;
nexp=2;
}

/* Perform expansion that occurs when ettringite is formed */
/* xexp, yexp and zexp are the coordinates of the last ettringite */
/* pixel to be added */
xexp=xnew;
yexp=ynew;
zexp=znew;
for(iexp=1;iexp<=nexp;iexp++){
newact=extettr(xexp,yexp,zexp);
/* update xexp, yexp and zexp */
switch (newact){
case 1:
xexp-=1;
if(xexp<0){xexp=(SYSIZE-1);}
break;
case 2:
xexp+=1;
if(xexp>=SYSIZE){xexp=0;}
break;
case 3:
yexp-=1;
if(yexp<0){yexp=(SYSIZE-1);}
break;
case 4:
yexp+=1;
if(yexp>=SYSIZE){yexp=0;}
break;
case 5:
zexp-=1;
if(zexp<0){zexp=(SYSIZE-1);}
break;
case 6:
zexp+=1;
if(zexp>=SYSIZE){zexp=0;}
break;
default:
break;
}
}

/* probabilistic-based expansion for last ettringite pixel */
pexp=ran1(seed);
if(pexp<=0.30){
newact=extettr(xexp,yexp,zexp);
}
}
/* check for reaction with diffusing or solid ettringite to form AFm */
/* reaction at solid ettringite only possible if ettringite is soluble */
/* and even then on a limited bases to avoid a great formation of AFm */
/* when ettringite first becomes soluble */
pgrow=ran1(seed);
if((check==DIFFETTR)||((check==ETTR)&&(soluble[ETTR]==1)&&(pgrow<=C3AETTR))){
/* convert diffusing or solid ettringite to AFm */
mic[xnew][ynew][znew]=AFM;
/* decrement counts of diffusing C3A and ettringite */
count[DIFFC3A]-=1;
if(check==DIFFETTR){
count[DIFFETTR]-=1;
}
action=0;

/* convert diffusing C3A to AFm or leave at diffusing C3A */
pexp=ran1(seed);
if(pexp<=0.2424){
mic[xcur][ycur][zcur]=AFM;
pafm=(-0.1);
}
else{
mic[xcur][ycur][zcur]=DIFFC3A;
count[DIFFC3A]+=1;
action=7;
pafm=(0.278-0.2424)/(1.-0.2424);
}

/* probabilistic-based expansion for new AFm pixel */
pexp=ran1(seed);
if(pexp<=pafm){
extafm(xnew,ynew,znew);
}
}
if((action!=0)&&(action!=7)){

/* if diffusion is possible, execute it */
if(check==POROSITY){
mic[xcur][ycur][zcur]=POROSITY;
mic[xnew][ynew][znew]=DIFFC3A;
}
else{
/* indicate that diffusing C3A remained */
/* at original location */
action=7;
}
}
}
return(action);
}

/* routine to oversee hydration by updating position of all */
/* remaining diffusing species */
/* Calls movech, movec3a, movefh3, moveettr, movecsh, and movegyp */
void hydrate (fincyc,stepmax,chpar1,chpar2,hgpar1,hgpar2,fhpar1,fhpar2)
int fincyc,stepmax;
float chpar1,chpar2,hgpar1,hgpar2,fhpar1,fhpar2;
{
int xpl,ypl,zpl,phpl,xpnew,ypnew,zpnew;
float chprob,c3ah6prob,fh3prob;
long int icnt,nleft,ntodo,ndale;
int istep,termflag,reactf;
float beterm;
struct ants *curant,*antgone;

termflag=0;

/* Perform diffusion until all reacted or max. # of diffusion steps reached */
for(istep=1;((istep<=stepmax)&&(nleft>0));istep++){
if((fincyc==1)&&(istep==stepmax)){termflag=1;}

nleft=0;
ndale=0;

/* determine probabilities for CH and C3AH6 nucleation */
beterm=exp(-(double)(count[DIFFCH])/chpar2);
chprob=chpar1*(1.-beterm);
beterm=exp(-(double)(count[DIFFC3A])/hgpar2);
c3ah6prob=hgpar1*(1.-beterm);
beterm=exp(-(double)(count[DIFFFH3])/fhpar2);
fh3prob=fhpar1*(1.-beterm);

/* Process each diffusing species in turn */
while(curant!=NULL){
ndale+=1;
xpl=curant->x;
ypl=curant->y;
zpl=curant->z;
phpl=curant->id;

/* based on ID, call appropriate routine to process diffusing species */
switch (phpl) {
case DIFFCSH:
reactf=movecsh(xpl,ypl,zpl,termflag);
break;
case DIFFCH:
reactf=movech(xpl,ypl,zpl,termflag,chprob);
break;
case DIFFFH3:
reactf=movefh3(xpl,ypl,zpl,termflag,fh3prob);
break;
case DIFFGYP:
reactf=movegyp(xpl,ypl,zpl,termflag);
break;
case DIFFC3A:
reactf=movec3a(xpl,ypl,zpl,termflag,c3ah6prob);
break;
case DIFFETTR:
reactf=moveettr(xpl,ypl,zpl,termflag);
break;
default:
printf("Error in ID of phase \n");
break;
}

/* if no reaction */
if(reactf!=0){
nleft+=1;
xpnew=xpl;
ypnew=ypl;
zpnew=zpl;

/* update location of diffusing species */
switch (reactf) {
case 1:
xpnew-=1;
if(xpnew<0){xpnew=(SYSIZE-1);}
break;
case 2:
xpnew+=1;
if(xpnew>=SYSIZE){xpnew=0;}
break;
case 3:
ypnew-=1;
if(ypnew<0){ypnew=(SYSIZE-1);}
break;
case 4:
ypnew+=1;
if(ypnew>=SYSIZE){ypnew=0;}
break;
case 5:
zpnew-=1;
if(zpnew<0){zpnew=(SYSIZE-1);}
break;
case 6:
zpnew+=1;
if(zpnew>=SYSIZE){zpnew=0;}
break;
default:
break;
}

/* store new location of diffusing species */
curant->x=xpnew;
curant->y=ypnew;
curant->z=zpnew;
curant->id=phpl;
curant=curant->nextant;
} /* end of reactf!=0 loop */
/* else remove ant from list */
else{
if(ndale==1){
}
else{
(curant->prevant)->nextant=curant->nextant;
}
if(curant->nextant!=NULL){
(curant->nextant)->prevant=curant->prevant;
}
else{
tailant=curant->prevant;
}
antgone=curant;
curant=curant->nextant;
free(antgone);
ngoing-=1;
}
} /* end of curant loop */
ntodo=nleft;
} /* end of istep loop */
}
```

Dale P Bentz
Fri Feb 21 08:44:14 EST 1997