//program file name <stream.cpp>
#include <iostream>
#include <fstream>
#include <cstring>
#include <cstdlib>
#include <cmath>

using namespace std;

//declare large arrays as global.
	float	v1[1202][1202];		//raster data 1 for input
	float	v2[1202][1202];		//raster data 2 for output

int main(int argc, char *argv[]) {
	int	i, j;
	char	rn[128];		//region name
	char	vn[128];		//variable name
	char	un[128];		//name of the unit
	int	dt;			//type of data (not in use)
	int	xas, yas;		//number of elements
	float	xs, ys;			//size of the region
	char	str[128];		//separator

	int	ii, jj;
	float	ux, uy, uxy;		//size of the cells - x, y, and diagonal
	float	gr, grmin;		//gradient and its minimum
	int	d;			//drain direction

	//initializing two raster data as zero.
	for (j=0; j<=1201; j++) {
		for (i=0; i<=1201; i++) {
			v1[i][j] = 0;
			v2[i][j] = 0;
		}
	}

	//exit if the number of arguments is not 1.
	if (argc != 2) {
		cerr << "Usage: stream <filename>\n";
		return 1;
	}
	
	//open the input file. exit if an error occurs.
	ifstream fin(argv[1]);
	if (!fin) {
		cerr << "Cannot open file.\n";
		return 1;
	}

	//read data.
	fin.getline(rn, 80);
	fin.getline(vn, 80);
	fin.getline(un, 80);
	fin >> dt >> xas >> yas;
	fin >> xs >> ys >> str;
	for (j=1; j<=yas; j++) {
		for (i=1; i<=xas; i++) {
			fin >> v1[i][j];
		}
	}
	fin.close();

	//make an extrapolation along the four edges of the data.
	for (j=1; j<=yas; j++) {
		v1[0][j] = v1[1][j] * 2 - v1[2][j];
		v1[xas + 1][j] = v1[xas][j] * 2 - v1[xas - 1][j];
	}

	for (i=1; i<=xas; i++) {
		v1[i][0] = v1[i][1] * 2 - v1[i][2];
		v1[i][yas + 1] = v1[i][yas] * 2 - v1[i][yas - 1];
	}

	v1[0][0] = v1[1][1] * 2 - v1[2][2];
	v1[xas + 1][0] = v1[xas][1] * 2 - v1[xas - 1][2];
	v1[0][yas + 1] = v1[1][yas] * 2 - v1[2][yas - 1];
	v1[xas + 1][yas + 1] = v1[xas][yas] * 2 - v1[xas - 1][yas - 1];

	//operation of data

	//calculating the cell size - x, y, and diagonal.
	ux = xs / xas;
	uy = ys / yas;
	uxy = sqrt(ux * ux + uy * uy);

	for (j=1; j<=yas; j++) {
		for (i=1; i<=xas; i++) {
			ii = i;
			jj = j;
			d = 0;
			//searching the flow route. if the route reaches the edge of the region, jump out of the loop.
			while (((ii >= 1) && (ii <= xas)) && ((jj >= 1) && (jj <= yas))) {
				//determining the steepest downhill direction. if cannot be determined, let d remain zero.
				d = 0;
				gr = 0;
				grmin = 0;
				gr = (v1[ii][jj-1] - v1[ii][jj]) / uy;
				if (gr < 0) {
					d = 1;	//north
					grmin = gr;
				}
				gr = (v1[ii+1][jj-1] - v1[ii][jj]) / uxy;
				if (gr < grmin) {
					grmin = gr;
					d = 2;	//northeast
				}
				gr = (v1[ii+1][jj] - v1[ii][jj]) / ux;
				if (gr < grmin) {
					grmin = gr;
					d = 3;	//east
				}
				gr = (v1[ii+1][jj+1] - v1[ii][jj]) / uxy;
				if (gr < grmin) {
					grmin = gr;
					d = 4;	//southeast
				}
				gr = (v1[ii][jj+1] - v1[ii][jj]) / uy;
				if (gr < grmin) {
					grmin = gr;
					d = 5;	//south
				}

				gr = (v1[ii-1][jj+1] - v1[ii][jj]) / uxy;
				if (gr < grmin) {
					grmin = gr;
					d = 6;	//southwest
				}
				gr = (v1[ii-1][jj] - v1[ii][jj]) / ux;
				if (gr < grmin) {
					grmin = gr;
					d = 7;	//west
				}
				gr = (v1[ii-1][jj-1] - v1[ii][jj]) / uxy;
				if (gr < grmin) {
					grmin = gr;
					d = 8;	//northwest
				}

				//move to the determined direction.
				if (d == 1) {	//to the north
					ii = ii;
					jj = jj - 1;
				} else if (d == 2) {	//to the northeast
					ii = ii + 1;
					jj = jj - 1;
				} else if (d == 3) {	//to the east
					ii = ii + 1;
					jj = jj;
				} else if (d == 4) {	//to the southeast
					ii = ii + 1;
					jj = jj + 1;
				} else if (d == 5) {	//to the south
					ii = ii;
					jj = jj + 1;
				} else if (d == 6) {	//to the southwest
					ii = ii - 1;
					jj = jj + 1;
				} else if (d == 7) {	//to the west
					ii = ii - 1;
					jj = jj;
				} else if (d == 8) {	//to the northwest
					ii = ii - 1;
					jj = jj - 1;
				} else break;	//if cannot be determined, jump out of the loop.

				v2[ii][jj]++;	//add up the times to be flow routes.
			}

		}
	}

	//change the variable and unit names.
	strcpy(vn, "catchment");
	strcpy(un, "cells");

	//output of the result
	//when you want to save them in a file, please redirect it.
	cout << rn << "\n" << vn << "\n" << un << "\n" ;
	cout << dt << "\n" << xas << "\n" << yas << "\n";
	cout << xs << "\n" << ys << "\n" << str << "\n";
	for (j=1; j<=yas; j++) {
		for (i=1; i<=xas; i++) {
			cout << v2[i][j] << "\n";
		}
	}
	return 0;
}	//end of the program