#include "PortModeSolver.h"
#include "../mesh/Mesh_Base.h"
#include "../material/Material_Base.h"
#include "../function/BF.h"
#include "../function/Gauss.h"
#include "../common/define.h"
#include "../nlohmann/json.hpp"
#include "../Eigen/Sparse"
#include "../Eigen/SparseLU"
#include "../Eigen/Eigenvalues"

#include <algorithm>
#include <array>
#include <cmath>
#include <complex>
#include <fstream>
#include <iostream>
#include <map>
#include <set>
#include <vector>

using namespace Eigen;
using json = nlohmann::json;

namespace {

struct PortBoundaryMesh
{
	MatrixXd nodes;           // nNodes x 3
	MatrixXi newFaces;        // nf x 3 (0-based local node)
	MatrixXi facesConn;       // nf x 2 (1-based tet, 1-based local face)
	VectorXi volDomain;       // nf (1-based volume domain for material)
	MatrixXi edgeOfFace;      // nf x 3 (0-based edge DOF)
	MatrixXi edgeNodes;       // nEdges x 2 (0-based local nodes)
	int nNodes = 0;
	int nEdges = 0;
	int nFaces = 0;
};

bool domainInList(int domain, const VectorXi& list)
{
	for (int i = 0; i < list.size(); ++i)
		if (list(i) == domain)
			return true;
	return false;
}

bool extractPortMesh(Mesh_3D* mesh, const VectorXi& domains, PortBoundaryMesh& out)
{
	std::vector<int> triIds;
	const int nTri = mesh->GetNbrTri();
	for (int t = 0; t < nTri; ++t)
	{
		if (domainInList(mesh->GetDomainOfTri(t), domains))
			triIds.push_back(t);
	}
	if (triIds.empty())
	{
		std::cerr << "[PortModeSolver] no port triangles for domains" << std::endl;
		return false;
	}

	const int nf = static_cast<int>(triIds.size());
	MatrixXi oldFaces(nf, 3);
	out.facesConn = MatrixXi::Zero(nf, 2);
	out.volDomain = VectorXi::Zero(nf);

	std::vector<int> globalNodes;
	globalNodes.reserve(nf * 3);
	for (int i = 0; i < nf; ++i)
	{
		const int tri = triIds[i];
		for (int k = 0; k < 3; ++k)
		{
			oldFaces(i, k) = mesh->GetTri(tri, k);
			globalNodes.push_back(oldFaces(i, k));
		}
		Vector3i conn;
		mesh->GetCoonOfTri(tri, conn);
		out.facesConn(i, 0) = conn(0) + 1;
		out.facesConn(i, 1) = conn(1) + 1;
		out.volDomain(i) = mesh->GetDomainOfTet(conn(0));
	}

	std::sort(globalNodes.begin(), globalNodes.end());
	globalNodes.erase(std::unique(globalNodes.begin(), globalNodes.end()), globalNodes.end());
	out.nNodes = static_cast<int>(globalNodes.size());
	out.nodes = MatrixXd::Zero(out.nNodes, 3);
	std::map<int, int> g2l;
	for (int i = 0; i < out.nNodes; ++i)
	{
		g2l[globalNodes[i]] = i;
		Vector3d v;
		mesh->GetVertex(globalNodes[i], v);
		out.nodes.row(i) = v.transpose();
	}

	out.newFaces = MatrixXi::Zero(nf, 3);
	for (int i = 0; i < nf; ++i)
		for (int k = 0; k < 3; ++k)
			out.newFaces(i, k) = g2l[oldFaces(i, k)];

	// Local edges like MATLAB: (n1,n2),(n1,n3),(n2,n3) without sorting pair; unique by lexicographic rows.
	std::vector<std::array<int, 2>> edgeRows;
	edgeRows.reserve(static_cast<size_t>(nf) * 3);
	for (int i = 0; i < nf; ++i)
	{
		const int a = out.newFaces(i, 0);
		const int b = out.newFaces(i, 1);
		const int c = out.newFaces(i, 2);
		edgeRows.push_back({ a, b });
		edgeRows.push_back({ a, c });
		edgeRows.push_back({ b, c });
	}
	std::vector<std::array<int, 2>> uniqueEdges = edgeRows;
	std::sort(uniqueEdges.begin(), uniqueEdges.end());
	uniqueEdges.erase(std::unique(uniqueEdges.begin(), uniqueEdges.end()), uniqueEdges.end());
	out.nEdges = static_cast<int>(uniqueEdges.size());
	std::map<std::array<int, 2>, int> e2i;
	for (int i = 0; i < out.nEdges; ++i)
		e2i[uniqueEdges[i]] = i;

	out.edgeOfFace = MatrixXi::Zero(nf, 3);
	out.nFaces = nf;
	for (int i = 0; i < nf; ++i)
	{
		for (int k = 0; k < 3; ++k)
			out.edgeOfFace(i, k) = e2i[edgeRows[static_cast<size_t>(i) * 3 + k]];
	}
	out.edgeNodes = MatrixXi::Zero(out.nEdges, 2);
	for (int i = 0; i < out.nEdges; ++i)
	{
		out.edgeNodes(i, 0) = uniqueEdges[static_cast<size_t>(i)][0];
		out.edgeNodes(i, 1) = uniqueEdges[static_cast<size_t>(i)][1];
	}
	return true;
}

void assemblePortEigen(
	const PortBoundaryMesh& pm,
	MaterialLib* matLib,
	double lam0,
	SparseMatrix<std::complex<double>, RowMajor>& Ared,
	SparseMatrix<std::complex<double>, RowMajor>& Bred,
	SparseMatrix<std::complex<double>, RowMajor>& P,
	int& dofFull)
{
	const double k0 = 2.0 * Pi / lam0;
	const int nNodes = pm.nNodes;
	const int nEdges = pm.nEdges;
	dofFull = nNodes + nEdges;

	Gauss gauss;
	const int nGP = gauss.GetNbrGaussPoints(TWODIM, TRIANGLE, BF_LINEFUNC * 2);
	std::vector<double> u(nGP), v(nGP), w(nGP), wght(nGP);
	gauss.GetGaussPoints(TWODIM, TRIANGLE, u.data(), v.data(), w.data(), wght.data());

	BF bfEt, bfCurlEt, bfEz, bfGradEz;
	bfEt.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);
	bfCurlEt.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_NEDELEC, BF_LINEFUNC);
	bfEz.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
	bfGradEz.GetNbrBF(TWODIM, TRIANGLE, BF_GRAD_LAGRANGE, BF_LINEFUNC);

	using Trip = Triplet<std::complex<double>>;
	std::vector<Trip> tripsA, tripsB;
	tripsA.reserve(static_cast<size_t>(pm.nFaces) * 9);
	tripsB.reserve(static_cast<size_t>(pm.nFaces) * 36);

	std::vector<int> edgeCount(static_cast<size_t>(nEdges), 0);
	for (int f = 0; f < pm.nFaces; ++f)
		for (int k = 0; k < 3; ++k)
			edgeCount[static_cast<size_t>(pm.edgeOfFace(f, k))]++;

	for (int n = 0; n < pm.nFaces; ++n)
	{
		double x[3], y[3], l[3];
		for (int i = 0; i < 3; ++i)
		{
			x[i] = pm.nodes(pm.newFaces(n, i), 0);
			y[i] = pm.nodes(pm.newFaces(n, i), 1);
		}
		l[0] = std::sqrt((x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]));
		l[1] = std::sqrt((x[0] - x[2]) * (x[0] - x[2]) + (y[0] - y[2]) * (y[0] - y[2]));
		l[2] = std::sqrt((x[2] - x[1]) * (x[2] - x[1]) + (y[2] - y[1]) * (y[2] - y[1]));

		Matrix3d Jac = Matrix3d::Zero();
		Jac(0, 0) = x[1] - x[0]; Jac(0, 1) = y[1] - y[0];
		Jac(1, 0) = x[2] - x[0]; Jac(1, 1) = y[2] - y[0];
		Jac(2, 2) = 1.0;
		const Matrix3d InvJac = Jac.inverse();
		const double DetJac = std::abs(Jac.determinant());
		const Matrix3d TJac = Jac.transpose() / Jac.determinant();

		std::vector<Vector3d> Et(3 * nGP), curlEt(3 * nGP), Ez(3 * nGP), gradEz(3 * nGP);
		auto at = [nGP](std::vector<Vector3d>& a, int j, int gp) -> Vector3d& {
			return a[static_cast<size_t>(j * nGP + gp)];
		};
		for (int gp = 0; gp < nGP; ++gp)
		{
			for (int j = 0; j < 3; ++j)
			{
				Vector3d temp;
				bfEt.GetValueBF(j + 1, u[gp], v[gp], 0.0, temp);
				at(Et, j, gp) = InvJac * temp * l[j];
				bfCurlEt.GetValueBF(j + 1, u[gp], v[gp], 0.0, temp);
				at(curlEt, j, gp) = TJac * temp * l[j];
				bfEz.GetValueBF(j + 1, u[gp], v[gp], 0.0, temp);
				at(Ez, j, gp) = temp;
				bfGradEz.GetValueBF(j + 1, u[gp], v[gp], 0.0, temp);
				at(gradEz, j, gp) = InvJac * temp;
			}
		}

		const int domain = pm.volDomain(n);
		Matrix3cd epsr = matLib->GetEpsr(domain);
		const Matrix3d sigma = matLib->GetSigma(domain);
		epsr = epsr - sigma * std::complex<double>(0.0, 1.0 / k0 * 120.0 * Pi);
		const Matrix3cd mur = matLib->GetMur(domain);
		const std::complex<double> mu = mur(0, 0);
		const std::complex<double> epsilon = epsr(0, 0);

		Matrix3cd St = Matrix3cd::Zero(), Tte = Matrix3cd::Zero();
		Matrix3cd Sz = Matrix3cd::Zero(), Tz = Matrix3cd::Zero();
		Matrix3cd G = Matrix3cd::Zero(), Ttu = Matrix3cd::Zero();
		for (int i = 0; i < 3; ++i)
		{
			for (int j = 0; j < 3; ++j)
			{
				for (int gp = 0; gp < nGP; ++gp)
				{
					const double wgt = wght[gp] * DetJac;
					St(i, j) += wgt / mu * at(curlEt, i, gp).dot(at(curlEt, j, gp));
					Tte(i, j) += wgt * k0 * k0 * epsilon * at(Et, i, gp).dot(at(Et, j, gp));
					Sz(i, j) += wgt / mu * at(gradEz, i, gp).dot(at(gradEz, j, gp));
					Tz(i, j) += wgt * k0 * k0 * epsilon * at(Ez, i, gp).dot(at(Ez, j, gp));
					G(i, j) += wgt / mu * at(Et, i, gp).dot(at(gradEz, j, gp));
					Ttu(i, j) += wgt / mu * at(Et, i, gp).dot(at(Et, j, gp));
				}
			}
		}
		const Matrix3cd Gt = G.transpose();

		for (int i = 0; i < 3; ++i)
		{
			for (int j = 0; j < 3; ++j)
			{
				const int si = pm.newFaces(n, i);
				const int sj = pm.newFaces(n, j);
				const int vi = pm.edgeOfFace(n, i) + nNodes;
				const int vj = pm.edgeOfFace(n, j) + nNodes;
				tripsA.emplace_back(vi, vj, St(i, j) - Tte(i, j));
				tripsB.emplace_back(vi, vj, Ttu(i, j));
				tripsB.emplace_back(si, sj, Sz(i, j) - Tz(i, j));
				tripsB.emplace_back(vi, sj, G(i, j));
				tripsB.emplace_back(si, vj, Gt(i, j));
			}
		}
	}

	SparseMatrix<std::complex<double>, RowMajor> Afull(dofFull, dofFull), Bfull(dofFull, dofFull);
	Afull.setFromTriplets(tripsA.begin(), tripsA.end());
	Bfull.setFromTriplets(tripsB.begin(), tripsB.end());

	// Auto PEC: boundary edges (count==1) and their nodes (MATLAB EigenModeMatrixAssembly)
	std::set<int> pec;
	for (int e = 0; e < nEdges; ++e)
	{
		if (edgeCount[static_cast<size_t>(e)] != 1)
			continue;
		pec.insert(pm.edgeNodes(e, 0));
		pec.insert(pm.edgeNodes(e, 1));
		pec.insert(e + nNodes);
	}

	std::vector<int> pecInd(pec.begin(), pec.end());
	std::sort(pecInd.begin(), pecInd.end());
	std::vector<char> isPec(static_cast<size_t>(dofFull), 0);
	for (int idx : pecInd)
		if (idx >= 0 && idx < dofFull)
			isPec[static_cast<size_t>(idx)] = 1;

	std::vector<int> freeOfFull;
	freeOfFull.reserve(static_cast<size_t>(dofFull));
	for (int i = 0; i < dofFull; ++i)
		if (!isPec[static_cast<size_t>(i)])
			freeOfFull.push_back(i);
	const int nFree = static_cast<int>(freeOfFull.size());

	std::vector<Trip> pTrips;
	pTrips.reserve(static_cast<size_t>(nFree));
	for (int j = 0; j < nFree; ++j)
		pTrips.emplace_back(freeOfFull[j], j, 1.0);
	P = SparseMatrix<std::complex<double>, RowMajor>(dofFull, nFree);
	P.setFromTriplets(pTrips.begin(), pTrips.end());

	Ared = SparseMatrix<std::complex<double>, RowMajor>(nFree, nFree);
	Bred = SparseMatrix<std::complex<double>, RowMajor>(nFree, nFree);
	Ared = SparseMatrix<std::complex<double>, RowMajor>(P.transpose() * Afull * P);
	Bred = SparseMatrix<std::complex<double>, RowMajor>(P.transpose() * Bfull * P);
}

