Added walker vector points output generation.

Removed confusing code, which seems to make things more complicted.
Added time series support to the r.sim.water spearfish example with vtk output.


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@45607 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/simwe/r.sim.sediment/main.c

@@ -161,6 +161,13 @@ int main(int argc, char *argv[])
     parm.maninval->description = _("Name of mannings n value");
     parm.maninval->guisection = _("Input");
 
+    parm.outwalk = G_define_standard_option(G_OPT_V_OUTPUT);
+    parm.outwalk->key = "outwalk";
+    parm.outwalk->required = NO;
+    parm.outwalk->description =
+	_("Base name of the output walkers vector points map");
+    parm.outwalk->guisection = _("Output_options");
+    
     /* needs to be updated to GRASS 6 vector format !! 
     parm.sfile = G_define_standard_option(G_OPT_V_INPUT);
     parm.sfile->key = "vector";
@@ -309,7 +316,7 @@ int main(int argc, char *argv[])
     conc = parm.conc->answer;
     flux = parm.flux->answer;
     erdep = parm.erdep->answer;
-/*    sfile = parm.sfile->answer; */
+    outwalk = parm.outwalk->answer; 
 
     /*      sscanf(parm.nwalk->answer, "%d", &maxwa); */
     sscanf(parm.niter->answer, "%d", &timesec);

raster/simwe/r.sim.water/main.c

@@ -212,14 +212,12 @@ int main(int argc, char *argv[])
     parm.err->description = _("Name for output simulation error raster map [m]");
     parm.err->guisection = _("Output");
 
-/*
     parm.outwalk = G_define_standard_option(G_OPT_V_OUTPUT);
     parm.outwalk->key = "outwalk";
     parm.outwalk->required = NO;
     parm.outwalk->description =
-	_("Name of the output walkers vector points map");
+	_("Base name of the output walkers vector points map");
     parm.outwalk->guisection = _("Output_options");
-*/
 
     parm.nwalk = G_define_option();
     parm.nwalk->key = "nwalk";
@@ -334,7 +332,7 @@ int main(int argc, char *argv[])
     depth = parm.depth->answer;
     disch = parm.disch->answer;
     err = parm.err->answer;
-/*    outwalk = parm.outwalk->answer; */
+    outwalk = parm.outwalk->answer; 
 /*    sfile = parm.sfile->answer; */
 
     sscanf(parm.niter->answer, "%d", &timesec);

raster/simwe/r.sim.water/spearfish.sh

@@ -3,43 +3,48 @@
 # sample script for Spearfish
 # written by
 #   andreas.philipp geo.uni-augsburg.de
-# modified by JH, HM, MN
+# modified by JH, HM, MN, SG
 
-dem=elevation.dem
-g.region rast=${dem} -p
+
+dem=elevation.10m
+output=simwe
+g.region rast=${output}
+g.region n=4920800 s=4917800 w=602500 e=606000
+g.region -p
 
 # Manning n
 man05=0.05
 
 # rainfall [mm/hr]
-rain01=5
+rain01=2.5
 
 # infitration [mm/hr]
 infil0=0.
 
 echo "Preparing input maps..."
