#include <iostream>

#include "include/praktable.hpp"
#include "include/prakphys.hpp"

using table = prak::table<f64>;
using f64p = prak::pvalue<f64>;
using f64v = std::vector<f64>;
using vecarg = const std::vector<f64> &;
using argvec = const std::vector<f64> &;

template <typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& vec) {
	os << "{ "; 
	for (const auto& val: vec) os << val << ", ";
	os << "} ";
	return os;
}

auto ex1(std::string file) {
	table t(file);
	t.column_width = 12;
	std::cout << t
		.apply(prak::diff<f64>, {"h1", "h2"}, "p")
		.fill_column("sp", 0.1 * std::sqrt(2.0))
		.apply(prak::log<f64>, {"p"}, "ln(p)")
		.add_column("s_ln(p)")
		.apply(prak::quot<f64>, {"sp", "p"}, "s_ln(p)")
		.apply(prak::add<f64, conv::C_to_K<f64>>, {"T"}, "T")
		.apply(prak::inv<f64>, {"T"}, "1/T")
		.apply(prak::avg<f64>, {"t1", "t2", "t3"}, "t")
		.apply(prak::stddev<f64>, {"t1", "t2", "t3"}, "st")
		.delete_cols({"t1", "t2", "t3"})
		.write_plot("plot_lnp.plot", "1/T", "ln(p)", "s_ln(p)");
	auto [A, B] = t.least_squares_linear("1/T", "ln(p)", "s_ln(p)", std::nullopt);
	std::cout << "Мнк ln(p) от (1/T): " << A << ' ' << B << std::endl;
	f64 a_tbl = .5451, mu = .09209, v_w = .00100177, v_s  = 57.8;
	f64p Q = A * (-prak::R<f64>);
	f64p q = Q/mu, p = 1333.2239023154 * (f64p){t["p", 4], t["sp", 4]};
	f64p a = mu*mu*(q - p*v_s)*v_w;
	std::cout << "\nmu = " << mu << ", v_воды = " << v_w << ", v_пар = " << v_s << ", p = " << p << " (темература 20 по методу наглых подгонов)"
		  << "\nQ = " << Q << ", q = " << q << ", a = " << a << ", a_теор = " << a_tbl << ", отклонение = " << a_tbl - a << std::endl;
	return t;
}

void ex2(prak::table<f64> t) {
	f64 beta = 6.1e-4, rho_0 = 1263.6, rho_1 = 8131, T_0 = 20 + conv::C_to_K<f64>, K = 114.4e-3, mu = .09209;
	t	.add_column("rho_2")
		.apply([=](vecarg a){ return rho_0 / (1 + beta * (a[0] - T_0)); }, {"T"}, "rho_2")
		.add_column("s_rho_2", 0.000001)
		.apply_function_n( [=](vecarg a){ return K * (rho_1 - a[0]) * a[1]; }, 
			{t.begin("rho_2"), t.begin("t")}, 
			{t.begin("s_rho_2"), t.begin("st")}, 
			t.rows, "nu", "snu")
		.add_column("ln(nu)", NAN)
		.add_column("s_ln(nu)", NAN)
		.apply_function_n(prak::log<f64>, {t.begin("nu")}, {t.begin("snu")}, t.rows, "ln(nu)", "s_ln(nu)")
		.add_column("1/nu", NAN)
		.add_column("s_1/nu", NAN)
		.add_column("1/rho_2", NAN)
		.apply(prak::inv<f64>, {"nu"}, "1/nu")
		.apply([](vecarg a){ return a[1] / a[0] / a[0]; }, {"nu", "snu"}, "s_1/nu")
		.apply(prak::inv<f64>, {"rho_2"}, "1/rho_2")
		.delete_cols({"h1", "h2", "p", "sp", "ln(p)", "t", "st"})
		.write_plot("plot_lnn.plot", "1/T", "ln(nu)", "s_ln(nu)")
		.write_plot("plot_nut.plot", "1/rho_2", "1/nu", "s_1/nu")
		.print();
	auto [A, B] = t.least_squares_linear("1/T", "ln(nu)", "s_ln(nu)", std::nullopt);
	f64p eps = prak::k<f64> * A;
	std::cout
		<< "Мнк ln(p) (1/T): (" << A << ")x + (" << B << ")\n"
		<< "энергия активации: " << eps * conv::J_to_eV<f64> << "эВ" << std::endl;
	auto [A1, B1] = t.least_squares_linear("1/rho_2", "1/nu", "s_1/nu", std::nullopt);
	f64p c = A / mu, b = -mu * B / A;
	std::cout 
		<< "Мнк 1/nu (1/rho_2): (" << A1 << ")x + (" << B1 << ")\n" 
		<< "b = " << b << "; c = " << c << std::endl;
}

int main() {
	ex2(ex1("data1"));
	return 0;
}