//  -*- Mode: C++;  -*-
//  File: ravel_class.cc
//  Author: Dino Bellugi (dino@Moray.Berkeley.EDU)
//  Copyright (C) Dino Bellugi, 1996
// *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// * FUNCTION:	Produces a grid containing the volume of soil moved by
// *		ravelling into each cell that is an accumulation point.
// *		Each cell contributing to an accumulation point will 
// *		contain the value of the area contributing up to that
// *		cell (uphill from it), identified by a negative value.
// *		Ravel is distributed in 8 directions proportionally
// *		to slope, when slope is greater than a specified (input)
// *		threshold.
// *
// * CLASSES:	Ravel_Class
// * 
// * RELATED PACKAGES: Grid
// *
// * HISTORY:
// * Created:	Thu Oct 24 12:48:27 1996 (dino)
// *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

#include <fstream.h>
#include <stdio.h>
#include "ravel_class.h"
#include "macros.h"

/* 
 * find the contributing area given the elevation grid
 * initilize to 0.0 then make the gridedges then find
 * contributing area
 */

#if (STREAMS)
Ravel_Class::Ravel_Class(Grid &inz, Grid &inc, double critical_slope_in,
	double ravel_constant_in) : Grid(inz) 
#else
Ravel_Class::Ravel_Class(Grid &inz, double critical_slope_in,
	double ravel_constant_in) : Grid(inz) 
#endif
{
	z = &inz;
#if (STREAMS)
	c = &inc;
#endif
	critical_slope = critical_slope_in;
	ravel_constant = ravel_constant_in;

#if (MESSAGES)
	cerr << "Initializing grid ...\n";
#endif

	for(int i = 0; i < xdim; i++) 
	{
		for(int j = 0; j < ydim; j++) 
		{
			at(i,j) = 0.0;
		}
	}

#if (MESSAGES)
	cerr << "Checking for sinks ...\n";
#endif

	checkForSinks(inz);

#if (MESSAGES)
	cerr << "Making weights ...\n";
#endif

	weights = new GridEdges(z->xdim, z->ydim);
	makeRavelWeights();

#if (MESSAGES)
	cerr << "Computing ravelling ...\n";
#endif

	initRavel();
}

Ravel_Class::~Ravel_Class() 
{
	delete weights;
}

/*
 * Produces a warning if there are sinks in the Grid.
 */

int Ravel_Class::checkForSinks(Grid &z) 
{
	int sinks = 0;

	for(int j = 1;j < z.ydim - 1; j++)  
	{ 
		for(int i = 1; i < z.xdim - 1; i++) 
		{
			if(z.isSink(i,j)) 
			{
				sinks = 1;
				break;
			}
		}
	}

#if (MESSAGES)
	if (sinks) 
		cerr << "WARNING: There are sinks in the elevation grid!!\n";
#endif
	return(sinks);
}

/*
 *  function initializes "weights". It sets the
 *  weight of each edge propotional to the slope
 *  between two neighboring cells, when this slope
 *  is greater than a specified threshold. The edge 
 *  weights are normalized so that the sum of all
 *  the weights leaving a cell = 1.0.
*/

void Ravel_Class::makeRavelWeights() 
{
#if (MESSAGES)
	unsigned long counter = 0L, ncells = z->xdim * z->ydim;
	int perc_comp = 0, last_comp = 0;

setbuf(stderr,NULL);

	cerr << "Completed with: 0 %\r";
#endif

	double slope[8] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
	double slptot = 0.0, slppct, val1, val2;
	int i, j, h, k;

	for (j = 0; j < z->ydim; j++) 
	{
    		for (i = 0; i < z->xdim; i++) 
    		{
#if (MESSAGES)
			// Let the user know how far we are
			perc_comp = int(((double)counter / (double)ncells) * 100.0);
			if (last_comp != perc_comp)
			{
				last_comp = perc_comp;
				cerr << "Completed with: " << perc_comp << " %\r";
			}
			++counter;
#endif

      			if ((val1 = z->at(i,j)) == z->nodata)
			{
				for (k = 0;  k< 8; ++k) 
	  				weights->set(i, j, (Direction)k, 0.0);
				continue;
			}

			for (slptot = 0.0, k = 0; k < 8; k++) 
			{
	  			if ((val2 = z->get(i,j,(Direction) k)) !=
					z->nodata) 
				{
	    				if ((slope[k] = (val1 - val2) /
 						(z->spacing *
						angleFactor((Direction) k))) >=
							critical_slope)
					{
	      					slptot += slope[k];
	    				}
	  			}
			}


 			for (k = 0; k < 8; k++) 
			{
	  			if ((val2 = z->get(i,j,(Direction)k)) !=
					z->nodata) 
				{
					if (slope[k] >= critical_slope) 
					{
						if (slptot != 0.0)
	      						slppct = slope[k] / slptot;
						else
							slppct = 0.0;
	      					weights->set(i,j,(Direction)k, -slppct);
	    				}
	  			}
			}

    		}
  	}
#if (MESSAGES)
	fprintf(stderr, "Done making weights!\n");
#endif
}

/*
 * Recursively add up all the contributing cells area to the
 * cell i,j.
 */