-r.slope.aspect --o elevation=$dem dx=${dem}_dx dy=${dem}_dy
-r.mapcalc "${dem}_rain =if(${dem},$rain01,null())"
-r.mapcalc "${dem}_manin=if(${dem},$man05,null())"
-r.mapcalc "${dem}_infil=if(${dem},$infil0,null())"
-
-echo "r.sim.water --o elevin=${dem} dxin=${dem}_dx dyin=${dem}_dy \
-      rain=${dem}_rain manin=${dem}_manin infil=${dem}_infil \
-      depth=${dem}_depth disch=${dem}_disch err=${dem}_err"
-  
-r.sim.water --o elevin=${dem} dxin=${dem}_dx dyin=${dem}_dy \
-  rain=${dem}_rain manin=${dem}_manin infil=${dem}_infil \
-  depth=${dem}_depth disch=${dem}_disch err=${dem}_err
-
-r.info -r ${dem}_depth
-r.info -r ${dem}_disch
-r.info -r ${dem}_err
-
-echo "Written:
- Output water depth raster file:     ${dem}_depth
- Output water discharge raster file: ${dem}_disch
- Output simulation error raster file: ${dem}_err
-"
+r.slope.aspect --o elevation=$dem dx=${output}_dx dy=${output}_dy
+r.mapcalc --o expr="${output}_rain =if(${dem},$rain01,null())"
+r.mapcalc --o expr="${output}_manin=if(${dem},$man05,null())"
+r.mapcalc --o expr="${output}_infil=if(${dem},$infil0,null())"
+
+echo "r.sim.water --o -t elevation=${dem} dx=${output}_dx dy=${output}_dy rain=${output}_rain man=${output}_manin infil=${output}_infil depth=${output}_depth disch=${output}_disch err=${output}_err outwalk=${output}_walker niter=5 outiter=1"
+      r.sim.water --o -t elevation=${dem} dx=${output}_dx dy=${output}_dy rain=${output}_rain man=${output}_manin infil=${output}_infil depth=${output}_depth disch=${output}_disch err=${output}_err outwalk=${output}_walker niter=5 outiter=1
+
+echo "Build topology and exporting walker vector points for each time step"
+for i in `g.mlist vect | grep ${output}` ; do
+	v.build $i
+	echo "v.out.vtk input=$i output=$i.vtk"
+	v.out.vtk input=$i output=$i.vtk
+done
+
+echo "Exporting the raster maps for each time step"
+for i in `g.mlist rast | grep ${output}` ; do
+	echo "r.out.vtk elevation=${dem} input=$i output=$i.vtk null=0.0"
+	r.out.vtk elevation=${dem} input=$i output=$i.vtk null=0.0
+done
+
+echo "Now open paraview from this command line and load the vtk files as time series using the File->Open menu entry"
+echo "Step throu the time steps and adjust the color tables"
 
 # cleanup
-g.remove --q rast=${dem}_dx,${dem}_dy,${dem}_rain,${dem}_manin,${dem}_infil
+g.remove --q rast=${output}_dx,${output}_dy,${output}_rain,${output}_manin,${output}_infil

raster/simwe/simlib/Makefile

@@ -1,7 +1,7 @@
 MODULE_TOPDIR = ../../..
 
-LIBES = $(GISLIB) $(DATETIMELIB) $(GMATHLIB) $(RASTERLIB)
-DEPENDENCIES = $(GISDEP) $(BITMAPDEP) $(DBMIDEP) $(GMATHDEP) $(LINKMDEP) $(XDRDEP) $(SITESDEP) $(VECTDEP) $(RASTERDEP)
+LIBES = $(GISLIB) $(DATETIMELIB) $(GMATHLIB) $(RASTERLIB) $(VECTLIB)
+DEPENDENCIES = $(GISDEP) $(BITMAPDEP) $(DBMIDEP) $(GMATHDEP) $(LINKMDEP) $(XDRDEP) $(SITESDEP) $(VECTDEP) $(RASTERDEP) $(VECTDEP)
 EXTRA_INC = $(VECT_INC)
 EXTRA_CFLAGS = $(VECT_CFLAGS)
 

raster/simwe/simlib/erod.c

@@ -19,7 +19,6 @@
 #include <stdlib.h>
 #include <math.h>
 #include <grass/gis.h>
-/* #include <grass/site.h> */
 #include <grass/bitmap.h>
 #include <grass/linkm.h>
 

raster/simwe/simlib/hydro.c

@@ -29,6 +29,21 @@
 struct options parm;
 struct flags flag;
 
+/*
+ * Soeren 8. Mar 2011 TODO:
+ * Put all these global variables into several meaningful structures and 
+ * document use and purpose.
+ * 
+ * Example:
+ * Put all file descriptors into a input_files struct and rename the variables:
+ * input_files.elev 
+ * input_files.dx 
+ * input_files.dy 
+ * input_files.drain
+ * ... 
+ * 
+ */
+
 FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps,
     *fdmanin, *fddepth, *fddisch, *fderr;
 FILE *fdwdepth, *fddetin, *fdtranin, *fdtauin, *fdtc, *fdet, *fdconc,
@@ -46,7 +61,7 @@ char *manin;
 char *depth;
 char *disch;
 char *err;
-/* char *outwalk; */
+char *outwalk; 
 char *mapset;
 char *mscale;
 char *tserie;
@@ -69,12 +84,6 @@ struct seed seed;
 
 struct Cell_head cellhd;
 
-/*
-struct Point *points;
-int npoints;
-int npoints_alloc;
-*/
-
 double xmin, ymin, xmax, ymax;
 double mayy, miyy, maxx, mixx;
 int mx, my;
