//  -*- Mode: C++;  -*-
//  File: ckravel.cc
//  Author: Dino Bellugi (dino@geomorph.berkeley.edu)
//  Copyright  Dino Bellugi, BlueG SoftWear, U.C. Berkeley, 1999 (C)
// *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// * FUNCTION: Checks watershed ravel values
// * 
// * RELATED PACKAGES: grid.cc stabfun.cc
// *
// * HISTORY:
// * Created:  Thu Mar 27 1997 (dino)
// * Modified: Mon May 10 1999 (dino)
// *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

// the following include file (notice.h) contains important copyright 
// information, as well as other legal statements. it must be present 
// in all files and in any distribution. removal or modification is
// strictly forbidden and violates u.s. and international copyright law.
#include "notice.h"
char *copyright_notice = "(C) 1999 Dino Bellugi, BlueG SoftWear, U.C. Berkeley";
//
#include <stdio.h>
#include <stdlib.h>
#include <fstream.h>
#include <string.h>
#include <math.h>
#include "grid.h"
#include "macros.h"

#define ELEVATION	0
#define WATERSHED	1

main(int argc, char **argv) 
{
	Grid zgrd, wgrd;
	int i, j, k, cells = 0;
	double slope, slopeval, val1, val2;

	if (argc != 4) 
	{
                // invalid arguments, print usage
                cerr << "usage:\n";
                cerr << "\t" << appName(argv) << " <elevation grid> <watershed grid> <critical slope>\n";
                cerr << "example:\n";
                cerr << "\t" << appName(argv) << " elev.grd watershed.grd 33\n";

    		exit(1);
  	}

	cerr << "Reading elevation grid ...\n";

	if (!zgrd.readGrid(argv[1]))
	{
		cerr << "can't read ravel grid " << argv[1] << " ...\n";
		exit(1);
	}

	cerr << "Reading watershed grid ...\n";

	if (!wgrd.readGrid(argv[2]))
	{
		cerr << "can't read watershed grid " << argv[2] << " ...\n";
		exit(1);
	}

	slope = tan(atof(argv[3]));

	for (i = 0; i < zgrd.xdim; ++i)
	{
		for (j = 0; j < zgrd.ydim; ++j)
		{
	  		if (wgrd.isValue(i, j) && (val1 = zgrd.at(i, j)) != zgrd.nodata) 
			{
 				for (k = 0; k < 8; k++) 
				{
	  				if (wgrd.isValue(i, j, (Direction)k) &&
						(val2 = zgrd.at(i, j,(Direction)k)) != zgrd.nodata) 
					{
						// check downhill
	  					if ((slopeval = (val1 - val2) / 
							(zgrd.spacing * angleFactor((Direction)k))) >= slope)
						{
							++cells; 
							break;	  	
						}
					}
				}
			}
		}
	}  
	
	slopeval = slope * 180 / VAL_PI;

	cerr << "Total ravel producing cells for " << slopeval << " degree slope: " << cells << "\n";

	exit(0);
}