double computePowerCoef(
	const PortBoundaryMesh& pm,
	MaterialLib* matLib,
	double lam0,
	const std::complex<double>& gamma,
	const VectorXcd& Ez,
	const VectorXcd& Et)
{
	const double k0 = 2.0 * Pi / lam0;
	Gauss gauss;
	const int nGP = gauss.GetNbrGaussPoints(TWODIM, TRIANGLE, BF_LINEFUNC * 2);
	std::vector<double> u(nGP), v(nGP), w(nGP), wght(nGP);
	gauss.GetGaussPoints(TWODIM, TRIANGLE, u.data(), v.data(), w.data(), wght.data());

	BF bfEt, bfCurlEt, bfEz, bfCurlEz;
	bfEt.GetNbrBF(TWODIM, TRIANGLE, BF_NEDELEC, BF_LINEFUNC);
	bfCurlEt.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_NEDELEC, BF_LINEFUNC);
	bfEz.GetNbrBF(TWODIM, TRIANGLE, BF_LAGRANGE, BF_LINEFUNC);
	bfCurlEz.GetNbrBF(TWODIM, TRIANGLE, BF_CURL_LAGRANGE, BF_LINEFUNC);

	double power = 0.0;
	for (int n = 0; n < pm.nFaces; ++n)
	{
		double x[3], y[3], l[3];
		for (int i = 0; i < 3; ++i)
		{
			x[i] = pm.nodes(pm.newFaces(n, i), 0);
			y[i] = pm.nodes(pm.newFaces(n, i), 1);
		}
		l[0] = std::sqrt((x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]));
		l[1] = std::sqrt((x[0] - x[2]) * (x[0] - x[2]) + (y[0] - y[2]) * (y[0] - y[2]));
		l[2] = std::sqrt((x[2] - x[1]) * (x[2] - x[1]) + (y[2] - y[1]) * (y[2] - y[1]));

		Matrix3d Jac = Matrix3d::Zero();
		Jac(0, 0) = x[1] - x[0]; Jac(0, 1) = y[1] - y[0];
		Jac(1, 0) = x[2] - x[0]; Jac(1, 1) = y[2] - y[0];
		Jac(2, 2) = 1.0;
		const Matrix3d InvJac = Jac.inverse();
		const double DetJac = std::abs(Jac.determinant());
		const Matrix3d TJac = Jac.transpose() / Jac.determinant();
		Matrix3d JacS = Matrix3d::Zero();
		JacS(0, 0) = InvJac(1, 1); JacS(0, 1) = -InvJac(1, 0);
		JacS(1, 0) = -InvJac(0, 1); JacS(1, 1) = InvJac(0, 0);
		JacS(2, 2) = 1.0;

		const int domain = pm.volDomain(n);
		const double mu = std::real(matLib->GetMur(domain)(0, 0));
		const std::complex<double> physicsIndex = std::complex<double>(0.0, 1.0) / (k0 * 120.0 * Pi) / mu;

		for (int gp = 0; gp < nGP; ++gp)
		{
			Vector3cd EE = Vector3cd::Zero();
			Vector3cd curlEE = Vector3cd::Zero();
			for (int i = 0; i < 3; ++i)
			{
				Vector3d temp;
				bfEt.GetValueBF(i + 1, u[gp], v[gp], 0.0, temp);
				const Vector3d et = InvJac * temp * l[i];
				bfEz.GetValueBF(i + 1, u[gp], v[gp], 0.0, temp);
				const Vector3d ezVec(0, 0, temp(2));
				bfCurlEt.GetValueBF(i + 1, u[gp], v[gp], 0.0, temp);
				const Vector3d curlEtVec = TJac * Vector3d(0, 0, temp(2)) * l[i];
				bfCurlEz.GetValueBF(i + 1, u[gp], v[gp], 0.0, temp);
				const Vector3d curlEzVec = JacS * Vector3d(temp(0), temp(1), 0.0);

				const int mapEt = pm.edgeOfFace(n, i);
				const int mapEz = pm.newFaces(n, i);
				EE += et.cast<std::complex<double>>() * Et(mapEt) + ezVec.cast<std::complex<double>>() * Ez(mapEz);
				curlEE += curlEtVec.cast<std::complex<double>>() * Et(mapEt)
					+ curlEzVec.cast<std::complex<double>>() * Ez(mapEz);
			}
			Vector3cd HH = curlEE;
			HH(0) += EE(1) * gamma;
			HH(1) -= EE(0) * gamma;
			HH *= physicsIndex;

			const double rEx = EE(0).real(), iEx = EE(0).imag();
			const double rEy = EE(1).real(), iEy = EE(1).imag();
			const double rHx = HH(0).real(), iHx = HH(0).imag();
			const double rHy = HH(1).real(), iHy = HH(1).imag();
			power += wght[gp] * 0.5 * DetJac * (-rEy * rHx - iEy * iHx + rEx * rHy + iEx * iHy);
		}
	}
	if (!(power > 0.0))
	{
		std::cerr << "[PortModeSolver] power<=0 (" << power << "), powerCoef=1" << std::endl;
		return 1.0;
	}
	return std::sqrt(1.0 / power);
}

