finished with bz tool

tools/bz/bz.cpp

@@ -195,6 +195,7 @@ BZDlg::BZDlg(QWidget* pParent) : QDialog{pParent},
 			(*cut)->setSingleStep(0.1);
 			(*cut)->setValue(0);
 
+			// signals
 			connect(*cut, static_cast<void (QDoubleSpinBox::*)(double)>
 					(&QDoubleSpinBox::valueChanged),
 				this, &BZDlg::CalcBZCut);
@@ -202,6 +203,11 @@ BZDlg::BZDlg(QWidget* pParent) : QDialog{pParent},
 		m_cutX->setValue(1);
 		m_cutNZ->setValue(1);
 
+		m_BZOrder = new QSpinBox(cutspanel);
+		m_BZOrder->setMinimum(0);
+		m_BZOrder->setMaximum(99);
+		m_BZOrder->setValue(4);
+
 		pGrid->addWidget(m_bzview, 0,0, 1,4);
 		pGrid->addWidget(new QLabel("In-Plane Vector:"), 1,0, 1,1);
 		pGrid->addWidget(m_cutX, 1,1, 1,1);
@@ -213,6 +219,13 @@ BZDlg::BZDlg(QWidget* pParent) : QDialog{pParent},
 		pGrid->addWidget(m_cutNZ, 2,3, 1,1);
 		pGrid->addWidget(new QLabel("Plane Offset:"), 3,0, 1,1);
 		pGrid->addWidget(m_cutD, 3,1, 1,1);
+		pGrid->addWidget(new QLabel("Peak Order:"), 3,2,1,1);
+		pGrid->addWidget(m_BZOrder, 3,3, 1,1);
+
+		// signals
+		connect(m_BZOrder,
+			static_cast<void (QSpinBox::*)(int)>(&QSpinBox::valueChanged),
+				this, [this]() { this->CalcBZCut(); });
 
 		tabs->addTab(cutspanel, "Cut");
 	}
@@ -392,3 +405,11 @@ void BZDlg::closeEvent(QCloseEvent *)
 			m_sett->setValue("geo_3dview", m_dlgPlot->saveGeometry());
 	}
 }
+
+
+void BZDlg::UpdateBZDescription()
+{
+	// brillouin zone description
+	std::string descr = m_descrBZ + "\n" + m_descrBZCut;
+	m_bz->setPlainText(descr.c_str());
+}

tools/bz/bz.h

@@ -101,11 +101,13 @@ protected:
 	QDoubleSpinBox *m_cutNY = nullptr;
 	QDoubleSpinBox *m_cutNZ = nullptr;
 	QDoubleSpinBox *m_cutD = nullptr;
+	QSpinBox *m_BZOrder = nullptr;
 
 	// brillouin zone panel
 	QPlainTextEdit *m_bz = nullptr;
 	QSpinBox *m_maxBZ = nullptr;
 	std::vector<std::vector<t_vec>> m_bz_polys;
+	std::string m_descrBZ, m_descrBZCut;
 
 	QMenu *m_pTabContextMenu = nullptr;        // menu in case a symop is selected
 	QMenu *m_pTabContextMenuNoItem = nullptr;  // menu if nothing is selected
@@ -149,6 +151,7 @@ protected:
 	void PlotShowCoordCross(bool show);
 	void PlotShowLabels(bool show);
 	void PlotShowPlane(bool show);
+	void UpdateBZDescription();
 
 	void PlotMouseDown(bool left, bool mid, bool right);
 	void PlotMouseUp(bool left, bool mid, bool right);

tools/bz/bz_calc.cpp

@@ -137,8 +137,7 @@ void BZDlg::CalcBZ(bool full_recalc)
 					idx000 = Qs_invA.size();
 
 				t_vec Q_invA = m_crystB * Q;
-				t_real Qabs_invA = tl2::norm(Q_invA);
-
+				//t_real Qabs_invA = tl2::norm(Q_invA);
 				Qs_invA.emplace_back(std::move(Q_invA));
 			}
 		}
@@ -215,11 +214,12 @@ void BZDlg::CalcBZ(bool full_recalc)
 	m_bz_polys = std::move(bz_triags);
 	PlotAddTriangles(bz_all_triags);
 
-	// brillouin zone description
-	m_bz->setPlainText(ostr.str().c_str());
+	m_descrBZ = ostr.str();
 
 	if(full_recalc)
 		CalcBZCut();