@@ -93,12 +102,10 @@ double **gama, **gammas, **si, **inf, **sigma;
 float **dc, **tau, **er, **ct, **trap;
 float **dif;
 
-/* double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3]; */
-double vavg[MAXW][2], w[MAXW][3];
+double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3]; 
 int iflag[MAXW];
 
 double hbeta;
-/* int ldemo; */
 double hhmax, sisum, vmean;
 double infsum, infmean;
 int maxw, maxwa, nwalk;
@@ -106,10 +113,9 @@ double rwalk, bresx, bresy, xrand, yrand;
 double stepx, stepy, xp0, yp0;
 double chmean, si0, deltap, deldif, cch, hhc, halpha;
 double eps;
-/* int maxwab, nstack; */
-int maxwab;
+int maxwab, nstack;
 int iterout, mx2o, my2o;
-int miter, nwalka, lwwfin;
+int miter, nwalka;
 double timec;
 int ts, timesec;
 
@@ -128,8 +134,7 @@ void main_loop(void)
 {
 
     int i, ii, l, k;
-/*    int icoub, lwout, nmult; */
-    int icoub, nmult;
+    int icoub, nmult; 
     int iw, iblock, lw;
     int itime, iter1;
     int nfiterh, nfiterw;
@@ -149,7 +154,7 @@ void main_loop(void)
     nblock = 1;
     icoub = 0;
     icfl = 0;
-/*    nstack = 0; */
+    nstack = 0; 
 
     if (maxwa > (MAXW - mx * my)) {
 	mitfac = maxwa / (MAXW - mx * my);
@@ -168,7 +173,6 @@ void main_loop(void)
 	sarea = bresx * bresy;
 	G_debug(2, " barea,sarea,rwalk,sisum: %f %f %f %f", barea, sarea,
 		rwalk, sisum);
-	lwwfin = 0;
 	/* write hh.walkers0 */
 
 	for (k = 0; k < my; k++) {
@@ -225,20 +229,11 @@ void main_loop(void)
 			    iflag[lw] = 1;
 			}
 
-/*
-			lwout = lw / ldemo;
-			lwout *= ldemo;
-
-			if (lwout == lw) {
-			    ++lwwfin;
-			}
-*/
 		    }
 		}		/*DEFined area */
 	    }
 	}
 	nwalk = lw;
-	G_debug(2, " number of written walkers: %d", lwwfin);
 	G_debug(2, " nwalk, maxw %d %d", nwalk, MAXW);
 	G_debug(2, " walkwe (walk weight),frac %f %f", walkwe, frac);
 
@@ -283,6 +278,7 @@ void main_loop(void)
 		addac = factor * .5;
 	    }
 	    nwalka = 0;
+	    nstack = 0;
 
 	    for (lw = 1; lw <= nwalk; lw++) {
 		if (w[lw][3] > EPS) {	/* check the walker weight */
@@ -302,9 +298,7 @@ void main_loop(void)
 		    }
 
 		    if (zz[k][l] != UNDEF) {
-
 			if (infil != NULL) {	/* infiltration part */
-
 			    if (inf[k][l] - si[k][l] > 0.) {
 
 				decr = pow(addac * w[lw][3], 3. / 5.);	/* decreasing factor in m */
@@ -317,12 +311,9 @@ void main_loop(void)
 				    inf[k][l] = 0.;
 
 				}
-
 			    }
-
 			}
 
-
 			gama[k][l] += (addac * w[lw][3]);	/* add walker weigh to water depth or conc. */
 
 			d1 = gama[k][l] * conn;
@@ -348,10 +339,8 @@ void main_loop(void)
 				velx = -0.1 * v1[k][l];	/* move it slightly back */
 				vely = -0.1 * v2[k][l];
 			    }
-
 			}
 
-
 			gaux = gasdev();
 			gauy = gasdev();
 
@@ -380,59 +369,49 @@ void main_loop(void)
 			w[lw][3] = 1e-10;	/* eliminate walker if it is out of area */
 		    }
 		}
