/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include "surfer_io.h"

void gctl::read_surfer6_grid(std::string filename, _2d_matrix &out_data, 
	int &xnum, int &ynum, double &xmin, double &xmax, double &ymin, double &ymax, 
	double &zmin, double &zmax, surfer_file_type_e file_format)
{
	std::ifstream infile;

	if (file_format == Surfer6Binary)
	{
		open_infile(infile, filename, ".grd", std::ios::in|std::ios::binary);

		char sym_id[5] = "NULL";
		for (int i = 0; i < 4; i++)
		{
			infile.read((char*)&sym_id[i], sizeof(char));
		}

		if (strcmp(sym_id, "DSBB"))
		{
			infile.close();
			throw runtime_error("Wrong file format. From read_surfer6_grid(...)");
		}

		float tmp_d;
		short int M, N;
		infile.read((char*)&N, sizeof(short int));
		infile.read((char*)&M, sizeof(short int));
		xnum = N;
		ynum = M;

		// resize space for output array
		out_data.resize(M, N, 0.0);

		infile.read((char*)&xmin, sizeof(double));
		infile.read((char*)&xmax, sizeof(double));
		infile.read((char*)&ymin, sizeof(double));
		infile.read((char*)&ymax, sizeof(double));
		infile.read((char*)&zmin, sizeof(double));
		infile.read((char*)&zmax, sizeof(double));
		for (int i = 0; i < M; i++)
		{
			for (int j = 0; j < N; j++)
			{
				infile.read((char*)&tmp_d, sizeof(float));
				out_data[i][j] = (double)tmp_d;
			}
		}
		infile.close();
	}
	else
	{
		open_infile(infile, filename, ".grd");

		std::string tmp_str;
		// read first line
		infile >> tmp_str;
		if (tmp_str != "DSAA")
		{
			infile.close();
			throw runtime_error("Wrong file format. From read_surfer6_grid(...)");
		}

		int M, N; // N for x number, M for y number
		// read dimensions
		infile >> N >> M;
		xnum = N;
		ynum = M;
		// resize space for output array
		out_data.resize(M, N, 0.0);

		// read grid's dimension
		infile >> xmin >> xmax;
		infile >> ymin >> ymax;
		infile >> zmin >> zmax;
		// read grid's data
		for (int i = 0; i < M; i++)
		{
			for (int j = 0; j < N; j++)
			{
				infile >> tmp_str;
				out_data[i][j] = str2double(tmp_str);
			}
		}
		infile.close();
	}
	return;
}

void gctl::read_surfer6_grid(std::string filename, array<double> &out_data, int &xnum, 
	int &ynum, double &xmin, double &xmax, double &ymin, double &ymax, double &zmin, 
	double &zmax, surfer_file_type_e file_format)
{
	std::ifstream infile;

	if (file_format == Surfer6Binary)
	{
		open_infile(infile, filename, ".grd", std::ios::in|std::ios::binary);

		char sym_id[5] = "NULL";
		for (int i = 0; i < 4; i++)
		{
			infile.read((char*)&sym_id[i], sizeof(char));
		}

		if (strcmp(sym_id, "DSBB"))
		{
			infile.close();
			throw runtime_error("Wrong file format. From read_surfer6_grid(...)");
		}

		float tmp_d;
		short int M, N;
		infile.read((char*)&N, sizeof(short int));
		infile.read((char*)&M, sizeof(short int));
		xnum = N;
		ynum = M;

		// resize space for output array
		out_data.resize(M*N, 0.0);

		infile.read((char*)&xmin, sizeof(double));
		infile.read((char*)&xmax, sizeof(double));
		infile.read((char*)&ymin, sizeof(double));
		infile.read((char*)&ymax, sizeof(double));
		infile.read((char*)&zmin, sizeof(double));
		infile.read((char*)&zmax, sizeof(double));
		for (int i = 0; i < M; i++)
		{
			for (int j = 0; j < N; j++)
			{
				infile.read((char*)&tmp_d, sizeof(float));
				out_data[N*i + j] = (double)tmp_d;
			}
		}
		infile.close();
	}
	else
	{
		open_infile(infile, filename, ".grd");

		std::string tmp_str;
		// read first line
		infile >> tmp_str;
		if (tmp_str != "DSAA")
		{
			infile.close();
			throw runtime_error("Wrong file format. From read_surfer6_grid(...)");
		}

		int M, N; // N for x number, M for y number
		// read dimensions
		infile >> N >> M;
		xnum = N;
		ynum = M;
		// resize space for output array
		out_data.resize(M*N, 0.0);

		// read grid's dimension
		infile >> xmin >> xmax;
		infile >> ymin >> ymax;
		infile >> zmin >> zmax;
		// read grid's data
		for (int i = 0; i < M; i++)
		{
			for (int j = 0; j < N; j++)
			{
				infile >> tmp_str;
				out_data[N*i + j] = str2double(tmp_str);
			}
		}
		infile.close();
	}
	return;
}

