/* erod.c (simlib), 20.nov.2002, JH */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <grass/gis.h>
#include <grass/site.h>
#include <grass/bitmap.h>
#include <grass/linkm.h>

#include <grass/waterglobs.h>

/* divergence computation from a given field */

void erod(double **hw)
{
/* hw = sigma or gamma */

  double dyp, dyn, dya, dxp, dxn, dxa;
  int k,l;
  int l1, lp, k1, kp, ln, kn, k2, l2;

            
        for (k = 0; k < my; k++) {
            for (l = 0; l < mx; l++) {
		lp = max(0,l-2);
		l1 = lp + 1;
                kp = max(0,k-2);
                k1 = kp + 1;
                ln = min(mx-1,l+1);
                l2 = ln - 1;
                kn = min(my-1,k+1);
                k2 = kn - 1;

        if (zz[k][l] != UNDEF || zz[k][ln] != UNDEF || zz[kp][l] != UNDEF || \
			zz[k][lp] != UNDEF || zz[k][l1] != UNDEF || \
			zz[k1][l] != UNDEF || zz[kn][l] != UNDEF) { /* jh fix */

		dxp = (v1[k][lp] * hw[k][lp] - v1[k][l1] * hw[k][l1]) / stepx;
                dxn = (v1[k][l2] * hw[k][l2] - v1[k][ln] * hw[k][ln]) / stepx;
		dxa = 0.5 * (dxp + dxn);


                dyp = (v2[kp][l] * hw[kp][l] - v2[k1][l] * hw[k1][l]) / stepy;
                dyn = (v2[k2][l] * hw[k2][l] - v2[kn][l] * hw[kn][l]) / stepy;
                dya = 0.5 * (dyp + dyn);

		er[k][l] = (dxa + dya) / deltap;

	} else er[k][l] = UNDEF;

		}
	}
}