bool solveOnePort(
	Mesh_3D* mesh,
	MaterialLib* matLib,
	const VectorXi& domains,
	double lambda0,
	double targetNeff,
	int modeNum,
	int portType,
	bool normalizePower,
	PortNumericMode& mode)
{
	PortBoundaryMesh pm;
	if (!extractPortMesh(mesh, domains, pm))
		return false;

	SparseMatrix<std::complex<double>, RowMajor> A, B, P;
	int dofFull = 0;
	assemblePortEigen(pm, matLib, lambda0, A, B, P, dofFull);

	const double k0 = 2.0 * Pi / lambda0;
	// MATLAB: sigma = -(k0*targetNeff)^2, A x = λ B x, gamma=sqrt(λ), neff=(-i*gamma)/k0
	const std::complex<double> sigma(-(k0 * targetNeff) * (k0 * targetNeff), 0.0);

	const int nFree = static_cast<int>(A.rows());
	if (nFree <= 0)
	{
		std::cerr << "[PortModeSolver] empty free DOF after PEC" << std::endl;
		return false;
	}
	const MatrixXcd Ad = MatrixXcd(A);
	const MatrixXcd Bd = MatrixXcd(B);
	const MatrixXcd As = Ad - sigma * Bd;
	Eigen::PartialPivLU<MatrixXcd> lu(As);
	if (lu.determinant() == std::complex<double>(0, 0))
	{
		std::cerr << "[PortModeSolver] (A-σB) singular" << std::endl;
		return false;
	}
	// M = (A-σB)^{-1} B ; eig(M)=ν=1/(λ-σ) ⇒ λ = σ + 1/ν
	const MatrixXcd M = lu.solve(Bd);
	Eigen::ComplexEigenSolver<MatrixXcd> ces(M);
	if (ces.info() != Eigen::Success)
	{
		std::cerr << "[PortModeSolver] ComplexEigenSolver failed" << std::endl;
		return false;
	}
	const VectorXcd nu = ces.eigenvalues();
	const MatrixXcd V = ces.eigenvectors();
	// Prefer largest |ν| (closest to σ), then pick among top few by Re(neff)~target
	std::vector<int> order(static_cast<size_t>(nu.size()));
	for (int i = 0; i < nu.size(); ++i)
		order[static_cast<size_t>(i)] = i;
	std::sort(order.begin(), order.end(), [&](int a, int b) {
		return std::abs(nu(a)) > std::abs(nu(b));
	});
	int best = order[0];
	double bestScore = 1e300;
	const int nCand = std::min(modeNum + 8, static_cast<int>(order.size()));
	for (int c = 0; c < nCand; ++c)
	{
		const int i = order[static_cast<size_t>(c)];
		if (std::abs(nu(i)) < 1e-30)
			continue;
		const std::complex<double> lam = sigma + std::complex<double>(1.0, 0.0) / nu(i);
		const std::complex<double> gam = std::sqrt(lam);
		const std::complex<double> ne = -std::complex<double>(0.0, 1.0) * gam / k0;
		const double score = std::abs(ne.real() - targetNeff) + 0.25 * std::abs(ne.imag());
		if (score < bestScore)
		{
			bestScore = score;
			best = i;
		}
	}
	const std::complex<double> lamBest = sigma + std::complex<double>(1.0, 0.0) / nu(best);
	std::complex<double> gamma = std::sqrt(lamBest);
	// Branch: prefer Im(gamma)>0 for decaying envelope convention used with neff≈-i*gamma/k0
	std::complex<double> neff = -std::complex<double>(0.0, 1.0) * gamma / k0;
	if (neff.real() < 0.0)
	{
		gamma = -gamma;
		neff = -std::complex<double>(0.0, 1.0) * gamma / k0;
	}

	VectorXcd xFree = V.col(best);
	VectorXcd xFull = P * xFree;

	std::cout << "[PortModeSolver] type=" << portType
		<< " neff=" << neff << " gamma=" << gamma
		<< " nFaces=" << pm.nFaces << " nNodes=" << pm.nNodes
		<< " nEdges=" << pm.nEdges << " nFree=" << nFree << std::endl;

	VectorXcd Ez = xFull.head(pm.nNodes);
	VectorXcd Et = xFull.segment(pm.nNodes, pm.nEdges);
	if (std::abs(gamma) > 0.0)
		Et /= gamma;

	// Phase gauge e^{iφ}=1: make the peak-|Et| DOF real and positive (COMSOL-like Mode phase 0).
	{
		Eigen::Index iPeak = 0;
		Et.cwiseAbs().maxCoeff(&iPeak);
		const std::complex<double> peak = Et(iPeak);
		if (std::abs(peak) > 0.0)
		{
			const std::complex<double> fix = std::abs(peak) / peak;
			Et *= fix;
			Ez *= fix;
			std::cout << "[PortModeSolver] phase gauge: arg(Et_peak)->0 (e^{iφ}=1)"
				<< " iPeak=" << iPeak << std::endl;
		}
	}

	double powerCoef = 1.0;
	if (normalizePower && portType == 1)
		powerCoef = computePowerCoef(pm, matLib, lambda0, gamma, Ez, Et);

	mode.type = portType;
	mode.gamma = gamma;
	mode.powerCoef = powerCoef;
	mode.normal = (portType == 1) ? Vector3d(0, 0, -1) : Vector3d(0, 0, 1);
	mode.facesConn = pm.facesConn;
	mode.portNewFaces = pm.newFaces;
	mode.portEdgeOfFace = pm.edgeOfFace;
	mode.portNodes = pm.nodes;
	mode.Ez = Ez;
	mode.Et = Et;

	std::cout << "[PortModeSolver] powerCoef=" << powerCoef
		<< " |Et|_max=" << Et.cwiseAbs().maxCoeff()
		<< " |Ez|_max=" << Ez.cwiseAbs().maxCoeff() << std::endl;
	return true;
}

} // namespace