void gctl::save_surfer6_grid(std::string filename, const _2d_matrix &out_data, 
	double xmin, double xmax, double ymin, double ymax, double zmin, double zmax, 
	surfer_file_type_e file_format)
{
	// 准备写出数据
	if (out_data.empty())
	{
		throw runtime_error("The output data is empty. From save_surfer6_grid(...)");
	}

	if (std::isnan(zmin))
	{
		zmin = GCTL_BDL_MAX;
		for (int i = 0; i < out_data.row_size(); i++)
		{
			for (int j = 0; j < out_data.col_size(); j++)
			{
				zmin = GCTL_MIN(zmin, out_data[i][j]);
			}
		}
	}

	if (std::isnan(zmax))
	{
		zmax = GCTL_BDL_MIN;
		for (int i = 0; i < out_data.row_size(); i++)
		{
			for (int j = 0; j < out_data.col_size(); j++)
			{
				zmax = GCTL_MAX(zmax, out_data[i][j]);
			}
		}
	}

	short int M = out_data.row_size(); // 占两个字节
	short int N = out_data.col_size();

	std::ofstream outfile;
	if (file_format == Surfer6Binary)
	{
		float tmp_d;
		char sym_id[] = "DSBB";

		open_outfile(outfile, filename, ".grd", std::ios::out|std::ios::binary);

		outfile.write((char*)&sym_id[0], sizeof(char)); // 注意不要直接输出sym_id 这里只能输出4个字节
		outfile.write((char*)&sym_id[1], sizeof(char));
		outfile.write((char*)&sym_id[2], sizeof(char));
		outfile.write((char*)&sym_id[3], sizeof(char));
		outfile.write((char*)&N, sizeof(short int));
		outfile.write((char*)&M, sizeof(short int));
		outfile.write((char*)&xmin, sizeof(double)); // 占8个字节
		outfile.write((char*)&xmax, sizeof(double));
		outfile.write((char*)&ymin, sizeof(double));
		outfile.write((char*)&ymax, sizeof(double));
		outfile.write((char*)&zmin, sizeof(double));
		outfile.write((char*)&zmax, sizeof(double));
		for (int i = 0; i < M; i++)
		{
			for (int j = 0; j < N; j++)
			{
				tmp_d = (float)out_data[i][j]; // 注意需要转换为float类型 占两个字节
				outfile.write((char*)&tmp_d, sizeof(float));
			}
		}
		outfile.close();
	}
	else
	{
		open_outfile(outfile, filename, ".grd");

		outfile << "DSAA" << std::endl;
		outfile << N << " " << M << std::endl;
		outfile << xmin << " " << xmax << std::endl;
		outfile << ymin << " " << ymax << std::endl;
		outfile << zmin << " " << zmax << std::endl;
		for (int i = 0; i < M; i++)
		{
			if (zmin <= out_data[i][0] && out_data[i][0] <= zmax)
				outfile << out_data[i][0];
			else outfile << "nan";

			for (int j = 1; j < N; j++)
			{
				if (zmin <= out_data[i][j] && out_data[i][j] <= zmax)
					outfile << " " << out_data[i][j];
				else outfile << " nan";
			}
			outfile << std::endl;
		}
		outfile.close();
	}
	return;
}

