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

using namespace std;

//３次元ベクトルの構造体
struct vect3d {
        float   x, y, z;
};

//３頂点がなす平面の法線ベクトルを算出する関数
vect3d normv(vect3d tri[3]) {

        float   vlen;   //ベクトルの長さ
        vect3d   nv;    //法線ベクトル

        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;
}


//大きな配列はグローバルとして宣言する。
	float	v1[1202][1202];		//ラスターデータ１
	float	v2[1202][1202];		//ラスターデータ２

int main(int argc, char *argv[]) {
	int	i, j;
	char	rn[128];		//地域名
	char	vn[128];		//変数名
	char	un[128];		//単位名
	int	dt;		//データのタイプ（未使用）
	int	xas, yas;		//要素数
	float	xs, ys;		//モデルの大きさ
	char	str[128];		//ヘッダーとデータを区切る文字列

        //以下は傾斜計算のために定義する変数
        vect3d  slopev;		//斜面の法線ベクトル
        vect3d  tri[3];		//三角形（３つの頂点）
        vect3d  nv[5];		//法線ベクトルを５つ用意する；この行に誤りがあり修正されました（2010/06/08）。
        float   vlen;		//ベクトルの長さ
        float   xus, yus;		//格子単位の大きさ



	//２つのラスターデータは、すべてゼロで初期化しておく。
	for (j=0; j<=1201; j++) {
		for (i=0; i<=1201; i++) {
			v1[i][j] = 0;
			v2[i][j] = 0;
		}
	}

	//コマンドライン引数が２つでない場合は、使用法を表示して終了する。
	if (argc != 3) {
		cerr << "Usage: slope <'-a' or '-g'> <filename>\n";
		return 1;
	}

	// コマンドラインの第１引数が '-a' でも '-g' でもない場合は，使用法を表示して終了する。
	if ((strcmp(argv[1], "-a") != 0) && (strcmp(argv[1], "-g") != 0)) {
		cerr << "Usage: slope <'-a' or '-g'> <filename>\n";
		return 1;
	}

	//ファイルが開けなかった場合は、その旨表示して終了する。
	ifstream fin(argv[2]);
	if (!fin) {
		cerr << "Cannot open file.\n";
		return 1;
	}

	//データの読み込み。
	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();

        //データの外挿
        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];

        //格子単位の計算
        
        xus = xs / xas * 1000;
        yus = ys / yas * 1000;

	//データの処理
	for (j=1; j<=yas; j++) {
		for (i=1; i<=xas; i++) {
			//斜面の法線ベクトルの計算
			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) {		//斜面方位の計算
				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 {					//傾斜の計算
				v2[i][j] = acos(- slopev.z) * 180 / 3.1416;
			}
		}
	}

	//変数名に処理内容を追加する。
	if (strcmp(argv[1], "-a") == 0) {		//斜面方位の場合
			strcat(vn, "−斜面方位");
			strcpy(un, "度");
	} else {					//傾斜の場合
			strcat(vn, "−傾斜");
			strcpy(un, "度");
	}

	//結果の出力。
	//標準出力（画面）に出力するので、ファイルに保存する場合は、リダイレクトする。
	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;
}