matrix.cpp

#include "matrix.h"

MatrixXd create_first_boundary_matrix_v1(int max_mode)
{
	// function creates the first boundary condition matrix U

	MatrixXd U = MatrixXd::Zero(max_mode, 3 * max_mode);

	for (int i = 0; i < max_mode; ++i)
	{
		for (int j = 0; j < 3 * max_mode; ++j)
		{
			if (j % max_mode == i)
				U(i, j) = 1.0;
		}
	}

	return U;
}

MatrixXd create_matrix_V(int max_mode)
{
	MatrixXd V(max_mode, max_mode);

	for (int p = 0; p < max_mode; ++p)
	{
		for (int q = 0; q < max_mode; ++q)
			V(p, q) = std::sqrt((2*p + 1) * (2*q + 1)) / 2;
	}

	return V;
}

vector<MatrixXd> create_matrix_block_SB_SL_SP_torque_v2(vector2d &L_1_0, valarray<double> &coeff_r_total_B,
	vector2d &L_1_a_B, vector2d &L_1_b_B_f, valarray<double> &coeff_r_total_L, vector2d &L_1_a_L, vector2d &L_1_b_L_f, 
	valarray<double> &coeff_r_total_P, vector2d &L_1_a_P, vector2d &L_1_b_P_f, vector<double> &w3j_all, int max_mode)
{
	MatrixXd S_B = MatrixXd::Zero(max_mode, max_mode);
	MatrixXd S_L = MatrixXd::Zero(max_mode, max_mode);
	MatrixXd S_P = MatrixXd::Zero(max_mode, max_mode);

	for (int p1 = 0; p1 < max_mode; ++p1)
	{
		for (int p2 = 0; p2 < max_mode; ++p2)
		{
			double prefactor_1 = std::pow(pi, 2) * std::sqrt((2 * p1 + 1)*(2 * p2 + 1));

			for (int n = 0; n < max_mode; ++n)
			{
				double prefactor_2 = prefactor_1 * (2 * n + 1);

				for (int k1 = std::abs(n - p1); k1 <= n + p1; ++k1)
				{
					double prefactor_3 = prefactor_2 * (2 * k1 + 1);

					for (int k2 = std::abs(n - p2); k2 <= n + p2; ++k2)
					{
						for (int r = 0; r <= std::min({ k1, k2, n }); ++r)
						{
							double common_part = prefactor_3 * (2 * k2 + 1) * L_1_0[k2][max_mode + r - 1] * 
								w3j_all[index(n, k1, p1, max_mode + r - 1, max_mode - 1)] * 
								w3j_all[index(n, k2, p2, max_mode + r - 1, max_mode - 1)];

							double temp = common_part * coeff_r_total_B[r] *
								L_1_b_B_f[k1][max_mode + r - 1] * L_1_a_B[n][max_mode + r - 1];
							
							S_B(p1, p2) = S_B(p1, p2) + common_part * coeff_r_total_B[r] *
								L_1_b_B_f[k1][max_mode + r - 1] * L_1_a_B[n][max_mode + r - 1];
							S_L(p1, p2) = S_L(p1, p2) + common_part * coeff_r_total_L[r] * 
								L_1_b_L_f[k1][max_mode + r - 1] * L_1_a_L[n][max_mode + r - 1];
							S_P(p1, p2) = S_P(p1, p2) + common_part * coeff_r_total_P[r] * 
								L_1_b_P_f[k1][max_mode + r - 1] * L_1_a_P[n][max_mode + r - 1];
						}
					}
				}
			}
		}
	}

	return vector<MatrixXd> {S_B, S_L, S_P};
}

MatrixXd create_transfer_matrix_torque_v2(Params &p, double mu_L, double mu_P, int max_mode, const vector<MatrixXd> &S_BLP)
{
	// function creates the main transfer matrix

	double exp_B = std::exp(-p.a_B - p.c_B);
	double exp_L = std::exp(-p.q * (mu_L - 2 * pi*p.tau*p.lk_L_0) - p.a_L - p.c_L);
	double exp_P = exp(-p.q * (mu_P - 2 * pi*p.tau*p.lk_P_0) - p.a_P - p.c_P);

	MatrixXd S(3 * max_mode, 3 * max_mode);

	MatrixXd S_11 = exp_B * S_BLP[0];
	MatrixXd S_12 = exp_B * S_BLP[0] * exp(-p.coop);
	MatrixXd S_13 = exp_B * S_BLP[0] * exp(-p.coop);

	MatrixXd S_21 = exp_L * S_BLP[1] * exp(-p.coop);
	MatrixXd S_22 = exp_L * S_BLP[1];
	MatrixXd S_23 = exp_L * S_BLP[1] * exp(-p.coop);

	MatrixXd S_31 = exp_P * S_BLP[2] * exp(-p.coop);
	MatrixXd S_32 = exp_P * S_BLP[2] * exp(-p.coop);
	MatrixXd S_33 = exp_P * S_BLP[2];

	S << S_11, S_12, S_13,
		S_21, S_22, S_23,
		S_31, S_32, S_33;

	return S;
}

MatrixXd quick_matrices_product_v3(MatrixXd &S, int N, double &log_coeff)
{
	string bin = std::bitset<32>(N).to_string();
	while (bin[0] == '0')
		bin.erase(0, 1);
	int bin_size = bin.size();

	int S_size = S.rows();
	MatrixXd A = MatrixXd::Identity(S_size, S_size);
	
	log_coeff = 0.0;
	double coeff_current = S.cwiseAbs().maxCoeff();
	MatrixXd A_current = S / coeff_current;
	double log_coeff_current = std::log(coeff_current);

	for (int i = 0; i < bin_size; ++i)
	{
		double coeff_A;

		switch (bin[bin_size - 1 - i])
		{
		case '1':
			A = A * A_current;
			coeff_A = A.cwiseAbs().maxCoeff();
			A = A / coeff_A;
			log_coeff = log_coeff + log_coeff_current + std::log(coeff_A);
			break;

		case '0':
			break;

		default:
			std::cerr << bin << " " << bin[bin_size - i] << "\n";
			std::cerr << "Matrices quick product error\n";
			exit(1);
		}

		if (i != bin_size - 1)
		{
			A_current = A_current * A_current;
			double coeff_current_new = A_current.cwiseAbs().maxCoeff();
			A_current = A_current / coeff_current_new;
			log_coeff_current = 2 * log_coeff_current + std::log(coeff_current_new);
		}
	}

	return A;
}

MatrixXd create_second_boundary_matrix_torque_v3(Params &p, double mu_L, double mu_P, double f, int max_mode, MatrixXd &V)
{
	// function creates the boundary condition block - matrix V_hat

	MatrixXd O(3 * max_mode, max_mode);

	MatrixXd O_1 = std::exp(p.b_B * f) * V;
	MatrixXd O_2 = std::exp(p.b_L * f - p.q * (mu_L - 2 * pi * p.tau * p.lk_L_0)) * V;
	MatrixXd O_3 = std::exp(p.b_P * f - p.q * (mu_P - 2 * pi * p.tau * p.lk_P_0)) * V;

	O << O_1, O_2, O_3;

	return O;
}