void gctl::save_surfer6_grid(std::string filename, const array<double> &out_data, 
	int xnum, int ynum, double xmin, double xmax, double ymin, double ymax, double zmin, 
	double zmax, surfer_file_type_e file_format)
{
	// 准备写出数据
	if (out_data.empty())
	{
		throw runtime_error("The output data is empty. From save_surfer6_grid(...)");
	}

	if (std::isnan(zmin))
	{
		zmin = GCTL_BDL_MAX;
		for (int i = 0; i < out_data.size(); i++)
		{
			zmin = GCTL_MIN(zmin, out_data[i]);
		}
	}

	if (std::isnan(zmax))
	{
		zmax = GCTL_BDL_MIN;
		for (int i = 0; i < out_data.size(); i++)
		{
			zmax = GCTL_MAX(zmax, out_data[i]);
		}
	}

	short int M = ynum; // 占两个字节
	short int N = xnum;

	std::ofstream outfile;
	if (file_format == Surfer6Binary)
	{
		float tmp_d;
		char sym_id[] = "DSBB";

		open_outfile(outfile, filename, ".grd", std::ios::out|std::ios::binary);

		outfile.write((char*)&sym_id[0], sizeof(char)); // 注意不要直接输出sym_id 这里只能输出4个字节
		outfile.write((char*)&sym_id[1], sizeof(char));
		outfile.write((char*)&sym_id[2], sizeof(char));
		outfile.write((char*)&sym_id[3], sizeof(char));
		outfile.write((char*)&N, sizeof(short int));
		outfile.write((char*)&M, sizeof(short int));
		outfile.write((char*)&xmin, sizeof(double)); // 占8个字节
		outfile.write((char*)&xmax, sizeof(double));
		outfile.write((char*)&ymin, sizeof(double));
		outfile.write((char*)&ymax, sizeof(double));
		outfile.write((char*)&zmin, sizeof(double));
		outfile.write((char*)&zmax, sizeof(double));
		for (int i = 0; i < M; i++)
		{
			for (int j = 0; j < N; j++)
			{
				tmp_d = (float)out_data[N*i+j]; // 注意需要转换为float类型 占两个字节
				outfile.write((char*)&tmp_d, sizeof(float));
			}
		}
		outfile.close();
	}
	else
	{
		open_outfile(outfile, filename, ".grd");

		outfile << "DSAA" << std::endl;
		outfile << N << " " << M << std::endl;
		outfile << xmin << " " << xmax << std::endl;
		outfile << ymin << " " << ymax << std::endl;
		outfile << zmin << " " << zmax << std::endl;
		for (int i = 0; i < M; i++)
		{
			if (zmin <= out_data[N*i] && out_data[N*i] <= zmax)
				outfile << out_data[N*i];
			else outfile << "nan";

			for (int j = 1; j < N; j++)
			{
				if (zmin <= out_data[N*i+j] && out_data[N*i+j] <= zmax)
					outfile << " " << out_data[N*i+j];
				else outfile << " nan";
			}
			outfile << std::endl;
		}
		outfile.close();
	}
	return;
}

void gctl::read_surfer7_grid(std::string filename, _2d_matrix &out_data, double &xmin, 
	double &ymin, double &dx, double &dy, double &zmin, double &zmax, double &blank_val)
{
	std::ifstream infile; 
	open_infile(infile, filename, ".grd", std::ios::in|std::ios::binary);

	int head_id, head_sec_size, head_version;
	infile.read((char*)&head_id, sizeof(int));
	infile.read((char*)&head_sec_size, sizeof(int));
	infile.read((char*)&head_version, sizeof(int));
	if (head_id != 0x42525344 || head_sec_size != 4)
	{
		infile.close();
		throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
	}

	int grid_id, grid_sec_size, M, N;
	double rotation;
	infile.read((char*)&grid_id, sizeof(int));
	infile.read((char*)&grid_sec_size, sizeof(int));
	infile.read((char*)&M, sizeof(int));
	infile.read((char*)&N, sizeof(int));
	if (grid_id != 0x44495247 || grid_sec_size != 72)
	{
		infile.close();
		throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
	}

	// resize space for output array
	out_data.resize(M, N, 0.0);

	infile.read((char*)&xmin, sizeof(double));
	infile.read((char*)&ymin, sizeof(double));
	infile.read((char*)&dx, sizeof(double));
	infile.read((char*)&dy, sizeof(double));
	infile.read((char*)&zmin, sizeof(double));
	infile.read((char*)&zmax, sizeof(double));
	infile.read((char*)&rotation, sizeof(double)); // 暂时无用
	infile.read((char*)&blank_val, sizeof(double));

	int data_id, data_sec_size;
	infile.read((char*)&data_id, sizeof(int));
	infile.read((char*)&data_sec_size, sizeof(int));
	if (data_id != 0x41544144 || data_sec_size != M*N*8)
	{
		infile.close();
		throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
	}

	for (int i = 0; i < M; i++)
	{
		for (int j = 0; j < N; j++)
		{
			infile.read((char*)&out_data[i][j], sizeof(double));
		}
	}
	infile.close();
	return;
}