bool ComputePortModesFromMesh(
	Mesh_3D* mesh,
	MaterialLib* matLib,
	const VectorXi& domainsInc,
	const VectorXi& domainsOut,
	double lambda0,
	double targetNeff,
	int modeNum,
	PortModeLibrary& outModes,
	bool normalizeInputPower)
{
	outModes.Clear();
	PortNumericMode inMode, outMode;
	if (!solveOnePort(mesh, matLib, domainsInc, lambda0, targetNeff, modeNum, 1, normalizeInputPower, inMode))
		return false;
	if (!solveOnePort(mesh, matLib, domainsOut, lambda0, targetNeff, modeNum, 2, false, outMode))
		return false;
	outModes.AddMode(inMode);
	outModes.AddMode(outMode);
	return true;
}

bool SavePortModesToFile(const PortModeLibrary& modes, const std::string& path)
{
	if (modes.GetNbrModes() < 2)
		return false;
	auto pack = [](const PortNumericMode& m) {
		json j;
		j["type"] = m.type;
		j["gamma"] = { m.gamma.real(), m.gamma.imag() };
		j["powerCoef"] = m.powerCoef;
		j["normal"] = { m.normal(0), m.normal(1), m.normal(2) };
		j["facesConn"] = json::array();
		for (int i = 0; i < m.facesConn.rows(); ++i)
			j["facesConn"].push_back({ m.facesConn(i, 0), m.facesConn(i, 1) });
		j["portNewFaces"] = json::array();
		for (int i = 0; i < m.portNewFaces.rows(); ++i)
			j["portNewFaces"].push_back({ m.portNewFaces(i, 0), m.portNewFaces(i, 1), m.portNewFaces(i, 2) });
		j["portEdgeOfFace"] = json::array();
		for (int i = 0; i < m.portEdgeOfFace.rows(); ++i)
			j["portEdgeOfFace"].push_back({ m.portEdgeOfFace(i, 0), m.portEdgeOfFace(i, 1), m.portEdgeOfFace(i, 2) });
		j["portNodes"] = json::array();
		for (int i = 0; i < m.portNodes.rows(); ++i)
			j["portNodes"].push_back({ m.portNodes(i, 0), m.portNodes(i, 1), m.portNodes(i, 2) });
		j["Ez_re"] = json::array();
		j["Ez_im"] = json::array();
		for (int i = 0; i < m.Ez.size(); ++i)
		{
			j["Ez_re"].push_back(m.Ez(i).real());
			j["Ez_im"].push_back(m.Ez(i).imag());
		}
		j["Et_re"] = json::array();
		j["Et_im"] = json::array();
		for (int i = 0; i < m.Et.size(); ++i)
		{
			j["Et_re"].push_back(m.Et(i).real());
			j["Et_im"].push_back(m.Et(i).imag());
		}
		return j;
	};
	json js;
	js["input"] = pack(modes.GetMode(0));
	js["output"] = pack(modes.GetMode(1));
	std::ofstream ofs(path);
	if (!ofs.is_open())
		return false;
	ofs << js.dump();
	return true;
}