• R/O
  • HTTP
  • SSH
  • HTTPS

fermisurfer: Commit

fermisurfer Git


Commit MetaInfo

Revisionf81636dc200b5619dd78d5555247ae0817cfe14f (tree)
Zeit2020-10-23 18:47:48
AutorMitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Backup

Ändern Zusammenfassung

Diff

--- /dev/null
+++ b/javascript/fermisurfer.js
@@ -0,0 +1,4414 @@
1+/*
2+The MIT License (MIT)
3+
4+Copyright (c) 2014 Mitsuaki KAWAMURA
5+
6+Permission is hereby granted, free of charge, to any person obtaining a copy
7+of this software and associated documentation files (the "Software"), to deal
8+in the Software without restriction, including without limitation the rights
9+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10+copies of the Software, and to permit persons to whom the Software is
11+furnished to do so, subject to the following conditions:
12+
13+The above copyright notice and this permission notice shall be included in
14+all copies or substantial portions of the Software.
15+
16+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
22+THE SOFTWARE.
23+*/
24+/**@file
25+@brief Mathematical operations used in various step
26+*/
27+/**
28+ @brief Work as Modulo function of fortran
29+ @return Modulated value
30+*/
31+function modulo(
32+ i = 0, //!< [in]
33+ n = 0 //!< [in]
34+) {
35+ let j = 0;
36+ j = (i + 100 * n) % n;
37+ return j;
38+};/*modulo(let i, let n)*/
39+/**
40+ @brief Solve linear system
41+ @return Determinant
42+*/
43+function solve3(
44+ a = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], //!< [in] Matix
45+ b = [0.0, 0.0, 0.0] //!< [in,out] Right hand side vector
46+)
47+{
48+ let i = 0;
49+ let det = 0.0, c = [0.0, 0.0, 0.0];
50+ /**/
51+ det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
52+ + a[0][1] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
53+ + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
54+ /**/
55+ c[0] = b[0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
56+ + b[1] * (a[0][2] * a[2][1] - a[0][1] * a[2][2])
57+ + b[2] * (a[0][1] * a[1][2] - a[0][2] * a[1][1]);
58+ /**/
59+ c[1] = b[0] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
60+ + b[1] * (a[0][0] * a[2][2] - a[0][2] * a[2][0])
61+ + b[2] * (a[0][2] * a[1][0] - a[0][0] * a[1][2]);
62+ /**/
63+ c[2] = b[0] * (a[1][0] * a[2][1] - a[1][1] * a[2][0])
64+ + b[1] * (a[0][1] * a[2][0] - a[0][0] * a[2][1])
65+ + b[2] * (a[0][0] * a[1][1] - a[0][1] * a[1][0]);
66+ /**/
67+ for (i = 0; i<3; ++i) b[i] = c[i] / det;
68+ return det;
69+}/*let solve3*/
70+/**
71+ @brief Simple sort
72+*/
73+function eigsort(
74+ n = 0, //!< [in] the number of components
75+ key = [0.0], //!< [in] Variables to be sorted [n].
76+ swap = [0] //!< [out] Order of index (sorted)
77+)
78+{
79+ let i = 0, j = 0, k = 0;
80+
81+ for (i = 0; i < n; ++i) swap[i] = i;
82+
83+ for (i = 0; i < n - 1; ++i) {
84+ for (j = i + 1; j < n; ++j) {
85+ if (key[swap[j]] < key[swap[i]]) {
86+ /*
87+ Swap
88+ */
89+ k = swap[j];
90+ swap[j] = swap[i];
91+ swap[i] = k;
92+ }/*if (sortee[j][0] < sortee[i][0])*/
93+ }/*for (j = i + 1; j < n; ++j)*/
94+ }/*for (i = 0; i < n - 1; ++i)*/
95+}/*eigsort*/
96+/**
97+ @brief Calculate normal vector from corners of triangle
98+*/
99+function normal_vec(
100+ in1 = [0.0, 0.0, 0.0], //!< [in] Corner 1
101+ in2 = [0.0, 0.0, 0.0], //!< [in] Corner 2
102+ in3 = [0.0, 0.0, 0.0], //!< [in] Corner 3
103+ out = [0.0, 0.0, 0.0] //!< [out] The normal vector
104+)
105+{
106+ let i = 0;
107+ let norm = 0.0;
108+ out[0] = in1[1] * in2[2] - in1[2] * in2[1]
109+ + in2[1] * in3[2] - in2[2] * in3[1]
110+ + in3[1] * in1[2] - in3[2] * in1[1];
111+ out[1] = in1[2] * in2[0] - in1[0] * in2[2]
112+ + in2[2] * in3[0] - in2[0] * in3[2]
113+ + in3[2] * in1[0] - in3[0] * in1[2];
114+ out[2] = in1[0] * in2[1] - in1[1] * in2[0]
115+ + in2[0] * in3[1] - in2[1] * in3[0]
116+ + in3[0] * in1[1] - in3[1] * in1[0];
117+ norm = Math.sqrt(out[0] * out[0] + out[1] * out[1] + out[2] * out[2]);
118+ for (i = 0; i<3; i++) out[i] = out[i] / norm;
119+} /* normal_vec */
120+/**
121+ @brief Judge wheser this line is the edge of 1st BZ
122+*/
123+function bragg_vert(
124+ bragg = [[0.0]],
125+ brnrm = [0.0],
126+ ibr = 0, //!< [in] Index of a Bragg plane
127+ jbr = 0, //!< [in] Index of a Bragg plane
128+ nbr = 0, //!< [in]
129+ vert = [0.0], //!< [in] start point of line
130+ vert2 = [0.0]//!< [in] end point of line
131+)
132+{
133+ let kbr = 0, i = 0, lbr = 0, nbr0 = 0;
134+ let bmat = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], rhs = [0.0, 0.0, 0.0], prod, thr, det;
135+ /**/
136+ nbr0 = nbr;
137+ /**/
138+ for (kbr = nbr0; kbr < 26; ++kbr) {
139+ /**/
140+ for (i = 0; i<3; ++i) bmat[0][i] = bragg[ibr][i];
141+ for (i = 0; i<3; ++i) bmat[1][i] = bragg[jbr][i];
142+ for (i = 0; i<3; ++i) bmat[2][i] = bragg[kbr][i];
143+ /**/
144+ rhs[0] = brnrm[ibr];
145+ rhs[1] = brnrm[jbr];
146+ rhs[2] = brnrm[kbr];
147+ thr = Math.sqrt(rhs[0] * rhs[1] * rhs[2]) * 0.001;
148+ /*
149+ if Bragg planes do not cross, roop next kbr
150+ */
151+ det = solve3(bmat, rhs);
152+ if (Math.abs(det) < thr) continue;
153+ /*
154+ if vert0 = vert1, roop next kbr
155+ */
156+ prod = (vert2[0] - rhs[0]) * (vert2[0] - rhs[0])
157+ + (vert2[1] - rhs[1]) * (vert2[1] - rhs[1])
158+ + (vert2[2] - rhs[2]) * (vert2[2] - rhs[2]);
159+ if (prod < thr) continue;
160+ /*
161+ is corner really in 1st BZ ?
162+ */
163+ i = 0;
164+ for (lbr = 0; lbr < 26; ++lbr) {
165+ prod = bragg[lbr][0] * rhs[0]
166+ + bragg[lbr][1] * rhs[1]
167+ + bragg[lbr][2] * rhs[2];
168+ /**/
169+ if (prod > brnrm[lbr] + thr) {
170+ i = 1;
171+ break;
172+ }
173+ }
174+ if (i == 1) {
175+ }
176+ else {
177+ for (i = 0; i<3; ++i) vert[i] = rhs[i];
178+ return kbr + 1;
179+ }
180+ }
181+ /*
182+ this line is not a BZ boundary
183+ */
184+ return 0;
185+}/* bragg_vert */
186+/**
187+ @brief Compute real number of Bragg plane at 1st BZ
188+*/
189+function check_bragg()
190+{
191+ let ibr = 0, ibzl = 0, ibzc = 0;
192+ let ii = 0, kk = 0, bzflag = 0, nbzcorner = 0, nn = 0;
193+ let thr = 0.0001, prod = 0.0, bzc = []; //[676][3];
194+ /*
195+ First, compute real number of corners of 1st BZ
196+ */
197+ nbzcorner = 0;
198+ for (ibzl = 0; ibzl < nbzl; ibzl++) {
199+ for (ii = 0; ii < 2; ii++) {
200+ bzflag = 0;
201+
202+ for (ibzc = 0; ibzc < nbzcorner; ibzc++) {
203+ prod = 0.0;
204+ for (kk = 0; kk < 3; kk++) prod += (bzl[ibzl][ii][kk] - bzc[ibzc][kk]) * (bzl[ibzl][ii][kk] - bzc[ibzc][kk]);
205+ if (prod < thr) bzflag = 1;
206+ }
207+
208+ if (bzflag == 0) {
209+ bzc.push([0.0, 0.0, 0.0]);
210+ for (kk = 0; kk < 3; kk++) bzc[nbzcorner][kk] = bzl[ibzl][ii][kk];
211+ nbzcorner += 1;
212+ }
213+
214+ }/*for (ii = 0; ii < 2; ii++)*/
215+ }/*for (ibzl = 0; ibzl < nbzl; ibzl++)*/
216+ terminal(" Number of corners of 1st BZ : " + toString(nbzcorner) + "\n");
217+ /**@brief
218+ Then, compute real number Bragg plane of 1st BZ (::nbragg),
219+ Re-order ::bragg and ::brnrm
220+ */
221+ nbragg = 0;
222+ for (ibr = 0; ibr < 26; ibr++) {
223+ nn = 0;
224+
225+ for (ibzc = 0; ibzc < nbzcorner; ibzc++) {
226+ prod = bragg[ibr][0] * bzc[ibzc][0] + bragg[ibr][1] * bzc[ibzc][1] + bragg[ibr][2] * bzc[ibzc][2];
227+ if (Math.abs(prod - brnrm[ibr]) < thr) nn += 1;
228+ }
229+ if (nn >= 3) {
230+ for (kk = 0; kk < 3; kk++) bragg[nbragg][kk] = bragg[ibr][kk];
231+ brnrm[nbragg] = brnrm[ibr];
232+ nbragg += 1;
233+ }
234+ }
235+ terminal(" Number of plane of 1st BZ : " + toString(nbragg) + "\n");
236+}/*function check_bragg*/
237+/**
238+ @brief Compute Brillouin zone boundariy lines
239+
240+ Modify : ::nbzl, ::bzl
241+*/
242+function bz_lines()
243+{
244+ /**/
245+ let ibr = 0, jbr = 0, nbr = 0, i = 0, j = 0, lvert = 0;
246+ let vert = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]];
247+ /**/
248+ nbzl = 0;
249+
250+ /**/
251+ for (ibr = 0; ibr < 26; ++ibr) {
252+ for (jbr = 0; jbr < 26; ++jbr) {
253+ /**/
254+ for (i = 0; i<3; ++i) vert[1][i] = 0.0;
255+ nbr = 0;
256+ lvert = bragg_vert(bragg, brnrm, ibr, jbr, nbr, vert[0], vert[1]);
257+ if (lvert == 0) continue;
258+ nbr = lvert;
259+ /**/
260+ lvert = bragg_vert(bragg, brnrm, ibr, jbr, nbr, vert[1], vert[0]);
261+ if (lvert == 0) continue;
262+ /**/
263+ for (i = 0; i < 2; ++i) for (j = 0; j < 3; ++j) bzl[nbzl][i][j] = vert[i][j];
264+ nbzl = nbzl + 1;
265+
266+ }/*for (jbr = 0; jbr < 26; ++jbr)*/
267+ }/*for (ibr = 0; ibr < 26; ++ibr)*/
268+
269+ check_bragg();
270+
271+}/*bz_lines*/
272+/**
273+@brief Compute node-line where \f$\Delta_{n k} = 0\f$
274+
275+ Modify : ::nnl, ::kvnl, ::kvnl_rot
276+
277+*/
278+function calc_nodeline()
279+{
280+ let ib = 0, itri = 0, i = 0, j = 0;
281+ let sw = [0, 0, 0];
282+ let a = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], matps = [0.0, 0.0, 0.0];
283+ let kvnl0 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]; // std::vector<std::vector<let> >
284+
285+ terminal(" band # of nodeline\n");
286+ for (ib = 0; ib < nb; ib++) {
287+
288+ kvnl.push([]);
289+
290+ for (itri = 0; itri < ntri[ib]; ++itri) {
291+
292+
293+ for (i = 0; i < 3; ++i) matps[i] = matp[ib][itri][i][0];
294+ eigsort(3, matps, sw);
295+
296+ for (i = 0; i < 3; ++i) {
297+ for (j = 0; j < 3; ++j) {
298+ a[i][j] = (0.0 - matp[ib][itri][sw[j]][0])
299+ / (matp[ib][itri][sw[i]][0] - matp[ib][itri][sw[j]][0]);
300+ }/*for (j = 0; j < 3; ++j)*/
301+ }/*for (i = 0; i < 3; ++i)*/
302+
303+ if ((matp[ib][itri][sw[0]][0] < 0.0 && 0.0 <= matp[ib][itri][sw[1]][0])
304+ || (matp[ib][itri][sw[0]][0] <= 0.0 && 0.0 < matp[ib][itri][sw[1]][0])) {
305+ for (i = 0; i < 3; ++i) {
306+ kvnl0[0][i] = kvp[ib][itri][sw[0]][i] * a[0][1] + kvp[ib][itri][sw[1]][i] * a[1][0];
307+ kvnl0[1][i] = kvp[ib][itri][sw[0]][i] * a[0][2] + kvp[ib][itri][sw[2]][i] * a[2][0];
308+ }/*for (i = 0; i < 3; ++i)*/
309+ kvnl[ib].push(kvnl0);
310+ }/*else if (matp[ib][itri][sw[0]] < 0.0 && 0.0 <= matp[ib][itri][sw[1]])*/
311+ else if ((matp[ib][itri][sw[1]][0] < 0.0 && 0.0 <= matp[ib][itri][sw[2]][0])
312+ || (matp[ib][itri][sw[1]][0] <= 0.0 && 0.0 < matp[ib][itri][sw[2]][0])) {
313+ for (i = 0; i < 3; ++i) {
314+ kvnl0[0][i] = kvp[ib][itri][sw[0]][i] * a[0][2] + kvp[ib][itri][sw[2]][i] * a[2][0];
315+ kvnl0[1][i] = kvp[ib][itri][sw[1]][i] * a[1][2] + kvp[ib][itri][sw[2]][i] * a[2][1];
316+ }/*for (i = 0; i < 3; ++i)*/
317+ kvnl[ib].push(kvnl0);
318+ }/*else if (matp[ib][itri][sw[1]] < 0.0 && 0.0 <= matp[ib][itri][sw[2]])*/
319+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
320+ /*
321+ Allocation of nodeline
322+ */
323+ nnl[ib] = kvnl[ib].size;
324+ terminal(" " + toString(ib + 1) + " " + toString(nnl[ib]) + "\n");
325+ }/*for (ib = 0; ib < nb; ib++)*/
326+}/*function calc_nodeline()*/
327+/**
328+ @brief Draw Fermi surfaces
329+
330+ Modify : ::nmlp_rot, ::kvp_rot, ::kvnl_rot
331+
332+ First, rotate ::nmlp and ::kvp with OpenMP then draw them.
333+ Also draw nodeline in the same way.
334+*/
335+function draw_fermi() {
336+ let ib = 0, a0 = 0, a1 = 0, a2 = 0, ia = 0;
337+ let kshift = [0.0, 0.0, 0.0];
338+ let i = 0, j = 0, l = 0, itri = 0;
339+ /*
340+ First, rotate k-vector and normal vector
341+ */
342+ for (a0 = 0; a0 < BZ_number[0]; a0++) {
343+ for (a1 = 0; a1 < BZ_number[1]; a1++) {
344+ for (a2 = 0; a2 < BZ_number[2]; a2++) {
345+ for (ia = 0; ia < 3; ia++) kshift[ia] = bvec[0][ia] * a0 + bvec[1][ia] * a1 + bvec[2][ia] * a2;
346+
347+ for (ib = 0; ib < nb; ib++) {
348+ if (draw_band[ib] == 1) {
349+ for (itri = 0; itri < ntri[ib]; ++itri) {
350+ for (i = 0; i < 3; ++i) {
351+ for (j = 0; j < 3; ++j) {
352+ kvp_rot[ib][j + 3 * i + 9 * itri]
353+ = rot[j][0] * (kvp[ib][itri][i][0] + kshift[0])
354+ + rot[j][1] * (kvp[ib][itri][i][1] + kshift[1])
355+ + rot[j][2] * (kvp[ib][itri][i][2] + kshift[2])
356+ + trans[j];
357+ nmlp_rot[ib][j + 3 * i + 9 * itri]
358+ = rot[j][0] * nmlp[ib][itri][i][0]
359+ + rot[j][1] * nmlp[ib][itri][i][1]
360+ + rot[j][2] * nmlp[ib][itri][i][2];
361+ nmlp_rot[ib][j + 3 * i + 9 * itri] *= side;
362+ for (l = 0; l < 2; ++l) {
363+ arw_rot[ib][j + 3 * l + 6 * i + 18 * itri]
364+ = rot[j][0] * (arw[ib][itri][i][l][0] + kshift[0])
365+ + rot[j][1] * (arw[ib][itri][i][l][1] + kshift[1])
366+ + rot[j][2] * (arw[ib][itri][i][l][2] + kshift[2])
367+ + trans[j];
368+ }
369+ }
370+ }/*for (i = 0; i < 3; ++i)*/
371+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
372+ }/*if (draw_band[ib] == 1)*/
373+ }/*for (ib = 0; ib < nb; ib++)*/
374+ /*
375+ Second, draw each triangle
376+ */
377+ glEnableClientState(GL_NORMAL_ARRAY);
378+ glEnableClientState(GL_COLOR_ARRAY);
379+ for (ib = 0; ib < nb; ib++) {
380+ if (draw_band[ib] == 1) {
381+ glVertexPointer(3, GL_FLOAT, 0, kvp_rot[ib]);
382+ glNormalPointer(GL_FLOAT, 0, nmlp_rot[ib]);
383+ glColorPointer(4, GL_FLOAT, 0, clr[ib]);
384+ glDrawArrays(GL_TRIANGLES, 0, ntri[ib] * 3);
385+ }/*if (draw_band[ib] == 1)*/
386+ }/*for (ib = 0; ib < nb; ib++)*/
387+ glDisableClientState(GL_NORMAL_ARRAY);
388+ glDisableClientState(GL_COLOR_ARRAY);
389+ /*
390+ Arrow for 3D value
391+ */
392+ glNormal3f(0.0, 0.0, 1.0);
393+ if (color_scale == 3) {
394+ for (ib = 0; ib < nb; ib++) {
395+ if (draw_band[ib] == 1) {
396+ glColor3f(1.0 - clr[ib][0], 1.0 - clr[ib][1], 1.0 - clr[ib][2]);
397+ glVertexPointer(3, GL_FLOAT, 0, arw_rot[ib]);
398+ glDrawArrays(GL_LINES, 0, ntri[ib] * 3 * 2);
399+ }
400+ }
401+ }
402+ /*
403+ Nodeline
404+ */
405+ if (nodeline == 1) {
406+ /*
407+ First, rotate k-vector
408+ */
409+ let i, j, itri;
410+ for (ib = 0; ib < nb; ib++) {
411+ if (draw_band[ib] == 1) {
412+ for (itri = 0; itri < nnl[ib]; ++itri) {
413+ for (i = 0; i < 2; ++i) {
414+ for (j = 0; j < 3; ++j)
415+ kvnl_rot[ib][j + 3 * i + 6 * itri]
416+ = rot[j][0] * (kvnl[ib][itri][i][0] + kshift[0])
417+ + rot[j][1] * (kvnl[ib][itri][i][1] + kshift[1])
418+ + rot[j][2] * (kvnl[ib][itri][i][2] + kshift[2])
419+ + trans[j];
420+ kvnl_rot[ib][2 + 3 * i + 6 * itri] += 0.001;
421+ }/*for (i = 0; i < 2; ++i)*/
422+ }/*for (itri = 0; itri < nnl[ib]; ++itri)*/
423+ }/*if (draw_band[ib] == 1)*/
424+ }/*for (ib = 0; ib < nb; ib++)*/
425+ /*
426+ Second, draw each lines
427+ */
428+ glColor3fv(black);
429+ glNormal3f(0.0, 0.0, 1.0);
430+ for (ib = 0; ib < nb; ib++) {
431+ if (draw_band[ib] == 1) {
432+ glVertexPointer(3, GL_FLOAT, 0, kvnl_rot[ib]);
433+ glDrawArrays(GL_LINES, 0, 2 * nnl[ib]);
434+ }/*if (draw_band[ib] == 1)*/
435+ }/* for (ib = 0; ib < nb; ib++)*/
436+ }/*if (nodeline == 1)*/
437+ /*
438+ Equator
439+ */
440+ if (lequator == 1) {
441+ /*
442+ First, rotate k-vector
443+ */
444+ for (ib = 0; ib < nb; ib++) {
445+ if (draw_band[ib] == 1) {
446+ for (itri = 0; itri < nequator[ib]; ++itri) {
447+ for (i = 0; i < 2; ++i) {
448+ for (j = 0; j < 3; ++j)
449+ kveq_rot[ib][j + 3 * i + 6 * itri]
450+ = rot[j][0] * (kveq[ib][itri][i][0] + kshift[0])
451+ + rot[j][1] * (kveq[ib][itri][i][1] + kshift[1])
452+ + rot[j][2] * (kveq[ib][itri][i][2] + kshift[2])
453+ + trans[j];
454+ kveq_rot[ib][2 + 3 * i + 6 * itri] += 0.001;
455+ }/*for (i = 0; i < 2; ++i)*/
456+ }/*for (itri = 0; itri < nequator[ib]; ++itri)*/
457+ }/*if (draw_band[ib] == 1)*/
458+ }/*for (ib = 0; ib < nb; ib++)*/
459+ /*
460+ Second, draw each lines
461+ */
462+ glColor3fv(black);
463+ glNormal3f(0.0, 0.0, 1.0);
464+ for (ib = 0; ib < nb; ib++) {
465+ if (draw_band[ib] == 1) {
466+ glVertexPointer(3, GL_FLOAT, 0, kveq_rot[ib]);
467+ glDrawArrays(GL_LINES, 0, 2 * nequator[ib]);
468+ }/*if (draw_band[ib] == 1)*/
469+ }/* for (ib = 0; ib < nb; ib++)*/
470+ }/*if (nodeline == 1)*/
471+ }
472+ }
473+ }
474+}/*function draw_fermi*/
475+/**
476+ @brief Draw lines of BZ boundaries
477+*/
478+function draw_bz_lines() {
479+ let ibzl = 0, i = 0, j = 0, a0 = 0, a1 = 0, a2 = 0, ia = 0;
480+ let bzl2 = [0.0, 0.0, 0.0], bvec2 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
481+ linecolor = [0.0, 0.0, 0.0, 0.0], secvec2 = [0.0, 0.0, 0.0], kshift = [0.0, 0.0, 0.0];
482+ let vertices = [];
483+ /*
484+ Line color is oposit of BG color
485+ */
486+ for (i = 0; i < 4; i++) linecolor[i] = LineColor[i];
487+ /**/
488+ glLineWidth(linewidth);
489+ for (a0 = 0; a0 < BZ_number[0]; a0++) {
490+ for (a1 = 0; a1 < BZ_number[1]; a1++) {
491+ for (a2 = 0; a2 < BZ_number[2]; a2++) {
492+ for (ia = 0; ia < 3; ia++) kshift[ia] = bvec[0][ia] * a0 + bvec[1][ia] * a1 + bvec[2][ia] * a2;
493+ /*
494+ First Brillouin zone mode
495+ */
496+ if (fbz == 1) {
497+ for (ibzl = 0; ibzl < nbzl; ++ibzl) {
498+ for (i = 0; i < 2; ++i) {
499+ for (j = 0; j < 3; ++j)
500+ bzl2[j] = rot[j][0] * (bzl[ibzl][i][0] + kshift[0])
501+ + rot[j][1] * (bzl[ibzl][i][1] + kshift[1])
502+ + rot[j][2] * (bzl[ibzl][i][2] + kshift[2])
503+ + trans[j];
504+ for (j = 0; j < 3; j++) vertices.push(bzl2[j]);
505+ }/*for (i = 0; i< 2; ++i)*/
506+ glColor3fv(linecolor);
507+ glNormal3f(0.0, 0.0, 1.0);
508+ glVertexPointer(3, GL_FLOAT, 0, vertices);
509+ glDrawArrays(GL_LINES, 0, 2);
510+ }/*for (ibzl = 0; ibzl < nbzl; ++ibzl)*/
511+ }/*if (fbz == 1)*/
512+ else {
513+ /*
514+ Premitive BZ mode
515+ */
516+ for (i = 0; i < 3; ++i) {
517+ for (j = 0; j < 3; ++j) {
518+ bvec2[i][j] = rot[j][0] * bvec[i][0]
519+ + rot[j][1] * bvec[i][1]
520+ + rot[j][2] * bvec[i][2];
521+ }/*for (j = 0; j < 3; ++j)*/
522+ }/*for (i = 0; i < 3; ++i)*/
523+ for (i = 0; i < 3; ++i) vertices.push(trans[i]);
524+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[0][i]);
525+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[0][i] + bvec2[1][i]);
526+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[0][i] + bvec2[1][i] + bvec2[2][i]);
527+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[1][i] + bvec2[2][i]);
528+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[1][i]);
529+ for (i = 0; i < 3; ++i) vertices.push(trans[i]);
530+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[2][i]);
531+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[0][i] + bvec2[2][i]);
532+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[0][i] + bvec2[1][i] + bvec2[2][i]);
533+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[0][i] + bvec2[2][i]);
534+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[0][i]);
535+ for (i = 0; i < 3; ++i) vertices.push(trans[i]);
536+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[2][i]);
537+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[1][i] + bvec2[2][i]);
538+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[1][i]);
539+ for (i = 0; i < 3; ++i) vertices.push(trans[i] + bvec2[0][i] + bvec2[1][i]);
540+ for (i = 0; i < 17; ++i) {
541+ for (j = 0; j < 3; ++j) {
542+ vertices[j + 3 * i] = vertices[j + 3 * i]
543+ + rot[j][0] * kshift[0]
544+ + rot[j][1] * kshift[1]
545+ + rot[j][2] * kshift[2];
546+ }/*for (j = 0; j < 3; ++j)*/
547+ }/*for (i = 0; i < 3; ++i)*/
548+ glColor3fv(linecolor);
549+ glNormal3f(0.0, 0.0, 1.0);
550+ glVertexPointer(3, GL_FLOAT, 0, vertices);
551+ glDrawArrays(GL_LINE_STRIP, 0, 17);
552+ }/*if (fbz != 1)*/
553+ }
554+ }
555+ }
556+ /*
557+ Section for the 2D Fermi line
558+ */
559+ if (lsection == 1 && fbz == 1) {
560+ for (j = 0; j < 3; ++j)
561+ secvec2[j] = rot[j][0] * secvec[0]
562+ + rot[j][1] * secvec[1]
563+ + rot[j][2] * secvec[2];
564+ for (ibzl = 0; ibzl < nbzl2d; ++ibzl) {
565+ for (j = 0; j < 3; ++j)
566+ bzl2[j] = rot[j][0] * bzl2d[ibzl][0]
567+ + rot[j][1] * bzl2d[ibzl][1]
568+ + rot[j][2] * bzl2d[ibzl][2]
569+ + trans[j];
570+ for (j = 0; j < 3; j++)vertices[j + 3 * ibzl] = bzl2[j];
571+ }/*for (ibzl = 0; ibzl < nbzl2d; ++ibzl)*/
572+ glColor3fv(gray);
573+ glNormal3fv(secvec2);
574+ glVertexPointer(3, GL_FLOAT, 0, vertices);
575+ glDrawArrays(GL_TRIANGLE_FAN, 0, nbzl2d);
576+ }/*if (lsection == 1)*/
577+}/*draw bz_lines */
578+/**
579+ @brief Draw color-bar or colr-circle (periodic)
580+ as a color scale
581+*/
582+function draw_colorbar()
583+{
584+ let i = 0, j = 0, k = 0;
585+ let mat2 = 0.0, vertices = [], colors = [], vector = [], vector_color = [],
586+ norm = 0.0;
587+
588+ glEnableClientState(GL_COLOR_ARRAY);
589+ /*
590+ Reciplocal lattice vectors
591+ */
592+ for (i = 0; i < 6; i++) {
593+ vector[3 * i] = -1.2;
594+ vector[3 * i + 1] = -1.05;
595+ vector[3 * i + 2] = 0.0;
596+ }
597+ for (i = 0; i < 3; i++) {
598+ norm = Math.sqrt(bvec[i][0] * bvec[i][0]+ bvec[i][1] * bvec[i][1]+ bvec[i][2] * bvec[i][2]);
599+ for (j = 0; j < 3; j++) {
600+ vector[j + 6 * i]
601+ += (rot[j][0] * bvec[i][0]
602+ + rot[j][1] * bvec[i][1]
603+ + rot[j][2] * bvec[i][2]) * 0.2 / norm;
604+ }
605+ for (j = 0; j < 4; j++) {
606+ vector_color[j + 8 * i] = BarColor[i * 2][j];
607+ vector_color[j + 8 * i + 4] = BarColor[i * 2][j];
608+ }
609+ }
610+ glNormal3f(0.0, 0.0, 1.0);
611+ glVertexPointer(3, GL_FLOAT, 0, vector);
612+ glColorPointer(4, GL_FLOAT, 0, vector_color);
613+ glDrawArrays(GL_LINES, 0, 6);
614+ /*
615+ Color bar/circle/cube
616+ */
617+ if (color_scale == 1 || color_scale == 4) {
618+ for (i = 0; i < 5; i++) {
619+ for (j = 0; j < 2; j++) {
620+ vertices[0 + j * 3 + i * 6] = -1.0 + 0.5*i;
621+ vertices[1 + j * 3 + i * 6] = -1.0 - 0.1*j;
622+ vertices[2 + j * 3 + i * 6] = 0.0;
623+ for (k = 0; k < 4; k++) colors[k + 4 * j + 8 * i] = BarColor[i][k];
624+ }
625+ }/*for (i = 0; i < 10; i++)*/
626+ glNormal3f(0.0, 0.0, 1.0);
627+ glVertexPointer(3, GL_FLOAT, 0, vertices);
628+ glColorPointer(4, GL_FLOAT, 0, colors);
629+ glDrawArrays(GL_TRIANGLE_STRIP, 0, 10);
630+ }/*if (color_scale == 1 || color_scale == 4)*/
631+ else if (color_scale == 2) {
632+ /*
633+ Periodic color scale
634+ */
635+ vertices[0] = 0.0;
636+ vertices[1] = -1.0;
637+ vertices[2] = 0.0;
638+ for (j = 0; j < 4; j++) colors[j] = 1.0 - BackGroundColor[j];
639+ /**/
640+ for (i = 0; i <= 60; i++) {
641+ /**/
642+ mat2 = i / 60.0 * 6.0;
643+ /**/
644+ if (mat2 <= 1.0) {
645+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = yellow[j] * mat2 + red[j] * (1.0 - mat2);
646+ }
647+ else if (mat2 <= 2.0) {
648+ mat2 = mat2 - 1.0;
649+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = green[j] * mat2 + yellow[j] * (1.0 - mat2);
650+ }
651+ else if (mat2 <= 3.0) {
652+ mat2 = mat2 - 2.0;
653+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = cyan[j] * mat2 + green[j] * (1.0 - mat2);
654+ }
655+ else if (mat2 <= 4.0) {
656+ mat2 = mat2 - 3.0;
657+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = blue[j] * mat2 + cyan[j] * (1.0 - mat2);
658+ }
659+ else if (mat2 <= 5.0) {
660+ mat2 = mat2 - 4.0;
661+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = magenta[j] * mat2 + blue[j] * (1.0 - mat2);
662+ }
663+ else {
664+ mat2 = mat2 - 5.0;
665+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = red[j] * mat2 + magenta[j] * (1.0 - mat2);
666+ }
667+ /**/
668+ vertices[0 + 3 * (i + 1)] = 0.2 * Math.cos(i / 60.0 * 6.283185307);
669+ vertices[1 + 3 * (i + 1)] = 0.2 * Math.sin(i / 60.0 * 6.283185307) - 1.0;
670+ vertices[2 + 3 * (i + 1)] = 0.0;
671+ }/*for (i = 0; i <= 60; i++)*/
672+ glNormal3f(0.0, 0.0, 1.0);
673+ glVertexPointer(3, GL_FLOAT, 0, vertices);
674+ glColorPointer(4, GL_FLOAT, 0, colors);
675+ glDrawArrays(GL_TRIANGLE_FAN, 0, 62);
676+ }/*else if (color_scale == 2)*/
677+ else if (color_scale == 6 || color_scale == 7) {
678+ for (i = 0; i < 2; i++) {
679+ for (j = 0; j < 2; j++) {
680+ vertices[0 + j * 3 + i * 6] = -1.0 + 2.0*i;
681+ vertices[1 + j * 3 + i * 6] = -1.0 - 0.1*j;
682+ vertices[2 + j * 3 + i * 6] = 0.0;
683+ if (i == 0) for (k = 0; k < 4; k++) colors[k + 4 * j + 8 * i] = bgray[k];
684+ else if (i == 1) for (k = 0; k < 4; k++) colors[k + 4 * j + 8 * i] = wgray[k];
685+ }
686+ }/*for (i = 0; i < 10; i++)*/
687+ glNormal3f(0.0, 0.0, 1.0);
688+ glVertexPointer(3, GL_FLOAT, 0, vertices);
689+ glColorPointer(4, GL_FLOAT, 0, colors);
690+ glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
691+ }/*if (color_scale == 6 || color_scale == 7)*/
692+ glDisableClientState(GL_COLOR_ARRAY);
693+}/*function draw_colorbar*/
694+/**
695+ @brief Draw eye-points for the stereogram
696+*/
697+function draw_circles(
698+ dx2d = 0.0 //!< [in] Translation used for the section-mode
699+) {
700+ let i = 0, j = 0;
701+ let r = 0.0, vertices = [];
702+ /**/
703+ r = 0.05;
704+ /**/
705+ vertices[0] = 0.7 - dx2d;
706+ vertices[1] = scl;
707+ vertices[2] = 0.0;
708+ for (i = 0; i <= 20; i++) {
709+ vertices[0 + (i + 1) * 3] = r * Math.cos(i / 20.0 * 6.283185307) + 0.7 - dx2d;
710+ vertices[1 + (i + 1) * 3] = r * Math.sin(i / 20.0 * 6.283185307) + scl;
711+ vertices[2 + (i + 1) * 3] = 0.0;
712+ }/*for (i = 0; i <= 20; i++)*/
713+ glNormal3f(0.0, 0.0, 1.0);
714+ glColor3fv(LineColor);
715+ glVertexPointer(3, GL_FLOAT, 0, vertices);
716+ glDrawArrays(GL_TRIANGLE_FAN, 0, 22);
717+ /**/
718+ for (i = 0; i < 22; i++) vertices[3 * i] += -1.4;
719+ glVertexPointer(3, GL_FLOAT, 0, vertices);
720+ glDrawArrays(GL_TRIANGLE_FAN, 0, 22);
721+}/*function draw_circles*/
722+/**
723+ @brief Draw 2D Fermi lines
724+*/
725+function draw_fermi_line() {
726+ let i = 0, ib = 0, ibzl = 0;
727+ let vertices = [];
728+ /*
729+ Draw 2D BZ lines
730+ */
731+ glLineWidth(linewidth);
732+ for (ibzl = 0; ibzl < nbzl2d; ++ibzl) {
733+ for (i = 0; i < 3; i++) vertices[i + 3 * ibzl] = bzl2d_proj[ibzl][i];
734+ }/*for (ibzl = 0; ibzl < nbzl2d; ++ibzl)*/
735+ glNormal3f(0.0, 0.0, 1.0);
736+ glColor3fv(LineColor);
737+ glVertexPointer(3, GL_FLOAT, 0, vertices);
738+ glDrawArrays(GL_LINE_LOOP, 0, nbzl2d);
739+ /*
740+ Draw Fermi lines
741+ */
742+ glLineWidth(linewidth);
743+ glEnableClientState(GL_COLOR_ARRAY);
744+ glNormal3f(0.0, 0.0, 1.0);
745+ for (ib = 0; ib < nb; ib++) {
746+ if (draw_band[ib] == 1) {
747+ glVertexPointer(3, GL_FLOAT, 0, kv2d[ib]);
748+ glColorPointer(4, GL_FLOAT, 0, clr2d[ib]);
749+ glDrawArrays(GL_LINES, 0, 2 * n2d[ib]);
750+ }/*if (draw_band[ib] == 1)*/
751+ }/* for (ib = 0; ib < nb; ib++)*/
752+ glDisableClientState(GL_COLOR_ARRAY);
753+}/*function draw_fermi_line*/
754+/**
755+ @brief Glut Display function
756+ called by glutDisplayFunc
757+*/
758+function draw()
759+{
760+ let pos = [1.0, 1.0, 1.0, 0.0];
761+ let amb = [0.2, 0.2, 0.2, 0.0];
762+ let dx = 0.0, dx2d = 0.0, theta = 0.0, posz = 0.0, phi = 0.0;
763+ let pos1 = [0.0, 0.0, 0.0, 0.0], pos2 = [0.0, 0.0, 0.0, 0.0];
764+
765+ if (draw_band == NULL) return;
766+
767+ if (lstereo == 2) {
768+ /*
769+ Parallel eyes
770+ */
771+ theta = 3.14159265359 / 180.0 * 2.5;
772+ posz = 5.0;
773+ dx = 0.7;
774+ phi = Math.atan(posz / dx) - theta;
775+ theta = (3.1416 * 0.5 - phi) / 3.1416 * 180.0;
776+ /**/
777+ pos1[0] = posz * Math.cos(phi) - dx;
778+ pos1[1] = 0.0;
779+ pos1[2] = posz * Math.sin(phi);
780+ pos1[3] = 1.0;
781+ /**/
782+ pos2[0] = -posz * Math.cos(phi) + dx;
783+ pos2[1] = 0.0;
784+ pos2[2] = posz * Math.sin(phi);
785+ pos2[3] = 1.0;
786+ }/*if (lstereo == 2)*/
787+ else if (lstereo == 3) {
788+ /*
789+ Cross eyes
790+ */
791+ theta = -3.1416 / 180.0 * 2.0;
792+ posz = 5.0;
793+ dx = -0.7;
794+ /**/
795+ pos1[0] = posz * Math.sin(theta) - dx;
796+ pos1[1] = 0.0;
797+ pos1[2] = posz * Math.cos(theta);
798+ pos1[3] = 1.0;
799+ /**/
800+ pos2[0] = -posz * Math.sin(theta) + dx;
801+ pos2[1] = 0.0;
802+ pos2[2] = posz * Math.cos(theta);
803+ pos2[3] = 1.0;
804+ }/*if (lstereo == 3)*/
805+ else {
806+ theta = 0.0;
807+ dx = 0.0;
808+ }/*lstero == 1*/
809+ if (lsection == 1 && fbz == 1) dx2d = 0.7;
810+ else dx2d = 0.0;
811+ /*
812+ Initialize
813+ */
814+ glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
815+ glLoadIdentity();
816+ /*
817+ Set view point & view line
818+ */
819+ glTranslatef(0.0, 0.0, -5.0);
820+ //gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
821+ /*
822+ Set position of light
823+ */
824+ if (lstereo == 1) {
825+ glLightfv(GL_LIGHT0, GL_POSITION, pos);
826+ glTranslatef(-dx2d, 0.0, 0.0);
827+ /*
828+ Draw color scale
829+ */
830+ if (lcolorbar == 1) draw_colorbar();
831+ }
832+ else {
833+ glLightfv(GL_LIGHT0, GL_POSITION, pos1);
834+ draw_circles(dx2d);
835+ glTranslatef(-dx-dx2d, 0.0, 0.0);
836+ glRotatef(theta, 0.0, 1.0, 0.0);
837+ }
838+ glLightfv(GL_LIGHT1, GL_AMBIENT, amb);
839+ /*
840+ Rotation & Zoom
841+ */
842+ glScalef(scl, scl, scl);
843+ /*
844+ Draw Brillouin zone boundaries
845+ */
846+ draw_bz_lines();
847+ /*
848+ Draw Fermi surfaces
849+ */
850+ draw_fermi();
851+ /*
852+ Draw the second object for stereogram
853+ */
854+ if (lstereo != 1) {
855+ glPushMatrix();
856+ glLoadIdentity();
857+ glTranslatef(0.0, 0.0, -5.0);
858+ //gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
859+ glLightfv(GL_LIGHT0, GL_POSITION, pos2);
860+ /**/
861+ glTranslatef(dx-dx2d, 0.0, 0.0);
862+ glRotatef(-theta, 0.0, 1.0, 0.0);
863+ /**/
864+ glScalef(scl, scl, scl);
865+ draw_bz_lines();
866+ draw_fermi();
867+ /**/
868+ glPopMatrix();
869+ }
870+ /*
871+ Draw 2D Fermi line
872+ */
873+ if (lsection == 1 && fbz == 1) {
874+ glPushMatrix();
875+ glLoadIdentity();
876+ glTranslatef(0.0, 0.0, -5.0);
877+ //gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
878+ glLightfv(GL_LIGHT0, GL_POSITION, pos);
879+ /**/
880+ if (lstereo == 1) glTranslatef(dx2d, 0.0, 0.0);
881+ else glTranslatef(2.0 * dx2d, 0.0, 0.0);
882+ /**/
883+ glScalef(scl, scl, scl);
884+ draw_fermi_line();
885+ /**/
886+ glPopMatrix();
887+ }/*if (lsection == 1)*/
888+ glFlush(); // Not really necessary: buffer swapping below implies glFlush()
889+ SwapBuffers();
890+}/*function display*/
891+/**
892+@brief Compute equator \f$\{\bf v}_{n k} \cdot {\bf k} = 0\f$
893+
894+ Modify : ::nequator, ::kveq, ::kveq_rot
895+*/
896+function equator() {
897+ let ib = 0, itri = 0, i = 0, j = 0, nequator0 = 0;
898+ let sw = [0,0,0];
899+ let a = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], prod = [0.0, 0.0, 0.0];
900+ let kveq_0 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]];
901+
902+ terminal(" band # of equator\n");
903+ for (ib = 0; ib < nb; ib++) {
904+ kveq.push([]);
905+
906+ for (itri = 0; itri < ntri[ib]; ++itri) {
907+ for (i = 0; i < 3; ++i) {
908+ prod[i] = 0.0;
909+ for (j = 0; j < 3; ++j) prod[i] += eqvec[j] * nmlp[ib][itri][i][j];
910+ }/*for (i = 0; i < 3; ++i)*/
911+
912+ eigsort(3, prod, sw);
913+ for (i = 0; i < 3; ++i) {
914+ for (j = 0; j < 3; ++j) {
915+ a[i][j] = (0.0 - prod[sw[j]]) / (prod[sw[i]] - prod[sw[j]]);
916+ }/*for (j = 0; j < 3; ++j)*/
917+ }/*for (i = 0; i < 3; ++i)*/
918+
919+ if ((prod[sw[0]] < 0.0 && 0.0 <= prod[sw[1]])
920+ || (prod[sw[0]] <= 0.0 && 0.0 < prod[sw[1]])) {
921+ for (i = 0; i < 3; ++i) {
922+ kveq_0[0][i] = kvp[ib][itri][sw[0]][i] * a[0][1] + kvp[ib][itri][sw[1]][i] * a[1][0];
923+ kveq_0[1][i] = kvp[ib][itri][sw[0]][i] * a[0][2] + kvp[ib][itri][sw[2]][i] * a[2][0];
924+ }/*for (i = 0; i < 3; ++i)*/
925+ kveq[ib].push(kveq_0);
926+ }/*else if (prod[sw[0]] < 0.0 && 0.0 <= prod[sw[1]])*/
927+ else if ((prod[sw[1]] < 0.0 && 0.0 <= prod[sw[2]])
928+ || (prod[sw[1]] <= 0.0 && 0.0 < prod[sw[2]])) {
929+ for (i = 0; i < 3; ++i) {
930+ kveq_0[0][i] = kvp[ib][itri][sw[0]][i] * a[0][2] + kvp[ib][itri][sw[2]][i] * a[2][0];
931+ kveq_0[1][i] = kvp[ib][itri][sw[1]][i] * a[1][2] + kvp[ib][itri][sw[2]][i] * a[2][1];
932+ }/*for (i = 0; i < 3; ++i)*/
933+ kveq[ib].push(kveq_0);
934+ }/*else if (prod[sw[1]] < 0.0 && 0.0 <= prod[sw[2]])*/
935+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
936+ /*
937+ Sum node-lines in all threads
938+ */
939+ nequator.push(kveq[ib].size);
940+ terminal(" " + toString(ib + 1) + " " + toString(nequator[ib])+ "\n");
941+ }/*for (ib = 0; ib < nb; ib++)*/
942+}/*function equator()*/
943+/**
944+ @brief Store triangle patch
945+
946+ Modify : ::matp, ::kvp
947+
948+ For the 1st BZ mode, this routine cuts triangle recursivly at the
949+ BZ boundary (Bragg plane).
950+
951+ - DO @f${\bf l}@f$ in Bragg vector
952+ - @f${p_i = {\bf l}\cdot{\bf k}}@f$
953+ - Sort : @f$p_0<p_1<p_2@f$
954+ - @f[
955+ a_{i j} \equiv \frac{-p_j}{p_i - p_j}
956+ @f]
957+ - if (@f$|{\bf l}| < p_0@f$)
958+ - This patch is not in the 1st BZ
959+ - if (@f$p_0 < |{\bf l}| < p_1@f$)
960+ - @f${\bf k}'_0 = {\bf k}_0@f$
961+ - @f${\bf k}'_1 = {\bf k}_0 a_{0 1} + {\bf k}_1 a_{1 0}@f$
962+ - @f${\bf k}'_2 = {\bf k}_0 a_{0 2} + {\bf k}_2 a_{2 0}@f$
963+ - if (@f$p_1 < |{\bf l}| < p_2@f$)
964+ - @f${\bf k}'_0 = {\bf k}_0@f$
965+ - @f${\bf k}'_1 = {\bf k}_1@f$
966+ - @f${\bf k}'_2 = {\bf k}_0 a_{0 2} + {\bf k}_2 a_{2 0}@f$
967+ - and
968+ - @f${\bf k}'_0 = {\bf k}_1 a_{1 2} + {\bf k}_2 a_{2 1}@f$
969+ - @f${\bf k}'_1 = {\bf k}_1@f$
970+ - @f${\bf k}'_2 = {\bf k}_0 a_{0 2} + {\bf k}_2 a_{2 0}@f$
971+ - if (@f$p_2 < |{\bf l}| < p_3@f$)
972+ - @f${\bf k}'_0 = {\bf k}_0@f$
973+ - @f${\bf k}'_1 = {\bf k}_1@f$
974+ - @f${\bf k}'_2 = {\bf k}_2@f$
975+ - END DO
976+*/
977+function triangle(
978+ nbr = 1, //!<[in] Bragg plane
979+ mat1 = [[0.0,0.0,0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], //!<[in] The matrix element
980+ kvec1 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], //!<[in] @f$k@f$-vector of corners
981+ vf1 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], //!<[in] @f$v_f@f$-vector of corners
982+ kvp_v =[],
983+ matp_v = [],
984+ nmlp_v = []
985+ )
986+{
987+ let ibr = 0, i = 0, j = 0, sw = [0, 0, 0];
988+ let prod = [0.0, 0.0, 0.0], thr = 0.0, thr2 = 0.001,
989+ mat2 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
990+ kvec2 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
991+ vf2 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
992+ a = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
993+ bshift = 0.0,
994+ vfave = [0.0, 0.0, 0.0], norm = [0.0, 0.0, 0.0],
995+ kvp_0 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
996+ matp_0 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
997+ nmlp_0 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]];
998+
999+ /*
1000+ If the area is nearly 0, it is ignored.
1001+ */
1002+ for (i = 0; i < 3; i++)norm[i] = 0.0;
1003+ for (i = 0; i < 3; i++) {
1004+ norm[0] += (kvec1[1][i] - kvec1[2][i])*(kvec1[1][i] - kvec1[2][i]);
1005+ norm[1] += (kvec1[2][i] - kvec1[0][i])*(kvec1[2][i] - kvec1[0][i]);
1006+ norm[2] += (kvec1[0][i] - kvec1[1][i])*(kvec1[0][i] - kvec1[1][i]);
1007+ }
1008+ for (i = 0; i < 3; i++) {
1009+ if (norm[i] < 1.0e-10*brnrm_min) return;
1010+ }
1011+ /*
1012+ For 1st BZ, it is cut at the BZ boundary.
1013+ */
1014+ if (fbz == 1) {
1015+ /**/
1016+ for (ibr = 0; ibr < nbragg; ++ibr) {
1017+
1018+ thr = brnrm[ibr] * 0.001;
1019+ /**/
1020+ for (i = 0; i < 3; ++i)
1021+ prod[i] = bragg[ibr][0] * kvec1[i][0]
1022+ + bragg[ibr][1] * kvec1[i][1]
1023+ + bragg[ibr][2] * kvec1[i][2];
1024+ eigsort(3, prod, sw);
1025+ for (i = 0; i < 3; ++i) {
1026+ for (j = 0; j < 3; ++j) {
1027+ a[i][j] = (brnrm[ibr] - prod[sw[j]]) / (prod[sw[i]] - prod[sw[j]]);
1028+ }/*for (j = 0; j < 3; ++j)*/
1029+ }/*for (i = 0; i < 3; ++i)*/
1030+ i = (0.5 * ((prod[sw[2]] / brnrm[ibr]) + 1.0));
1031+ bshift = -2.0 *i;
1032+
1033+ if (brnrm[ibr] + thr > prod[sw[2]]) continue;
1034+
1035+ if (brnrm[ibr] < prod[sw[0]]) {
1036+ /*
1037+ All corners are outside of the Bragg plane
1038+ */
1039+ for (i = 0; i < 3; ++i) {
1040+ for (j = 0; j < 3; ++j) {
1041+ kvec2[i][j] = kvec1[sw[i]][j] + bshift * bragg[ibr][j];
1042+ mat2[i][j] = mat1[sw[i]][j];
1043+ vf2[i][j] = vf1[sw[i]][j];
1044+ }
1045+ }
1046+ triangle(ibr + 1, mat2, kvec2, vf2, kvp_v, matp_v, nmlp_v);
1047+ return;
1048+ }
1049+ else if (brnrm[ibr] < prod[sw[1]]) {
1050+ /*
1051+ Single corner (#0) is inside of the Bragg plane
1052+ */
1053+ for (i = 0; i < 3; ++i) {
1054+ kvec2[0][i] = kvec1[sw[0]][i];
1055+ kvec2[1][i] = kvec1[sw[0]][i] * a[0][1] + kvec1[sw[1]][i] * a[1][0];
1056+ kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
1057+
1058+ mat2[0][i] = mat1[sw[0]][i];
1059+ mat2[1][i] = mat1[sw[0]][i] * a[0][1] + mat1[sw[1]][i] * a[1][0];
1060+ mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
1061+
1062+ vf2[0][i] = vf1[sw[0]][i];
1063+ vf2[1][i] = vf1[sw[0]][i] * a[0][1] + vf1[sw[1]][i] * a[1][0];
1064+ vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
1065+ }/*for (i = 0; i < 3; ++i)*/
1066+ triangle(ibr + 1, mat2, kvec2, vf2, kvp_v, matp_v, nmlp_v);
1067+
1068+ for (i = 0; i < 3; ++i) {
1069+ kvec2[0][i] = kvec1[sw[0]][i] * a[0][1] + kvec1[sw[1]][i] * a[1][0];
1070+ kvec2[1][i] = kvec1[sw[1]][i];
1071+ kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
1072+
1073+ mat2[0][i] = mat1[sw[0]][i] * a[0][1] + mat1[sw[1]][i] * a[1][0];
1074+ mat2[1][i] = mat1[sw[1]][i];
1075+ mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
1076+
1077+ vf2[0][i] = vf1[sw[0]][i] * a[0][1] + vf1[sw[1]][i] * a[1][0];
1078+ vf2[1][i] = vf1[sw[1]][i];
1079+ vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
1080+ }/*for (i = 0; i < 3; ++i)*/
1081+ for (i = 0; i < 3; ++i) for (j = 0; j < 3; ++j)
1082+ kvec2[i][j] += bshift * bragg[ibr][j];
1083+ triangle(ibr + 1, mat2, kvec2, vf2, kvp_v, matp_v, nmlp_v);
1084+
1085+ for (i = 0; i < 3; ++i) {
1086+ kvec2[0][i] = kvec1[sw[2]][i];
1087+ kvec2[1][i] = kvec1[sw[1]][i];
1088+ kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
1089+
1090+ mat2[0][i] = mat1[sw[2]][i];
1091+ mat2[1][i] = mat1[sw[1]][i];
1092+ mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
1093+
1094+ vf2[0][i] = vf1[sw[2]][i];
1095+ vf2[1][i] = vf1[sw[1]][i];
1096+ vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
1097+ }/*for (i = 0; i < 3; ++i)*/
1098+ for (i = 0; i < 3; ++i) for (j = 0; j < 3; ++j)
1099+ kvec2[i][j] += bshift * bragg[ibr][j];
1100+ triangle(ibr + 1, mat2, kvec2, vf2, kvp_v, matp_v, nmlp_v);
1101+ return;
1102+ }
1103+ else if (brnrm[ibr] < prod[sw[2]]) {
1104+ /*
1105+ Two corners (#0, #1) are inside of the Bragg plane
1106+ */
1107+ for (i = 0; i < 3; ++i) {
1108+ kvec2[0][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
1109+ kvec2[1][i] = kvec1[sw[1]][i] * a[1][2] + kvec1[sw[2]][i] * a[2][1];
1110+ kvec2[2][i] = kvec1[sw[2]][i];
1111+
1112+ mat2[0][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
1113+ mat2[1][i] = mat1[sw[1]][i] * a[1][2] + mat1[sw[2]][i] * a[2][1];
1114+ mat2[2][i] = mat1[sw[2]][i];
1115+
1116+ vf2[0][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
1117+ vf2[1][i] = vf1[sw[1]][i] * a[1][2] + vf1[sw[2]][i] * a[2][1];
1118+ vf2[2][i] = vf1[sw[2]][i];
1119+ }/*for (i = 0; i < 3; ++i)*/
1120+ for (i = 0; i < 3; ++i) for (j = 0; j < 3; ++j)
1121+ kvec2[i][j] += bshift * bragg[ibr][j];
1122+ triangle(ibr + 1, mat2, kvec2, vf2, kvp_v, matp_v, nmlp_v);
1123+
1124+ for (i = 0; i < 3; ++i) {
1125+ kvec2[0][i] = kvec1[sw[0]][i];
1126+ kvec2[1][i] = kvec1[sw[1]][i];
1127+ kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
1128+
1129+ mat2[0][i] = mat1[sw[0]][i];
1130+ mat2[1][i] = mat1[sw[1]][i];
1131+ mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
1132+
1133+ vf2[0][i] = vf1[sw[0]][i];
1134+ vf2[1][i] = vf1[sw[1]][i];
1135+ vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
1136+ }/*for (i = 0; i < 3; ++i)*/
1137+ triangle(ibr + 1, mat2, kvec2, vf2, kvp_v, matp_v, nmlp_v);
1138+
1139+ for (i = 0; i < 3; ++i) {
1140+ kvec2[0][i] = kvec1[sw[1]][i] * a[1][2] + kvec1[sw[2]][i] * a[2][1];
1141+ kvec2[1][i] = kvec1[sw[1]][i];
1142+ kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
1143+
1144+ mat2[0][i] = mat1[sw[1]][i] * a[1][2] + mat1[sw[2]][i] * a[2][1];
1145+ mat2[1][i] = mat1[sw[1]][i];
1146+ mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
1147+
1148+ vf2[0][i] = vf1[sw[1]][i] * a[1][2] + vf1[sw[2]][i] * a[2][1];
1149+ vf2[1][i] = vf1[sw[1]][i];
1150+ vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
1151+ }/*for (i = 0; i < 3; ++i)*/
1152+ triangle(ibr + 1, mat2, kvec2, vf2, kvp_v, matp_v, nmlp_v);
1153+ return;
1154+ }
1155+ else {
1156+ /*
1157+ All corners are inside of the Bragg plane
1158+ */
1159+ } /* brnrm[ibr] + thr < prod */
1160+ } /* for ibr = 1; ibr < nbragg*/
1161+ } /* if fbz == 1 */
1162+ /*
1163+ Store patch
1164+ */
1165+ normal_vec(kvec1[0], kvec1[1], kvec1[2], norm);
1166+ for (i = 0; i < 3; ++i) {
1167+ vfave[i] = 0.0;
1168+ for (j = 0; j < 3; ++j) vfave[i] += vf1[j][i];
1169+ }
1170+ prod[0] = 0.0;
1171+ for (i = 0; i < 3; ++i) prod[0] += vfave[i] * norm[i];
1172+
1173+ if (prod[0] < 0.0) {
1174+ for (i = 0; i < 3; ++i) {
1175+ for (j = 0; j < 3; ++j) {
1176+ kvp_0[i][j] = kvec1[2 - i][j];
1177+ matp_0[i][j] = mat1[2 - i][j];
1178+ nmlp_0[i][j] = vf1[2 - i][j];
1179+ }
1180+ }/*for (i = 0; i < 3; ++i)*/
1181+ }
1182+ else {
1183+ for (i = 0; i < 3; ++i) {
1184+ for (j = 0; j < 3; ++j) {
1185+ kvp_0[i][j] = kvec1[i][j];
1186+ matp_0[i][j] = mat1[i][j];
1187+ nmlp_0[i][j] = vf1[i][j];
1188+ }
1189+ }/*for (i = 0; i < 3; ++i)*/
1190+ }
1191+ kvp_v.push(kvp_0);
1192+ matp_v.push(matp_0);
1193+ nmlp_v.push(nmlp_0);
1194+}/* triangle */
1195+/**
1196+@brief Cut triangle patch with the tetrahedron method.
1197+
1198+ - Sort : @f$\varepsilon_0<\varepsilon_1<\varepsilon_2<\varepsilon_3@f$
1199+ - @f[
1200+ a_{i j} \equiv \frac{-\varepsilon_j}{\varepsilon_i - \varepsilon_j}
1201+ @f]
1202+ - if (@f$\varepsilon_0 < 0 < \varepsilon_1@f$)
1203+ - @f${\bf k}'_0 = {\bf k}_0 a_{0 1} + {\bf k}_1 a_{1 0}@f$
1204+ - @f${\bf k}'_1 = {\bf k}_0 a_{0 2} + {\bf k}_2 a_{2 0}@f$
1205+ - @f${\bf k}'_2 = {\bf k}_0 a_{0 3} + {\bf k}_3 a_{3 0}@f$
1206+ - if (@f$\varepsilon_1 < 0 < \varepsilon_2@f$)
1207+ - @f${\bf k}'_0 = {\bf k}_0 a_{0 2} + {\bf k}_2 a_{2 0}@f$
1208+ - @f${\bf k}'_1 = {\bf k}_0 a_{0 3} + {\bf k}_3 a_{3 0}@f$
1209+ - @f${\bf k}'_2 = {\bf k}_1 a_{1 2} + {\bf k}_2 a_{2 1}@f$
1210+ - and
1211+ - @f${\bf k}'_0 = {\bf k}_1 a_{1 3} + {\bf k}_3 a_{3 1}@f$
1212+ - @f${\bf k}'_1 = {\bf k}_0 a_{0 3} + {\bf k}_3 a_{3 0}@f$
1213+ - @f${\bf k}'_2 = {\bf k}_1 a_{1 2} + {\bf k}_2 a_{2 1}@f$
1214+ - if (@f$\varepsilon_2 < 0 < \varepsilon_3@f$)
1215+ - @f${\bf k}'_0 = {\bf k}_3 a_{3 0} + {\bf k}_0 a_{0 3}@f$
1216+ - @f${\bf k}'_1 = {\bf k}_3 a_{3 1} + {\bf k}_1 a_{1 3}@f$
1217+ - @f${\bf k}'_2 = {\bf k}_3 a_{3 2} + {\bf k}_2 a_{2 3}@f$
1218+*/
1219+function tetrahedron(
1220+ eig1 = [], //!< [in] Orbital energies @f$\varepsilon_{n k}@f$
1221+ mat1 = [], //!< [in] Matrix elements @f$\Delta_{n k}@f$
1222+ kvec1 = [], //!< [in] @f$k@f$-vectors
1223+ vf1 = [], //!< [in] @f$v_f@f$-vectors
1224+ kvp_v = [],
1225+ matp_v = [],
1226+ nmlp_v = []
1227+)
1228+{
1229+ let it = 0, i = 0, j = 0, sw = [0, 0, 0, 0];
1230+ let eig2 = [0.0, 0.0, 0.0, 0.0],
1231+ mat2 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
1232+ kvec2 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
1233+ vf2 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
1234+ a = [[0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0]],
1235+ kvec3 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
1236+ mat3 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
1237+ vf3 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
1238+ vol = 0.0, thr = 0.000;
1239+
1240+ for (it = 0; it < 6; ++it) {
1241+ /*
1242+ Define corners of the tetrahedron
1243+ */
1244+ for (i = 0; i < 4; ++i) {
1245+ eig2[i] = eig1[corner[it][i]];
1246+ for (j = 0; j < 3; ++j) {
1247+ mat2[i][j] = mat1[corner[it][i]][j];
1248+ vf2[i][j] = vf1[corner[it][i]][j];
1249+ }
1250+ /*
1251+ Fractional -> Cartecian
1252+ */
1253+ for (j = 0; j < 3; ++j)
1254+ kvec2[i][j] = bvec[0][j] * kvec1[corner[it][i]][0]
1255+ + bvec[1][j] * kvec1[corner[it][i]][1]
1256+ + bvec[2][j] * kvec1[corner[it][i]][2];
1257+ }/*for (i = 0; i < 4; ++i)*/
1258+ eigsort(4, eig2, sw);
1259+
1260+ for (i = 0; i < 4; ++i) {
1261+ for (j = 0; j < 4; ++j) {
1262+ a[i][j] = (0.0 - eig2[sw[j]]) / (eig2[sw[i]] - eig2[sw[j]]);
1263+ }/*for (j = 0; j < 4; ++j)*/
1264+ }/*for (i = 0; i < 4; ++i)*/
1265+ /*
1266+ Draw triangle in each cases
1267+ */
1268+ if (eig2[sw[0]] <= 0.0 && 0.0 < eig2[sw[1]]) {
1269+ for (i = 0; i < 3; ++i) {
1270+ kvec3[0][i] = kvec2[sw[0]][i] * a[0][1] + kvec2[sw[1]][i] * a[1][0];
1271+ kvec3[1][i] = kvec2[sw[0]][i] * a[0][2] + kvec2[sw[2]][i] * a[2][0];
1272+ kvec3[2][i] = kvec2[sw[0]][i] * a[0][3] + kvec2[sw[3]][i] * a[3][0];
1273+
1274+ mat3[0][i] = mat2[sw[0]][i] * a[0][1] + mat2[sw[1]][i] * a[1][0];
1275+ mat3[1][i] = mat2[sw[0]][i] * a[0][2] + mat2[sw[2]][i] * a[2][0];
1276+ mat3[2][i] = mat2[sw[0]][i] * a[0][3] + mat2[sw[3]][i] * a[3][0];
1277+
1278+ vf3[0][i] = vf2[sw[0]][i] * a[0][1] + vf2[sw[1]][i] * a[1][0];
1279+ vf3[1][i] = vf2[sw[0]][i] * a[0][2] + vf2[sw[2]][i] * a[2][0];
1280+ vf3[2][i] = vf2[sw[0]][i] * a[0][3] + vf2[sw[3]][i] * a[3][0];
1281+ }
1282+
1283+ vol = a[1][0] * a[2][0] * a[3][0];
1284+ triangle(0, mat3, kvec3, vf3, kvp_v, matp_v, nmlp_v);
1285+ }
1286+ else if (eig2[sw[1]] <= 0.0 && 0.0 < eig2[sw[2]]) {
1287+ for (i = 0; i < 3; ++i) {
1288+ kvec3[0][i] = kvec2[sw[0]][i] * a[0][2] + kvec2[sw[2]][i] * a[2][0];
1289+ kvec3[1][i] = kvec2[sw[0]][i] * a[0][3] + kvec2[sw[3]][i] * a[3][0];
1290+ kvec3[2][i] = kvec2[sw[1]][i] * a[1][2] + kvec2[sw[2]][i] * a[2][1];
1291+
1292+ mat3[0][i] = mat2[sw[0]][i] * a[0][2] + mat2[sw[2]][i] * a[2][0];
1293+ mat3[1][i] = mat2[sw[0]][i] * a[0][3] + mat2[sw[3]][i] * a[3][0];
1294+ mat3[2][i] = mat2[sw[1]][i] * a[1][2] + mat2[sw[2]][i] * a[2][1];
1295+
1296+ vf3[0][i] = vf2[sw[0]][i] * a[0][2] + vf2[sw[2]][i] * a[2][0];
1297+ vf3[1][i] = vf2[sw[0]][i] * a[0][3] + vf2[sw[3]][i] * a[3][0];
1298+ vf3[2][i] = vf2[sw[1]][i] * a[1][2] + vf2[sw[2]][i] * a[2][1];
1299+ }
1300+
1301+ vol = a[1][2] * a[2][0] * a[3][0];
1302+ triangle(0, mat3, kvec3, vf3, kvp_v, matp_v, nmlp_v);
1303+ /**/
1304+ for (i = 0; i < 3; ++i) {
1305+ kvec3[0][i] = kvec2[sw[1]][i] * a[1][3] + kvec2[sw[3]][i] * a[3][1];
1306+ kvec3[1][i] = kvec2[sw[0]][i] * a[0][3] + kvec2[sw[3]][i] * a[3][0];
1307+ kvec3[2][i] = kvec2[sw[1]][i] * a[1][2] + kvec2[sw[2]][i] * a[2][1];
1308+
1309+ mat3[0][i] = mat2[sw[1]][i] * a[1][3] + mat2[sw[3]][i] * a[3][1];
1310+ mat3[1][i] = mat2[sw[0]][i] * a[0][3] + mat2[sw[3]][i] * a[3][0];
1311+ mat3[2][i] = mat2[sw[1]][i] * a[1][2] + mat2[sw[2]][i] * a[2][1];
1312+
1313+ vf3[0][i] = vf2[sw[1]][i] * a[1][3] + vf2[sw[3]][i] * a[3][1];
1314+ vf3[1][i] = vf2[sw[0]][i] * a[0][3] + vf2[sw[3]][i] * a[3][0];
1315+ vf3[2][i] = vf2[sw[1]][i] * a[1][2] + vf2[sw[2]][i] * a[2][1];
1316+ }
1317+
1318+ vol = a[1][3] * a[3][0] * a[2][1];
1319+ triangle(0, mat3, kvec3, vf3, kvp_v, matp_v, nmlp_v);
1320+ }
1321+ else if (eig2[sw[2]] <= 0.0 && 0.0 < eig2[sw[3]]) {
1322+ for (i = 0; i < 3; ++i) {
1323+ kvec3[0][i] = kvec2[sw[3]][i] * a[3][0] + kvec2[sw[0]][i] * a[0][3];
1324+ kvec3[1][i] = kvec2[sw[3]][i] * a[3][1] + kvec2[sw[1]][i] * a[1][3];
1325+ kvec3[2][i] = kvec2[sw[3]][i] * a[3][2] + kvec2[sw[2]][i] * a[2][3];
1326+
1327+ mat3[0][i] = mat2[sw[3]][i] * a[3][0] + mat2[sw[0]][i] * a[0][3];
1328+ mat3[1][i] = mat2[sw[3]][i] * a[3][1] + mat2[sw[1]][i] * a[1][3];
1329+ mat3[2][i] = mat2[sw[3]][i] * a[3][2] + mat2[sw[2]][i] * a[2][3];
1330+
1331+ vf3[0][i] = vf2[sw[3]][i] * a[3][0] + vf2[sw[0]][i] * a[0][3];
1332+ vf3[1][i] = vf2[sw[3]][i] * a[3][1] + vf2[sw[1]][i] * a[1][3];
1333+ vf3[2][i] = vf2[sw[3]][i] * a[3][2] + vf2[sw[2]][i] * a[2][3];
1334+ }
1335+
1336+ vol = a[0][3] * a[1][3] * a[2][3];
1337+ triangle(0, mat3, kvec3, vf3, kvp_v, matp_v, nmlp_v);
1338+ }
1339+ else {
1340+ }
1341+ }/*for (it = 0; it < 6; ++it)*/
1342+}/* tetrahedron */
1343+/**
1344+ @brief Compute patches for Fermi surfaces
1345+
1346+ Modify : ::ntri, nmlp, ::matp, ::kvp, ::clr, ::nmlp_rot, ::kvp_rot
1347+*/
1348+function fermi_patch()
1349+{
1350+ let ntri0 = 0, ib = 0, i0 = 0, i1 = 0, j0 = 0, start = [0, 0, 0], last = [0, 0, 0];
1351+ let i = 0, j = 0, i2 = 0, j1 = 0, j2 = 0, ii0 = 0, ii1 = 0, ii2 = 0;
1352+ let kvec1 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0],
1353+ [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
1354+ mat1 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0],
1355+ [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
1356+ eig1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1357+ vf1 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0],
1358+ [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]];
1359+
1360+ if (fbz == 1) {
1361+ terminal("\n");
1362+ terminal(" ## First Brillouin zone mode #######\n");
1363+ for (i0 = 0; i0 < 3; ++i0) {
1364+ start[i0] = ng[i0] / 2 - ng[i0];
1365+ last[i0] = ng[i0] / 2;
1366+ }
1367+ }
1368+ else {
1369+ terminal("\n");
1370+ terminal(" ## Premitive Brillouin zone mode #######\n");
1371+ for (i0 = 0; i0 < 3; ++i0) {
1372+ start[i0] = 0;
1373+ last[i0] = ng[i0];
1374+ }
1375+ }
1376+ terminal(" Computing patch ...\n");
1377+ terminal(" band # of patchs\n");
1378+ /**/
1379+ for (ib = 0; ib < nb; ++ib) {
1380+ for (j0 = start[0]; j0 < last[0]; ++j0) {
1381+ for (j1 = start[1]; j1 < last[1]; ++j1) {
1382+ for (j2 = start[2]; j2 < last[2]; ++j2) {
1383+
1384+ i0 = j0;
1385+ i1 = j1;
1386+ i2 = j2;
1387+ ii0 = j0 + 1;
1388+ ii1 = j1 + 1;
1389+ ii2 = j2 + 1;
1390+ /**/
1391+ kvec1[0][0] = i0 / ng[0];
1392+ kvec1[1][0] = i0 / ng[0];
1393+ kvec1[2][0] = i0 / ng[0];
1394+ kvec1[3][0] = i0 / ng[0];
1395+ kvec1[4][0] = ii0 / ng[0];
1396+ kvec1[5][0] = ii0 / ng[0];
1397+ kvec1[6][0] = ii0 / ng[0];
1398+ kvec1[7][0] = ii0 / ng[0];
1399+ /**/
1400+ kvec1[0][1] = i1 / ng[1];
1401+ kvec1[1][1] = i1 / ng[1];
1402+ kvec1[2][1] = ii1 / ng[1];
1403+ kvec1[3][1] = ii1 / ng[1];
1404+ kvec1[4][1] = i1 / ng[1];
1405+ kvec1[5][1] = i1 / ng[1];
1406+ kvec1[6][1] = ii1 / ng[1];
1407+ kvec1[7][1] = ii1 / ng[1];
1408+ /**/
1409+ kvec1[0][2] = i2 / ng[2];
1410+ kvec1[1][2] = ii2 / ng[2];
1411+ kvec1[2][2] = i2 / ng[2];
1412+ kvec1[3][2] = ii2 / ng[2];
1413+ kvec1[4][2] = i2 / ng[2];
1414+ kvec1[5][2] = ii2 / ng[2];
1415+ kvec1[6][2] = i2 / ng[2];
1416+ kvec1[7][2] = ii2 / ng[2];
1417+ /**/
1418+ for (i = 0; i < 8; i++)
1419+ for (j = 0; j < 3; j++)
1420+ kvec1[i][j] = kvec1[i][j] + shiftk[j] / (2.0 * ng0[j]);
1421+ /**/
1422+ i0 = modulo(i0, ng[0]);
1423+ i1 = modulo(i1, ng[1]);
1424+ i2 = modulo(i2, ng[2]);
1425+ ii0 = modulo(ii0, ng[0]);
1426+ ii1 = modulo(ii1, ng[1]);
1427+ ii2 = modulo(ii2, ng[2]);
1428+ /**/
1429+ eig1[0] = eig[ib][i0][i1][i2] - EF;
1430+ eig1[1] = eig[ib][i0][i1][ii2] - EF;
1431+ eig1[2] = eig[ib][i0][ii1][i2] - EF;
1432+ eig1[3] = eig[ib][i0][ii1][ii2] - EF;
1433+ eig1[4] = eig[ib][ii0][i1][i2] - EF;
1434+ eig1[5] = eig[ib][ii0][i1][ii2] - EF;
1435+ eig1[6] = eig[ib][ii0][ii1][i2] - EF;
1436+ eig1[7] = eig[ib][ii0][ii1][ii2] - EF;
1437+ /**/
1438+ for (j = 0; j < 3; j++) {
1439+ mat1[0][j] = mat[ib][i0][i1][i2][j];
1440+ mat1[1][j] = mat[ib][i0][i1][ii2][j];
1441+ mat1[2][j] = mat[ib][i0][ii1][i2][j];
1442+ mat1[3][j] = mat[ib][i0][ii1][ii2][j];
1443+ mat1[4][j] = mat[ib][ii0][i1][i2][j];
1444+ mat1[5][j] = mat[ib][ii0][i1][ii2][j];
1445+ mat1[6][j] = mat[ib][ii0][ii1][i2][j];
1446+ mat1[7][j] = mat[ib][ii0][ii1][ii2][j];
1447+ /**/
1448+ vf1[0][j] = vf[ib][i0][i1][i2][j];
1449+ vf1[1][j] = vf[ib][i0][i1][ii2][j];
1450+ vf1[2][j] = vf[ib][i0][ii1][i2][j];
1451+ vf1[3][j] = vf[ib][i0][ii1][ii2][j];
1452+ vf1[4][j] = vf[ib][ii0][i1][i2][j];
1453+ vf1[5][j] = vf[ib][ii0][i1][ii2][j];
1454+ vf1[6][j] = vf[ib][ii0][ii1][i2][j];
1455+ vf1[7][j] = vf[ib][ii0][ii1][ii2][j];
1456+ }/*for (j = 0; j < 3; j++)*/
1457+ /**/
1458+ tetrahedron(eig1, mat1, kvec1, vf1,
1459+ kvp[ib], matp[ib], nmlp[ib]);
1460+ }/*for (j0 = start[0]; j0 < ng[0]; ++j0)*/
1461+ }/*for (j1 = start[1]; j1 < ng[1]; ++j1)*/
1462+ }/*for (j0 = start[0]; j0 < ng[0]; ++j0)*/
1463+ }/*for (ib = 0; ib < nb; ++ib)*/
1464+ terminal(" ... Done\n");
1465+} /* fermi_patch */
1466+/**@file
1467+ @brief Main routine
1468+*/
1469+/**@mainpage FermiSurfer Main Page
1470+
1471+Fermisurfer displays Fermi surfaces
1472+with a color-plot of the arbitraly matrix element
1473+
1474+@section Notation
1475+
1476+- @f$\varepsilon_{n k}@f$ : Energy
1477+- @f$\Delta_{n k}@f$ : Any @f$(n, k)@f$-dependent value for the color-plot.
1478+- @f$N_b@f$ : Number of bands
1479+
1480+@section sec_routine Important routines
1481+- Main routine: main()
1482+- Right click menu : FS_CreateMenu(), FS_ModifyMenu(let status)
1483+
1484+@section sec_file Important files
1485+- Main routine : fermisurfer.cpp
1486+- Global valiable : variable.hpp
1487+
1488+@section sec_flow Flow
1489+
1490+- main() :
1491+ - read_file() : Read .frmsf file
1492+ - compute_patch_segment() : Compute patch and segment
1493+ - fermi_patch() : Compute patches constructing Fermi surface
1494+ - display() : Display figures with OpenGL
1495+
1496+*/
1497+/*
1498+ Input variables
1499+*/
1500+let ng0 = [0, 0, 0]; //!< @f$k@f$-point grid in the input file
1501+let shiftk = [0, 0, 0]; //!< Wherether @f$k@f$-grid is shifted or not
1502+let nb = 0; //!< The number of Bands
1503+let avec = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]; //!< Direct lattice vector
1504+let bvec = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]; //!< Reciprocal lattice vector
1505+let eig0 = []; //!< Eigenvalues @f$\varepsilon_{n k}@f$[::nb][::ng0[0]][::ng0[1]][::ng0[2]]
1506+let mat0 = []; //!< Matrix element [::nb][::ng0[0]][::ng0[1]][::ng0[2]][3]
1507+/*
1508+ Interpolation
1509+*/
1510+let ng = [0, 0, 0]; //!< @b Interpolated @f$k@f$-grids
1511+let eig = []; //!< Eigenvalues @f$\varepsilon_{n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]]
1512+let mat =[]; //!< Matrix element @f$\delta_{n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]][3]
1513+let vf =[]; //!< Matrix element @f$\{\bf v}_{{\rm f} n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]][3]
1514+let interpol = 1; //!< Ratio of interpolation
1515+/*
1516+ Switch for some modes
1517+*/
1518+let color_scale = 1; //!< Switch for full color scale mode
1519+let fbz = 1; //!< Switch for 1st Brillouin zone mode
1520+let nodeline = 0; //!< Switch for node lines
1521+let lcolorbar = 1; //!< Switch for colorbar
1522+let lstereo = 1; //!< Switch for the stereogram
1523+let lmouse = 1; //!< Switch for the mouse function
1524+let lsection = 0; //!< Switch for the 2D Fermi lines
1525+let lequator = 0; //!< Switch for equator
1526+let BZ_number = [1, 1, 1];
1527+/*
1528+ Variables for Brillouin zone boundaries
1529+*/
1530+let nbzl = 1; //!< The number of Lines of 1st Brillouin zone
1531+let bzl = []; //!< Lines of 1st BZ [nbzl(max:26*26=676)][2][3]
1532+let bragg = []; //!< Bragg plane vectors
1533+let brnrm = []; //!< Norms of Bragg plane vectors
1534+let brnrm_min = 0.0; //!< Minimum scale of the reciplocal space
1535+let nbragg = 0; //!< Number of Bragg plane og 1st BZ
1536+/*
1537+ Variables for patchs
1538+*/
1539+let ntri = []; //!< The number of triangle patch [::nb]
1540+let draw_band = []; //!< Switch for drawn bands [::nb]
1541+let nmlp = []; //!< Normal vector of patchs [::nb][::ntri][3][3]
1542+let kvp = []; //!< @f$k@f$-vectors of points [::nb][::ntri][3][3]
1543+let arw = [];
1544+let nmlp_rot = []; //!< Normal vector of patchs [::nb][::ntri*3*3]
1545+let kvp_rot = []; //!< @f$k@f$-vectors of points [::nb][::ntri*3*3]
1546+let arw_rot = [];
1547+let matp = []; //!< Matrix elements of points [::nb][::ntri][3][3]
1548+let clr = []; //!< Colors of points [::nb][::ntri*3*4]
1549+let itet = 0; //!< Counter for tetrahedron
1550+let side = 1.0; //!< Which side is lighted
1551+let patch_max = 0.0; //!< Max value across patch
1552+let patch_min = 0.0; //!< Max value across patch
1553+/*
1554+ Variables for nodeline
1555+*/
1556+let nnl = []; //!< The number of nodeline
1557+let kvnl = []; //!< @f$k@f$-vector of nodeline [::nb][::nnl][2][3]
1558+let kvnl_rot = []; //!< @f$k@f$-vector of nodeline [::nb][::nnl*2*3]
1559+/*
1560+ 2D Fermi line
1561+*/
1562+let secvec = [0.0, 0.0, 0.0]; //!< @f$k@f$-vector to define section
1563+let secvec_fr = [0.0, 0.0, 0.0]; //!< @f$k@f$-vector to define section
1564+let secscale; //!< 0.0 (across @f$\Gamma@f$) or 1.0
1565+let axis2d = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]; //!< @f$k@f$-vector to define section
1566+let n2d = []; //!< Number of line segment
1567+let kv2d = []; //!< @f$k@f$-vector for 2D plot [::nb][::n2d*2*3]
1568+let clr2d = []; //!< Matrix element for 2D plot [::nb][::n2d*2*4]
1569+let nbzl2d = 0; //!< The number of Lines of 1st Brillouin zone
1570+let bzl2d = []; //!< Lines of 1st BZ [::nbzl2d (max:26)][3]
1571+let bzl2d_proj = []; //!< Lines of 1st BZ [::nbzl2d (max:26)][3], projected into 2D plane
1572+/*
1573+ Equator
1574+*/
1575+let eqvec = [0.0, 0.0, 0.0]; //!< @f$k@f$-vector for equator
1576+let eqvec_fr = [0.0, 0.0, 0.0]; //!< @f$k@f$-vector for equator
1577+let nequator = []; //!< The number of equator
1578+let kveq = []; //!< @f$k@f$-vector of equator [::nb][::nequator][2][3]
1579+let kveq_rot = []; //!< @f$k@f$-vector of equator [::nb][::nequator*2*3]
1580+/*
1581+ Variables for mouse & cursorkey
1582+*/
1583+let sx = 0.0; //!< Scale of mouse movement
1584+let sy = 0.0; //!< Scale of mouse movement
1585+let cx = 0.0; //!< Starting point of drug
1586+let cy = 0.0; //!< Starting point of drug
1587+let scl = 1.0; //!< Initial scale
1588+let trans = [0.0, 0.0, 0.0]; //!< Translation
1589+let rot = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]; //!< Rotation matrix
1590+let thetax = 0.0; //!< Rotation angle
1591+let thetay = 0.0; //!< Rotation angle
1592+let thetaz = 0.0; //!< Rotation angle
1593+let linewidth = 3.0; //!< BZ/nodal-line/Fermiline width
1594+/*
1595+ Colors
1596+*/
1597+let black = [0.0, 0.0, 0.0, 1.0]; //!< Black color code
1598+let gray = [0.5, 0.5, 0.5, 1.0]; //!< Gray color code
1599+let wgray = [0.9, 0.9, 0.9, 1.0]; //!< Gray color code
1600+let bgray = [0.1, 0.1, 0.1, 1.0]; //!< Gray color code
1601+let white = [1.0, 1.0, 1.0, 1.0]; //!< White color code
1602+let cyan = [0.0, 1.0, 1.0, 1.0]; //!< Cyan color code
1603+let magenta = [1.0, 0.0, 1.0, 1.0]; //!< Magenta color code
1604+let yellow = [1.0, 1.0, 0.0, 1.0]; //!< Yellow color code
1605+let red = [1.0, 0.0, 0.0, 1.0]; //!< Red color code
1606+let green = [0.0, 1.0, 0.0, 1.0]; //!< Green color code
1607+let blue = [0.0, 0.0, 1.0, 1.0]; //!< Blue color code
1608+let BackGroundColor = [0.0, 0.0, 0.0, 1.0];//!< BackGround color code
1609+let LineColor = [1.0, 1.0, 1.0, 1.0];//!< Line color code
1610+let BarColor = [[0.0, 0.0, 1.0, 1.0], [0.0, 1.0, 1.0, 1.0],
1611+[0.0, 1.0, 0.0, 1.0], [1.0, 1.0, 0.0, 1.0], [1.0, 0.0, 0.0, 1.0]];
1612+/*
1613+ Others
1614+*/
1615+let corner = [[0.0, 0.0, 0.0, 0.0],
1616+ [0.0, 0.0, 0.0, 0.0],
1617+ [0.0, 0.0, 0.0, 0.0],
1618+ [0.0, 0.0, 0.0, 0.0],
1619+ [0.0, 0.0, 0.0, 0.0],
1620+ [0.0, 0.0, 0.0, 0.0]]; //!< Corners of tetrahedron
1621+let EF = 0.0; //!< Fermi energy
1622+/*
1623+Batch mode
1624+*/
1625+let batch_name;
1626+let frmsf_file_name;
1627+let lbatch = 0;
1628+
1629+wxTextCtrl* terminal;
1630+let refresh_interpol = 0;
1631+let refresh_patch = 1;
1632+let refresh_color = 1;
1633+let refresh_nodeline = 1;
1634+let refresh_equator = 1;
1635+let refresh_section = 1;
1636+let skip_minmax = 0;
1637+let windowx = 1100;
1638+let windowy = 850;
1639+/**
1640+ @brief Glut Display function
1641+ called by glutDisplayFunc
1642+*/
1643+function batch_draw()
1644+{
1645+ let iminmax = 0;
1646+ let minmax = [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]];
1647+
1648+ printf("\n Batch mode.\n");
1649+
1650+ read_batch();
1651+ refresh_patch_segment();
1652+}
1653+/**
1654+ @brief Main routine of FermiSurfer
1655+
1656+*/
1657+function OnInit()
1658+{
1659+ let ierr;
1660+
1661+ myf = new MyFrame(NULL, argv[1], wxDefaultPosition, wxSize(windowx, windowy));
1662+
1663+ terminal("\n");
1664+ terminal("##### Welocome to FermiSurfer ver. ") <<
1665+ wxT(VERSION) << wxT(" #####\n");
1666+ terminal("\n");
1667+ if (argc < 2) {
1668+ printf("\n");
1669+ printf(" Input file is not specified !\n");
1670+ printf(" Press any key to exit.\n");
1671+ ierr = getchar();
1672+ exit(-1);
1673+ }
1674+ /**/
1675+ terminal(" Initialize variables ...\n");
1676+ terminal("\n");
1677+ /*
1678+ Input from BXSF or FRMSF file
1679+ */
1680+ if (frmsf_file_name.AfterLast(wxUniChar('.')).CmpNoCase(wxT("bxsf")) == 0) {
1681+ read_bxsf();
1682+ }
1683+ else {
1684+ color_scale = read_file();
1685+ if (color_scale == 0)color_scale = 4;
1686+ }
1687+ /**/
1688+ interpol_energy();
1689+ init_corner();
1690+ bragg_vector();
1691+ modify_band();
1692+ /*
1693+ Brillouin zone
1694+ */
1695+ bz_lines();
1696+ calc_2dbz();
1697+ /**/
1698+ max_and_min_bz();
1699+ /**/
1700+ compute_patch_segment();
1701+ /*
1702+ Description
1703+ */
1704+ terminal("\n");
1705+ terminal(" ## How to handle ###################\n");
1706+ terminal("\n");
1707+ terminal(" mouse drag : Rotate objects\n");
1708+ terminal(" mousewheel : Resize objects\n");
1709+ terminal(" cursorkey or w,a,s,d : Move objects\n");
1710+ terminal("\n");
1711+ /**/
1712+ if (lbatch == 1) {
1713+ batch_draw();
1714+ }
1715+ return true;
1716+} /* main */
1717+
1718+function OnInitCmdLine(parser)
1719+{
1720+ OnInitCmdLine(parser);
1721+
1722+ parser.AddParam("FRMSF file to plot.",
1723+ wxCMD_LINE_VAL_STRING, wxCMD_LINE_OPTION_MANDATORY);
1724+ parser.AddParam("Batch file",
1725+ wxCMD_LINE_VAL_STRING, wxCMD_LINE_PARAM_OPTIONAL);
1726+ parser.AddParam("Window Size x",
1727+ wxCMD_LINE_VAL_STRING, wxCMD_LINE_PARAM_OPTIONAL);
1728+ parser.AddParam("Window Size y",
1729+ wxCMD_LINE_VAL_STRING, wxCMD_LINE_PARAM_OPTIONAL);
1730+}
1731+
1732+function OnCmdLineParsed(parser)
1733+{
1734+ if (parser.GetParamCount() > 0) {
1735+ frmsf_file_name = parser.GetParam(0);
1736+ if (parser.GetParamCount() > 1) {
1737+ batch_name = parser.GetParam(1);
1738+ lbatch = 1;
1739+ if (parser.GetParamCount() > 2) {
1740+ windowx = wxAtoi(parser.GetParam(2));
1741+ if (parser.GetParamCount() > 3) {
1742+ windowy = wxAtoi(parser.GetParam(3));
1743+ }
1744+ }
1745+ }
1746+ }
1747+
1748+ return OnCmdLineParsed(parser);
1749+}
1750+/**
1751+ @brief Free variables for patch before new patch is computed
1752+
1753+ Free : ::nmlp, ::matp, ::clr, ::kvp, ::nmlp_rot, ::kvp_rot,
1754+ ::kvnl, ::kvnl_rot, ::kv2d, ::clr2d
1755+*/
1756+function free_patch()
1757+{
1758+ let ib, i0, i1, i2;
1759+ /*
1760+ Fermi patch
1761+ */
1762+ if (refresh_patch == 1) {
1763+ for (ib = 0; ib < nb; ++ib) {
1764+ for (i0 = 0; i0 < ntri[ib]; ++i0) {
1765+ for (i1 = 0; i1 < 3; ++i1) {
1766+ for (i2 = 0; i2 < 2; ++i2)
1767+ arw[ib][i0][i1][i2] = [];
1768+ nmlp[ib][i0][i1] = [];
1769+ matp[ib][i0][i1] = [];
1770+ kvp[ib][i0][i1] = [];
1771+ arw[ib][i0][i1] = [];
1772+ }
1773+ nmlp[ib][i0] = [];
1774+ matp[ib][i0] = [];
1775+ kvp[ib][i0] = [];
1776+ arw[ib][i0] = [];
1777+ }
1778+ nmlp[ib] = [];
1779+ matp[ib] = [];
1780+ clr[ib] = [];
1781+ kvp[ib] = [];
1782+ arw[ib] = [];
1783+ nmlp_rot[ib] = [];
1784+ kvp_rot[ib] = [];
1785+ arw_rot[ib] = [];
1786+ }
1787+ nmlp = [];
1788+ matp = [];
1789+ clr = [];
1790+ kvp = [];
1791+ arw = [];
1792+ nmlp_rot = [];
1793+ kvp_rot = [];
1794+ arw_rot = [];
1795+ }/*if (refresh_patch == 1)*/
1796+ /*
1797+ Nodal line
1798+ */
1799+ if (refresh_nodeline == 1) {
1800+ for (ib = 0; ib < nb; ++ib) {
1801+ for (i0 = 0; i0 < nnl[ib]; ++i0) {
1802+ for (i1 = 0; i1 < 2; ++i1) {
1803+ kvnl[ib][i0][i1] = [];
1804+ }/*for (i1 = 0; i1 < 2; ++i1)*/
1805+ kvnl[ib][i0] = [];
1806+ }/*for (i0 = 0; i0 < nnl[ib]; ++i0)*/
1807+ kvnl[ib] = [];
1808+ kvnl_rot[ib] = [];
1809+ }/*for (ib = 0; ib < nb; ++ib)*/
1810+ kvnl = [];
1811+ kvnl_rot = [];
1812+ }/*if (refresh_nodeline == 1)*/
1813+ /*
1814+ 2D Fermi line
1815+ */
1816+ if (refresh_section == 1) {
1817+ for (ib = 0; ib < nb; ++ib) {
1818+ kv2d[ib] = [];
1819+ clr2d[ib] = [];
1820+ }/*for (ib = 0; ib < nb; ++ib)*/
1821+ kv2d = [];
1822+ clr2d = [];
1823+ }/*if (refresh_section == 1)*/
1824+ /*
1825+ equator
1826+ */
1827+ if (refresh_equator == 1) {
1828+ for (ib = 0; ib < nb; ++ib) {
1829+ for (i0 = 0; i0 < nequator[ib]; ++i0) {
1830+ for (i1 = 0; i1 < 2; ++i1) {
1831+ kveq[ib][i0][i1] = [];
1832+ }/*for (i1 = 0; i1 < 2; ++i1)*/
1833+ kveq[ib][i0] = [];
1834+ }/*for (i0 = 0; i0 < nequator[ib]; ++i0)*/
1835+ kveq[ib] = [];
1836+ kveq_rot[ib] = [];
1837+ }/*for (ib = 0; ib < nb; ++ib)*/
1838+ kveq = [];
1839+ kveq_rot = [];
1840+ }/*if (refresh_equator == 1)*/
1841+}/*function free_patch()*/
1842+/**
1843+ @brief Compute Max. & Min. of matrix elements.
1844+ Compute color of each patch
1845+
1846+ Modify : ::clr
1847+*/
1848+function max_and_min()
1849+{
1850+ let itri = 0;
1851+ let i = 0, ib = 0;
1852+ let abs = 0.0;
1853+
1854+ terminal("\n");
1855+ if (color_scale == 1) terminal(" ## Color Scale as Input Quantity (Real) #############\n");
1856+ else if (color_scale == 2) terminal(" ## Color Scale as Input Quantity (Complex) #############\n");
1857+ else if (color_scale == 3) terminal(" ## Color Scale as Fermi Velocity #############\n");
1858+ else if (color_scale == 4) terminal(" ## Color Scale as Band Index #############\n");
1859+ else if (color_scale == 5) terminal(" ## Gray Scale as Input Quantity (Real) #############\n");
1860+ else if (color_scale == 6) terminal(" ## Gray Scale as Fermi Velocity #############\n");
1861+ terminal("\n");
1862+ /*
1863+ Search max and min.
1864+ */
1865+ if (color_scale == 1 || color_scale == 6) {
1866+
1867+ patch_max = -1.0e10;
1868+ patch_min = 1.0e10;
1869+
1870+ for (ib = 0; ib < nb; ib++) {
1871+ for (itri = 0; itri < ntri[ib]; ++itri) {
1872+ for (i = 0; i < 3; ++i) {
1873+ if (matp[ib][itri][i][0] > patch_max) patch_max = matp[ib][itri][i][0];
1874+ if (matp[ib][itri][i][0] < patch_min) patch_min = matp[ib][itri][i][0];
1875+ }
1876+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
1877+ }/*for (ib = 0; ib < nb; ib++)*/
1878+ }/*if (color_scale == 0 || color_scale == 4)*/
1879+ else if (color_scale == 2) {
1880+
1881+ patch_max = -1.0e10;
1882+ patch_min = 1.0e10;
1883+
1884+ for (ib = 0; ib < nb; ib++) {
1885+ for (itri = 0; itri < ntri[ib]; ++itri) {
1886+ for (i = 0; i < 3; ++i) {
1887+ abs = Math.sqrt(matp[ib][itri][i][0] * matp[ib][itri][i][0]
1888+ + matp[ib][itri][i][1] * matp[ib][itri][i][1]);
1889+ if (abs > patch_max) patch_max = abs;
1890+ if (abs < patch_min) patch_min = abs;
1891+ }
1892+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
1893+ }/*for (ib = 0; ib < nb; ib++)*/
1894+ }/*if (color_scale == 2)*/
1895+ else if (color_scale == 3) {
1896+
1897+ patch_min = 1.0e10;
1898+ patch_max = -1.0e10;
1899+
1900+ for (ib = 0; ib < nb; ib++) {
1901+ for (itri = 0; itri < ntri[ib]; ++itri) {
1902+ for (i = 0; i < 3; ++i) {
1903+ abs = Math.sqrt(matp[ib][itri][i][0] * matp[ib][itri][i][0]
1904+ + matp[ib][itri][i][1] * matp[ib][itri][i][1]
1905+ + matp[ib][itri][i][2] * matp[ib][itri][i][2]);
1906+ if (abs > patch_max) patch_max = abs;
1907+ if (abs < patch_min) patch_min = abs;
1908+ }
1909+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
1910+ }/*for (ib = 0; ib < nb; ib++)*/
1911+ }/*if (color_scale == 3)*/
1912+ else if (color_scale == 4 || color_scale == 7) {
1913+
1914+ patch_max = -1.0e10;
1915+ patch_min = 1.0e10;
1916+
1917+ for (ib = 0; ib < nb; ib++) {
1918+ for (itri = 0; itri < ntri[ib]; ++itri) {
1919+ for (i = 0; i < 3; ++i) {
1920+ norm = 0.0;
1921+ for (j = 0; j < 3; ++j) norm += nmlp[ib][itri][i][j]*nmlp[ib][itri][i][j];
1922+ norm = Math.sqrt(norm);
1923+
1924+ if (norm > patch_max) patch_max = norm;
1925+ if (norm < patch_min) patch_min = norm;
1926+ }
1927+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
1928+ }/*for (ib = 0; ib < nb; ib++)*/
1929+ }/*if (color_scale == 5 || color_scale == 6)*/
1930+
1931+ myf->textbox_min->ChangeValue(wxString::Format(wxT("%f"), patch_min));
1932+ myf->textbox_max->ChangeValue(wxString::Format(wxT("%f"), patch_max));
1933+}/* max_and_min */
1934+ /**
1935+ @brief Compute Max. & Min. of matrix elements.
1936+ Compute color of each patch
1937+
1938+ Modify : ::clr
1939+ */
1940+function paint()
1941+{
1942+ let itri, j;
1943+ let origin = [0.0, 0.0, 0.0, 0.0];
1944+ let i, ib;
1945+ let mat2;
1946+
1947+ if (color_scale == 1) {
1948+
1949+ for (ib = 0; ib < nb; ib++) {
1950+
1951+ for (itri = 0; itri < ntri[ib]; ++itri) {
1952+ for (i = 0; i < 3; ++i) {
1953+ /**/
1954+ mat2 = (matp[ib][itri][i][0] - patch_min) / (patch_max - patch_min);
1955+ mat2 = mat2 * 4.0;
1956+ /**/
1957+ if (mat2 <= 1.0) {
1958+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
1959+ = BarColor[1][j] * mat2 + BarColor[0][j] * (1.0 - mat2);
1960+ }
1961+ else if (mat2 <= 2.0) {
1962+ mat2 = mat2 - 1.0;
1963+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
1964+ = BarColor[2][j] * mat2 + BarColor[1][j] * (1.0 - mat2);
1965+ }
1966+ else if (mat2 <= 3.0) {
1967+ mat2 = mat2 - 2.0;
1968+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
1969+ = BarColor[3][j] * mat2 + BarColor[2][j] * (1.0 - mat2);
1970+ }
1971+ else {
1972+ mat2 = mat2 - 3.0;
1973+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
1974+ = BarColor[4][j] * mat2 + BarColor[3][j] * (1.0 - mat2);
1975+ }
1976+ }/*for (i = 0; i < 3; ++i)*/
1977+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
1978+ }/*for (ib = 0; ib < nb; ib++)*/
1979+ }/*if (color_scale == 1 || color_scale == 2)*/
1980+ else if (color_scale == 2) {
1981+
1982+ for (j = 0; j < 4; ++j) origin[j] = 1.0 - BackGroundColor[j];
1983+
1984+ let i, ib;
1985+ let theta, abs, theta2;
1986+
1987+ for (ib = 0; ib < nb; ib++) {
1988+ for (itri = 0; itri < ntri[ib]; ++itri) {
1989+ for (i = 0; i < 3; ++i) {
1990+ /**/
1991+ abs = Math.sqrt(matp[ib][itri][i][0] * matp[ib][itri][i][0]
1992+ + matp[ib][itri][i][1] * matp[ib][itri][i][1]);
1993+ if (abs / patch_max < 0.00001) theta = 0.0;
1994+ else if (matp[ib][itri][i][1] > 0.0) theta = Math.acos(matp[ib][itri][i][0] / abs);
1995+ else theta = -Math.acos(matp[ib][itri][i][0] / abs);
1996+ abs /= patch_max;
1997+ theta = 3.0 * theta / Math.acos(-1.0);
1998+ /**/
1999+ if (-3.0 <= theta && theta < -2.0) {
2000+ theta2 = theta + 3.0;
2001+ for (j = 0; j < 4; ++j)
2002+ clr[ib][j + 4 * i + 12 * itri] = blue[j] * theta2 + cyan[j] * (1.0 - theta2);
2003+ }
2004+ else if (-2.0 <= theta && theta < -1.0) {
2005+ theta2 = theta + 2.0;
2006+ for (j = 0; j < 4; ++j)
2007+ clr[ib][j + 4 * i + 12 * itri] = magenta[j] * theta2 + blue[j] * (1.0 - theta2);
2008+ }
2009+ else if (-1.0 <= theta && theta < 0.0) {
2010+ theta2 = theta + 1.0;
2011+ for (j = 0; j < 4; ++j)
2012+ clr[ib][j + 4 * i + 12 * itri] = red[j] * theta2 + magenta[j] * (1.0 - theta2);
2013+ }
2014+ else if (0.0 <= theta && theta < 1.0) {
2015+ theta2 = theta;
2016+ for (j = 0; j < 4; ++j)
2017+ clr[ib][j + 4 * i + 12 * itri] = yellow[j] * theta2 + red[j] * (1.0 - theta2);
2018+ }
2019+ else if (1.0 <= theta && theta < 2.0) {
2020+ theta2 = theta - 1.0;
2021+ for (j = 0; j < 4; ++j)
2022+ clr[ib][j + 4 * i + 12 * itri] = green[j] * theta2 + yellow[j] * (1.0 - theta2);
2023+ }
2024+ else {
2025+ theta2 = theta - 2.0;
2026+ for (j = 0; j < 4; ++j)
2027+ clr[ib][j + 4 * i + 12 * itri] = cyan[j] * theta2 + green[j] * (1.0 - theta2);
2028+ }
2029+ clr[ib][j + 4 * i + 12 * itri] = clr[ib][j + 4 * i + 12 * itri] * abs + origin[j] * (1.0 - abs);
2030+
2031+ }/*for (i = 0; i < 3; ++i)*/
2032+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
2033+ }/*for (ib = 0; ib < nb; ib++)*/
2034+ }/*if (color_scale == 2)*/
2035+ else if (color_scale == 4) {
2036+ let i, ib;
2037+ let mat2;
2038+
2039+ for (ib = 0; ib < nb; ib++) {
2040+
2041+ for (itri = 0; itri < ntri[ib]; ++itri) {
2042+ for (i = 0; i < 3; ++i) {
2043+ /**/
2044+ mat2 = 0.0;
2045+ for (j = 0; j < 3; ++j) mat2 += nmlp[ib][itri][i][j] * nmlp[ib][itri][i][j];
2046+ mat2 = Math.sqrt(mat2);
2047+ mat2 = (mat2 - patch_min) / (patch_max - patch_min);
2048+ mat2 = mat2 * 4.0;
2049+ /**/
2050+ if (mat2 <= 1.0) {
2051+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
2052+ = BarColor[1][j] * mat2 + BarColor[0][j] * (1.0 - mat2);
2053+ }
2054+ else if (mat2 <= 2.0) {
2055+ mat2 = mat2 - 1.0;
2056+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
2057+ = BarColor[2][j] * mat2 + BarColor[1][j] * (1.0 - mat2);
2058+ }
2059+ else if (mat2 <= 3.0) {
2060+ mat2 = mat2 - 2.0;
2061+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
2062+ = BarColor[3][j] * mat2 + BarColor[2][j] * (1.0 - mat2);
2063+ }
2064+ else {
2065+ mat2 = mat2 - 3.0;
2066+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
2067+ = BarColor[4][j] * mat2 + BarColor[3][j] * (1.0 - mat2);
2068+ }
2069+ }/*for (i = 0; i < 3; ++i)*/
2070+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
2071+ }/*for (ib = 0; ib < nb; ib++)*/
2072+ }/*if (color_scale == 4)*/
2073+ else if (color_scale == 3 || color_scale == 5) {
2074+ let i, ib;
2075+ let mat2;
2076+
2077+ for (ib = 0; ib < nb; ib++) {
2078+ /**/
2079+ if (nb == 1) mat2 = 0.5;
2080+ else mat2 = 1.0 / (nb - 1) * ib;
2081+ mat2 *= 4.0;
2082+ /**/
2083+ if (mat2 <= 1.0) {
2084+
2085+ for (itri = 0; itri < ntri[ib]; ++itri) {
2086+ for (i = 0; i < 3; ++i) {
2087+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
2088+ = BarColor[1][j] * mat2 + BarColor[0][j] * (1.0 - mat2);
2089+ }
2090+ }
2091+ }
2092+ else if (mat2 <= 2.0) {
2093+ mat2 = mat2 - 1.0;
2094+
2095+ for (itri = 0; itri < ntri[ib]; ++itri) {
2096+ for (i = 0; i < 3; ++i) {
2097+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
2098+ = BarColor[2][j] * mat2 + BarColor[1][j] * (1.0 - mat2);
2099+ }
2100+ }
2101+ }
2102+ else if (mat2 <= 3.0) {
2103+ mat2 = mat2 - 2.0;
2104+
2105+ for (itri = 0; itri < ntri[ib]; ++itri) {
2106+ for (i = 0; i < 3; ++i) {
2107+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
2108+ = BarColor[3][j] * mat2 + BarColor[2][j] * (1.0 - mat2);
2109+ }
2110+ }
2111+ }
2112+ else {
2113+ mat2 = mat2 - 3.0;
2114+
2115+ for (itri = 0; itri < ntri[ib]; ++itri) {
2116+ for (i = 0; i < 3; ++i) {
2117+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri]
2118+ = BarColor[4][j] * mat2 + BarColor[3][j] * (1.0 - mat2);
2119+ }
2120+ }
2121+ }
2122+ }/*for (ib = 0; ib < nb; ib++*/
2123+
2124+ if (color_scale == 3) {
2125+ for (ib = 0; ib < nb; ib++) {
2126+
2127+ for (itri = 0; itri < ntri[ib]; ++itri) {
2128+ for (i = 0; i < 3; ++i) {
2129+ for (j = 0; j < 3; ++j) {
2130+ arw[ib][itri][i][0][j] = kvp[ib][itri][i][j]
2131+ - 0.5 * matp[ib][itri][i][j] / patch_max;
2132+ arw[ib][itri][i][1][j] = kvp[ib][itri][i][j]
2133+ + 0.5 * matp[ib][itri][i][j] / patch_max;
2134+ }/*for (j = 0; j < 3; ++j)*/
2135+ }/*for (i = 0; i < 3; ++i)*/
2136+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
2137+ }/*for (ib = 0; ib < nb; ib++)*/
2138+ }/*if (color_scale == 3)*/
2139+ }/*if (color_scale == 5)*/
2140+ else if (color_scale == 6) {
2141+ let i, ib;
2142+ let mat2;
2143+
2144+ for (ib = 0; ib < nb; ib++) {
2145+
2146+ for (itri = 0; itri < ntri[ib]; ++itri) {
2147+ for (i = 0; i < 3; ++i) {
2148+ /**/
2149+ mat2 = (matp[ib][itri][i][0] - patch_min) / (patch_max - patch_min);
2150+ /**/
2151+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = wgray[j] * mat2 + bgray[j] * (1.0 - mat2);
2152+ }/*for (i = 0; i < 3; ++i)*/
2153+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
2154+ }/*for (ib = 0; ib < nb; ib++)*/
2155+ }/*if (color_scale == 6)*/
2156+ else if (color_scale == 7) {
2157+ let i, ib;
2158+ let mat2;
2159+
2160+ for (ib = 0; ib < nb; ib++) {
2161+
2162+ for (itri = 0; itri < ntri[ib]; ++itri) {
2163+ for (i = 0; i < 3; ++i) {
2164+ /**/
2165+ mat2 = 0.0;
2166+ for (j = 0; j < 3; ++j) mat2 += nmlp[ib][itri][i][j] * nmlp[ib][itri][i][j];
2167+ mat2 = Math.sqrt(mat2);
2168+ mat2 = (mat2 - patch_min) / (patch_max - patch_min);
2169+ /**/
2170+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = wgray[j] * mat2 + bgray[j] * (1.0 - mat2);
2171+ }/*for (i = 0; i < 3; ++i)*/
2172+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
2173+ }/*for (ib = 0; ib < nb; ib++)*/
2174+ }/*if (color_scale == 7)*/
2175+}/* paint */
2176+/**
2177+ @brief Specify corners of tetrahedron
2178+
2179+ Modify : ::corner
2180+*/
2181+function init_corner()
2182+{
2183+ let i = 0, j = 0;
2184+ let corner1 = [
2185+ /*
2186+ [0] min = 0-7
2187+ */
2188+ [ [ 0, 1, 3, 7 ],
2189+ [ 0, 1, 5, 7 ],
2190+ [ 0, 2, 3, 7 ],
2191+ [ 0, 2, 6, 7 ],
2192+ [ 0, 4, 5, 7 ],
2193+ [ 0, 4, 6, 7 ] ],
2194+ /*
2195+ [1] min = 1-6
2196+ */
2197+ [ [ 1, 0, 2, 6 ],
2198+ [ 1, 0, 4, 6 ],
2199+ [ 1, 3, 2, 6 ],
2200+ [ 1, 3, 7, 6 ],
2201+ [ 1, 5, 4, 6 ],
2202+ [ 1, 5, 7, 6 ] ],
2203+ /*
2204+ [2] min = 2-5
2205+ */
2206+ [ [ 2, 0, 1, 5 ],
2207+ [ 2, 0, 4, 5 ],
2208+ [ 2, 3, 1, 5 ],
2209+ [ 2, 3, 7, 5 ],
2210+ [ 2, 6, 4, 5 ],
2211+ [ 2, 6, 7, 5 ] ],
2212+ /*
2213+ [3] min = 3-4
2214+ */
2215+ [ [ 4, 0, 1, 3 ],
2216+ [ 4, 0, 2, 3 ],
2217+ [ 4, 5, 1, 3 ],
2218+ [ 4, 5, 7, 3 ],
2219+ [ 4, 6, 2, 3 ],
2220+ [ 4, 6, 7, 3 ] ],
2221+ /*
2222+ [4] min = 0-7, max = 3-4
2223+ */
2224+ [ [ 0, 4, 5, 6 ],
2225+ [ 1, 2, 3, 7 ],
2226+ [ 0, 7, 2, 6 ],
2227+ [ 0, 7, 1, 5 ],
2228+ [ 0, 7, 1, 2 ],
2229+ [ 0, 7, 6, 5 ] ],
2230+ /*
2231+ [5] min = 1-6, max = 3-4
2232+ */
2233+ [ [ 0, 4, 5, 6 ],
2234+ [ 1, 2, 3, 7 ],
2235+ [ 1, 6, 5, 7 ],
2236+ [ 1, 6, 7, 2 ],
2237+ [ 1, 6, 2, 0 ],
2238+ [ 1, 6, 0, 5 ] ],
2239+ /*
2240+ [6] min = 2-5, max = 3-4
2241+ */
2242+ [ [ 0, 4, 5, 6 ],
2243+ [ 1, 2, 3, 7 ],
2244+ [ 2, 5, 7, 6 ],
2245+ [ 2, 5, 6, 0 ],
2246+ [ 2, 5, 0, 1 ],
2247+ [ 2, 5, 1, 7 ] ],
2248+ /*
2249+ [7] min = 3-4, max = 0-7
2250+ */
2251+ [ [ 0, 1, 2, 4 ],
2252+ [ 7, 3, 5, 6 ],
2253+ [ 3, 4, 1, 5 ],
2254+ [ 3, 4, 5, 6 ],
2255+ [ 3, 4, 6, 2 ],
2256+ [ 3, 4, 2, 1 ] ],
2257+ /*
2258+ [8] min = 2-5, max = 0-7
2259+ */
2260+ [ [ 0, 1, 2, 4 ],
2261+ [ 7, 3, 5, 6 ],
2262+ [ 2, 5, 1, 3 ],
2263+ [ 2, 5, 3, 6 ],
2264+ [ 2, 5, 6, 4 ],
2265+ [ 2, 5, 4, 1 ] ],
2266+ /*
2267+ [9] min = 1-6, max = 0-7
2268+ */
2269+ [ [ 0, 1, 2, 4 ],
2270+ [ 7, 3, 5, 6 ],
2271+ [ 1, 6, 2, 3 ],
2272+ [ 1, 6, 3, 5 ],
2273+ [ 1, 6, 5, 4 ],
2274+ [ 1, 6, 4, 2 ] ],
2275+ /*
2276+ [10] min = 0-7, max = 1-6
2277+ */
2278+ [ [ 1, 0, 3, 5 ],
2279+ [ 6, 2, 4, 7 ],
2280+ [ 0, 7, 2, 3 ],
2281+ [ 0, 7, 3, 5 ],
2282+ [ 0, 7, 5, 4 ],
2283+ [ 0, 7, 4, 2 ] ],
2284+ /*
2285+ [11] min = 3-4, max = 1-6
2286+ */
2287+ [ [ 1, 0, 3, 5 ],
2288+ [ 6, 2, 4, 7 ],
2289+ [ 3, 4, 0, 2 ],
2290+ [ 3, 4, 2, 7 ],
2291+ [ 3, 4, 7, 5 ],
2292+ [ 3, 4, 5, 0 ] ],
2293+ /*
2294+ [12] min = 2-5, max = 1-6
2295+ */
2296+ [ [ 1, 0, 3, 5 ],
2297+ [ 6, 2, 4, 7 ],
2298+ [ 2, 5, 0, 3 ],
2299+ [ 2, 5, 3, 7 ],
2300+ [ 2, 5, 7, 4 ],
2301+ [ 2, 5, 4, 0 ] ],
2302+ /*
2303+ [13] min = 0-7, max = 2-5
2304+ */
2305+ [ [ 2, 0, 3, 6 ],
2306+ [ 5, 1, 4, 7 ],
2307+ [ 0, 7, 1, 3 ],
2308+ [ 0, 7, 3, 6 ],
2309+ [ 0, 7, 6, 4 ],
2310+ [ 0, 7, 4, 1 ] ],
2311+ /*
2312+ [14] min = 1-6, max = 2-5
2313+ */
2314+ [ [ 2, 0, 3, 6 ],
2315+ [ 5, 1, 4, 7 ],
2316+ [ 1, 6, 0, 3 ],
2317+ [ 1, 6, 3, 7 ],
2318+ [ 1, 6, 7, 4 ],
2319+ [ 1, 6, 4, 0 ] ],
2320+ /*
2321+ [15] min = 3-4, max = 2-5
2322+ */
2323+ [ [ 2, 0, 3, 6 ],
2324+ [ 5, 1, 4, 7 ],
2325+ [ 3, 4, 0, 1 ],
2326+ [ 3, 4, 1, 7 ],
2327+ [ 3, 4, 7, 6 ],
2328+ [ 3, 4, 6, 0 ] ],
2329+ ];
2330+ /**/
2331+ for (i = 0; i < 6; ++i) {
2332+ for (j = 0; j < 4; ++j) {
2333+ corner[i][j] = corner1[itet][i][j];
2334+ }
2335+ }
2336+}
2337+/**
2338+ @brief Compute Bragg vetor
2339+
2340+ Modify : ::bragg, ::brnrm
2341+*/
2342+function bragg_vector()
2343+{
2344+ let i0, i1, i2, i, ibr;
2345+ /**/
2346+ ibr = 0;
2347+ /**/
2348+ for (i0 = -1; i0 <= 1; ++i0) {
2349+ for (i1 = -1; i1 <= 1; ++i1) {
2350+ for (i2 = -1; i2 <= 1; ++i2) {
2351+ /*
2352+ Excepte Gamma points
2353+ */
2354+ if (i0 == 0 && i1 == 0 && i2 == 0) continue;
2355+ /*
2356+ Fractional -> Cartecian
2357+ */
2358+ for (i = 0; i < 3; ++i)
2359+ bragg[ibr][i] = (i0 * bvec[0][i]
2360+ + i1 * bvec[1][i]
2361+ + i2 * bvec[2][i]) * 0.5;
2362+ /*
2363+ And its norm
2364+ */
2365+ brnrm[ibr] = bragg[ibr][0] * bragg[ibr][0]
2366+ + bragg[ibr][1] * bragg[ibr][1]
2367+ + bragg[ibr][2] * bragg[ibr][2];
2368+ /**/
2369+ ibr = ibr + 1;
2370+ }/*for (i2 = -1; i2 <= 1; ++i2)*/
2371+ }/*for (i1 = -1; i1 <= 1; ++i1)*/
2372+ }/*for (i0 = -1; i0 <= 1; ++i0)*/
2373+ /*
2374+ Search min. of brnrm
2375+ */
2376+ brnrm_min = brnrm[0];
2377+ for (ibr = 1; ibr < 26; ibr++) {
2378+ if (brnrm_min > brnrm[ibr]) brnrm_min = brnrm[ibr];
2379+ }
2380+ terminal(" Minimum Bragg norm : " + toString(brnrm_min) + "\n");
2381+}/* bragg_vector */
2382+/**
2383+ @brief Print max and minimum @f$\varepsilon_{n k}, \Delta_{n k}@f$
2384+ in the whole Brillouine zone
2385+*/
2386+function max_and_min_bz()
2387+{
2388+ let ib, i0, i1, i2;
2389+ let eigmin, eigmax, matmin, matmax;
2390+ /**/
2391+ terminal("\n");
2392+ terminal(" ## Max. and Min. of each bands #######################\n");
2393+ terminal("\n");
2394+ terminal(" Band Eig_Min. Eig_Max Mat_Min Mat_Max\n");
2395+ for (ib = 0; ib < nb; ib++) {
2396+ eigmax = eig0[0][0][0][0];
2397+ eigmin = eig0[0][0][0][0];
2398+ matmax = mat0[0][0][0][0][0];
2399+ matmin = mat0[0][0][0][0][0];
2400+ for (i0 = 0; i0 < ng0[0]; ++i0) {
2401+ for (i1 = 0; i1 < ng0[1]; ++i1) {
2402+ for (i2 = 0; i2 < ng0[2]; ++i2) {
2403+ if (eig0[ib][i0][i1][i2] > eigmax) eigmax = eig0[ib][i0][i1][i2];
2404+ if (eig0[ib][i0][i1][i2] < eigmin) eigmin = eig0[ib][i0][i1][i2];
2405+ if (mat0[ib][i0][i1][i2][0] > matmax) matmax = mat0[ib][i0][i1][i2][0];
2406+ if (mat0[ib][i0][i1][i2][0] < matmin) matmin = mat0[ib][i0][i1][i2][0];
2407+ }/*for (i2 = 0; i2 < ng0[2]; ++i2)*/
2408+ }/*for (i1 = 0; i1 < ng0[1]; ++i1)*/
2409+ }/*for (i0 = 0; i0 < ng0[0]; ++i0)*/
2410+ terminal(" " + toString(ib + 1) + " " + toString(eigmin) + " " +
2411+ toString(eigmax) + " " + toString(matmin) + " " + toString(matmax) + "\n");
2412+ }/*for (ib = 0; ib < nb; ib++)*/
2413+}/* max_and_min_bz */
2414+/**
2415+ @brief Compute coefficient for the French-curve (Kumo) interpolation
2416+ @f[
2417+ A^{\rm intp} = \sum_{i = 1}^4 C_i A_i^{\rm orig}
2418+ @f]
2419+*/
2420+function kumo_coef(
2421+ j = 0, //!< [in] Interpolated grid index
2422+ coef = [0.0, 0.0, 0.0, 0.0],//!< [out] Coefficient of interpolation @f$C_i@f$
2423+ interpol = 0
2424+) {
2425+ let x = 0.0, mx = 0.0;
2426+ x = j / interpol;
2427+ mx = 1.0 - x;
2428+ coef[0] = -0.5 * x * mx * mx;
2429+ coef[1] = mx * (mx*mx + 3.0* x*mx + 0.5* x* x);
2430+ coef[2] = x * ( x* x + 3.0*mx* x + 0.5*mx*mx);
2431+ coef[3] = -0.5 * x * x * mx;
2432+}
2433+/**
2434+ @brief Interpolation of energy and matrix
2435+ with the French-curve (Kumo) interpolation.
2436+
2437+ Modify : ::eig, ::mat
2438+*/
2439+function interpol_energy()
2440+{
2441+ let ib, i0, i1, i2, ii;
2442+
2443+ terminal(" Interpolating ... ");
2444+ /*
2445+ Reallocate
2446+ */
2447+ for (ib = 0; ib < nb; ib++) {
2448+ for (i0 = 0; i0 < ng[0]; i0++) {
2449+ for (i1 = 0; i1 < ng[1]; i1++) {
2450+ for (i2 = 0; i2 < ng[2]; i2++) {
2451+ vf[ib][i0][i1][i2] = [];
2452+ mat[ib][i0][i1][i2] = [];
2453+ }
2454+ eig[ib][i0][i1] = [];
2455+ mat[ib][i0][i1] = [];
2456+ vf[ib][i0][i1] = [];
2457+ }/*for (i1 = 0; i1 < ng[1]; i1++)*/
2458+ eig[ib][i0] = [];
2459+ mat[ib][i0] = [];
2460+ vf[ib][i0] = [];
2461+ }/*for (i0 = 0; i0 < ng[0]; i0++)*/
2462+ eig[ib] = [];
2463+ mat[ib] = [];
2464+ vf[ib] = [];
2465+ }/*for (ib = 0; ib < nb; ib++)*/
2466+ eig = [];
2467+ mat = [];
2468+ vf = [];
2469+ for (ii = 0; ii < 3; ii++)ng[ii] = ng0[ii] * interpol;
2470+ /**/
2471+ for (ib = 0; ib < nb; ib++) {
2472+ eig[ib] = new let**[ng[0]];
2473+ mat[ib] = new let***[ng[0]];
2474+ vf[ib] = new let***[ng[0]];
2475+ for (i0 = 0; i0 < ng[0]; i0++) {
2476+ eig[ib][i0] = new let*[ng[1]];
2477+ mat[ib][i0] = new let**[ng[1]];
2478+ vf[ib][i0] = new let**[ng[1]];
2479+ for (i1 = 0; i1 < ng[1]; i1++) {
2480+ eig[ib][i0][i1] = new let[ng[2]];
2481+ mat[ib][i0][i1] = new let*[ng[2]];
2482+ vf[ib][i0][i1] = new let*[ng[2]];
2483+ for (i2 = 0; i2 < ng[2]; i2++) {
2484+ mat[ib][i0][i1][i2] = new let[3];
2485+ vf[ib][i0][i1][i2] = new let[3];
2486+ }
2487+ }/*for (i1 = 0; i1 < ng[1]; i1++)*/
2488+ }/*for (i0 = 0; i0 < ng[0]; i0++)*/
2489+ }/*for (ib = 0; ib < nb; ib++)*/
2490+ /*
2491+ 3rd order - three dimensional Kumo interpolation
2492+ */
2493+ let j0, j1, j2, jj;
2494+ let coef[4],
2495+ mat1[4][4][4][3], mat2[4][4][3], mat3[4][3],
2496+ eig1[4][4][4], eig2[4][4], eig3[4];
2497+
2498+ for (ib = 0; ib < nb; ib++) {
2499+ for (i0 = 0; i0 < ng0[0]; i0++) {
2500+ //if (ith == 1) continue;
2501+ for (i1 = 0; i1 < ng0[1]; i1++) {
2502+ for (i2 = 0; i2 < ng0[2]; i2++) {
2503+ for (j0 = 0; j0 < 4; j0++) {
2504+ for (j1 = 0; j1 < 4; j1++) {
2505+ for (j2 = 0; j2 < 4; j2++) {
2506+ eig1[j0][j1][j2] = eig0[ib][modulo(i0 + j0 - 1, ng0[0])]
2507+ [modulo(i1 + j1 - 1, ng0[1])]
2508+ [modulo(i2 + j2 - 1, ng0[2])];
2509+ for (jj = 0; jj < 3; jj++) {
2510+ mat1[j0][j1][j2][jj] = mat0[ib][modulo(i0 + j0 - 1, ng0[0])]
2511+ [modulo(i1 + j1 - 1, ng0[1])]
2512+ [modulo(i2 + j2 - 1, ng0[2])][jj];
2513+ }
2514+ }/*for (j2 = 0; j2 < 4; j2++)*/
2515+ }/*for (j1 = 0; j1 < 4; j1++)*/
2516+ }/*for (i2 = 0; i2 < ng0[2]; i2++)*/
2517+ for (j0 = 0; j0 < interpol; j0++) {
2518+ kumo_coef(j0, &coef[0], interpol);
2519+ for (j1 = 0; j1 < 4; j1++) {
2520+ for (j2 = 0; j2 < 4; j2++) {
2521+ eig2[j1][j2] = 0.0;
2522+ for (jj = 0; jj < 3; jj++) mat2[j1][j2][jj] = 0.0;
2523+ for (ii = 0; ii < 4; ii++) {
2524+ eig2[j1][j2] += coef[ii] * eig1[ii][j1][j2];
2525+ for (jj = 0; jj < 3; jj++)
2526+ mat2[j1][j2][jj] += coef[ii] * mat1[ii][j1][j2][jj];
2527+ }/*for (ii = 0; ii < 4; ii++)*/
2528+ }/*for (j2 = 0; j2 < 4; j2++)*/
2529+ }/*for (j1 = 0; j1 < 4; j1++)*/
2530+ for (j1 = 0; j1 < interpol; j1++) {
2531+ kumo_coef(j1, &coef[0], interpol);
2532+ for (j2 = 0; j2 < 4; j2++) {
2533+ eig3[j2] = 0.0;
2534+ for (jj = 0; jj < 3; jj++) mat3[j2][jj] = 0.0;
2535+ for (ii = 0; ii < 4; ii++) {
2536+ eig3[j2] += coef[ii] * eig2[ii][j2];
2537+ for (jj = 0; jj < 3; jj++)
2538+ mat3[j2][jj] += coef[ii] * mat2[ii][j2][jj];
2539+ }/*for (ii = 0; ii < 4; ii++)*/
2540+ }/*for (j2 = 0; j2 < 4; j2++)*/
2541+ for (j2 = 0; j2 < interpol; j2++) {
2542+ kumo_coef(j2, &coef[0], interpol);
2543+ eig[ib][i0*interpol + j0]
2544+ [i1*interpol + j1]
2545+ [i2*interpol + j2] = 0.0;
2546+ for (jj = 0; jj < 3; jj++)
2547+ mat[ib][i0*interpol + j0]
2548+ [i1*interpol + j1]
2549+ [i2*interpol + j2][jj] = 0.0;
2550+ for (ii = 0; ii < 4; ii++) {
2551+ eig[ib][i0*interpol + j0]
2552+ [i1*interpol + j1]
2553+ [i2*interpol + j2] += coef[ii] * eig3[ii];
2554+ for (jj = 0; jj < 3; jj++)
2555+ mat[ib][i0*interpol + j0]
2556+ [i1*interpol + j1]
2557+ [i2*interpol + j2][jj] += coef[ii] * mat3[ii][jj];
2558+ }/*for (ii = 0; ii < 4; ii++)*/
2559+ }/*for (j2 = 0; j2 < interpol; j2++)*/
2560+ }/*for (j1 = 0; j1 < interpol; j1++)*/
2561+ }/*for (j0 = 0; j0 < interpol; j0++)*/
2562+ }/*for (i2 = 0; i2 < ng0[2]; i2++)*/
2563+ }/*for (i1 = 0; i1 < ng0[1]; i1++)*/
2564+ }/*for (i0 = 0; i0 < ng0[0]; i0++)*/
2565+ }/*for (ib = 0; ib < nb; ib++)*/
2566+ /*
2567+ Fermi velocity
2568+ */
2569+#pragma omp parallel default(none) \
2570+ shared(nb,ng,eig,vf,avec) \
2571+ private (ib,i0,i1,i2,ii)
2572+ {
2573+ let i0p, i0m, i1p, i1m, i2p, i2m;
2574+ let de[3];
2575+
2576+ for (ib = 0; ib < nb; ib++) {
2577+ for (i0 = 0; i0 < ng[0]; i0++) {
2578+ i0p = modulo(i0 + 1, ng[0]);
2579+ i0m = modulo(i0 - 1, ng[0]);
2580+ for (i1 = 0; i1 < ng[1]; i1++) {
2581+ i1p= modulo(i1 + 1, ng[1]);
2582+ i1m = modulo(i1 - 1, ng[1]);
2583+ for (i2 = 0; i2 < ng[2]; i2++) {
2584+ i2p = modulo(i2 + 1, ng[2]);
2585+ i2m = modulo(i2 - 1, ng[2]);
2586+
2587+ de[0] = eig[ib][i0p][i1][i2] - eig[ib][i0m][i1][i2];
2588+ de[1] = eig[ib][i0][i1p][i2] - eig[ib][i0][i1m][i2];
2589+ de[2] = eig[ib][i0][i1][i2p] - eig[ib][i0][i1][i2m];
2590+ for (ii = 0; ii < 3; ii++)de[ii] *= 0.5f * ng[ii];
2591+ for (ii = 0; ii < 3; ii++) vf[ib][i0][i1][i2][ii] =
2592+ avec[0][ii] * de[0] + avec[1][ii] * de[1] + avec[2][ii] * de[2];
2593+
2594+ }/*for (i2 = 0; i2 < ng[2]; i2++)*/
2595+ }/*for (i1 = 0; i1 < ng[1]; i1++)*/
2596+ }/*for (i0 = 0; i0 < ng[0]; i0++)*/
2597+ }/*for (ib = 0; ib < nb; ib++)*/
2598+ }/*End of parallel region*/
2599+ terminal("Done\n\n");
2600+}/*function interpol_energy() */
2601+
2602+function compute_patch_segment() {
2603+ if (refresh_interpol == 1){
2604+ interpol_energy();
2605+ refresh_patch = 1;
2606+ refresh_interpol = 0;
2607+ }
2608+ if (refresh_patch == 1) {
2609+ fermi_patch();
2610+ refresh_color = 1;
2611+ refresh_section = 1;
2612+ refresh_equator = 1;
2613+ refresh_nodeline = 1;
2614+ refresh_patch = 0;
2615+ }
2616+ if (refresh_color == 1) {
2617+ if (skip_minmax == 1) skip_minmax = 0;
2618+ else max_and_min();
2619+ paint();
2620+ refresh_section = 1;
2621+ refresh_color = 0;
2622+ }
2623+ if (refresh_nodeline == 1) {
2624+ calc_nodeline();
2625+ refresh_nodeline = 0;
2626+ }
2627+ if (refresh_section == 1) {
2628+ calc_2dbz();
2629+ calc_section();
2630+ refresh_section = 0;
2631+ }
2632+ if (refresh_equator == 1) {
2633+ equator();
2634+ refresh_equator = 0;
2635+ }
2636+}
2637+
2638+function refresh_patch_segment() {
2639+ free_patch();
2640+ compute_patch_segment();
2641+}
2642+
2643+enum
2644+{
2645+ ibutton_reflesh,
2646+ ibutton_compute,
2647+ iradio_background,
2648+ iradio_brillouinzone,
2649+ iradio_lighting,
2650+ iradio_mouse,
2651+ iradio_stereo,
2652+ iradio_tetra,
2653+ iradio_colorscale,
2654+ itext_colorscalemin,
2655+ itext_colorscalemax,
2656+ icheck_section,
2657+ icheck_gamma,
2658+ itext_sectionx,
2659+ itext_sectiony,
2660+ itext_sectionz,
2661+ icheck_nodeline,
2662+ icheck_colorbar,
2663+ icheck_equator,
2664+ itext_equatorx,
2665+ itext_equatory,
2666+ itext_equatorz,
2667+ itext_line,
2668+ itext_shift,
2669+ itext_interpol,
2670+ itext_scale,
2671+ itext_positionx,
2672+ itext_positiony,
2673+ itext_rotx,
2674+ itext_roty,
2675+ itext_rotz,
2676+ ibutton_rotate,
2677+ icheck_band,
2678+ itext_BackGroundR,
2679+ itext_BackGroundG,
2680+ itext_BackGroundB,
2681+ itext_LineColorR,
2682+ itext_LineColorG,
2683+ itext_LineColorB,
2684+ iradio_BarColor,
2685+ ibutton_section,
2686+ itext_BZ_number0,
2687+ itext_BZ_number1,
2688+ itext_BZ_number2
2689+};
2690+
2691+function MyFrame::button_refresh(
2692+ wxCommandEvent& event//!<[in] Selected menu
2693+) {
2694+ Refresh(false);
2695+}
2696+
2697+function MyFrame::button_compute(
2698+ wxCommandEvent& event//!<[in] Selected menu
2699+) {
2700+ free_patch();
2701+ compute_patch_segment();
2702+ Refresh(false);
2703+}
2704+function MyFrame::button_section(
2705+ wxCommandEvent& event//!<[in] Selected menu
2706+) {
2707+ let ib, i2d, i, j;
2708+ FILE *fp;
2709+ fp = fopen("fermi_line.dat", "w");
2710+ for (ib = 0; ib < nb; ib++) {
2711+ if (draw_band[ib] == 1) {
2712+ for (i2d = 0; i2d < n2d[ib]; i2d++) {
2713+ for (i = 0; i < 2; i++) {
2714+ fprintf(fp, "%15.5e %15.5e\n",
2715+ kv2d[ib][i * 3 + 6 * i2d],
2716+ kv2d[ib][1 + i * 3 + 6 * i2d]);
2717+ }
2718+ fprintf(fp, "\n\n");
2719+ }
2720+ }/*if (draw_band[ib] == 1)*/
2721+ }/* for (ib = 0; ib < nb; ib++)*/
2722+ fclose(fp);
2723+ terminal(" fermi_line.dat was written.\n");
2724+
2725+ fp = fopen("bz_line.dat", "w");
2726+ for (i2d = 0; i2d < nbzl2d; i2d++) {
2727+ fprintf(fp, "%15.5e %15.5e\n",
2728+ bzl2d_proj[i2d][0], bzl2d_proj[i2d][1]);
2729+ }
2730+ fprintf(fp, "%15.5e %15.5e\n",
2731+ bzl2d_proj[0][0], bzl2d_proj[0][1]);
2732+ fclose(fp);
2733+ terminal(" bz_line.dat was written.\n");
2734+}
2735+/**
2736+@brief Change Line color color (::blackback)
2737+*/
2738+function MyFrame::textctrl_LineColor(
2739+ wxCommandEvent& event //!<[in] Selected menu
2740+)
2741+{
2742+ let ierr;
2743+ double dvalue;
2744+
2745+ if (event.GetId() == itext_LineColorR) {
2746+ if (event.GetString().ToDouble(&dvalue)) {
2747+ LineColor[0] = dvalue;
2748+ Refresh(false);
2749+ }
2750+ }
2751+ else if (event.GetId() == itext_LineColorG) {
2752+ if (event.GetString().ToDouble(&dvalue)) {
2753+ LineColor[1] = dvalue;
2754+ Refresh(false);
2755+ }
2756+ }
2757+ else if (event.GetId() == itext_LineColorB) {
2758+ if (event.GetString().ToDouble(&dvalue)) {
2759+ LineColor[2] = dvalue;
2760+ Refresh(false);
2761+ }
2762+ }
2763+}
2764+/**
2765+@brief Change background color (::blackback)
2766+*/
2767+function MyFrame::textctrl_BackGround(
2768+ wxCommandEvent& event //!<[in] Selected menu
2769+)
2770+{
2771+ let ierr;
2772+ double dvalue;
2773+
2774+ if (event.GetId() == itext_BackGroundR) {
2775+ if (event.GetString().ToDouble(&dvalue)) {
2776+ BackGroundColor[0] = dvalue;
2777+ glClearColor(BackGroundColor[0], BackGroundColor[1],
2778+ BackGroundColor[2], BackGroundColor[3]);
2779+ Refresh(false);
2780+ }
2781+ }
2782+ else if (event.GetId() == itext_BackGroundG) {
2783+ if (event.GetString().ToDouble(&dvalue)) {
2784+ BackGroundColor[1] = dvalue;
2785+ glClearColor(BackGroundColor[0], BackGroundColor[1],
2786+ BackGroundColor[2], BackGroundColor[3]);
2787+ Refresh(false);
2788+ }
2789+ }
2790+ else if (event.GetId() == itext_BackGroundB) {
2791+ if (event.GetString().ToDouble(&dvalue)) {
2792+ BackGroundColor[2] = dvalue;
2793+ glClearColor(BackGroundColor[0], BackGroundColor[1],
2794+ BackGroundColor[2], BackGroundColor[3]);
2795+ Refresh(false);
2796+ }
2797+ }
2798+}
2799+/**
2800+@brief Change Number of Brillouin zone
2801+*/
2802+function MyFrame::textctrl_BZ_number(
2803+ wxCommandEvent& event //!<[in] Selected menu
2804+)
2805+{
2806+ let ierr;
2807+ double dvalue;
2808+
2809+ if (event.GetId() == itext_BZ_number0) {
2810+ if (event.GetString().ToDouble(&dvalue)) {
2811+ BZ_number[0] = dvalue;
2812+ Refresh(false);
2813+ }
2814+ }
2815+ else if (event.GetId() == itext_BZ_number1) {
2816+ if (event.GetString().ToDouble(&dvalue)) {
2817+ BZ_number[1] = dvalue;
2818+ Refresh(false);
2819+ }
2820+ }
2821+ else if (event.GetId() == itext_BZ_number2) {
2822+ if (event.GetString().ToDouble(&dvalue)) {
2823+ BZ_number[2] = dvalue;
2824+ Refresh(false);
2825+ }
2826+ }
2827+}
2828+ /**
2829+ @brief Toggle the appearance of each band (::draw_band)
2830+*/
2831+function MyFrame::check_band(
2832+ wxCommandEvent& event //!<[in] Selected menu
2833+)
2834+{
2835+ let ib = event.GetId() - icheck_band;
2836+ if (draw_band[ib] == 0) {
2837+ draw_band[ib] = 1;
2838+ }
2839+ else {
2840+ draw_band[ib] = 0;
2841+ }
2842+ Refresh(false);
2843+} /* menu_band */
2844+/**
2845+ @brief Change Brillouin zone (::fbz)
2846+*/
2847+function MyFrame::radio_brillouinzone(
2848+ wxCommandEvent& event //!<[in] Selected menu
2849+)
2850+{
2851+ if (event.GetString().Cmp(wxT("First Brillouin zone")) == 0) {
2852+ fbz = 1;
2853+ }
2854+ else if (event.GetString().Cmp(wxT("Primitive Brillouin zone")) == 0) {
2855+ fbz = -1;
2856+ }
2857+ refresh_patch = 1;
2858+} /* menu_brillouinzone */
2859+/**
2860+ @brief Change Brillouin zone (::fbz)
2861+*/
2862+function MyFrame::radio_BarColor(
2863+ wxCommandEvent& event //!<[in] Selected menu
2864+)
2865+{
2866+ let ii;
2867+ if (event.GetString().Cmp(wxT("BGR")) == 0) {
2868+ for (ii = 0; ii < 4; ii++) {
2869+ BarColor[0][ii] = blue[ii];
2870+ BarColor[1][ii] = cyan[ii];
2871+ BarColor[2][ii] = green[ii];
2872+ BarColor[3][ii] = yellow[ii];
2873+ BarColor[4][ii] = red[ii];
2874+ }
2875+ }
2876+ else if (event.GetString().Cmp(wxT("CMY")) == 0) {
2877+ for (ii = 0; ii < 4; ii++) {
2878+ BarColor[0][ii] = cyan[ii];
2879+ BarColor[1][ii] = blue[ii];
2880+ BarColor[2][ii] = magenta[ii];
2881+ BarColor[3][ii] = red[ii];
2882+ BarColor[4][ii] = yellow[ii];
2883+ }
2884+ }
2885+ else if (event.GetString().Cmp(wxT("MCY")) == 0) {
2886+ for (ii = 0; ii < 4; ii++) {
2887+ BarColor[0][ii] = magenta[ii];
2888+ BarColor[1][ii] = blue[ii];
2889+ BarColor[2][ii] = cyan[ii];
2890+ BarColor[3][ii] = green[ii];
2891+ BarColor[4][ii] = yellow[ii];
2892+ }
2893+ }
2894+ paint();
2895+ Refresh(false);
2896+} /* menu_brillouinzone */
2897+/**
2898+ @brief Toggle Colorbar (::lcolorbar)
2899+*/
2900+function MyFrame::check_colorbar(
2901+ wxCommandEvent& event //!<[in] Selected menu
2902+)
2903+{
2904+ if (lcolorbar != 1) lcolorbar = 1;
2905+ else lcolorbar = 0;
2906+ Refresh(false);
2907+} /* menu_colorbar */
2908+/**
2909+ @brief Change color scale mode (::color_scale)
2910+*/
2911+function MyFrame::radiovalue_colorscale(
2912+ wxCommandEvent& event //!<[in] Selected menu
2913+)
2914+{
2915+ let ierr, ii;
2916+ double dminmax;
2917+
2918+ if (event.GetId() == itext_colorscalemin) {
2919+ if (event.GetString().ToDouble(&dminmax))
2920+ patch_min = dminmax;
2921+ skip_minmax = 1;
2922+ }
2923+ else if (event.GetId() == itext_colorscalemax) {
2924+ if (event.GetString().ToDouble(&dminmax))
2925+ patch_max = dminmax;
2926+ skip_minmax = 1;
2927+ }
2928+ else if (event.GetString().Cmp(wxT("Input (1D)")) == 0)
2929+ color_scale = 1;
2930+ else if (event.GetString().Cmp(wxT("Input (2D)")) == 0)
2931+ color_scale = 2;
2932+ else if (event.GetString().Cmp(wxT("Input (3D)")) == 0)
2933+ color_scale = 3;
2934+ else if (event.GetString().Cmp(wxT("Fermi Velocity")) == 0)
2935+ color_scale = 4;
2936+ else if (event.GetString().Cmp(wxT("Band Index")) == 0)
2937+ color_scale = 5;
2938+ else if (event.GetString().Cmp(wxT("Input (1D, Gray Scale)")) == 0)
2939+ color_scale = 6;
2940+ else if (event.GetString().Cmp(wxT("Fermi Velocity (Gray Scale)")) == 0)
2941+ color_scale = 7;
2942+ refresh_color = 1;
2943+} /* menu_colorscale */
2944+/**
2945+ @brief Modify and toggle appearance of equator (::lequator)
2946+*/
2947+function MyFrame::checkvalue_equator(
2948+ wxCommandEvent& event //!<[in] Selected menu
2949+)
2950+{
2951+ let ierr, ii, jj, ib;
2952+ double deqvec;
2953+
2954+ if (event.GetId() == icheck_equator) {
2955+ if (lequator != 1) lequator = 1;
2956+ else lequator = 0;
2957+ Refresh(false);
2958+ }/*if (event.GetId() == 1)*/
2959+ else {
2960+ if (event.GetId() == itext_equatorx) {
2961+ if (event.GetString().ToDouble(&deqvec)) eqvec_fr[0] = deqvec;
2962+ }
2963+ else if (event.GetId() == itext_equatory) {
2964+ if (event.GetString().ToDouble(&deqvec)) eqvec_fr[1] = deqvec;
2965+ }
2966+ else if (event.GetId() == itext_equatorz) {
2967+ if (event.GetString().ToDouble(&deqvec)) eqvec_fr[2] = deqvec;
2968+ }
2969+ /*
2970+ Fractional -> Cartecian
2971+ */
2972+ for (ii = 0; ii < 3; ii++) {
2973+ eqvec[ii] = 0.0;
2974+ for (jj = 0; jj < 3; jj++) {
2975+ eqvec[ii] += eqvec_fr[jj] * bvec[jj][ii];
2976+ }/*for (jj = 0; jj < 3; jj++)*/
2977+ }/*for (ii = 0; ii < 3; ii++)*/
2978+ refresh_equator = 1;
2979+ }/*else if (event.GetId() > 1)*/
2980+} /*function menu_equator*/
2981+/**
2982+ @brief Modify interpolation ratio
2983+
2984+ This routine modify interpolation ratio (::interpol)
2985+ then compute Fermi surfaces, etc.
2986+*/
2987+function MyFrame::textctrl_interpol(
2988+ wxCommandEvent& event //!<[in] Selected menu
2989+)
2990+{
2991+ let ierr;
2992+ long let long_interpol;
2993+
2994+ if (event.GetString().ToLong(&long_interpol)) {
2995+ interpol = long_interpol;
2996+ refresh_interpol = 1;
2997+ }
2998+}/*function menu_interpol*/
2999+/**
3000+ @brief Toggle Lighting
3001+*/
3002+function MyFrame::radio_lighting(
3003+ wxCommandEvent& event //!<[in] Selected menu
3004+)
3005+{
3006+ if (event.GetString().Cmp(wxT("Both")) == 0) {
3007+ side = 1.0;
3008+ glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
3009+ }
3010+ if (event.GetString().Cmp(wxT("Unoccupy")) == 0) {
3011+ side = 1.0;
3012+ glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, GL_FALSE);
3013+ }
3014+ if (event.GetString().Cmp(wxT("Occupy")) == 0) {
3015+ side = -1.0;
3016+ glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, GL_FALSE);
3017+ }
3018+ Refresh(false);
3019+} /* menu_lighting */
3020+/**
3021+ @brief Line width
3022+*/
3023+function MyFrame::textctrl_line(
3024+ wxCommandEvent& event //!<[in] Selected menu
3025+)
3026+{
3027+ let ierr;
3028+ double dlinewidth;
3029+
3030+ if (event.GetString().ToDouble(&dlinewidth)) linewidth = dlinewidth;
3031+ Refresh(false);
3032+} /* menu_line */
3033+/**
3034+ @brief Change the function associated to the mouse movement(::lmouse)
3035+*/
3036+function MyFrame::radio_mouse(
3037+ wxCommandEvent& event //!<[in] Selected menu
3038+)
3039+{
3040+ if (event.GetString().Cmp(wxT("Rotate")) == 0) lmouse = 1;
3041+ else if (event.GetString().Cmp(wxT("Scale")) == 0) lmouse = 2;
3042+ else if (event.GetString().Cmp(wxT("Translate")) == 0) lmouse = 3;
3043+ Refresh(false);
3044+} /* menu_mouse */
3045+/**
3046+ @brief Toggle apearance of nodale-line
3047+*/
3048+function MyFrame::check_nodeline(
3049+ wxCommandEvent& event //!<[in] Selected menu
3050+)
3051+{
3052+ if (nodeline != 1) nodeline = 1;
3053+ else nodeline = 0;
3054+ Refresh(false);
3055+}/*menu_nodeline*/
3056+/**
3057+ @brief Modify and toggle appearance of 2D Fermi lines (::lsection)
3058+*/
3059+function MyFrame::radiovalue_section(
3060+ wxCommandEvent& event //!<[in] Selected menu
3061+)
3062+{
3063+ let ierr, ii, jj, ib;
3064+ double dsecvec;
3065+
3066+ if (event.GetId() == icheck_section) {
3067+ if (lsection == 0)lsection = 1;
3068+ else lsection = 0;
3069+ Refresh(false);
3070+ }
3071+ else if(event.GetId() == icheck_gamma) {
3072+ if(secscale > 0.5) secscale = 0.0;
3073+ else secscale = 1.0;
3074+ refresh_section = 1;
3075+ }
3076+ else {
3077+ if (event.GetId() == itext_sectionx) {
3078+ if (event.GetString().ToDouble(&dsecvec)) secvec_fr[0] = dsecvec;
3079+ }
3080+ else if (event.GetId() == itext_sectiony) {
3081+ if (event.GetString().ToDouble(&dsecvec)) secvec_fr[1] = dsecvec;
3082+ }
3083+ else if (event.GetId() == itext_sectionz) {
3084+ if (event.GetString().ToDouble(&dsecvec)) secvec_fr[2] = dsecvec;
3085+ }
3086+ /*
3087+ Fractional -> Cartecian
3088+ */
3089+ for (ii = 0; ii < 3; ii++) {
3090+ secvec[ii] = 0.0;
3091+ for (jj = 0; jj < 3; jj++) {
3092+ secvec[ii] += secvec_fr[jj] * bvec[jj][ii];
3093+ }/*for (jj = 0; jj < 3; jj++)*/
3094+ }/*for (ii = 0; ii < 3; ii++)*/
3095+ refresh_section = 1;
3096+ }/*else if (event.GetId() > 1)*/
3097+} /*function menu_section*/
3098+/**
3099+ @brief Shift Fermi energy
3100+*/
3101+function MyFrame::textctrl_shift(
3102+ wxCommandEvent& event //!<[in] Selected menu
3103+)
3104+{
3105+ double dEF;
3106+ if (event.GetString().ToDouble(&dEF)) {
3107+ EF = dEF;
3108+ refresh_patch = 1;
3109+ }
3110+} /* menu_shift */
3111+/**
3112+ @brief Tern stereogram (::lstereo)
3113+*/
3114+function MyFrame::radio_stereo(
3115+ wxCommandEvent& event //!<[in] Selected menu
3116+) {
3117+ if (event.GetString().Cmp(wxT("None")) == 0) lstereo = 1;
3118+ else if (event.GetString().Cmp(wxT("Parallel")) == 0) lstereo = 2;
3119+ else if (event.GetString().Cmp(wxT("Cross")) == 0) lstereo = 3;
3120+ Refresh(false);
3121+} /* menu_stereo */
3122+/**
3123+ @brief Change tetrahedron (::itet)
3124+*/
3125+function MyFrame::radio_tetra(
3126+ wxCommandEvent& event //!<[in] Selected menu
3127+)
3128+{
3129+ itet = wxAtoi(event.GetString()) - 1;
3130+ init_corner();
3131+ refresh_patch = 1;
3132+}/*menu_tetra*/
3133+ /**
3134+ @brief Setting of view
3135+
3136+ This modify scale (::scl) & tarnslation (::trans) &
3137+ rotation (::thetax, ::thetay, ::thetaz, ::rot),
3138+ */
3139+function MyFrame::textctrl_view(
3140+ wxCommandEvent& event //!<[in] Selected menu
3141+)
3142+{
3143+ let ierr;
3144+ double dvalue;
3145+
3146+ if (event.GetId() == itext_scale) {
3147+ if (event.GetString().ToDouble(&dvalue)) {
3148+ linewidth *= dvalue/scl;
3149+ scl = dvalue;
3150+ textbox_linewidth->ChangeValue(wxString::Format(wxT("%f"), linewidth));
3151+ Refresh(false);
3152+ }
3153+ }
3154+ else if (event.GetId() == itext_positionx) {
3155+ if (event.GetString().ToDouble(&dvalue)) {
3156+ trans[0] = dvalue;
3157+ Refresh(false);
3158+ }
3159+ }
3160+ else if (event.GetId() == itext_positiony) {
3161+ if (event.GetString().ToDouble(&dvalue)) {
3162+ trans[1] = dvalue;
3163+ Refresh(false);
3164+ }
3165+ }
3166+ else if (event.GetId() == itext_rotx ||
3167+ event.GetId() == itext_roty ||
3168+ event.GetId() == itext_rotz) {
3169+ /*
3170+ thetay = Math.asin(rot[0][2]);
3171+ if (Math.cos(thetay) != 0.0) {
3172+ if (-rot[1][2] / Math.cos(thetay) >= 0.0) thetax = Math.acos(rot[2][2] / Math.cos(thetay));
3173+ else thetax = 6.283185307f - Math.acos(rot[2][2] / Math.cos(thetay));
3174+ if (-rot[0][1] / Math.cos(thetay) >= 0.0) thetaz = Math.acos(rot[0][0] / Math.cos(thetay));
3175+ else thetaz = 6.283185307f - Math.acos(rot[0][0] / Math.cos(thetay));
3176+ }
3177+ else {
3178+ thetax = 0.0;
3179+ if (rot[1][0] >= 0.0) thetaz = Math.acos(rot[1][1]);
3180+ else thetaz = 6.283185307f - Math.acos(rot[1][1]);
3181+ }
3182+ thetax = 180.0 / 3.14159265f * thetax;
3183+ thetay = 180.0 / 3.14159265f * thetay;
3184+ thetaz = 180.0 / 3.14159265f * thetaz;
3185+ printf(" Current Rotation (theta_x theta_y teta_z) in degree : %f %f %f\n", thetax, thetay, thetaz);
3186+ printf(" New Rotation (theta_x theta_y teta_z) in degree : ");
3187+ */
3188+ if (event.GetId() == itext_rotx) {
3189+ if (event.GetString().ToDouble(&dvalue))
3190+ thetax = (dvalue) / 180.0 * 3.14159265;
3191+ }
3192+ else if (event.GetId() == itext_roty) {
3193+ if (event.GetString().ToDouble(&dvalue))
3194+ thetay = (dvalue) / 180.0 * 3.14159265;
3195+ }
3196+ else if (event.GetId() == itext_rotz) {
3197+ if (event.GetString().ToDouble(&dvalue))
3198+ thetaz = (dvalue) / 180.0 * 3.14159265;
3199+ }
3200+ }
3201+ else if(event.GetId() == ibutton_rotate){
3202+ rot[0][0] = Math.cos(thetay)* Math.cos(thetaz);
3203+ rot[0][1] = -Math.cos(thetay)* Math.sin(thetaz);
3204+ rot[0][2] = Math.sin(thetay);
3205+ rot[1][0] = Math.cos(thetaz)* Math.sin(thetax)* Math.sin(thetay) + Math.cos(thetax)* Math.sin(thetaz);
3206+ rot[1][1] = Math.cos(thetax) * Math.cos(thetaz) - Math.sin(thetax)* Math.sin(thetay)* Math.sin(thetaz);
3207+ rot[1][2] = -Math.cos(thetay)* Math.sin(thetax);
3208+ rot[2][0] = -Math.cos(thetax)* Math.cos(thetaz)* Math.sin(thetay) + Math.sin(thetax)* Math.sin(thetaz);
3209+ rot[2][1] = Math.cos(thetaz)* Math.sin(thetax) + Math.cos(thetax)* Math.sin(thetay)* Math.sin(thetaz);
3210+ rot[2][2] = Math.cos(thetax)* Math.cos(thetay);
3211+ Refresh(false);
3212+ }
3213+}
3214+
3215+MyFrame::MyFrame(wxFrame* frame, const wxString& title, const wxPoint& pos,
3216+ const wxSize& size, long style)
3217+ : wxFrame(frame, wxID_ANY, title, pos, size, style),
3218+ m_canvas(NULL)
3219+{
3220+ let ib, itet;
3221+ char menuname[8];
3222+
3223+ SetIcon(wxICON(fermisurfer));
3224+
3225+ // Make a menubar
3226+ //wxMenu* fileMenu = new wxMenu;
3227+ //wxMenuBar* menuBar = new wxMenuBar;
3228+ //menuBar->Append(fileMenu, wxT("File"));
3229+ //SetMenuBar(menuBar);
3230+ // JACS
3231+#ifdef __WXMSW__
3232+ let* gl_attrib = NULL;
3233+#else
3234+ let gl_attrib[20] =
3235+ { WX_GL_RGBA, WX_GL_MIN_RED, 1, WX_GL_MIN_GREEN, 1,
3236+ WX_GL_MIN_BLUE, 1, WX_GL_DEPTH_SIZE, 1,
3237+ WX_GL_DOUBLEBUFFER,
3238+# if defined(__WXMAC__) || defined(__WXCOCOA__)
3239+ GL_NONE };
3240+# else
3241+ None
3242+};
3243+# endif
3244+#endif
3245+
3246+ wxBoxSizer* sizermain = new wxBoxSizer(wxVERTICAL);
3247+
3248+ splitterV = new wxSplitterWindow(this, wxID_ANY);
3249+ splitterV->SetSashGravity(1.0);
3250+ splitterV->SetMinimumPaneSize(0);
3251+ sizermain->Add(splitterV, 1, wxEXPAND, 0);
3252+
3253+ panel = new wxPanel(splitterV);
3254+
3255+ gbsizer = new wxGridBagSizer();
3256+
3257+ splitterH = new wxSplitterWindow(splitterV, wxID_ANY);
3258+ splitterH->SetSashGravity(1.0);
3259+ splitterH->SetMinimumPaneSize(0);
3260+
3261+ m_canvas = new TestGLCanvas(splitterH, wxID_ANY, gl_attrib);
3262+
3263+ terminal = new wxTextCtrl(splitterH, wxID_ANY, wxT(""),
3264+ wxPoint(0, 250), wxSize(100, 0), wxTE_MULTILINE);
3265+ if (lbatch == 1) splitterH->SplitHorizontally(m_canvas, terminal, -1);
3266+ else splitterH->SplitHorizontally(m_canvas, terminal, -200);
3267+ panel->SetSizer(gbsizer);
3268+
3269+ if (lbatch == 1)splitterV->SplitVertically(splitterH, panel, -1);
3270+ else splitterV->SplitVertically(splitterH, panel, 0);
3271+
3272+ Bind(wxEVT_COMMAND_BUTTON_CLICKED, &MyFrame::button_compute, this, ibutton_compute);
3273+ gbsizer->Add(new wxButton(panel, ibutton_compute, wxT("Update")), wxGBPosition(0,0), wxGBSpan(1,1));
3274+
3275+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Line width : ")),
3276+ wxGBPosition(0, 1), wxGBSpan(1, 1), wxALIGN_RIGHT);
3277+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_line, this, itext_line);
3278+ textbox_linewidth = new wxTextCtrl(panel, itext_line, wxT(""));
3279+ gbsizer->Add(textbox_linewidth, wxGBPosition(0, 2), wxGBSpan(1, 1));
3280+ textbox_linewidth->ChangeValue(wxString::Format(wxT("%f"), linewidth));
3281+
3282+ Bind(wxEVT_COMMAND_CHECKBOX_CLICKED, &MyFrame::radiovalue_section, this, icheck_gamma);
3283+ wxCheckBox *check_ongamma = new wxCheckBox(panel, icheck_gamma, wxT("On Gamma"));
3284+ gbsizer->Add(check_ongamma, wxGBPosition(0, 3), wxGBSpan(1, 1));
3285+ check_ongamma->SetValue(true);
3286+
3287+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Section-v : ")),
3288+ wxGBPosition(1, 0), wxGBSpan(1, 1), wxALIGN_RIGHT);
3289+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::radiovalue_section, this, itext_sectionx);
3290+ wxTextCtrl* textbox_sectionx = new wxTextCtrl(panel, itext_sectionx, wxT(""));
3291+ gbsizer->Add(textbox_sectionx, wxGBPosition(1, 1), wxGBSpan(1, 1));
3292+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::radiovalue_section, this, itext_sectiony);
3293+ wxTextCtrl* textbox_sectiony = new wxTextCtrl(panel, itext_sectiony, wxT(""));
3294+ gbsizer->Add(textbox_sectiony, wxGBPosition(1, 2), wxGBSpan(1, 1));
3295+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::radiovalue_section, this, itext_sectionz);
3296+ wxTextCtrl* textbox_sectionz = new wxTextCtrl(panel, itext_sectionz, wxT(""));
3297+ gbsizer->Add(textbox_sectionz, wxGBPosition(1, 3), wxGBSpan(1, 1));
3298+ textbox_sectionx->ChangeValue(wxT("0"));
3299+ textbox_sectiony->ChangeValue(wxT("0"));
3300+ textbox_sectionz->ChangeValue(wxT("1"));
3301+
3302+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Equator-v : ")),
3303+ wxGBPosition(2, 0), wxGBSpan(1, 1), wxALIGN_RIGHT);
3304+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::checkvalue_equator, this, itext_equatorx);
3305+ wxTextCtrl* textbox_equatorx = new wxTextCtrl(panel, itext_equatorx, wxT(""));
3306+ gbsizer->Add(textbox_equatorx, wxGBPosition(2, 1), wxGBSpan(1, 1));
3307+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::checkvalue_equator, this, itext_equatory);
3308+ wxTextCtrl* textbox_equatory = new wxTextCtrl(panel, itext_equatory, wxT(""));
3309+ gbsizer->Add(textbox_equatory, wxGBPosition(2, 2), wxGBSpan(1, 1));
3310+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::checkvalue_equator, this, itext_equatorz);
3311+ wxTextCtrl* textbox_equatorz = new wxTextCtrl(panel, itext_equatorz, wxT(""));
3312+ gbsizer->Add(textbox_equatorz, wxGBPosition(2, 3), wxGBSpan(1, 1));
3313+ textbox_equatorx->ChangeValue(wxT("0"));
3314+ textbox_equatory->ChangeValue(wxT("0"));
3315+ textbox_equatorz->ChangeValue(wxT("1"));
3316+
3317+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Interpol ratio : ")),
3318+ wxGBPosition(3, 0), wxGBSpan(1, 1), wxALIGN_RIGHT);
3319+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_interpol, this, itext_interpol);
3320+ wxTextCtrl* textbox_interpol = new wxTextCtrl(panel, itext_interpol, wxT(""));
3321+ gbsizer->Add(textbox_interpol, wxGBPosition(3, 1), wxGBSpan(1, 1));
3322+ textbox_interpol->ChangeValue(wxT("1"));
3323+
3324+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Fermi energy : ")),
3325+ wxGBPosition(4, 0), wxGBSpan(1, 1), wxALIGN_RIGHT);
3326+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_shift, this, itext_shift);
3327+ wxTextCtrl* textbox_shift = new wxTextCtrl(panel, itext_shift, wxT(""));
3328+ gbsizer->Add(textbox_shift, wxGBPosition(4, 1), wxGBSpan(1, 1));
3329+ textbox_shift->ChangeValue(wxT("0"));
3330+
3331+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Min of Scale : ")),
3332+ wxGBPosition(5, 0), wxGBSpan(1, 1), wxALIGN_RIGHT);
3333+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::radiovalue_colorscale, this, itext_colorscalemin);
3334+ textbox_min = new wxTextCtrl(panel, itext_colorscalemin, wxT(""));
3335+ gbsizer->Add(textbox_min, wxGBPosition(5, 1), wxGBSpan(1, 1));
3336+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Max of Scale : ")),
3337+ wxGBPosition(6, 0), wxGBSpan(1, 1), wxALIGN_RIGHT);
3338+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::radiovalue_colorscale, this, itext_colorscalemax);
3339+ textbox_max = new wxTextCtrl(panel, itext_colorscalemax, wxT(""));
3340+ gbsizer->Add(textbox_max, wxGBPosition(6, 1), wxGBSpan(1, 1));
3341+
3342+ wxString choices_tetra[] = { wxT("1"), wxT("2"), wxT("3"), wxT("4"), wxT("5"), wxT("6"), wxT("7") ,
3343+wxT("8"), wxT("9"), wxT("10"), wxT("11"), wxT("12"), wxT("13"), wxT("14"),
3344+ wxT("15"), wxT("16") };
3345+ Bind(wxEVT_COMMAND_RADIOBOX_SELECTED, &MyFrame::radio_tetra, this, iradio_tetra);
3346+ gbsizer->Add(new wxRadioBox(panel, iradio_tetra, wxT("Tetrahedron"),
3347+ wxDefaultPosition, wxDefaultSize,
3348+ WXSIZEOF(choices_tetra), choices_tetra,
3349+ 4, wxRA_SPECIFY_COLS), wxGBPosition(3,2), wxGBSpan(5, 2));
3350+
3351+ wxString choices_colorscale[] = { wxT("Input (1D)"), wxT("Input (2D)"),
3352+ wxT("Input (3D)"), wxT("Fermi Velocity"), wxT("Band Index"),
3353+ wxT("Input (1D, Gray Scale)"), wxT("Fermi Velocity (Gray Scale)") };
3354+ radiobox_color = new wxRadioBox(panel, iradio_colorscale,
3355+ wxT("Color scale mode"),
3356+ wxDefaultPosition, wxDefaultSize,
3357+ WXSIZEOF(choices_colorscale), choices_colorscale,
3358+ 1, wxRA_SPECIFY_COLS);
3359+ Bind(wxEVT_COMMAND_RADIOBOX_SELECTED, &MyFrame::radiovalue_colorscale, this, iradio_colorscale);
3360+ gbsizer->Add(radiobox_color, wxGBPosition(7, 0), wxGBSpan(3, 2));
3361+
3362+ wxString choices_bz[] = { wxT("First Brillouin zone"), wxT("Primitive Brillouin zone") };
3363+ Bind(wxEVT_COMMAND_RADIOBOX_SELECTED, &MyFrame::radio_brillouinzone, this, iradio_brillouinzone);
3364+ gbsizer->Add(new wxRadioBox(panel, iradio_brillouinzone, wxT("Brillouin zone"),
3365+ wxDefaultPosition, wxDefaultSize,
3366+ WXSIZEOF(choices_bz), choices_bz,
3367+ 1, wxRA_SPECIFY_COLS), wxGBPosition(8, 2), wxGBSpan(1, 2));
3368+
3369+ wxString choices_stereo[] = { wxT("None"), wxT("Parallel"), wxT("Cross") };
3370+ Bind(wxEVT_COMMAND_RADIOBOX_SELECTED, &MyFrame::radio_stereo, this, iradio_stereo);
3371+ gbsizer->Add(new wxRadioBox(panel, iradio_stereo, wxT("Stereogram"),
3372+ wxDefaultPosition, wxDefaultSize,
3373+ WXSIZEOF(choices_stereo), choices_stereo,
3374+ 1, wxRA_SPECIFY_COLS), wxGBPosition(9, 2), wxGBSpan(1, 1));
3375+
3376+ wxString choices_mouse[] = { wxT("Rotate"), wxT("Scale"), wxT("Translate") };
3377+ Bind(wxEVT_COMMAND_RADIOBOX_SELECTED, &MyFrame::radio_mouse, this, iradio_mouse);
3378+ gbsizer->Add(new wxRadioBox(panel, iradio_mouse, wxT("Mouse Drag"),
3379+ wxDefaultPosition, wxDefaultSize,
3380+ WXSIZEOF(choices_mouse), choices_mouse,
3381+ 1, wxRA_SPECIFY_COLS), wxGBPosition(9, 3), wxGBSpan(1, 1));
3382+
3383+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("BZ number : ")),
3384+ wxGBPosition(10, 0), wxGBSpan(1, 1), wxALIGN_RIGHT);
3385+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_BZ_number, this, itext_BZ_number0);
3386+ textbox_BZ_number0 = new wxTextCtrl(panel, itext_BZ_number0, wxT(""));
3387+ gbsizer->Add(textbox_BZ_number0, wxGBPosition(10, 1), wxGBSpan(1, 1));
3388+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_BZ_number, this, itext_BZ_number1);
3389+ textbox_BZ_number1 = new wxTextCtrl(panel, itext_BZ_number1, wxT(""));
3390+ gbsizer->Add(textbox_BZ_number1, wxGBPosition(10, 2), wxGBSpan(1, 1));
3391+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_BZ_number, this, itext_BZ_number2);
3392+ textbox_BZ_number2 = new wxTextCtrl(panel, itext_BZ_number2, wxT(""));
3393+ gbsizer->Add(textbox_BZ_number2, wxGBPosition(10, 3), wxGBSpan(1, 1));
3394+ textbox_BZ_number0->ChangeValue(wxT("1"));
3395+ textbox_BZ_number1->ChangeValue(wxT("1"));
3396+ textbox_BZ_number2->ChangeValue(wxT("1"));
3397+
3398+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Background (RGB) : ")),
3399+ wxGBPosition(11, 0), wxGBSpan(1, 1), wxALIGN_RIGHT);
3400+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_BackGround, this, itext_BackGroundR);
3401+ textbox_BackGroundR = new wxTextCtrl(panel, itext_BackGroundR, wxT(""));
3402+ gbsizer->Add(textbox_BackGroundR, wxGBPosition(11, 1), wxGBSpan(1, 1));
3403+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_BackGround, this, itext_BackGroundG);
3404+ textbox_BackGroundG = new wxTextCtrl(panel, itext_BackGroundG, wxT(""));
3405+ gbsizer->Add(textbox_BackGroundG, wxGBPosition(11, 2), wxGBSpan(1, 1));
3406+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_BackGround, this, itext_BackGroundB);
3407+ textbox_BackGroundB = new wxTextCtrl(panel, itext_BackGroundB, wxT(""));
3408+ gbsizer->Add(textbox_BackGroundB, wxGBPosition(11, 3), wxGBSpan(1, 1));
3409+ textbox_BackGroundR->ChangeValue(wxT("0"));
3410+ textbox_BackGroundG->ChangeValue(wxT("0"));
3411+ textbox_BackGroundB->ChangeValue(wxT("0"));
3412+
3413+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Line color (RGB) : ")),
3414+ wxGBPosition(12, 0), wxGBSpan(1, 1), wxALIGN_RIGHT);
3415+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_LineColor, this, itext_LineColorR);
3416+ textbox_LineColorR = new wxTextCtrl(panel, itext_LineColorR, wxT(""));
3417+ gbsizer->Add(textbox_LineColorR, wxGBPosition(12, 1), wxGBSpan(1, 1));
3418+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_LineColor, this, itext_LineColorG);
3419+ textbox_LineColorG = new wxTextCtrl(panel, itext_LineColorG, wxT(""));
3420+ gbsizer->Add(textbox_LineColorG, wxGBPosition(12, 2), wxGBSpan(1, 1));
3421+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_LineColor, this, itext_LineColorB);
3422+ textbox_LineColorB = new wxTextCtrl(panel, itext_LineColorB, wxT(""));
3423+ gbsizer->Add(textbox_LineColorB, wxGBPosition(12, 3), wxGBSpan(1, 1));
3424+ textbox_LineColorR->ChangeValue(wxT("1"));
3425+ textbox_LineColorG->ChangeValue(wxT("1"));
3426+ textbox_LineColorB->ChangeValue(wxT("1"));
3427+
3428+ Bind(wxEVT_COMMAND_BUTTON_CLICKED, &MyFrame::textctrl_view, this, ibutton_rotate);
3429+ gbsizer->Add(new wxButton(panel, ibutton_rotate, wxT("Rotate")),
3430+ wxGBPosition(13, 0), wxGBSpan(1, 1));
3431+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_view, this, itext_rotx);
3432+ textbox_rotatex = new wxTextCtrl(panel, itext_rotx, wxT(""));
3433+ gbsizer->Add(textbox_rotatex, wxGBPosition(13, 1), wxGBSpan(1, 1));
3434+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_view, this, itext_roty);
3435+ textbox_rotatey = new wxTextCtrl(panel, itext_roty, wxT(""));
3436+ gbsizer->Add(textbox_rotatey, wxGBPosition(13, 2), wxGBSpan(1, 1));
3437+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_view, this, itext_rotz);
3438+ textbox_rotatez = new wxTextCtrl(panel, itext_rotz, wxT(""));
3439+ gbsizer->Add(textbox_rotatez, wxGBPosition(13, 3), wxGBSpan(1, 1));
3440+ textbox_rotatex->ChangeValue(wxT("0"));
3441+ textbox_rotatey->ChangeValue(wxT("0"));
3442+ textbox_rotatez->ChangeValue(wxT("0"));
3443+
3444+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Position : ")),
3445+ wxGBPosition(14, 1), wxGBSpan(1, 1), wxALIGN_RIGHT);
3446+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_view, this, itext_positionx);
3447+ textbox_positionx = new wxTextCtrl(panel, itext_positionx, wxT(""));
3448+ gbsizer->Add(textbox_positionx, wxGBPosition(14, 2), wxGBSpan(1, 1));
3449+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_view, this, itext_positiony);
3450+ textbox_positiony = new wxTextCtrl(panel, itext_positiony, wxT(""));
3451+ gbsizer->Add(textbox_positiony, wxGBPosition(14, 3), wxGBSpan(1, 1));
3452+ textbox_positionx->ChangeValue(wxT("0"));
3453+ textbox_positiony->ChangeValue(wxT("0"));
3454+
3455+ gbsizer->Add(new wxStaticText(panel, wxID_ANY, wxT("Scale : ")),
3456+ wxGBPosition(15, 2), wxGBSpan(1, 1), wxALIGN_RIGHT);
3457+ Bind(wxEVT_COMMAND_TEXT_UPDATED, &MyFrame::textctrl_view, this, itext_scale);
3458+ textbox_scale = new wxTextCtrl(panel, itext_scale, wxT(""));
3459+ gbsizer->Add(textbox_scale, wxGBPosition(15, 3), wxGBSpan(1, 1));
3460+ textbox_scale->ChangeValue(wxString::Format(wxT("%f"), scl));
3461+
3462+ Bind(wxEVT_COMMAND_CHECKBOX_CLICKED, &MyFrame::check_colorbar, this, icheck_colorbar);
3463+ wxCheckBox* check = new wxCheckBox(panel, icheck_colorbar, wxT("Color bar"));
3464+ gbsizer->Add(check, wxGBPosition(15, 1), wxGBSpan(1, 1));
3465+ check->SetValue(true);
3466+ // debug fileMenu->Check(menu_colorbar_check, true);
3467+
3468+ Bind(wxEVT_COMMAND_CHECKBOX_CLICKED, &MyFrame::checkvalue_equator, this, icheck_equator);
3469+ gbsizer->Add(new wxCheckBox(panel, icheck_equator, wxT("Equator")), wxGBPosition(16, 1), wxGBSpan(1, 1));
3470+
3471+ Bind(wxEVT_COMMAND_CHECKBOX_CLICKED, &MyFrame::check_nodeline, this, icheck_nodeline);
3472+ gbsizer->Add(new wxCheckBox(panel, icheck_nodeline, wxT("Nodal line")), wxGBPosition(17, 1), wxGBSpan(1, 1));
3473+
3474+ Bind(wxEVT_COMMAND_CHECKBOX_CLICKED, &MyFrame::radiovalue_section, this, icheck_section);
3475+ gbsizer->Add(new wxCheckBox(panel, icheck_section, wxT("Section")), wxGBPosition(18, 1), wxGBSpan(1, 1));
3476+
3477+ Bind(wxEVT_COMMAND_BUTTON_CLICKED, &MyFrame::button_section, this, ibutton_section);
3478+ gbsizer->Add(new wxButton(panel, ibutton_section, wxT("Section file")),
3479+ wxGBPosition(19, 1), wxGBSpan(1, 1));
3480+
3481+ wxString choices_light[] = { wxT("Both"), wxT("Unoccupy"), wxT("Occupy") };
3482+ Bind(wxEVT_COMMAND_RADIOBOX_SELECTED, &MyFrame::radio_lighting, this, iradio_lighting);
3483+ gbsizer->Add(new wxRadioBox(panel, iradio_lighting, wxT("Lighting"),
3484+ wxDefaultPosition, wxDefaultSize,
3485+ WXSIZEOF(choices_light), choices_light,
3486+ 1, wxRA_SPECIFY_COLS), wxGBPosition(16, 2), wxGBSpan(4, 1));
3487+
3488+ wxString choices_BarColor[] = { wxT("BGR"), wxT("CMY"), wxT("MCY")};
3489+ Bind(wxEVT_COMMAND_RADIOBOX_SELECTED, &MyFrame::radio_BarColor, this, iradio_BarColor);
3490+ gbsizer->Add(new wxRadioBox(panel, iradio_BarColor, wxT("Bar Color"),
3491+ wxDefaultPosition, wxDefaultSize,
3492+ WXSIZEOF(choices_BarColor), choices_BarColor,
3493+ 1, wxRA_SPECIFY_COLS), wxGBPosition(16, 3), wxGBSpan(4, 1));
3494+
3495+ SetSizer(sizermain);
3496+ SetAutoLayout(true);
3497+
3498+ // Show the frame
3499+ Show(true);
3500+ Raise();
3501+
3502+ m_canvas->InitGL();
3503+}
3504+
3505+MyFrame::~MyFrame()
3506+{
3507+ delete m_canvas;
3508+}
3509+
3510+// Intercept menu commands
3511+function MyFrame::OnExit(wxCommandEvent& WXUNUSED(event))
3512+{
3513+ // true is to force the frame to close
3514+ Close(true);
3515+}
3516+
3517+function MyFrame::modify_band() {
3518+ let ib, width, height;
3519+ wxCheckBox** check;
3520+
3521+ radiobox_color->SetSelection(color_scale - 1);
3522+
3523+ check = new wxCheckBox * [nb];
3524+
3525+ for (ib = 0; ib < nb; ib++) {
3526+ Bind(wxEVT_COMMAND_CHECKBOX_CLICKED, &MyFrame::check_band, this, icheck_band + ib);
3527+ check[ib] = new wxCheckBox(panel, icheck_band + ib,
3528+ wxString::Format(wxT("Band " + toString() + ""), ib));
3529+ gbsizer->Add(check[ib], wxGBPosition(14 + ib, 0), wxGBSpan(1, 1));
3530+ check[ib]->SetValue(true);
3531+ }
3532+ gbsizer->Layout();
3533+ if (lbatch == 1) {
3534+ splitterH->Unsplit();
3535+ splitterV->Unsplit();
3536+ }
3537+ else {
3538+ splitterV->SetSashPosition(-gbsizer->CalcMin().GetX());
3539+ }
3540+ splitterV->Layout();
3541+ splitterH->Layout();
3542+ Refresh(false);
3543+ Raise();
3544+}
3545+/**
3546+ @brief Window resize
3547+
3548+ Modify : ::sx, ::sy
3549+*/
3550+function TestGLCanvas::OnSize(wxSizeEvent& event)
3551+{
3552+ if (!IsShownOnScreen())
3553+ return;
3554+ // This is normally only necessary if there is more than one wxGLCanvas
3555+ // or more than one wxGLContext in the application.
3556+ SetCurrent(*m_glRC);
3557+
3558+ /*
3559+ Scale of translation of mousepointer
3560+ */
3561+ sx = 1.0 / event.GetSize().x;
3562+ sy = 1.0 / event.GetSize().y;
3563+ // It's up to the application code to update the OpenGL viewport settings.
3564+ // This is OK here only because there is only one canvas that uses the
3565+ // context. See the cube sample for that case that multiple canvases are
3566+ // made current with one context.
3567+ glViewport(0, 0, event.GetSize().x, event.GetSize().y);
3568+ /**/
3569+ glMatrixMode(GL_PROJECTION);
3570+ /**/
3571+ glLoadIdentity();
3572+ gluPerspective(30.0, event.GetSize().x / event.GetSize().y, 1.0, 100.0);
3573+ /**/
3574+ glMatrixMode(GL_MODELVIEW);
3575+ Refresh(false);
3576+ //myf->Show(true);
3577+}
3578+/**
3579+ @brief Glut mouse function
3580+
3581+ Modify : ::cx, ::cy, ::scl
3582+*/
3583+function TestGLCanvas::OnMouseEvent(wxMouseEvent& event)
3584+{
3585+ static let dragging = 0;
3586+ static float last_x, last_y;
3587+ let i, j, wheel;
3588+ let dx, dy, a, rot0[3][3], rot1[3][3], ax, ay;
3589+
3590+ // Allow default processing to happen, or else the canvas cannot gain focus
3591+ // (for key events).
3592+ event.Skip();
3593+
3594+ if (event.LeftIsDown())
3595+ {
3596+ if (!dragging)
3597+ {
3598+ dragging = 1;
3599+ }
3600+ else
3601+ {
3602+ /*
3603+ Translation of mousepointer from starting point
3604+ */
3605+ dx = (event.GetX() - last_x) * sx;
3606+ dy = (event.GetY() - last_y) * sy;
3607+ /*
3608+ Distanse from starting point
3609+ */
3610+ a = Math.sqrt(dx * dx + dy * dy);
3611+ /**/
3612+ if (lmouse == 1) {
3613+ /**/
3614+ if (a != 0.0) {
3615+ /*
3616+ Compute rotational matrix from translation of mousepointer
3617+ */
3618+ ax = -dy;
3619+ ay = dx;
3620+ /**/
3621+ a = a * 10.0;
3622+ /**/
3623+ rot0[0][0] = (ax * ax + ay * ay * Math.cos(a)) / (ax * ax + ay * ay);
3624+ rot0[0][1] = ax * ay * (Math.cos(a) - 1.0) / (ax * ax + ay * ay);
3625+ rot0[0][2] = ay * Math.sin(a) / Math.sqrt(ax * ax + ay * ay);
3626+ rot0[1][0] = ax * ay * (Math.cos(a) - 1.0) / (ax * ax + ay * ay);
3627+ rot0[1][1] = (ax * ax * Math.cos(a) + ay * ay) / (ax * ax + ay * ay);
3628+ rot0[1][2] = ax * Math.sin(a) / Math.sqrt(ax * ax + ay * ay);
3629+ rot0[2][0] = -ay * Math.sin(a) / Math.sqrt(ax * ax + ay * ay);
3630+ rot0[2][1] = -ax * Math.sin(a) / Math.sqrt(ax * ax + ay * ay);
3631+ rot0[2][2] = Math.cos(a);
3632+ /**/
3633+ for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) rot1[i][j] = rot[i][j];
3634+ /**/
3635+ for (i = 0; i < 3; i++) {
3636+ for (j = 0; j < 3; j++) {
3637+ rot[i][j] = rot0[i][0] * rot1[0][j]
3638+ + rot0[i][1] * rot1[1][j]
3639+ + rot0[i][2] * rot1[2][j];
3640+ }
3641+ }
3642+ /*
3643+ Print angle to text Box
3644+ */
3645+ thetay = Math.asin(rot[0][2]);
3646+ if (Math.cos(thetay) != 0.0) {
3647+ if (-rot[1][2] / Math.cos(thetay) >= 0.0) thetax = Math.acos(rot[2][2] / Math.cos(thetay));
3648+ else thetax = 6.283185307f - Math.acos(rot[2][2] / Math.cos(thetay));
3649+ if (-rot[0][1] / Math.cos(thetay) >= 0.0) thetaz = Math.acos(rot[0][0] / Math.cos(thetay));
3650+ else thetaz = 6.283185307f - Math.acos(rot[0][0] / Math.cos(thetay));
3651+ }
3652+ else {
3653+ thetax = 0.0;
3654+ if (rot[1][0] >= 0.0) thetaz = Math.acos(rot[1][1]);
3655+ else thetaz = 6.283185307f - Math.acos(rot[1][1]);
3656+ }
3657+ thetax *= 180.0 / 3.14159265;
3658+ thetay *= 180.0 / 3.14159265;
3659+ thetaz *= 180.0 / 3.14159265;
3660+ myf->textbox_rotatex->ChangeValue(wxString::Format(wxT("%f"), thetax));
3661+ myf->textbox_rotatey->ChangeValue(wxString::Format(wxT("%f"), thetay));
3662+ myf->textbox_rotatez->ChangeValue(wxString::Format(wxT("%f"), thetaz));
3663+ myf->Show(true);
3664+ }
3665+ }
3666+ else if (lmouse == 2) {
3667+ scl = scl * expf(-dy);
3668+ myf->textbox_scale->ChangeValue(wxString::Format(wxT("%f"), scl));
3669+ linewidth *= expf(-dy);
3670+ myf->textbox_linewidth->ChangeValue(wxString::Format(wxT("%f"), linewidth));
3671+ myf->Show(true);
3672+ }
3673+ else {
3674+ trans[0] = trans[0] + dx;
3675+ trans[1] = trans[1] - dy;
3676+ myf->textbox_positionx->ChangeValue(wxString::Format(wxT("%f"), trans[0]));
3677+ myf->textbox_positiony->ChangeValue(wxString::Format(wxT("%f"), trans[1]));
3678+ myf->Show(true);
3679+ }
3680+ Refresh(false);
3681+ }
3682+ last_x = event.GetX();
3683+ last_y = event.GetY();
3684+ }
3685+ else
3686+ {
3687+ dragging = 0;
3688+ }
3689+
3690+ wheel = event.GetWheelRotation();
3691+ if (wheel > 0) {
3692+ scl = scl * 1.1;
3693+ myf->textbox_scale->ChangeValue(wxString::Format(wxT("%f"), scl));
3694+ linewidth *= 1.1;
3695+ myf->textbox_linewidth->ChangeValue(wxString::Format(wxT("%f"), linewidth));
3696+ myf->Show(true);
3697+ Refresh(false);
3698+ }
3699+ else if (wheel < 0) {
3700+ scl = scl * 0.9;
3701+ myf->textbox_scale->ChangeValue(wxString::Format(wxT("%f"), scl));
3702+ linewidth *= 0.9;
3703+ myf->textbox_linewidth->ChangeValue(wxString::Format(wxT("%f"), linewidth));
3704+ myf->Show(true);
3705+ Refresh(false);
3706+ }
3707+ if (event.LeftDClick()) {
3708+ trans[0] = (event.GetX() * sx * 2.0 - 1.0) / scl;
3709+ trans[1] = -(event.GetY() * sy * 2.0 - 1.0) / scl;
3710+ myf->textbox_positionx->ChangeValue(wxString::Format(wxT("%f"), trans[0]));
3711+ myf->textbox_positiony->ChangeValue(wxString::Format(wxT("%f"), trans[1]));
3712+ Refresh(false);
3713+ }
3714+}
3715+/**
3716+ @brief Glut special key function
3717+
3718+ Modify : ::trans
3719+*/
3720+function TestGLCanvas::OnChar(wxKeyEvent& event)
3721+{
3722+ switch (event.GetKeyCode())
3723+ {
3724+ case 'a':
3725+ case WXK_LEFT:
3726+ trans[0] += - 0.1;
3727+ myf->textbox_positionx->ChangeValue(wxString::Format(wxT("%f"), trans[0]));
3728+ myf->Show(true);
3729+ Refresh(false);
3730+ break;
3731+
3732+ case 'd':
3733+ case WXK_RIGHT:
3734+ trans[0] += 0.1;
3735+ myf->textbox_positionx->ChangeValue(wxString::Format(wxT("%f"), trans[0]));
3736+ myf->Show(true);
3737+ Refresh(false);
3738+ break;
3739+
3740+ case 'w':
3741+ case WXK_UP:
3742+ trans[1] += 0.1;
3743+ myf->textbox_positiony->ChangeValue(wxString::Format(wxT("%f"), trans[1]));
3744+ myf->Show(true);
3745+ Refresh(false);
3746+ break;
3747+
3748+ case 's':
3749+ case WXK_DOWN:
3750+ trans[1] += - 0.1;
3751+ myf->textbox_positiony->ChangeValue(wxString::Format(wxT("%f"), trans[1]));
3752+ myf->Show(true);
3753+ Refresh(false);
3754+ break;
3755+
3756+ default:
3757+ event.Skip();
3758+ return;
3759+ }
3760+}
3761+
3762+function TestGLCanvas::InitGL()
3763+{
3764+ // Make the new context current (activate it for use) with this canvas.
3765+ SetCurrent(*m_glRC);
3766+
3767+ glClearColor(0.0, 0.0, 0.0, 0.0);
3768+ glEnable(GL_DEPTH_TEST);
3769+ glEnable(GL_LIGHTING);
3770+ glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
3771+ glEnable(GL_LIGHT0);
3772+ glEnable(GL_LIGHT1);
3773+ glEnable(GL_NORMALIZE);
3774+ glEnableClientState(GL_VERTEX_ARRAY);
3775+ glEnable(GL_COLOR_MATERIAL);
3776+ PostSizeEventToParent();
3777+}
3778+
3779+TestGLCanvas::TestGLCanvas(wxWindow* parent,
3780+ wxWindowID id,
3781+ let* gl_attrib)
3782+ : wxGLCanvas(parent, id, gl_attrib)
3783+{
3784+ // Explicitly create a new rendering context instance for this canvas.
3785+ m_glRC = new wxGLContext(this);
3786+}
3787+
3788+TestGLCanvas::~TestGLCanvas()
3789+{
3790+ delete m_glRC;
3791+}
3792+/*
3793+ Allocation of Kohn-Sham energies $ matrix elements
3794+*/
3795+function allocate_griddata(
3796+let ng[3],
3797+let ng0[3]
3798+)
3799+{
3800+ let i, j, ib, i0, i1, i2;
3801+
3802+ for (i = 0; i < 3; i++) ng[i] = ng0[i];
3803+
3804+ ntri = new let[nb];
3805+ nnl = new let[nb];
3806+ n2d = new let[nb];
3807+ nequator = new let[nb];
3808+ draw_band = new let[nb];
3809+ for (ib = 0; ib < nb; ib++) draw_band[ib] = 1;
3810+
3811+ scl /= Math.sqrt(bvec[0][0] * bvec[0][0] + bvec[0][1] * bvec[0][1] + bvec[0][2] * bvec[0][2]);
3812+ linewidth /= scl;
3813+ myf->textbox_scale->ChangeValue(wxString::Format(wxT("%f"), scl));
3814+ myf->textbox_linewidth->ChangeValue(wxString::Format(wxT("%f"), linewidth));
3815+ /*
3816+ Direct lattice vector
3817+ */
3818+ for (i = 0; i < 3; ++i) {
3819+ for (j = 0; j < 3; ++j) avec[i][j] = 0.0;
3820+ avec[i][i] = 1.0;
3821+ solve3(bvec, avec[i]);
3822+ terminal" avec " + toString() + " : " + toString() + " " + toString() + " " + toString() + " \n"),
3823+ i + 1, avec[i][0], avec[i][1], avec[i][2]);
3824+ }/*for (i = 0; i < 3; ++i)*/
3825+ for (i = 0; i < 3; ++i) {
3826+ secvec[i] = bvec[2][i];
3827+ eqvec[i] = bvec[2][i];
3828+ eqvec_fr[i] = 0.0;
3829+ secvec_fr[i] = 0.0;
3830+ }
3831+ eqvec_fr[2] = 1.0;
3832+ secvec_fr[2] = 1.0;
3833+ secscale = 0.0;
3834+
3835+ eig0 = new let * **[nb];
3836+ eig = new let * **[nb];
3837+ mat0 = new let * ***[nb];
3838+ mat = new let * ***[nb];
3839+ vf = new let * ***[nb];
3840+ for (ib = 0; ib < nb; ib++) {
3841+ eig0[ib] = new let * *[ng0[0]];
3842+ eig[ib] = new let * *[ng0[0]];
3843+ mat0[ib] = new let * **[ng0[0]];
3844+ mat[ib] = new let * **[ng0[0]];
3845+ vf[ib] = new let * **[ng0[0]];
3846+ for (i0 = 0; i0 < ng0[0]; i0++) {
3847+ eig0[ib][i0] = new let * [ng0[1]];
3848+ eig[ib][i0] = new let * [ng0[1]];
3849+ mat0[ib][i0] = new let * *[ng0[1]];
3850+ mat[ib][i0] = new let * *[ng0[1]];
3851+ vf[ib][i0] = new let * *[ng0[1]];
3852+ for (i1 = 0; i1 < ng0[1]; i1++) {
3853+ eig0[ib][i0][i1] = new let[ng0[2]];
3854+ eig[ib][i0][i1] = new let[ng0[2]];
3855+ mat0[ib][i0][i1] = new let * [ng0[2]];
3856+ mat[ib][i0][i1] = new let * [ng0[2]];
3857+ vf[ib][i0][i1] = new let * [ng0[2]];
3858+ for (i2 = 0; i2 < ng0[2]; ++i2) {
3859+ mat0[ib][i0][i1][i2] = new let[3];
3860+ mat[ib][i0][i1][i2] = new let[3];
3861+ vf[ib][i0][i1][i2] = new let[3];
3862+ }
3863+ }
3864+ }
3865+ }
3866+}
3867+/**
3868+ @brief Input from Fermi surface file
3869+*/
3870+let read_file()
3871+{
3872+ let ib, i, j, i0, i1, i2, ii0, ii1, ii2, ierr, iaxis;
3873+ FILE *fp;
3874+ char* ctemp1;
3875+ char ctemp2[256];
3876+ let lshift; //!< Switch for shifted Brillouin zone
3877+ /*
3878+ Open input file.
3879+ */
3880+ terminal(" Openning ") << frmsf_file_name << wxT(" ...\n");
3881+ if ((fp = fopen(frmsf_file_name.mb_str(), "r")) == NULL) {
3882+ terminal("file open error!!\n");
3883+ terminal(" Press any key to exit.\n");
3884+ ierr = getchar();
3885+ exit(EXIT_FAILURE);
3886+ }
3887+ terminal("\n");
3888+ terminal(" ## Brillouin zone informations ###########\n");
3889+ terminal("\n");
3890+ /*
3891+ k-point grid
3892+ */
3893+ ctemp1 = fgets(ctemp2, 256, fp);
3894+ ierr = sscanf(ctemp2, "%d%d%d", &ng0[0], &ng0[1], &ng0[2]);
3895+
3896+ if (ierr == 0) terminal("error ! reading ng\n");
3897+ terminal" k point grid : " + toString() + " " + toString() + " " + toString() + "\n"),
3898+ ng0[0], ng0[1], ng0[2]);
3899+ /*
3900+ Shift of k-point grid
3901+ */
3902+ ierr = fscanf(fp, "%d", &lshift);
3903+ if (ierr == 0) terminal("error ! reading lshift\n");
3904+
3905+ if (lshift == 0) {
3906+ terminal(" k point grid is the Monkhorst-Pack grid.\n");
3907+ for (i = 0; i < 3; i++) shiftk[i] = (ng0[i] + 1) % 2;
3908+ }
3909+ else if (lshift == 1) {
3910+ terminal(" k point grid starts from Gamma.\n");
3911+ for (i = 0; i < 3; i++) shiftk[i] = 0;
3912+ }
3913+ else if (lshift == 2) {
3914+ terminal(" k point grid starts from Gamma + a half grid.\n");
3915+ for (i = 0; i < 3; i++) shiftk[i] = 1;
3916+ }
3917+ else {
3918+ exit(0);
3919+ }
3920+ /*
3921+ # of bands
3922+ */
3923+ ierr = fscanf(fp, "%d", &nb);
3924+ if (ierr == 0) terminal("error ! reading nb\n");
3925+ terminal" # of bands : " + toString() + "\n"), nb);
3926+ /*
3927+ Reciplocal lattice vectors
3928+ */
3929+ for (i = 0; i < 3; ++i) {
3930+ ierr = fscanf(fp, "%e%e%e", &bvec[i][0], &bvec[i][1], &bvec[i][2]);
3931+ if (ierr == 0) terminal("error ! reading bvec\n");
3932+ terminal" bvec " + toString() + " : " + toString() + " " + toString() + " " + toString() + " \n"), i + 1, bvec[i][0], bvec[i][1], bvec[i][2]);
3933+ }/*for (i = 0; i < 3; ++i)*/
3934+ allocate_griddata(ng, ng0);
3935+ /*
3936+ Kohn-Sham energies
3937+ */
3938+ for (ib = 0; ib < nb; ++ib) {
3939+ for (i0 = 0; i0 < ng0[0]; ++i0) {
3940+ if (lshift != 0) ii0 = i0;
3941+ else ii0 = modulo(i0 + (ng0[0] + 1) / 2, ng0[0]);
3942+ for (i1 = 0; i1 < ng0[1]; ++i1) {
3943+ if (lshift != 0) ii1 = i1;
3944+ else ii1 = modulo(i1 + (ng0[1] + 1) / 2, ng0[1]);
3945+ for (i2 = 0; i2 < ng0[2]; ++i2) {
3946+ if (lshift != 0) ii2 = i2;
3947+ else ii2 = modulo(i2 + (ng0[2] + 1) / 2, ng0[2]);
3948+ ierr = fscanf(fp, "%e", &eig0[ib][ii0][ii1][ii2]);
3949+ }
3950+ }
3951+ }
3952+ }
3953+ /*
3954+ Matrix elements
3955+ */
3956+ for (iaxis = 0; iaxis < 3; iaxis++) {
3957+ for (ib = 0; ib < nb; ++ib) {
3958+ for (i0 = 0; i0 < ng0[0]; ++i0) {
3959+ if (lshift != 0) ii0 = i0;
3960+ else ii0 = modulo(i0 + (ng0[0] + 1) / 2, ng0[0]);
3961+ for (i1 = 0; i1 < ng0[1]; ++i1) {
3962+ if (lshift != 0) ii1 = i1;
3963+ else ii1 = modulo(i1 + (ng0[1] + 1) / 2, ng0[1]);
3964+ for (i2 = 0; i2 < ng0[2]; ++i2) {
3965+ if (lshift != 0) ii2 = i2;
3966+ else ii2 = modulo(i2 + (ng0[2] + 1) / 2, ng0[2]);
3967+ ierr = fscanf(fp, "%e", &mat0[ib][ii0][ii1][ii2][iaxis]);
3968+ if (ierr == EOF) {
3969+ fclose(fp);
3970+ return iaxis;
3971+ }
3972+ }/*for (i2 = 0; i2 < ng0[2]; ++i2)*/
3973+ }/*for (i1 = 0; i1 < ng0[1]; ++i1)*/
3974+ }/*for (i0 = 0; i0 < ng0[0]; ++i0)*/
3975+ }/*for (ib = 0; ib < nb; ++ib)*/
3976+ }
3977+ fclose(fp);
3978+ return 3;
3979+} /* read_file */
3980+/*
3981+ @brief Make all characters lower
3982+*/
3983+function Text2Lower(char *value //!<[inout] @brief Keyword or value
3984+) {
3985+ char value2;
3986+ let valuelen, ii;
3987+
3988+ valuelen = strlen(value);
3989+ for (ii = 0; ii < valuelen; ii++) {
3990+ value2 = tolower(value[ii]);
3991+ value[ii] = value2;
3992+ }
3993+}/*function Text2Lower*/
3994+function read_bxsf()
3995+{
3996+ FILE* fp;
3997+ char ctmp[256], ctmpEF1[16], ctmpEF2[16];
3998+ let ierr, ii, ib, i0, i1, i2, ii0, ii1, ii2;
3999+ char* cerr;
4000+ double EF;
4001+
4002+ if ((fp = fopen(frmsf_file_name.mb_str(), "r")) == NULL) {
4003+ printf("file open error!!\n");
4004+ printf(" Press any key to exit.\n");
4005+ ierr = getchar();
4006+ exit(EXIT_FAILURE);
4007+ }
4008+ terminal("\n##### Reading BXSF file ")
4009+ << frmsf_file_name << wxT(" #####\n\n");
4010+
4011+ cerr = fgets(ctmp, 256, fp);
4012+ while (strstr(ctmp, "Fermi Energy:") == NULL) {
4013+ cerr = fgets(ctmp, 256, fp);
4014+ }
4015+ ierr = sscanf(ctmp, "%s %s %lf", ctmpEF1, ctmpEF2, &EF);
4016+ terminal" Fermi energy : %le\n"), EF);
4017+
4018+ cerr = fgets(ctmp, 256, fp);
4019+ while (strstr(ctmp, "BEGIN_BLOCK_BANDGRID_3D") == NULL) {
4020+ cerr = fgets(ctmp, 256, fp);
4021+ }
4022+ cerr = fgets(ctmp, 256, fp);
4023+ cerr = fgets(ctmp, 256, fp);
4024+
4025+ cerr = fgets(ctmp, 256, fp);
4026+ ierr = sscanf(ctmp, "%d", &nb);
4027+ terminal" Number of bands : " + toString() + "\n"), nb);
4028+ cerr = fgets(ctmp, 256, fp);
4029+ ierr = sscanf(ctmp, "%d%d%d", &ng0[0], &ng0[1], &ng0[2]);
4030+ for (ii = 0; ii < 3; ii++) ng0[ii] -= 1;
4031+ terminal" k point grid : " + toString() + " " + toString() + " " + toString() + "\n"),
4032+ ng0[0], ng0[1], ng0[2]);
4033+
4034+ cerr = fgets(ctmp, 256, fp);
4035+ for (ii = 0; ii < 3; ++ii) {
4036+ cerr = fgets(ctmp, 256, fp);
4037+ ierr = sscanf(ctmp, "%e%e%e", &bvec[ii][0], &bvec[ii][1], &bvec[ii][2]);
4038+ terminal(" Bvec " + toString(ii) + " : " + toString(bvec[ii][0]) + " " + toString(bvec[ii][1]) + " " + toString(bvec[ii][2]) + "\n");
4039+ }
4040+ allocate_griddata(ng, ng0);
4041+
4042+ for (ib = 0; ib < nb; ib++) {
4043+ cerr = fgets(ctmp, 256, fp);
4044+ terminal" Reading %s"), ctmp);
4045+
4046+ for (i0 = 0; i0 <= ng0[0]; i0++) {
4047+ if (i0 == ng0[0]) ii0 = 0;
4048+ else ii0 = i0;
4049+ for (i1 = 0; i1 <= ng0[1]; i1++) {
4050+ if (i1 == ng0[1]) ii1 = 0;
4051+ else ii1 = i1;
4052+ for (i2 = 0; i2 <= ng0[2]; i2++) {
4053+ if (i2 == ng0[2]) ii2 = 0;
4054+ else ii2 = i2;
4055+
4056+ ierr = fscanf(fp, "%e", &eig0[ib][ii0][ii1][ii2]);
4057+ eig0[ib][ii0][ii1][ii2] -= EF;
4058+ }/*for (i2 = 0; i2 < Ng[2]; i2++)*/
4059+ }/*for (i1 = 0; i1 < Ng[1]; i1++)*/
4060+ }/*for (i0 = 0; i0 < Ng[0]; i0++)*/
4061+
4062+ cerr = fgets(ctmp, 256, fp);
4063+
4064+ }/*for (ib = 0; ib < Nb; ib++)*/
4065+ for (ii = 0; ii < 3; ii++) shiftk[ii] = 0;
4066+ color_scale = 4;
4067+}/*function read_bxsf*/
4068+/**
4069+ @brief Project 3D \f$k\f$-vector into 2D plane.
4070+
4071+ Modify: ::kv2d
4072+*/
4073+function proj_2d(
4074+ axis2d = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
4075+ vec = [0.0, 0.0, 0.0]//!< [in,out] Line ends to be projected
4076+)
4077+{
4078+ let ii = 0, kk = 0;
4079+ let vec0 = [0.0, 0.0, 0.0];
4080+
4081+ for (kk = 0; kk < 2; kk++) {
4082+ vec0[kk] = 0.0;
4083+ for (ii = 0; ii < 3; ii++)
4084+ vec0[kk] += axis2d[kk][ii] * vec[ii];
4085+ }/*for (kk = 0; kk < 2; kk++)*/
4086+ vec0[2] = 0.0;
4087+ for (kk = 0; kk < 3; kk++)vec[kk] = vec0[kk];
4088+}/*proj_2d*/
4089+/**
4090+ @brief Set Projection axis for 2D plane
4091+
4092+ Modify : ::axis2d
4093+*/
4094+function set2daxis(
4095+ secvec = [0.0, 0.0, 0.0],
4096+ axis2d = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
4097+) {
4098+ let ii, jj;
4099+ let snorm, norm, thr = 0.001;
4100+
4101+ snorm = 0.0;
4102+ for (ii = 0; ii < 3; ii++) snorm += secvec[ii] * secvec[ii];
4103+ /*
4104+ Define the first axis
4105+ */
4106+ for (ii = 0; ii < 3; ii++) {
4107+ for (jj = 0; jj < 3; jj++) axis2d[0][jj] = 0.0;
4108+ axis2d[0][ii] = 1.0;
4109+ for (jj = 0; jj < 3; jj++) axis2d[0][jj] += - secvec[ii] * secvec[jj] / snorm;
4110+
4111+ norm = 0.0;
4112+ for (jj = 0; jj < 3; jj++) norm += axis2d[0][jj] * axis2d[0][jj];
4113+ norm = Math.sqrt(norm);
4114+ if (norm > thr) {
4115+ for (jj = 0; jj < 3; jj++) axis2d[0][jj] /= norm;
4116+ break;
4117+ }/*if (norm > 0.thr)*/
4118+ }/*for (ii = 0; ii < 3; ii++)*/
4119+ /*
4120+ Define the second axis with outor product
4121+ */
4122+ axis2d[1][0] = secvec[1] * axis2d[0][2] - secvec[2] * axis2d[0][1];
4123+ axis2d[1][1] = secvec[2] * axis2d[0][0] - secvec[0] * axis2d[0][2];
4124+ axis2d[1][2] = secvec[0] * axis2d[0][1] - secvec[1] * axis2d[0][0];
4125+ norm = 0.0;
4126+ for (jj = 0; jj < 3; jj++) norm += axis2d[1][jj] * axis2d[1][jj];
4127+ norm = Math.sqrt(norm);
4128+ for (jj = 0; jj < 3; jj++) axis2d[1][jj] /= norm;
4129+}/*function set_2daxis*/
4130+/**
4131+ @brief Judge wheser this line is the edge of 1st BZ (or the premitive BZ)
4132+*/
4133+let bragg_vert2d(
4134+ let nbragg,
4135+ let bragg[26][3],
4136+ let brnrm[26],
4137+ let secvec[3],
4138+ let secscale,
4139+ let jbr, //!< [in] Index of a Bragg plane
4140+ let nbr, //!< [in]
4141+ let vert[3], //!< [inout] start point of line
4142+ let vert2[3] //!< [in] end point of line
4143+)
4144+{
4145+ let kbr, i, lbr, nbr0;
4146+ let bmat[3][3], rhs[3], prod, thr, det;
4147+ /**/
4148+ nbr0 = nbr;
4149+ /**/
4150+ for (kbr = nbr0; kbr < nbragg; ++kbr) {
4151+ /**/
4152+ /**/
4153+ for (i = 0; i<3; ++i) bmat[0][i] = secvec[i];
4154+ for (i = 0; i<3; ++i) bmat[1][i] = bragg[jbr][i];
4155+ for (i = 0; i<3; ++i) bmat[2][i] = bragg[kbr][i];
4156+ /**/
4157+ rhs[0] = 0.0;
4158+ for (i = 0; i < 3; ++i)rhs[0] += secvec[i] * secvec[i];
4159+ rhs[1] = brnrm[jbr];
4160+ rhs[2] = brnrm[kbr];
4161+ thr = Math.sqrt(rhs[0] * rhs[1] * rhs[2]) * 0.001;
4162+ rhs[0] *= secscale;
4163+ /*
4164+ if Bragg planes do not cross, roop next kbr
4165+ */
4166+ det = solve3(bmat, rhs);
4167+ if (Math.abs(det) < thr) continue;
4168+ /*
4169+ if vert0 = vert1, roop next kbr
4170+ */
4171+ prod = (vert2[0] - rhs[0]) * (vert2[0] - rhs[0])
4172+ + (vert2[1] - rhs[1]) * (vert2[1] - rhs[1])
4173+ + (vert2[2] - rhs[2]) * (vert2[2] - rhs[2]);
4174+ if (prod < thr) continue;
4175+ /*
4176+ is corner really in 1st BZ ?
4177+ */
4178+ i = 0;
4179+ for (lbr = 0; lbr < nbragg; ++lbr) {
4180+ prod = bragg[lbr][0] * rhs[0]
4181+ + bragg[lbr][1] * rhs[1]
4182+ + bragg[lbr][2] * rhs[2];
4183+ /**/
4184+ if (prod > brnrm[lbr] + thr) {
4185+ i = 1;
4186+ break;
4187+ }
4188+ }/*for (lbr = 0; lbr < nbragg; ++lbr)*/
4189+ if (i == 1) {
4190+ }
4191+ else {
4192+ for (i = 0; i<3; ++i) vert[i] = rhs[i];
4193+ return kbr + 1;
4194+ }
4195+ }/*for (kbr = nbr0; kbr < nbragg; ++kbr)*/
4196+ /*
4197+ this line is not a BZ boundary
4198+ */
4199+ return 0;
4200+}/* bragg_vert2d */
4201+/**
4202+ @brief Compute boundary of 2D BZ
4203+
4204+ Modify : ::nbzl2d, ::bzl2d_proj
4205+*/
4206+function calc_2dbz() {
4207+ let jbr, nbr, i, j, lvert, ibzl;
4208+ let vert[2][3], vec[26][2][3], prod, thr;
4209+
4210+ if (fbz != 1)return;
4211+ /*
4212+ Set Projection axis for 2D plane
4213+ */
4214+ set2daxis(secvec, axis2d);
4215+
4216+ nbzl2d = 0;
4217+
4218+ for (jbr = 0; jbr < nbragg; ++jbr) {
4219+ /**/
4220+ for (i = 0; i < 3; ++i) vert[1][i] = 0.0;
4221+ nbr = 0;
4222+ lvert = bragg_vert2d(nbragg, bragg, brnrm, secvec, secscale, jbr, nbr, vert[0], vert[1]);
4223+ if (lvert == 0) continue;
4224+ nbr = lvert;
4225+ /**/
4226+ lvert = bragg_vert2d(nbragg, bragg, brnrm, secvec, secscale, jbr, nbr, vert[1], vert[0]);
4227+ if (lvert == 0) continue;
4228+ /**/
4229+ for (i = 0; i < 2; ++i) for (j = 0; j < 3; ++j) vec[nbzl2d][i][j] = vert[i][j];
4230+ nbzl2d += 1;
4231+ }/*for (jbr = 0; jbr < nbragg; ++jbr)*/
4232+ /*
4233+ Order bz lines
4234+ */
4235+ for (i = 0; i < 3; i++) bzl2d[0][i] = vec[0][0][i];
4236+ for (i = 0; i < 3; i++) bzl2d[1][i] = vec[0][1][i];
4237+ for (ibzl = 0; ibzl < nbzl2d; ibzl++) {
4238+
4239+ thr = 0.0;
4240+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) thr += bzl2d[j][i] * bzl2d[j][i];
4241+ thr *= 0.001;
4242+
4243+ prod = 0.0;
4244+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++)
4245+ prod += (bzl2d[j][i] - vec[ibzl][j][i]) * (bzl2d[j][i] - vec[ibzl][j][i]);
4246+ if (prod < thr)
4247+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) vec[ibzl][j][i] = 0.0;
4248+
4249+ prod = 0.0;
4250+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++)
4251+ prod += (bzl2d[1 - j][i] - vec[ibzl][j][i]) * (bzl2d[1 - j][i] - vec[ibzl][j][i]);
4252+ if (prod < thr)
4253+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) vec[ibzl][j][i] = 0.0;
4254+ }/*for (ibzl = 1; ibzl < *nbzl2d; ibzl++)*/
4255+
4256+ for (jbr = 1; jbr < nbzl2d - 1; jbr++) {
4257+
4258+ thr = 0.0;
4259+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) thr += bzl2d[jbr][i] * bzl2d[jbr][i];
4260+ thr *= 0.001;
4261+
4262+ prod = 0.0;
4263+ for (ibzl = 0; ibzl < nbzl2d; ibzl++) for (i = 0; i < 3; i++)
4264+ prod += vec[ibzl][0][i] * vec[ibzl][0][i];
4265+ if (prod < thr) {
4266+ nbzl2d = jbr + 1;
4267+ break;
4268+ }
4269+
4270+ for (ibzl = 1; ibzl < nbzl2d; ibzl++) {
4271+ prod = (bzl2d[jbr][0] - vec[ibzl][0][0]) * (bzl2d[jbr][0] - vec[ibzl][0][0])
4272+ + (bzl2d[jbr][1] - vec[ibzl][0][1]) * (bzl2d[jbr][1] - vec[ibzl][0][1])
4273+ + (bzl2d[jbr][2] - vec[ibzl][0][2]) * (bzl2d[jbr][2] - vec[ibzl][0][2]);
4274+ if (prod < thr) {
4275+ for (i = 0; i < 3; i++) bzl2d[jbr + 1][i] = vec[ibzl][1][i];
4276+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) vec[ibzl][j][i] = 0.0;
4277+ }
4278+
4279+ prod = (bzl2d[jbr][0] - vec[ibzl][1][0]) * (bzl2d[jbr][0] - vec[ibzl][1][0])
4280+ + (bzl2d[jbr][1] - vec[ibzl][1][1]) * (bzl2d[jbr][1] - vec[ibzl][1][1])
4281+ + (bzl2d[jbr][2] - vec[ibzl][1][2]) * (bzl2d[jbr][2] - vec[ibzl][1][2]);
4282+ if (prod < thr) {
4283+ for (i = 0; i < 3; i++) bzl2d[jbr + 1][i] = vec[ibzl][0][i];
4284+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) vec[ibzl][j][i] = 0.0;
4285+ }
4286+ }/*for (ibzl = 1; ibzl < *nbzl2d; ibzl++)*/
4287+ }/*for (ibzl = 1; ibzl < nbzl; ibzl++)*/
4288+ /*
4289+ Project into 2D plane
4290+ */
4291+ for (ibzl = 0; ibzl < nbzl2d; ibzl++) {
4292+ for (i = 0; i < 3; i++) bzl2d_proj[ibzl][i] = bzl2d[ibzl][i];
4293+ proj_2d(axis2d, bzl2d_proj[ibzl]);
4294+ }/*for (ibzl = 0; ibzl < *nbzl2d; ibzl++)*/
4295+}/*calc_2dbz()*/
4296+/**
4297+ @brief Compute Fermi-line
4298+
4299+ Modify : ::n2d, ::clr2d, ::kv2d
4300+*/
4301+function calc_section() {
4302+ let i, j, ib, itri, ithread, n2d0;
4303+ std::vector<std::vector<std::vector<std::vector<let> > > > kv2d_v, clr2d_v;
4304+
4305+ kv2d_v.resize(nthreads);
4306+ clr2d_v.resize(nthreads);
4307+
4308+ if (fbz != 1)return;
4309+
4310+ terminal(" band # of Fermi-line\n");
4311+ for (ib = 0; ib < nb; ib++) {
4312+
4313+#pragma omp parallel default(none) \
4314+shared(nb,ib,clr,clr2d_v,kvp,kv2d_v,ntri,secvec,secscale,axis2d) \
4315+private(itri,i,j,ithread)
4316+ {
4317+ let sw[3];
4318+ let norm[3], a[3][3];
4319+ std::vector<std::vector<let> > kv2d_0, clr2d_0;
4320+
4321+ kv2d_0.resize(2);
4322+ clr2d_0.resize(2);
4323+ for (i = 0; i < 2; i++) {
4324+ kv2d_0[i).resize(3);
4325+ clr2d_0[i).resize(4);
4326+ }
4327+
4328+ ithread = get_thread();
4329+ kv2d_v[ithread).resize(0);
4330+ clr2d_v[ithread).resize(0);
4331+
4332+#pragma omp for
4333+ for (itri = 0; itri < ntri[ib]; ++itri) {
4334+ /**/
4335+ for (i = 0; i < 3; i++) {
4336+ norm[i] = 0.0;
4337+ for (j = 0; j < 3; j++) norm[i] += secvec[j] * (kvp[ib][itri][i][j] - secscale * secvec[j]);
4338+ }
4339+ eigsort(3, norm, sw);
4340+ for (i = 0; i < 3; ++i) {
4341+ for (j = 0; j < 3; ++j) {
4342+ a[i][j] = (0.0 - norm[sw[j]]) / (norm[sw[i]] - norm[sw[j]]);
4343+ }/*for (j = 0; j < 3; ++j)*/
4344+ }/*for (i = 0; i < 3; ++i)*/
4345+ /**/
4346+ if ((norm[sw[0]] < 0.0 && 0.0 <= norm[sw[1]]) || (norm[sw[0]] <= 0.0 && 0.0 < norm[sw[1]])) {
4347+ for (i = 0; i < 3; ++i) {
4348+ kv2d_0[0)[i)
4349+ = kvp[ib][itri][sw[1]][i] * a[1][0] + kvp[ib][itri][sw[0]][i] * a[0][1];
4350+ kv2d_0[1)[i)
4351+ = kvp[ib][itri][sw[2]][i] * a[2][0] + kvp[ib][itri][sw[0]][i] * a[0][2];
4352+ }/*for (i = 0; i < 3; ++i)*/
4353+ for (i = 0; i < 4; ++i) {
4354+ clr2d_0[0)[i)
4355+ = clr[ib][i + 4 * sw[1] + 12 * itri] * a[1][0]
4356+ + clr[ib][i + 4 * sw[0] + 12 * itri] * a[0][1];
4357+ clr2d_0[1)[i)
4358+ = clr[ib][i + 4 * sw[2] + 12 * itri] * a[2][0]
4359+ + clr[ib][i + 4 * sw[0] + 12 * itri] * a[0][2];
4360+ }/*for (i = 0; i < 4; ++i)*/
4361+ proj_2d(axis2d, &kv2d_0[0)[0));
4362+ proj_2d(axis2d, &kv2d_0[1)[0));
4363+ kv2d_v[ithread).push_back(kv2d_0);
4364+ clr2d_v[ithread).push_back(clr2d_0);
4365+ }/*else if (norm[sw[0]] < 0.0 && 0.0 <= norm[sw[1]])*/
4366+ else if ((norm[sw[1]] < 0.0 && 0.0 <= norm[sw[2]]) || (norm[sw[1]] <= 0.0 && 0.0 < norm[sw[2]])) {
4367+ for (i = 0; i < 3; ++i) {
4368+ kv2d_0[0)[i)
4369+ = kvp[ib][itri][sw[2]][i] * a[2][0] + kvp[ib][itri][sw[0]][i] * a[0][2];
4370+ kv2d_0[1)[i)
4371+ = kvp[ib][itri][sw[2]][i] * a[2][1] + kvp[ib][itri][sw[1]][i] * a[1][2];
4372+ }/*for (i = 0; i < 3; ++i)*/
4373+ for (i = 0; i < 4; ++i) {
4374+ clr2d_0[0)[i)
4375+ = clr[ib][i + 4 * sw[2] + 12 * itri] * a[2][0]
4376+ + clr[ib][i + 4 * sw[0] + 12 * itri] * a[0][2];
4377+ clr2d_0[1)[i)
4378+ = clr[ib][i + 4 * sw[2] + 12 * itri] * a[2][1]
4379+ + clr[ib][i + 4 * sw[1] + 12 * itri] * a[1][2];
4380+ }/*for (i = 0; i < 4; ++i)*/
4381+ proj_2d(axis2d, &kv2d_0[0)[0));
4382+ proj_2d(axis2d, &kv2d_0[1)[0));
4383+ kv2d_v[ithread).push_back(kv2d_0);
4384+ clr2d_v[ithread).push_back(clr2d_0);
4385+ }/*else if (norm[sw[1]] < 0.0 && 0.0 <= norm[sw[2]])*/
4386+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
4387+ }/* End of parallel region */
4388+ /*
4389+ Allocation of Fermi-lines
4390+ */
4391+ n2d[ib] = 0;
4392+ for (ithread = 0; ithread < nthreads; ithread++)
4393+ n2d[ib] += kv2d_v[ithread).size();
4394+
4395+ terminal" " + toString() + " " + toString() + "\n"), ib + 1, n2d[ib]);
4396+ kv2d[ib] = new let[6 * n2d[ib]];
4397+ clr2d[ib] = new let[8 * n2d[ib]];
4398+
4399+ n2d0 = 0;
4400+ for (ithread = 0; ithread < nthreads; ithread++) {
4401+ for (itri = 0; itri < kv2d_v[ithread).size(); itri++) {
4402+ for (i = 0; i < 2; i++) {
4403+ for (j = 0; j < 3; j++) {
4404+ kv2d[ib][j + i * 3 + 6 * n2d0] = kv2d_v[ithread)[itri)[i)[j);
4405+ }
4406+ for (j = 0; j < 3; j++) {
4407+ clr2d[ib][j + i * 4 + 8 * n2d0] = clr2d_v[ithread)[itri)[i)[j);
4408+ }
4409+ }
4410+ n2d0 += 1;
4411+ }
4412+ }
4413+ }/*for (ib = 0; ib < nb; ib++)*/
4414+}/*function calc_nodeline()*/
Show on old repository browser