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

using namespace std;

//structure for 3 diemnsional vector.
struct vect3d {
        float   x, y, z;
};

//function calculating the normal vector from a triangle.
vect3d normv(vect3d tri[3]) {

        float   vlen;	//length of a vector
        vect3d   nv;	//normal vector

        nv.x = (tri[0].y - tri[1].y) * (tri[0].z + tri[1].z) + (tri[1].y - tri[2].y) * (tri[1].z + tri[2].z) + (tri[2].y - tri[0].y) * (tri[2].z + tri[0].z);
        nv.y = (tri[0].z - tri[1].z) * (tri[0].x + tri[1].x) + (tri[1].z - tri[2].z) * (tri[1].x + tri[2].x) + (tri[2].z - tri[0].z) * (tri[2].x + tri[0].x);
        nv.z = (tri[0].x - tri[1].x) * (tri[0].y + tri[1].y) + (tri[1].x - tri[2].x) * (tri[1].y + tri[2].y) + (tri[2].x - tri[0].x) * (tri[2].y + tri[0].y);
        vlen = sqrt(nv.x * nv.x + nv.y * nv.y + nv.z * nv.z);
        nv.x =  nv.x / vlen;
        nv.y =  nv.y / vlen;
        nv.z =  nv.z / vlen;

        return nv;
}


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

        //the following variables are for calculating slope parameters.
        vect3d  slopev;		//normal vector at the sample point on a DEM
        vect3d  tri[3];		//triangle
        vect3d  nv[5];		//normal vectors; This line was corrected on June 8, 2010.
        float   vlen;		//length of a vector
        float   xus, yus;		//unit sizes of 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: slope <'-a' or '-g'> <filename>\n";
		return 1;
	}

	//if 1st argument is neither '-a' nor '-g', show an error message and exit.
	if ((strcmp(argv[1], "-a") != 0) && (strcmp(argv[1], "-g") != 0)) {
		cerr << "Usage: slope <'-a' or '-g'> <filename>\n";
		return 1;
	}

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

	//unit sizes of raster grid.
	xus = xs / xas * 1000;
	yus = ys / yas * 1000;

	//manipulation of data. for calculation of slope parameters here.
	for (j=1; j<=yas; j++) {
		for (i=1; i<=xas; i++) {
			//calculate normal vector for the slope at a grid cell.
			tri[0].x = 0;
			tri[0].y = 0;
			tri[0].z = v1[i][j];

			tri[1].x = 0;
			tri[1].y = -yus;
			tri[1].z = v1[i][j - 1];
			tri[2].x = -xus;
			tri[2].y = 0;
			tri[2].z = v1[i - 1][j];
			nv[1] = normv(tri);

			tri[1].x = xus;
			tri[1].y = 0;
			tri[1].z = v1[i + 1][j];
			tri[2].x = 0;
			tri[2].y = -yus;
			tri[2].z = v1[i][j - 1];
			nv[2] = normv(tri);

			tri[1].x = 0;
			tri[1].y = yus;
			tri[1].z = v1[i][j + 1];
			tri[2].x = xus;
			tri[2].y = 0;
			tri[2].z = v1[i + 1][j];
			nv[3] = normv(tri);

			tri[1].x = -xus;
			tri[1].y = 0;
			tri[1].z = v1[i - 1][j];
			tri[2].x = 0;
			tri[2].y = yus;
			tri[2].z = v1[i][j + 1];
			nv[4] = normv(tri);

			nv[0].x = nv[1].x + nv[2].x + nv[3].x + nv[4].x;
			nv[0].y = nv[1].y + nv[2].y + nv[3].y + nv[4].y;
			nv[0].z = nv[1].z + nv[2].z + nv[3].z + nv[4].z;

			vlen = sqrt(nv[0].x * nv[0].x + nv[0].y * nv[0].y + nv[0].z * nv[0].z);

			if (vlen == 0) {
				slopev.x = 0;
				slopev.y = 0;
				slopev.z = -1;
			} else {
				slopev.x = nv[0].x / vlen;
				slopev.y = nv[0].y / vlen;
				slopev.z = nv[0].z / vlen;
			}

			if (strcmp(argv[1], "-a") == 0) {		//calculate the aspect angle.
				if (slopev.x == 0) {
					if (slopev.y >= 0) {
						v2[i][j] = 0;
					} else {
						v2[i][j] = 180;
					}
				} else if (slopev.x > 0) {
					v2[i][j] = atan(slopev.y / slopev.x) * 180 / 3.1416 + 270;
				} else {
					v2[i][j] = atan(slopev.y / slopev.x) * 180 / 3.1416 + 90;
				}
			} else {					//calculate the gradient.
				v2[i][j] = acos(- slopev.z) * 180 / 3.1416;
			}
		}
	}

	//add manipulation to the variable name.
	if (strcmp(argv[1], "-a") == 0) {		//in the case of slope angle
			strcat(vn, "-aspect");
			strcpy(un, "degree");
	} else {					//in the case of gradient
			strcat(vn, "-gradient");
			strcpy(un, "degree");
	}

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