-
-/* output walkers commented out */
-		/*        write the walkers which reach the box */
-/*
-		if (w[lw][1] >= xmin && w[lw][2] >= ymin && w[lw][1]
-		    <= xmax && w[lw][2] <= ymax && iflag[lw] == 0) {
-
-		    lwout = lw / ldemo;
-		    lwout *= ldemo;
-		    if (lwout == lw && (i == miter || i == iter1)) {	
-			stack[nstack][1] = mixx / conv + w[lw][1] / conv;
-			stack[nstack][2] = miyy / conv + w[lw][2] / conv;
-			stack[nstack][3] = w[lw][3];
-
-			nstack++;
-*/
-			/*                  iflag[lw] = 1; */
-/*
-		    }
-		}
-*/
-
-		/* walker is leaving the box */
-
-		/*              if (w[lw][1] <= mixx && w[lw][2] <= bymi && w[lw][1] 
-		   >= bxma && w[lw][2] >= byma && iflag[lw] == 1) {
-		   iflag[lw] = 0;
-		   } */
-
-		/* every iterout iteration write out the selected ldemo walkers */
-		/* and water depth and time */
-
-
-	    }			/* lw loop */
-
-	    if (i == iter1) {
-
-		/* call output for iteration output */
-
-		if (ts == 1) {
-		    if (erdep != NULL)
-			erod(gama);	/* divergence of gama field */
-
-		    conn = (double)nblock / (double)iblock;
-		    itime = (int)(i * deltap * timec);
-		    ii = output_data(itime, conn);
-/*		    nstack = 0; */
-		    if (ii != 1)
-			G_fatal_error(_("Unable to write raster maps"));
-		}
-
+            } /* lw loop */
+            
+            /* Changes made by Soeren 8. Mar 2011 to replace the site walker output implementation */
+            /* Save all walkers located within the computational region and with valid 
+               z coordinates */
+            if ((i == miter || i == iter1)) {	
+                nstack = 0;
+                
+                for (lw = 1; lw <= nwalk; lw++) {
+                    /* Compute the  elevation raster map index */
+                    l = (int)((w[lw][1] + stxm) / stepx) - mx - 1;
+                    k = (int)((w[lw][2] + stym) / stepy) - my - 1;
+                    
+		    /* Check for correct elevation raster map index */
+		    if(l < 0 || l >= mx || k < 0 || k >= my)
+			 continue;
+
+                    if (w[lw][3] > EPS && zz[k][l] != UNDEF) {
+
+                        /* Save the 3d position of the walker */
+                        stack[nstack][1] = mixx / conv + w[lw][1] / conv;
+                        stack[nstack][2] = miyy / conv + w[lw][2] / conv;
+                        stack[nstack][3] = zz[k][l];
+
+                        nstack++;
+                    }
+                }
+            } /* lw loop */
+
+	    if (i == iter1 && ts == 1) {
+            /* call output for iteration output */
+                if (erdep != NULL)
+                    erod(gama);	/* divergence of gama field */
+
+                conn = (double)nblock / (double)iblock;
+                itime = (int)(i * deltap * timec);
+                ii = output_data(itime, conn);
+                if (ii != 1)
+                    G_fatal_error(_("Unable to write raster maps"));
 	    }
 
+            /* Soeren 8. Mar 2011 TODO:
+             *  This hould be replaced by vector functionality and sql database storage */
 /* ascii data site file output for gamma  - hydrograph or sediment*/
 /* cchez incl. sqrt(sinsl) */
 /* sediment */
@@ -469,6 +448,7 @@ void main_loop(void)
 	}			/* miter */
 
       L_800:
+      /* Soeren 8. Mar 2011: Why is this commented out?*/
 	/*        if (iwrib != nblock) {
 	   icount = icoub / iwrib;
 

raster/simwe/simlib/input.c

@@ -18,7 +18,7 @@
 
 /*!
  * \brief allocate memory, read input rasters, assign UNDEF to NODATA
- *
+ * 
  *  \return int
  * sites related input/output commented out - needs update to vect, HM nov 2008
  */
@@ -32,7 +32,7 @@ int input_data(void)
     DCELL *dxin_cell, *dyin_cell, *rain_cell, *infil_cell, *wdepth_cell;
     int elevin_fd, dxin_fd, dyin_fd, rain_fd, infil_fd, traps_fd, manin_fd, row, row_rev;
     int detin_fd, trainin_fd, tauin_fd, wdepth_fd;
-    int l, j;
+    int j;
 /*    int nn, cc, ii, dd; */
     double unitconv = 0.0000002;	/* mm/hr to m/s */
     const char *mapset;