double Ravel_Class::getRavelArea(int i,int j) 
{
	double ravelArea, weight, rflag = 0;

	if ((ravelArea = at(i, j)) == 0.0) 
	{
		ravelArea = 0.0;

		if ((weight = weights->get(i,j,L)) > 0.0) 
    		{
      			ravelArea += (weight * getRavelArea(i-1,j));
			rflag = 1;
    		}

    		if ((weight = weights->get(i,j,R)) > 0.0) 
    		{
      			ravelArea += (weight * getRavelArea(i+1,j));
			rflag = 1;
    		}

    		if ((weight = weights->get(i,j,U)) > 0.0) 
    		{
    			ravelArea += (weight * getRavelArea(i,j+1));
 			rflag = 1;
   		}

    		if ((weight = weights->get(i,j,D)) > 0.0)
    		{
    			ravelArea += (weight * getRavelArea(i,j-1));
			rflag = 1;
    		}

    		if ((weight = weights->get(i,j,UL)) > 0.0) 
    		{

      			ravelArea += (weight * getRavelArea(i-1,j+1));
			rflag = 1;
    		}

    		if ((weight = weights->get(i,j,UR)) > 0.0) 
    		{
      			ravelArea += (weight * getRavelArea(i+1,j+1));
			rflag = 1;
    		}

    		if ((weight = weights->get(i,j,DL)) > 0.0) 
    		{
    			ravelArea += (weight * getRavelArea(i-1,j-1));
			rflag = 1;
    		}

    		if ((weight = weights->get(i,j,DR)) > 0.0) 
    		{
      			ravelArea += (weight * getRavelArea(i+1,j-1));
			rflag = 1;
    		}

		// this eliminates ravelling caused by one cell only
		// something to think about ...
		if (rflag)
			ravelArea += sqr(z->spacing);

		// flip the sign to identify a trasitory value
		// and to keep track of contributing area
    		at(i,j) = -ravelArea;
  	}

	else if (ravelArea < 0.0) 
	  	return(-ravelArea);

  	return(ravelArea);
}

void Ravel_Class::initRavel() 
{
#if (MESSAGES)
	unsigned long counter = 0L, ncells = z->xdim * z->ydim;
	int perc_comp = 0, last_comp = 0;

	cerr << "Completed with: 0 %\r";
#endif

	int i, j, k;

	// Initialize ravel to each cell with no downhill neighbor 
	// towards which slope >= critical_slope

	for (j = 0; j < z->ydim; ++j) 
  	{ 
    		for (i = 0; i < z->xdim; ++i) 
    		{
#if (MESSAGES)
			// Let the user know how far we are
			perc_comp = int(((double)counter / (double)ncells) * 100.0);
			if (last_comp != perc_comp)
			{
				last_comp = perc_comp;
				cerr << "Completed with: " << perc_comp << " %\r";
			}
			++counter;
#endif
			if (isAccPoint(i, j))
			{
#if (MESSAGES == 2)
				cerr << "Accumulation point at (" << i << ", " << j << ")\n";
#endif
				at(i, j) = getRavelArea(i,j) * ravel_constant;
			}
    		}
  	}

#if (MESSAGES)
	cerr << "Setting nodata and zero values ...\n";
#endif

  	for (j = 0; j < z->ydim; ++j) 
  	{ 
    		for (i = 0; i < z->xdim; ++i) 
    		{
			// nodata values
      			if (z->at(i,j) == z->nodata) 
				at(i,j) = nodata;

#if (!AREA)
			//transitory value have a negative sign
      			if (at(i,j) < 0.0) 
				at(i,j) = 0.0;
#endif
    		}
  	}
}

int Ravel_Class::isAccPoint(int i, int j)
{
	int k, flg = 0;
	double val1, val2, slope;
 
	if ((val1 = z->at(i, j)) == z->nodata) 
		return(0);

 	for (k = 0; k < 8; k++) 
	{
	  	if ((val2 = z->get(i, j,(Direction)k)) != z->nodata) 
		{
			// check downhill
	    		if ((slope = (val1 - val2) / 
				(z->spacing *
				angleFactor((Direction)k))) >=
				critical_slope) 
				return(0);	  	

			// check uphill
	    		if (slope <= -(critical_slope))
				flg = 1; 
		}
	}

	return(flg);
}

#if (STREAMS)

int Ravel_Class::isUpStream(int i, int j, int k) 
{
	double elev, elev0;

	elev = z->at(i, j);
	c->goTo(i, j);
	c->goTo(i, j, (Direction)k);
	c->followEnd();
	elev0 = z->at(c->xcur, c->ycur);
	c->goTo(i, j);
	return(elev0 > elev ? 1 : 0);
}

int Ravel_Class::isDownStream(int i, int j, int k) 
{
	double elev, elev0;

 	elev = z->at(i, j);
	c->goTo(i, j);
	c->goTo(i, j, (Direction)k);
	c->followEnd();
	elev0 = z->at(c->xcur, c->ycur);
	c->goTo(i, j);
	return(elev > elev0 ? 1 : 0);
}

int Ravel_Class::isDownStream2(int i, int j, int k) 
{
	double elev, elev0;

 	elev = z->at(i, j);
	c->goTo(i, j);
	c->goTo(i, j, (Direction)k);

	for (int i = 0; i < CKCELLS; ++i)
		if (!c->followNext())
			break;

	elev0 = z->at(c->xcur, c->ycur);
	c->goTo(i, j);
	return(elev > elev0 ? 1 : 0);
}

#endif