//program file name <stream2.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
	long	v2[1202][1202];		//raster data 2 for output
	int		v3[1202][1202];		//raster data 3 for the stream footprint

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, k;
	float	ux, uy, uxy;	//size of the cells - x, y, and diagonal
	float	gr[9];			//array for gradients
	int		grf[9];			//array for the flags of steepest downhill gradient
	int		fp[9];			//array for the flags of stream footprint
	float	grmin, grmin2;	//gradient and its minimum
	int		d, dd;		//drain direction

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

	//exit if the number of arguments is not 1.
	if (argc != 2) {
		cerr << "Usage: stream2 <filename>\n";
		cerr << "During the processing, the line under process will be shown on the monitor.\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 * 1000;
	uy = ys / yas * 1000;
	uxy = sqrt(ux * ux + uy * uy) * 1000;

	for (j=1; j<=yas; j++) {
		for (i=1; i<=xas; i++) {
			//initialize the array for footprint as zero.
			for (jj=1; jj<=yas; jj++) {
				for (ii=1; ii<=xas; ii++) {
					v3[ii][jj] = 0;
				}
			}

			ii = i;
			jj = j;
			dd = 1;
			//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))) {
				d = 0;
				grmin = 1000;		//minimum of the gradient
				grmin2 = 1000;	//minimum of the gradient in the case of avoiding footprint
				
				for (k=0; k<=8; k++) {	//set the array for footprint as zero.
					fp[k] = 0;
				}
				//calculating the gradient
				gr[1] = (v1[ii][jj-1] - v1[ii][jj]) / uy;	//north

				gr[2] = (v1[ii+1][jj-1] - v1[ii][jj]) / uxy;	//northeast

				gr[3] = (v1[ii+1][jj] - v1[ii][jj]) / ux;		//east

				gr[4] = (v1[ii+1][jj+1] - v1[ii][jj]) / uxy;	//southeast

				gr[5] = (v1[ii][jj+1] - v1[ii][jj]) / uy;		//south

				gr[6] = (v1[ii-1][jj+1] - v1[ii][jj]) / uxy;	//southwest

				gr[7] = (v1[ii-1][jj] - v1[ii][jj]) / ux;		//west

				gr[8] = (v1[ii-1][jj-1] - v1[ii][jj]) / uxy;	//northwest

				//set a flag in the direction where the flow route left footprints.
				if (v3[ii][jj-1] == 1) fp[1] = 1;	//north
				
				if (v3[ii+1][jj-1] == 1) fp[2] = 1;	//northeast
				
				if (v3[ii+1][jj] == 1) fp[3] = 1;	//east
				
				if (v3[ii+1][jj+1] == 1) fp[4] = 1;	//southeast
				
				if (v3[ii][jj+1] == 1) fp[5] = 1;	//south
				
				if (v3[ii-1][jj+1] == 1) fp[6] = 1;	//southwest
				
				if (v3[ii-1][jj] == 1) fp[7] = 1;	//west
				
				if (v3[ii-1][jj-1] == 1) fp[8] = 1;	//northwest
				

			//specify the minimum value of the gradients
				for (k=1; k<=8; k++) {
					if (gr[k] < grmin) grmin = gr[k];
					if ((gr[k] < grmin2) && (fp[k] != 1)) grmin2 = gr[k];	//avoid the footprints
				}
				
			//if the minimum gradient is grater than the prescribed value, jump out of the loop.
				if (grmin2 > 0.01) break;

			//if the difference of grmin and grmin2 is greater than the prescribed valude, jump out of the loop.
				if ((grmin2 - grmin) > 0.01) break;
				
				for (k=0; k<=8; k++) {	//set the array for flags of the steepest downhill direction as zero.
					grf[k] = 0;
				}

				for (k=1; k<=8; k++) {	//set a flag in the direction of the steepest downhill gradient.
					if (gr[k] == grmin2) grf[k] = 1;
				}

			//determining the flow direction
			//determine the flow direction considering the previous flow direction. choose the steepest downhill direction avoiding footprints.
				if ((grf[dd] == 1) && (fp[dd] != 1)) {
					d = dd;
					goto jump;
				}
				
				dd = dd + 1;
				if (dd > 8) dd = dd - 8;
				if (dd < 1) dd = dd + 8;
				if ((grf[dd] == 1) && (fp[dd] != 1)) {
					d = dd;
					goto jump;
				}
				
				dd = dd - 2;
				if (dd > 8) dd = dd - 8;
				if (dd < 1) dd = dd + 8;
				if ((grf[dd] == 1) && (fp[dd] != 1)) {
					d = dd;
					goto jump;
				}

				dd = dd + 3;
				if (dd > 8) dd = dd - 8;
				if (dd < 1) dd = dd + 8;
				if ((grf[dd] == 1) && (fp[dd] != 1)) {
					d = dd;
					goto jump;
				}
				
				dd = dd - 4;
				if (dd > 8) dd = dd - 8;
				if (dd < 1) dd = dd + 8;
				if ((grf[dd] == 1) && (fp[dd] != 1)) {
					d = dd;
					goto jump;
				}
				
				dd = dd + 5;
				if (dd > 8) dd = dd - 8;
				if (dd < 1) dd = dd + 8;
				if ((grf[dd] == 1) && (fp[dd] != 1)) {
					d = dd;
					goto jump;
				}
				
				dd = dd - 6;
				if (dd > 8) dd = dd - 8;
				if (dd < 1) dd = dd + 8;
				if ((grf[dd] == 1) && (fp[dd] != 1)) {
					d = dd;
					goto jump;
				}

				break;	//if the flow direction cannot be determined, jump out of the loop.

				jump:		//jump here.

				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;

				v3[ii][jj] = 1;		//record the footprint.

				dd = d;		//record the present flow direction.
			}


			//add up the times to be the stream routes (footprints).
			for (jj=1; jj<=yas; jj++) {
				for (ii=1; ii<=xas; ii++) {
					v2[ii][jj] += v3[ii][jj];
				}
			}

		}
		cerr << "line: " << j << '\n';	//show the status of processing on the monitor.
	}

	//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