@@ -116,7 +116,6 @@ int input_data(void)
 	wdepth_cell = Rast_allocate_d_buf();
 
     /* Allocate some double dimension arrays for each input */
-    /* with length of matrix Y */
     zz = G_alloc_fmatrix(my, mx);
     v1 = G_alloc_matrix(my, mx);
     v2 = G_alloc_matrix(my, mx);
@@ -143,7 +142,7 @@ int input_data(void)
     if (wdepth != NULL)
 	gama = G_alloc_matrix(my, mx);
 
-    G_debug(3, "Running JAN 2010 version, started modifications on 20080211");
+    G_debug(3, "Running MAR 2011 version, started modifications on 20080211");
 
     /* Check if data available in mapsets
      * if found, then open the files */

raster/simwe/simlib/output.c

@@ -5,13 +5,70 @@
 #include <math.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
-/* #include <grass/site.h> */
 #include <grass/bitmap.h>
 #include <grass/linkm.h>
+#include <grass/vector.h>
 
 #include <grass/waterglobs.h>
 #include <grass/glocale.h>
 
+
+static void output_walker_as_vector(int tt, int ndigit);
+
+
+/* This function was added by Soeren 8. Mar 2011 
+ * It replaces the site walker output implementation */
+void output_walker_as_vector(int tt, int ndigit)
+{
+    char buf[256];
+    char *outwalk_time = NULL; 
+    double x, y, z;
+    struct Map_info Out;
+    struct line_pnts *Points;
+    struct line_cats *Cats;
+    int i;
+
+    if (outwalk != NULL) {
+
+	/* In case of time series we extent the output name with the time value */
+	if (ts == 1) {
+	    sprintf(buf, "%s_%.*d", outwalk, ndigit, tt);
+	    outwalk_time = G_store(buf);
+	    Vect_open_new(&Out, outwalk_time, WITH_Z);
+	    G_message("Writing %i walker into vector file %s", nstack, outwalk_time);
+	}
+	else {
+	    Vect_open_new(&Out, outwalk, WITH_Z);
+	    G_message("Writing %i walker into vector file %s", nstack, outwalk);
+	}
+
+	Points = Vect_new_line_struct();
+	Cats = Vect_new_cats_struct();
+
+	for (i = 0; i < nstack; i++) {
+	    x = (float)stack[i][1];
+	    y = (float)stack[i][2];
+	    z = (float)stack[i][3];
+
+	    Vect_reset_line(Points);
+	    Vect_reset_cats(Cats);
+			    
+	    Vect_cat_set(Cats, 1, i + 1);
+	    Vect_append_point(Points, x, y, z);
+	    Vect_write_line(&Out, GV_POINT, Points, Cats);
+	}
+	/* Close vector file */
+        Vect_close(&Out);
+                
+        Vect_destroy_line_struct(Points);
+        Vect_destroy_cats_struct(Cats);
+    }
+    
+    return;
+}
+
+/* Soeren 8. Mar 2011 TODO: 
+ * This function needs to be refractured and splittet into smaller parts */
 int output_data(int tt, double ft)
 {
 
@@ -25,7 +82,6 @@ int output_data(int tt, double ft)
     struct History hist, hist1;	/* hist2, hist3, hist4, hist5 */
     char *depth0 = NULL, *disch0 = NULL, *err0 = NULL;
     char *conc0 = NULL, *flux0 = NULL;
-/*    char *erdep0 = NULL, *outwalk0 = NULL; */
     char *erdep0 = NULL;
     const char *mapst = NULL;
     char *type;
@@ -33,9 +89,6 @@ int output_data(int tt, double ft)
     int ndigit;
     FCELL dat1, dat2;
     float a1, a2;
-/*    Site_head walkershead;
-    Site *sd;
-*/
 
     ndigit = 2;
     if (timesec >= 10)
@@ -46,50 +99,11 @@ int output_data(int tt, double ft)
 	ndigit = 5;
     if (timesec >= 10000)
 	ndigit = 6;
+    
+    /* Write the output walkers */
+    output_walker_as_vector(tt, ndigit);
 
     Rast_set_window(&cellhd);
-
-/*
-    if (outwalk) {
-	if (ts == 1) {
-	    sprintf(buf, "%s%.*d", outwalk, ndigit, tt);
-	    outwalk0 = G_store(buf);
-	    fdoutwalk = G_fopen_sites_new(outwalk0);
-	}
-	else
-	    fdoutwalk = G_fopen_sites_new(outwalk);
-
-	if (fdoutwalk == NULL)
-	    G_fatal_error("Cannot open %s", outwalk);
-	else {
-	    char buf[GNAME_MAX + 40];
-
-	    if (NULL == (sd = G_site_new_struct(-1, 2, 0, 1)))
-		G_fatal_error("memory allocation failed for site");
-
-	    if (ts == 1)
-		walkershead.name = outwalk0;
-	    else
-		walkershead.name = outwalk;
-
-	    sprintf(buf, "output walkers of %s [raster]", depth);
-	    walkershead.desc = G_store(buf);
-	    walkershead.time = NULL;
-	    walkershead.stime = NULL;
-	    walkershead.labels = NULL;
-	    walkershead.form = NULL;
-
-	    G_site_put_head(fdoutwalk, &walkershead);
-
-	    for (i = 0; i < nstack; i++) {
-		sd->east = (float)stack[i][1];
-		sd->north = (float)stack[i][2];
-		sd->fcat = (float)stack[i][3];
-		G_site_put(fdoutwalk, sd);
-	    }
-	}
-    }
-*/
     if (depth) {
 	depth_cell = Rast_allocate_f_buf();
 	if (ts == 1) {
@@ -477,8 +491,8 @@ int output_data(int tt, double ft)
 	    &hist, "duration (sec.)=%d, time-serie iteration=%d",
 	    timesec, tt);
 	Rast_append_format_history(
-	    &hist, "written walkers=%d, deltap=%f, mean vel.=%f",
-	    lwwfin, deltap, vmean);
+	    &hist, "written deltap=%f, mean vel.=%f",
+	    deltap, vmean);
 	Rast_append_format_history(
 	    &hist, "mean source (si)=%e, mean infil=%e",
 	    si0, infmean);
@@ -517,8 +531,8 @@ int output_data(int tt, double ft)
 	    &hist, "duration (sec.)=%d, time-serie iteration=%d",
 	    timesec, tt);
 	Rast_append_format_history(
-	    &hist, "written walkers=%d, deltap=%f, mean vel.=%f",
-	    lwwfin, deltap, vmean);
+	    &hist, "written deltap=%f, mean vel.=%f",
+	    deltap, vmean);
 	Rast_append_format_history(
 	    &hist, "mean source (si)=%e, mean infil=%e",
 	    si0, infmean);
