//program file name <quad.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
	float	v2[1202][1202];		//raster data 2

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		opt;		//option of the operation
	float	a, b, c, d, e, f, g, h, k;	//coefficients of the quadratic polynomial
	float	xus;		//unit size of the x direction raster grid

	//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 2.
	if (argc != 3) {
		cerr << "Usage: quad <option> <filename>\n";
		cerr << '\t' << "Options are as follows.\n";
		cerr << '\t' << "-s  : slope angle\n";
		cerr << '\t' << "-a  : aspect angle\n";
		cerr << '\t' << "-pr : profile curvature\n";
		cerr << '\t' << "-pl : plan curvature\n";
		return 1;
	}

	//according to the 1st argument, specify the option of the operation.
	if (strcmp(argv[1], "-s") == 0) opt = 1;
	else if (strcmp(argv[1], "-a") == 0) opt = 2;
	else if (strcmp(argv[1], "-pr") == 0) opt = 3;
	else if (strcmp(argv[1], "-pl") == 0) opt = 4;
	else {
		//if the 1st argument is not correct, show an error message and exit.
		cerr << "Usage: quad <option> <filename>\n";
		cerr << '\t' << "Options are as follows.\n";
		cerr << '\t' << "-s  : slope angle\n";
		cerr << '\t' << "-a  : aspect angle\n";
		cerr << '\t' << "-pr : profile curvature\n";
		cerr << '\t' << "-pl : plan curvature\n";
		return 1;
	}

	//open the input file. exit if an error occurs.
	ifstream fin(argv[2]);
	if (!fin) {
		cerr << "Cannot open data 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();

	xus = xs / xas * 1000;	//unit size of the x direction raster grid.

	//make 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];

	//manipulation of data.
	for (j=1; j<=yas; j++) {
		for (i=1; i<=xas; i++) {
			//calculate the coefficients of the quadratic polynomial
			//a = ((v1[i-1][j-1]+v1[i+1][j-1]+v1[i-1][j+1]+v1[i+1][j+1])/4 - (v1[i][j-1] + v1[i-1][j] + v1[i+1][j] + v1[i][j+1])/2 + v1[i][j]) / pow(xus, 4);
			//b = ((v1[i-1][j-1]+v1[i+1][j-1]-v1[i-1][j+1]-v1[i+1][j+1])/4 - (v1[i][j-1] - v1[i][j+1])/2) / pow(xus, 3);
			//c = ((-v1[i-1][j-1]+v1[i+1][j-1]-v1[i-1][j+1]+v1[i+1][j+1])/4 + (v1[i-1][j] - v1[i+1][j])/2) / pow(xus, 3);
			d = ((v1[i-1][j] + v1[i+1][j])/2 - v1[i][j]) / pow(xus, 2);
			e = ((v1[i][j-1] + v1[i][j+1])/2 - v1[i][j]) / pow(xus, 2);
			f = (-v1[i-1][j-1]+v1[i+1][j-1]+v1[i-1][j+1]-v1[i+1][j+1]) / (4*xus*xus);
			g = (-v1[i-1][j]+v1[i+1][j]) / (2*xus);
			h = (v1[i][j-1] - v1[i][j+1]) / (2*xus);
			//k = v1[i][j];

			//calculate the slope angle.
			if (opt == 1) {
				v2[i][j] = sqrt(g*g + h*h);
			}

			//calculate the aspect angle.
			else if (opt == 2) {
				if (g != 0) {
					v2[i][j] = atan(-h / -g);
				} else {
					v2[i][j] = 0;
				}
			}
			
			//calculate the profile curvature.
			else if (opt == 3) {
				if ((g*g + h*h) != 0) {
					v2[i][j] = 2 * (d*g*g + e*h*h + f*g*h) / (g*g + h*h);
				} else {
					v2[i][j] = 0;
				}
			}

			//calculate the plan curvature.
			else if (opt == 4) {
				if ((g*g + h*h) != 0) {
					v2[i][j] = -2 * (d*h*h + e*g*g - f*g*h) / (g*g + h*h);
				} else {
					v2[i][j] = 0;
				}
			} else {
				v2[i][j] = 0;
			}
		}
	}

	//add manipulation to the variable name.
	if (opt == 1) strcat(vn, "-slope angle");
	else if (opt == 2) strcat(vn, "-aspect angle");
	else if (opt == 3) strcat(vn, "-profile curvature");
	else if (opt == 4) strcat(vn, "-plan curvature");
	else strcpy (vn, vn);

	//make name of the unit null string.
	strcpy(un, "");

	//output the results of calculation.
	//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

//---------------------------------------------------------------------------