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

using namespace std;

//大きな配列はグローバルとして宣言する。
	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];		//ヘッダーとデータを区切る文字列

	int		opt;		//操作オプション
	float	a, b, c, d, e, f, g, h, k;	//２次式のパラメータ
	float	xus;		//格子単位の大きさ（ｘ軸方向）

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

	//コマンドラインの第1引数が正しくない場合は、使用法を表示して終了する。
	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;
	}

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

	//ファイルが開けなかった場合は、その旨表示して終了する。
	ifstream fin(argv[2]);
	if (!fin) {
		cerr << "Cannot open data 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();

	xus = xs / xas * 1000;	//格子単位の大きさ（ｘ軸方向）の計算

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

	//データの処理
	for (j=1; j<=yas; j++) {
		for (i=1; i<=xas; i++) {
			//２次多項式のパラメータを計算する。
			//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];

			//傾斜を算出する。
			if (opt == 1) {
				v2[i][j] = sqrt(g*g + h*h);
			}

			//斜面方位を算出する。
			else if (opt == 2) {
				if (g != 0) {
					v2[i][j] = atan(-h / -g);
				} else {
					v2[i][j] = 0;
				}
			}
			
			//断面曲率を算出する。
			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;
				}
			}

			//平面曲率を算出する。
			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;
			}
		}
	}

	//変数名に処理内容を適宜追加する。
	if (opt == 1) strcat(vn, "・傾斜");
	else if (opt == 2) strcat(vn, "・斜面方位");
	else if (opt == 3) strcat(vn, "・断面曲率");
	else if (opt == 4) strcat(vn, "・平面曲率");
	else strcpy (vn, 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;
}	//プログラムの終わり。

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