@@ -557,8 +571,8 @@ int output_data(int tt, double ft)
 	    &hist, "duration (sec.)=%d, time-serie iteration=%d",
 	    timesec, tt);
 	Rast_append_format_history(
-	    &hist, "written walkers=%d, deltap=%f, mean vel.=%f",
-	    lwwfin, deltap, vmean);
+	    &hist, "written deltap=%f, mean vel.=%f",
+	    deltap, vmean);
 	Rast_append_format_history(
 	    &hist, "mean source (si)=%f", si0);
 

raster/simwe/simlib/waterglobs.h

@@ -28,7 +28,7 @@ extern char *manin;
 extern char *depth;
 extern char *disch;
 extern char *err;
-/* extern char *outwalk; */
+extern char *outwalk;
 extern char *mapset;
 extern char *mscale;
 extern char *tserie;
@@ -120,12 +120,10 @@ extern double **gama, **gammas, **si, **inf, **sigma;
 extern float **dc, **tau, **er, **ct, **trap;
 extern float **dif;
 
-/* extern double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3]; */
-extern double vavg[MAXW][2], w[MAXW][3];
+extern double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3]; 
 extern int iflag[MAXW];
 
 extern double hbeta;
-/* extern int ldemo; */
 extern double hhmax, sisum, vmean;
 extern double infsum, infmean;
 extern int maxw, maxwa, nwalk;
@@ -133,10 +131,9 @@ extern double rwalk, bresx, bresy, xrand, yrand;
 extern double stepx, stepy, xp0, yp0;
 extern double chmean, si0, deltap, deldif, cch, hhc, halpha;
 extern double eps;
-/* extern int maxwab, nstack; */
-extern int maxwab;
+extern int maxwab, nstack; 
 extern int iterout, mx2o, my2o;
-extern int miter, nwalka, lwwfin;
+extern int miter, nwalka;
 extern double timec;
 extern int ts, timesec;