void gctl::read_surfer7_grid(std::string filename, array<double> &out_data, int &xnum, 
	int &ynum, double &xmin, double &ymin, double &dx, double &dy, double &zmin, 
	double &zmax, double &blank_val)
{
	std::ifstream infile; 
	open_infile(infile, filename, ".grd", std::ios::in|std::ios::binary);

	int head_id, head_sec_size, head_version;
	infile.read((char*)&head_id, sizeof(int));
	infile.read((char*)&head_sec_size, sizeof(int));
	infile.read((char*)&head_version, sizeof(int));
	if (head_id != 0x42525344 || head_sec_size != 4)
	{
		infile.close();
		throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
	}

	int grid_id, grid_sec_size, M, N;
	double rotation;
	infile.read((char*)&grid_id, sizeof(int));
	infile.read((char*)&grid_sec_size, sizeof(int));
	infile.read((char*)&M, sizeof(int));
	infile.read((char*)&N, sizeof(int));
	if (grid_id != 0x44495247 || grid_sec_size != 72)
	{
		infile.close();
		throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
	}
	xnum = N;
	ynum = M;

	// resize space for output array
	out_data.resize(M*N, 0.0);

	infile.read((char*)&xmin, sizeof(double));
	infile.read((char*)&ymin, sizeof(double));
	infile.read((char*)&dx, sizeof(double));
	infile.read((char*)&dy, sizeof(double));
	infile.read((char*)&zmin, sizeof(double));
	infile.read((char*)&zmax, sizeof(double));
	infile.read((char*)&rotation, sizeof(double)); // 暂时无用
	infile.read((char*)&blank_val, sizeof(double));

	int data_id, data_sec_size;
	infile.read((char*)&data_id, sizeof(int));
	infile.read((char*)&data_sec_size, sizeof(int));
	if (data_id != 0x41544144 || data_sec_size != M*N*8)
	{
		infile.close();
		throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
	}

	for (int i = 0; i < M; i++)
	{
		for (int j = 0; j < N; j++)
		{
			infile.read((char*)&out_data[N*i+j], sizeof(double));
		}
	}
	infile.close();
	return;
}

