#include <iostream>

#include "include/praktable.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> &;

f64p (ex1) (std::string file) {
	table t(file);
	t	.add_columns({"sI", "sU", "sR", "sW"})
		.multiply_column("I", 1e-3)
		.apply(prak::mul<f64, 0.0025>, "U", "sU")
		.apply(prak::add<f64, 0.02>, "sU", "sU")
		.apply(prak::mul<f64, 0.01>, "I", "sI")
		.apply_function([](argvec a){ return a[0]/a[1]; }, {"U", "I"}, {"sU", "sI"}, "R", "sR")
		.apply_function(prak::prod<f64>, {"U", "I"}, {"sU", "sI"}, "W", "sW")
		.write_plot_4("R(W).plot", "W", "R", "sW", "sR")
		.print();
	auto [A, B] = t.least_squares_linear("W", "R", "sR", std::nullopt);
	std::cout << "МНК R(W): A = " << A << "; B = " << B << '\n';
	f64p alpha = {3.5e-3, 5e-5},
	     L = {0.25, 0.005},
	     r_2 = {0.006, 5e-6}, r_1 = {0.075e-3, 5e-6},
	     k = prak::function<f64>([](vecarg a){ return a[0] * a[1] / (2 * prak::PI * a[2] * a[3]) * std::log(a[4]/a[5]);}, {alpha, B, L, A, r_2, r_1});
	std::cout var(alpha) var (L) var(r_1) var(r_2) << "k = " << k << '\n';
	f64p eps = {0.35, 0.005}, sgm = {5.67e-8, 0.01e-8}, 
	     S = prak::function<f64>([](vecarg a){ return 2 * prak::PI * a[0] * a[1]; }, {r_1, L}),
	     R = {t.col_avg("R"), 0}, dT = prak::function<f64>([](vecarg a){return std::abs((a[0] - a[1])/(a[1]*a[2]));}, {R, B, alpha}),
	     T_0 = {295, 5}, T_1 = T_0 + dT, W_rad = eps * sgm * (T_1*T_1*T_1*T_1 - T_0*T_0*T_0*T_0)*S, 
	     g = {9.8, 0.05}, beta = 1.0 / T_0, nu = {1.6e-5, 0.1e-5}, Gr = g * beta * L * L * L * dT / nu / nu,
	     alpha2 = 0.75 * k / L * ((0.7 * Gr) ^ 0.25), W_con = alpha2 * dT * S;
	std::cout var(eps) var(sgm) var(R) var(dT) var(T_0) var(T_1) var(S) var(W_rad) var(beta) var(nu) var(Gr) var(alpha2) var(W_con);
	return k;
}

f64p (ex2) (const std::string& s) {
	table t(s);
	t	.multiply_column("I_V", 1e-6)
		.apply(prak::mul<f64, 1333.224>, "DP", "DP'")
		.write_plot("I_V(DP).plot", "DP'", "I_V", std::nullopt)
		.print();
	// gnuplot дал оценку лучше
	f64p A = {1.24834e-09, 5.18e-11},
	     B = {1.57878e-07, 2.787e-07};

	f64p R = {1e-3, 5e-5}, l = {0.23, 0.005}, eta = prak::PI * (R^4.0) / (8.0 * l * A),
	     Re = 2.0 * f64p{1.2041, 0.02} * t.col_max("I_V") / (eta * prak::PI * R);
	std::cout var(A) var(B) var(eta) var(Re) << std::endl;
	return eta;
}

int (main) (void) {
	std::cout << ex1("data") / ex2("data2") / 2.0  
		  << " = C_V\nпо табличным значениям: C_V = " << 0.0259 / 0.001004 / 2
		  << "\nПо здравому смыслу:5/2R = " << 5.0/2 * 8.31  << std::endl;

	return 0;
}