+	else
+		UpdateBZDescription();
 }
 
 
@@ -231,6 +231,9 @@ void BZDlg::CalcBZCut()
 	if(m_ignoreCalc || !m_bz_polys.size())
 		return;
 
+	std::ostringstream ostr;
+	ostr.precision(g_prec);
+
 	t_real x = m_cutX->value();
 	t_real y = m_cutY->value();
 	t_real z = m_cutZ->value();
@@ -251,24 +254,68 @@ void BZDlg::CalcBZCut()
 
 	t_mat matPlane = tl2::create<t_mat, t_vec>({ vec1, vec2, norm }, true);
 
-	std::vector<std::pair<t_vec, t_vec>> lines;
-	for(const auto& bz_poly : m_bz_polys)
-	{
-		auto vecs = tl2::intersect_plane_poly<t_vec>(
-			norm, d, bz_poly, g_eps);
-		vecs = tl2::remove_duplicates(vecs, g_eps);
+	std::vector<std::pair<t_vec, t_vec>> cut_lines, cut_lines000;
+
+	const auto order = m_BZOrder->value();
+	const auto ops = GetSymOps(true);
 
-		if(vecs.size() >= 2)
+	for(t_real h=-order; h<=order; ++h)
+	{
+		for(t_real k=-order; k<=order; ++k)
 		{
-			t_vec pt1 = matPlane * vecs[0];
-			t_vec pt2 = matPlane * vecs[1];
+			for(t_real l=-order; l<=order; ++l)
+			{
+				t_vec Q = tl2::create<t_vec>({ h, k, l });
 
-			lines.emplace_back(std::make_pair(pt1, pt2));
+				if(!is_reflection_allowed<t_mat, t_vec, t_real>(
+					Q, ops, g_eps).first)
+					continue;
+
+				// (000) peak?
+				bool is_000 = tl2::equals_0(Q, g_eps);
+				t_vec Q_invA = m_crystB * Q;
+
+				for(const auto& _bz_poly : m_bz_polys)
+				{
+					// centre bz around bragg peak
+					auto bz_poly = _bz_poly;
+					for(t_vec& vec : bz_poly)
+						vec += Q_invA;
+
+					auto vecs = tl2::intersect_plane_poly<t_vec>(
+						norm, d, bz_poly, g_eps);
+					vecs = tl2::remove_duplicates(vecs, g_eps);
+
+					if(vecs.size() >= 2)
+					{
+						t_vec pt1 = matPlane * vecs[0];
+						t_vec pt2 = matPlane * vecs[1];
+						tl2::set_eps_0(pt1, g_eps);
+						tl2::set_eps_0(pt2, g_eps);
+
+						cut_lines.emplace_back(std::make_pair(pt1, pt2));
+						if(is_000)
+							cut_lines000.emplace_back(std::make_pair(pt1, pt2));
+					}
+				}
+			}
 		}
 	}
 
+
 	m_bzscene->clear();
-	m_bzscene->AddCut(lines);
+	m_bzscene->AddCut(cut_lines);
+
+	ostr << "# Brillouin zone cut" << std::endl;
+	for(std::size_t i=0; i<cut_lines000.size(); ++i)
+	{
+		const auto& line = cut_lines000[i];
+
+		ostr << "line " << i << ":\n\tvertex 0: " << line.first
+			<< "\n\tvertex 1: " << line.second << std::endl;
+	}
+	m_descrBZCut = ostr.str();
 
 	PlotSetPlane(norm, d);
+	UpdateBZDescription();
 }

tools/bz/plot_cut.cpp

@@ -26,6 +26,7 @@
  */
 
 #include "plot_cut.h"
+#include <QtWidgets/QApplication>
 
 
 // --------------------------------------------------------------------------------
@@ -43,9 +44,14 @@ void BZCutScene::AddCut(const std::vector<std::pair<t_vec, t_vec>>& lines)
 {
 	for(const auto& line : lines)
 	{
+		QPen pen;
+		pen.setCosmetic(true);
+		pen.setColor(qApp->palette().color(QPalette::WindowText));
+
 		addLine(QLineF(
 			line.first[0]*m_scale, line.first[1]*m_scale,
-			line.second[0]*m_scale, line.second[1]*m_scale));
+			line.second[0]*m_scale, line.second[1]*m_scale),
+			pen);
 	}
 }
 // --------------------------------------------------------------------------------