void gctl::save_surfer7_grid(std::string filename, const _2d_matrix &out_data, double xmin, 
	double ymin, double dx, double dy, double zmin, double zmax, double blank_val)
{
	// 准备写出数据
	if (out_data.empty())
	{
		throw runtime_error("The output data is empty. From save_surfer7_grid(...)");
	}

	if (std::isnan(zmin))
	{
		zmin = GCTL_BDL_MAX;
		for (int i = 0; i < out_data.row_size(); i++)
		{
			for (int j = 0; j < out_data.col_size(); j++)
			{
				zmin = GCTL_MIN(zmin, out_data[i][j]);
			}
		}
	}

	if (std::isnan(zmax))
	{
		zmax = GCTL_BDL_MIN;
		for (int i = 0; i < out_data.row_size(); i++)
		{
			for (int j = 0; j < out_data.col_size(); j++)
			{
				zmax = GCTL_MAX(zmax, out_data[i][j]);
			}
		}
	}

	std::ofstream outfile;
	open_outfile(outfile, filename, ".grd", std::ios::out|std::ios::binary);

	int head_id = 0x42525344; // 固定值 查看官方说明文档
	int head_sec_size = sizeof(int);
	int head_version = 1;
	outfile.write((char*)&head_id, sizeof(int));
	outfile.write((char*)&head_sec_size, sizeof(int));
	outfile.write((char*)&head_version, sizeof(int));

	int grid_id = 0x44495247;
	int grid_sec_size = 8*sizeof(double) + 2*sizeof(int);
	int M = out_data.row_size();
	int N = out_data.col_size();
	double rotation = 0.0;
	outfile.write((char*)&grid_id, sizeof(int));
	outfile.write((char*)&grid_sec_size, sizeof(int));
	outfile.write((char*)&M, sizeof(int));
	outfile.write((char*)&N, sizeof(int));
	outfile.write((char*)&xmin, sizeof(double));
	outfile.write((char*)&ymin, sizeof(double));
	outfile.write((char*)&dx, sizeof(double));
	outfile.write((char*)&dy, sizeof(double));
	outfile.write((char*)&zmin, sizeof(double));
	outfile.write((char*)&zmax, sizeof(double));
	outfile.write((char*)&rotation, sizeof(double));
	outfile.write((char*)&blank_val, sizeof(double));

	int data_id = 0x41544144;
	int data_sec_size = M * N * 8;
	double tmp_d;
	outfile.write((char*)&data_id, sizeof(int));
	outfile.write((char*)&data_sec_size, sizeof(int));
	for (int i = 0; i < M; i++)
	{
		for (int j = 0; j < N; j++)
		{
			tmp_d = out_data[i][j];
			outfile.write((char*)&tmp_d, sizeof(double));
		}
	}
	outfile.close();
	return;
}

void gctl::save_surfer7_grid(std::string filename, const array<double> &out_data, 
	int xnum, int ynum, double xmin, double ymin, double dx, double dy, 
	double zmin, double zmax, double blank_val)
{
	// 准备写出数据
	if (out_data.empty())
	{
		throw runtime_error("The output data is empty. From save_surfer7_grid(...)");
	}

	if (std::isnan(zmin))
	{
		zmin = GCTL_BDL_MAX;
		for (int i = 0; i < out_data.size(); i++)
		{
			zmin = GCTL_MIN(zmin, out_data[i]);
		}
	}

	if (std::isnan(zmax))
	{
		zmax = GCTL_BDL_MIN;
		for (int i = 0; i < out_data.size(); i++)
		{
			zmax = GCTL_MAX(zmax, out_data[i]);
		}
	}

	std::ofstream outfile;
	open_outfile(outfile, filename, ".grd", std::ios::out|std::ios::binary);

	int head_id = 0x42525344; // 固定值 查看官方说明文档
	int head_sec_size = sizeof(int);
	int head_version = 1;
	outfile.write((char*)&head_id, sizeof(int));
	outfile.write((char*)&head_sec_size, sizeof(int));
	outfile.write((char*)&head_version, sizeof(int));

	int grid_id = 0x44495247;
	int grid_sec_size = 8*sizeof(double) + 2*sizeof(int);
	int M = ynum;
	int N = xnum;
	double rotation = 0.0;
	outfile.write((char*)&grid_id, sizeof(int));
	outfile.write((char*)&grid_sec_size, sizeof(int));
	outfile.write((char*)&M, sizeof(int));
	outfile.write((char*)&N, sizeof(int));
	outfile.write((char*)&xmin, sizeof(double));
	outfile.write((char*)&ymin, sizeof(double));
	outfile.write((char*)&dx, sizeof(double));
	outfile.write((char*)&dy, sizeof(double));
	outfile.write((char*)&zmin, sizeof(double));
	outfile.write((char*)&zmax, sizeof(double));
	outfile.write((char*)&rotation, sizeof(double));
	outfile.write((char*)&blank_val, sizeof(double));

	int data_id = 0x41544144;
	int data_sec_size = M * N * 8;
	double tmp_d;
	outfile.write((char*)&data_id, sizeof(int));
	outfile.write((char*)&data_sec_size, sizeof(int));
	for (int i = 0; i < M; i++)
	{
		for (int j = 0; j < N; j++)
		{
			tmp_d = out_data[N*i+j];
			outfile.write((char*)&tmp_d, sizeof(double));
		}
	}
	outfile.close();
	return;
}