• R/O
  • HTTP
  • SSH
  • HTTPS

android: Commit

FermiSurfer for android


Commit MetaInfo

Revisionf55d49029873285aca290c1407102ea3e81a3637 (tree)
Zeit2019-08-31 18:28:42
Autormitsuaki1987 <kawamitsuaki@gmai...>
Commitermitsuaki1987

Log Message

Backup

Ändern Zusammenfassung

Diff

--- a/Android1.NativeActivity/Android1.NativeActivity.vcxproj
+++ b/Android1.NativeActivity/Android1.NativeActivity.vcxproj
@@ -149,7 +149,7 @@
149149 <CompileAs>CompileAsCpp</CompileAs>
150150 </ClCompile>
151151 <Link>
152- <LibraryDependencies>%(LibraryDependencies);GLESv1_CM;EGL;</LibraryDependencies>
152+ <LibraryDependencies>%(LibraryDependencies);GLESv1_CM;EGL;m</LibraryDependencies>
153153 </Link>
154154 </ItemDefinitionGroup>
155155 <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|ARM'">
@@ -204,10 +204,35 @@
204204 </ItemDefinitionGroup>
205205 <ItemGroup>
206206 <ClInclude Include="android_native_app_glue.h" />
207+ <ClInclude Include="basic_math.hpp" />
208+ <ClInclude Include="bz_lines.hpp" />
209+ <ClInclude Include="calc_nodeline.hpp" />
210+ <ClInclude Include="draw.hpp" />
211+ <ClInclude Include="equator.hpp" />
212+ <ClInclude Include="fermi_patch.hpp" />
213+ <ClInclude Include="free_patch.hpp" />
214+ <ClInclude Include="initialize.hpp" />
215+ <ClInclude Include="kumo.hpp" />
216+ <ClInclude Include="menu.hpp" />
217+ <ClInclude Include="read_file.hpp" />
218+ <ClInclude Include="section.hpp" />
219+ <ClInclude Include="variable.hpp" />
207220 </ItemGroup>
208221 <ItemGroup>
209222 <ClCompile Include="android_native_app_glue.c" />
223+ <ClCompile Include="basic_math.cpp" />
224+ <ClCompile Include="bz_lines.cpp" />
225+ <ClCompile Include="calc_nodeline.cpp" />
226+ <ClCompile Include="draw.cpp" />
227+ <ClCompile Include="equator.cpp" />
228+ <ClCompile Include="fermi_patch.cpp" />
229+ <ClCompile Include="free_patch.cpp" />
230+ <ClCompile Include="initialize.cpp" />
231+ <ClCompile Include="kumo.cpp" />
210232 <ClCompile Include="main.cpp" />
233+ <ClCompile Include="menu.cpp" />
234+ <ClCompile Include="read_file.cpp" />
235+ <ClCompile Include="section.cpp" />
211236 </ItemGroup>
212237 <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
213238 <ImportGroup Label="ExtensionTargets">
--- a/Android1.NativeActivity/Android1.NativeActivity.vcxproj.filters
+++ b/Android1.NativeActivity/Android1.NativeActivity.vcxproj.filters
@@ -2,9 +2,34 @@
22 <Project ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
33 <ItemGroup>
44 <ClInclude Include="android_native_app_glue.h" />
5+ <ClInclude Include="read_file.hpp" />
6+ <ClInclude Include="variable.hpp" />
7+ <ClInclude Include="basic_math.hpp" />
8+ <ClInclude Include="kumo.hpp" />
9+ <ClInclude Include="bz_lines.hpp" />
10+ <ClInclude Include="calc_nodeline.hpp" />
11+ <ClInclude Include="draw.hpp" />
12+ <ClInclude Include="equator.hpp" />
13+ <ClInclude Include="fermi_patch.hpp" />
14+ <ClInclude Include="free_patch.hpp" />
15+ <ClInclude Include="initialize.hpp" />
16+ <ClInclude Include="menu.hpp" />
17+ <ClInclude Include="section.hpp" />
518 </ItemGroup>
619 <ItemGroup>
720 <ClCompile Include="android_native_app_glue.c" />
821 <ClCompile Include="main.cpp" />
22+ <ClCompile Include="read_file.cpp" />
23+ <ClCompile Include="basic_math.cpp" />
24+ <ClCompile Include="kumo.cpp" />
25+ <ClCompile Include="bz_lines.cpp" />
26+ <ClCompile Include="calc_nodeline.cpp" />
27+ <ClCompile Include="draw.cpp" />
28+ <ClCompile Include="equator.cpp" />
29+ <ClCompile Include="fermi_patch.cpp" />
30+ <ClCompile Include="free_patch.cpp" />
31+ <ClCompile Include="initialize.cpp" />
32+ <ClCompile Include="menu.cpp" />
33+ <ClCompile Include="section.cpp" />
934 </ItemGroup>
1035 </Project>
\ No newline at end of file
--- /dev/null
+++ b/Android1.NativeActivity/basic_math.cpp
@@ -0,0 +1,139 @@
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+#include <GLES/gl.h>
28+#include <math.h>
29+
30+#if defined(_OPENMP)
31+#include <omp.h>
32+#endif
33+
34+/**
35+ @brief Work as Modulo function of fortran
36+ @return Modulated value
37+*/
38+int modulo(
39+ int i, //!< [in]
40+ int n //!< [in]
41+) {
42+ int j;
43+ j = (i + 100 * n) % n;
44+ return j;
45+}/*modulo(int i, int n)*/
46+/**
47+ @brief Solve linear system
48+ @return Determinant
49+*/
50+GLfloat solve3(
51+ GLfloat a[3][3], //!< [in] Matix
52+ GLfloat b[3] //!< [in,out] Right hand side vector
53+)
54+{
55+ int i;
56+ GLfloat det, c[3];
57+ /**/
58+ det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
59+ + a[0][1] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
60+ + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
61+ /**/
62+ c[0] = b[0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
63+ + b[1] * (a[0][2] * a[2][1] - a[0][1] * a[2][2])
64+ + b[2] * (a[0][1] * a[1][2] - a[0][2] * a[1][1]);
65+ /**/
66+ c[1] = b[0] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
67+ + b[1] * (a[0][0] * a[2][2] - a[0][2] * a[2][0])
68+ + b[2] * (a[0][2] * a[1][0] - a[0][0] * a[1][2]);
69+ /**/
70+ c[2] = b[0] * (a[1][0] * a[2][1] - a[1][1] * a[2][0])
71+ + b[1] * (a[0][1] * a[2][0] - a[0][0] * a[2][1])
72+ + b[2] * (a[0][0] * a[1][1] - a[0][1] * a[1][0]);
73+ /**/
74+ for (i = 0; i<3; ++i) b[i] = c[i] / det;
75+ return det;
76+}/*GLfloat solve3*/
77+/**
78+ @brief Simple sort
79+*/
80+void eigsort(
81+ int n, //!< [in] the number of components
82+ GLfloat *key, //!< [in] Variables to be sorted [n].
83+ int *swap //!< [out] Order of index (sorted)
84+)
85+{
86+ int i, j, k;
87+
88+ for (i = 0; i < n; ++i) swap[i] = i;
89+
90+ for (i = 0; i < n - 1; ++i) {
91+ for (j = i + 1; j < n; ++j) {
92+ if (key[swap[j]] < key[swap[i]]) {
93+ /*
94+ Swap
95+ */
96+ k = swap[j];
97+ swap[j] = swap[i];
98+ swap[i] = k;
99+ }/*if (sortee[j][0] < sortee[i][0])*/
100+ }/*for (j = i + 1; j < n; ++j)*/
101+ }/*for (i = 0; i < n - 1; ++i)*/
102+}/*eigsort*/
103+/**
104+ @brief Calculate normal vector from corners of triangle
105+*/
106+void normal_vec(
107+ GLfloat in1[3], //!< [in] Corner 1
108+ GLfloat in2[3], //!< [in] Corner 2
109+ GLfloat in3[3], //!< [in] Corner 3
110+ GLfloat out[3] //!< [out] The normal vector
111+)
112+{
113+ int i;
114+ GLfloat norm;
115+ out[0] = in1[1] * in2[2] - in1[2] * in2[1]
116+ + in2[1] * in3[2] - in2[2] * in3[1]
117+ + in3[1] * in1[2] - in3[2] * in1[1];
118+ out[1] = in1[2] * in2[0] - in1[0] * in2[2]
119+ + in2[2] * in3[0] - in2[0] * in3[2]
120+ + in3[2] * in1[0] - in3[0] * in1[2];
121+ out[2] = in1[0] * in2[1] - in1[1] * in2[0]
122+ + in2[0] * in3[1] - in2[1] * in3[0]
123+ + in3[0] * in1[1] - in3[1] * in1[0];
124+ norm = sqrtf(out[0] * out[0] + out[1] * out[1] + out[2] * out[2]);
125+ for (i = 0; i<3; i++) out[i] = out[i] / norm;
126+} /* normal_vec */
127+/**
128+ @brief OpenMP wrapper, get the number of threads
129+ @return the number of threads
130+*/
131+int get_thread() {
132+ int ithread;
133+#if defined(_OPENMP)
134+ ithread = omp_get_thread_num();
135+#else
136+ ithread = 0;
137+#endif
138+ return ithread;
139+}
--- /dev/null
+++ b/Android1.NativeActivity/basic_math.hpp
@@ -0,0 +1,30 @@
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+#include <GLES/gl.h>
25+
26+int modulo(int i, int n);
27+GLfloat solve3(GLfloat a[3][3], GLfloat b[3]);
28+void eigsort(int n, GLfloat *key, int *swap);
29+void normal_vec(GLfloat in1[3], GLfloat in2[3], GLfloat in3[3], GLfloat out[3]);
30+int get_thread();
--- /dev/null
+++ b/Android1.NativeActivity/bz_lines.cpp
@@ -0,0 +1,178 @@
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 Compute lines of BZ boundary
26+*/
27+#include <GLES/gl.h>
28+#include <stdlib.h>
29+#include <stdio.h>
30+#include <math.h>
31+#include "basic_math.hpp"
32+#include "variable.hpp"
33+
34+/**
35+ @brief Judge wheser this line is the edge of 1st BZ
36+*/
37+static int bragg_vert(
38+ int ibr, //!< [in] Index of a Bragg plane
39+ int jbr, //!< [in] Index of a Bragg plane
40+ int nbr, //!< [in]
41+ GLfloat vert[3], //!< [in] start point of line
42+ GLfloat vert2[3] //!< [in] end point of line
43+)
44+{
45+ int kbr, i, lbr, nbr0;
46+ GLfloat bmat[3][3], rhs[3], prod, thr, det;
47+ /**/
48+ nbr0 = nbr;
49+ /**/
50+ for (kbr = nbr0; kbr < 26; ++kbr) {
51+ /**/
52+ for (i = 0; i<3; ++i) bmat[0][i] = bragg[ibr][i];
53+ for (i = 0; i<3; ++i) bmat[1][i] = bragg[jbr][i];
54+ for (i = 0; i<3; ++i) bmat[2][i] = bragg[kbr][i];
55+ /**/
56+ rhs[0] = brnrm[ibr];
57+ rhs[1] = brnrm[jbr];
58+ rhs[2] = brnrm[kbr];
59+ thr = sqrtf(rhs[0] * rhs[1] * rhs[2]) * 0.001f;
60+ /*
61+ if Bragg planes do not cross, roop next kbr
62+ */
63+ det = solve3(bmat, rhs);
64+ if (fabsf(det) < thr) continue;
65+ /*
66+ if vert0 = vert1, roop next kbr
67+ */
68+ prod = (vert2[0] - rhs[0]) * (vert2[0] - rhs[0])
69+ + (vert2[1] - rhs[1]) * (vert2[1] - rhs[1])
70+ + (vert2[2] - rhs[2]) * (vert2[2] - rhs[2]);
71+ if (prod < thr) continue;
72+ /*
73+ is corner really in 1st BZ ?
74+ */
75+ i = 0;
76+ for (lbr = 0; lbr < 26; ++lbr) {
77+ prod = bragg[lbr][0] * rhs[0]
78+ + bragg[lbr][1] * rhs[1]
79+ + bragg[lbr][2] * rhs[2];
80+ /**/
81+ if (prod > brnrm[lbr] + thr) {
82+ i = 1;
83+ break;
84+ }
85+ }
86+ if (i == 1) {
87+ }
88+ else {
89+ for (i = 0; i<3; ++i) vert[i] = rhs[i];
90+ return kbr + 1;
91+ }
92+ }
93+ /*
94+ this line is not a BZ boundary
95+ */
96+ return 0;
97+}/* bragg_vert */
98+/**
99+ @brief Compute real number of Bragg plane at 1st BZ
100+*/
101+static void check_bragg() {
102+ int ibr, ibzl, ibzc;
103+ int ii, kk, bzflag, nbzcorner, nn;
104+ GLfloat thr = (GLfloat)0.0001, prod, bzc[676][3];
105+ /*
106+ First, compute real number of corners of 1st BZ
107+ */
108+ nbzcorner = 0;
109+ for (ibzl = 0; ibzl < nbzl; ibzl++) {
110+ for (ii = 0; ii < 2; ii++) {
111+ bzflag = 0;
112+
113+ for (ibzc = 0; ibzc < nbzcorner; ibzc++) {
114+ prod = 0.0f;
115+ for (kk = 0; kk < 3; kk++) prod += (bzl[ibzl][ii][kk] - bzc[ibzc][kk]) * (bzl[ibzl][ii][kk] - bzc[ibzc][kk]);
116+ if (prod < thr) bzflag = 1;
117+ }
118+
119+ if (bzflag == 0) {
120+ for (kk = 0; kk < 3; kk++) bzc[nbzcorner][kk] = bzl[ibzl][ii][kk];
121+ nbzcorner += 1;
122+ }
123+
124+ }/*for (ii = 0; ii < 2; ii++)*/
125+ }/*for (ibzl = 0; ibzl < nbzl; ibzl++)*/
126+ /**@brief
127+ Then, compute real number Bragg plane of 1st BZ (::nbragg),
128+ Re-order ::bragg and ::brnrm
129+ */
130+ nbragg = 0;
131+ for (ibr = 0; ibr < 26; ibr++) {
132+ nn = 0;
133+
134+ for (ibzc = 0; ibzc < nbzcorner; ibzc++) {
135+ prod = bragg[ibr][0] * bzc[ibzc][0] + bragg[ibr][1] * bzc[ibzc][1] + bragg[ibr][2] * bzc[ibzc][2];
136+ if (fabsf(prod - brnrm[ibr]) < thr) nn += 1;
137+ }
138+ if (nn >= 3) {
139+ for (kk = 0; kk < 3; kk++) bragg[nbragg][kk] = bragg[ibr][kk];
140+ brnrm[nbragg] = brnrm[ibr];
141+ nbragg += 1;
142+ }
143+ }
144+
145+}/*static void check_bragg*/
146+/**
147+ @brief Compute Brillouin zone boundariy lines
148+
149+ Modify : ::nbzl, ::bzl
150+*/
151+void bz_lines() {
152+ /**/
153+ int ibr, jbr, nbr, i, j, lvert;
154+ GLfloat vert[2][3];
155+ /**/
156+ nbzl = 0;
157+ /**/
158+ for (ibr = 0; ibr < 26; ++ibr) {
159+ for (jbr = 0; jbr < 26; ++jbr) {
160+ /**/
161+ for (i = 0; i<3; ++i) vert[1][i] = 0.0;
162+ nbr = 0;
163+ lvert = bragg_vert(ibr, jbr, nbr, vert[0], vert[1]);
164+ if (lvert == 0) continue;
165+ nbr = lvert;
166+ /**/
167+ lvert = bragg_vert(ibr, jbr, nbr, vert[1], vert[0]);
168+ if (lvert == 0) continue;
169+ /**/
170+ for (i = 0; i < 2; ++i) for (j = 0; j < 3; ++j) bzl[nbzl][i][j] = vert[i][j];
171+ nbzl = nbzl + 1;
172+
173+ }/*for (jbr = 0; jbr < 26; ++jbr)*/
174+ }/*for (ibr = 0; ibr < 26; ++ibr)*/
175+
176+ check_bragg();
177+
178+}/*bz_lines*/
--- /dev/null
+++ b/Android1.NativeActivity/bz_lines.hpp
@@ -0,0 +1,25 @@
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+
25+void bz_lines();
--- /dev/null
+++ b/Android1.NativeActivity/calc_nodeline.cpp
@@ -0,0 +1,127 @@
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 Compute nodal lines
26+*/
27+#include <GLES/gl.h>
28+
29+#include <stdio.h>
30+#include <stdlib.h>
31+#include <math.h>
32+#include "basic_math.hpp"
33+#include "variable.hpp"
34+
35+/**
36+@brief Compute node-line where \f$\Delta_{n k} = 0\f$
37+
38+ Modify : ::ntri_th, ::nnl, ::kvnl, ::kvnl_rot
39+
40+ If ::query = 1, this routine only compute the number of
41+ line segmants and malloc variables.
42+*/
43+void calc_nodeline() {
44+ int ib, itri, i, j, ithread;
45+
46+#pragma omp parallel default(none) \
47+ shared(nb,nnl,matp,kvnl,kvp,ntri,ntri_th,query) \
48+ private(ib,itri,i,j,ithread)
49+ {
50+ int sw[3], nnl0;
51+ GLfloat a[3][3], matps[3];
52+
53+ ithread = get_thread();
54+ for (ib = 0; ib < nb; ib++) {
55+ if(query == 1) nnl0 = 0;
56+ else nnl0 = ntri_th[ib][ithread];
57+#pragma omp for
58+ for (itri = 0; itri < ntri[ib]; ++itri) {
59+
60+ for (i = 0; i < 3; ++i) matps[i] = matp[ib][itri][i][0];
61+ eigsort(3, matps, sw);
62+
63+ for (i = 0; i < 3; ++i) {
64+ for (j = 0; j < 3; ++j) {
65+ a[i][j] = (0.f - matp[ib][itri][sw[j]][0])
66+ / (matp[ib][itri][sw[i]][0] - matp[ib][itri][sw[j]][0]);
67+ }/*for (j = 0; j < 3; ++j)*/
68+ }/*for (i = 0; i < 3; ++i)*/
69+
70+ if (( matp[ib][itri][sw[0]][0] < 0.0 && 0.0 <= matp[ib][itri][sw[1]][0])
71+ || (matp[ib][itri][sw[0]][0] <= 0.0 && 0.0 < matp[ib][itri][sw[1]][0])) {
72+ if (query == 0) {
73+ for (i = 0; i < 3; ++i) {
74+ kvnl[ib][nnl0][0][i] = kvp[ib][itri][sw[0]][i] * a[0][1] + kvp[ib][itri][sw[1]][i] * a[1][0];
75+ kvnl[ib][nnl0][1][i] = kvp[ib][itri][sw[0]][i] * a[0][2] + kvp[ib][itri][sw[2]][i] * a[2][0];
76+ }/*for (i = 0; i < 3; ++i)*/
77+ }/*if (query == 0)*/
78+ nnl0 += 1;
79+ }/*else if (matp[ib][itri][sw[0]] < 0.0 && 0.0 <= matp[ib][itri][sw[1]])*/
80+ else if ((matp[ib][itri][sw[1]][0] < 0.0 && 0.0 <= matp[ib][itri][sw[2]][0])
81+ || (matp[ib][itri][sw[1]][0] <= 0.0 && 0.0 < matp[ib][itri][sw[2]][0])) {
82+ if (query == 0) {
83+ for (i = 0; i < 3; ++i) {
84+ kvnl[ib][nnl0][0][i] = kvp[ib][itri][sw[0]][i] * a[0][2] + kvp[ib][itri][sw[2]][i] * a[2][0];
85+ kvnl[ib][nnl0][1][i] = kvp[ib][itri][sw[1]][i] * a[1][2] + kvp[ib][itri][sw[2]][i] * a[2][1];
86+ }/*for (i = 0; i < 3; ++i)*/
87+ }/*if (query == 0)*/
88+ nnl0 += 1;
89+ }/*else if (matp[ib][itri][sw[1]] < 0.0 && 0.0 <= matp[ib][itri][sw[2]])*/
90+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
91+ if (query == 1) ntri_th[ib][ithread] = nnl0;
92+ }/*for (ib = 0; ib < nb; ib++)*/
93+ }/* End of parallel region */
94+ /*
95+ Allocation of nodeline
96+ */
97+ if (query == 1) {
98+ /*
99+ Sum node-lines in all threads
100+ */
101+ for (ib = 0; ib < nb; ib++) {
102+ for (ithread = 1; ithread < nthreads; ithread++) {
103+ ntri_th[ib][ithread] += ntri_th[ib][ithread - 1];
104+ }
105+ nnl[ib] = ntri_th[ib][nthreads - 1];
106+ for (ithread = nthreads - 1; ithread > 0; ithread--) {
107+ ntri_th[ib][ithread] = ntri_th[ib][ithread - 1];
108+ }
109+ ntri_th[ib][0] = 0;
110+ }/*for (ib = 0; ib < nb; ib++)*/
111+ /*
112+ Allocation of nodeline
113+ */
114+ kvnl = new GLfloat***[nb];
115+ kvnl_rot = new GLfloat*[nb];
116+ for (ib = 0; ib < nb; ++ib) {
117+ kvnl[ib] = new GLfloat**[nnl[ib]];
118+ kvnl_rot[ib] = new GLfloat[6*nnl[ib]];
119+ for (itri = 0; itri < nnl[ib]; ++itri) {
120+ kvnl[ib][itri] = new GLfloat * [2];
121+ for (i = 0; i < 2; ++i) {
122+ kvnl[ib][itri][i] = new GLfloat[3];
123+ }/*for (j = 0; j < 2; ++j)*/
124+ }/*for (i = 0; i < nnl[ib]; ++i)*/
125+ }/*for (ib = 0; ib < nb; ++ib)*/
126+ }/*if (query == 1)*/
127+}/*void calc_nodeline()*/
--- /dev/null
+++ b/Android1.NativeActivity/calc_nodeline.hpp
@@ -0,0 +1,25 @@
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+
25+void calc_nodeline();
--- /dev/null
+++ b/Android1.NativeActivity/draw.cpp
@@ -0,0 +1,577 @@
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 Functions for actual displaying
26+*/
27+#include <GLES/gl.h>
28+#include <math.h>
29+#include <stdio.h>
30+#include <stdlib.h>
31+#include "variable.hpp"
32+/**
33+ @brief Draw Fermi surfaces
34+
35+ Modify : ::nmlp_rot, ::kvp_rot, ::kvnl_rot
36+
37+ First, rotate ::nmlp and ::kvp with OpenMP then draw them.
38+ Also draw nodeline in the same way.
39+*/
40+static void draw_fermi() {
41+ int ib;
42+ /*
43+ First, rotate k-vector and normal vector
44+ */
45+#pragma omp parallel default(none) \
46+ shared(nb,draw_band,ntri,rot,nmlp,arw,nmlp_rot,kvp,kvp_rot,arw_rot,trans,side) \
47+ private(ib)
48+ {
49+ int i, j, l, itri;
50+
51+ for (ib = 0; ib < nb; ib++) {
52+ if (draw_band[ib] == 1) {
53+#pragma omp for nowait
54+ for (itri = 0; itri < ntri[ib]; ++itri) {
55+ for (i = 0; i < 3; ++i) {
56+ for (j = 0; j < 3; ++j) {
57+ kvp_rot[ib][j + 3 * i + 9 * itri]
58+ = rot[j][0] * kvp[ib][itri][i][0]
59+ + rot[j][1] * kvp[ib][itri][i][1]
60+ + rot[j][2] * kvp[ib][itri][i][2]
61+ + trans[j];
62+ nmlp_rot[ib][j + 3 * i + 9 * itri]
63+ = rot[j][0] * nmlp[ib][itri][i][0]
64+ + rot[j][1] * nmlp[ib][itri][i][1]
65+ + rot[j][2] * nmlp[ib][itri][i][2];
66+ nmlp_rot[ib][j + 3 * i + 9 * itri] *= side;
67+ for (l = 0; l < 2; ++l) {
68+ arw_rot[ib][j + 3*l + 6 * i + 18 * itri]
69+ = rot[j][0] * arw[ib][itri][i][l][0]
70+ + rot[j][1] * arw[ib][itri][i][l][1]
71+ + rot[j][2] * arw[ib][itri][i][l][2];
72+ }
73+ }
74+ }/*for (i = 0; i < 3; ++i)*/
75+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
76+ }/*if (draw_band[ib] == 1)*/
77+ }/*for (ib = 0; ib < nb; ib++)*/
78+ }/*End of parallel region*/
79+ /*
80+ Second, draw each triangle
81+ */
82+ glEnableClientState(GL_NORMAL_ARRAY);
83+ glEnableClientState(GL_COLOR_ARRAY);
84+ for (ib = 0; ib < nb; ib++) {
85+ if (draw_band[ib] == 1) {
86+ glVertexPointer(3, GL_FLOAT, 0, kvp_rot[ib]);
87+ glNormalPointer(GL_FLOAT, 0, nmlp_rot[ib]);
88+ glColorPointer(4, GL_FLOAT, 0, clr[ib]);
89+ glDrawArrays(GL_TRIANGLES, 0, ntri[ib] * 3);
90+ }/*if (draw_band[ib] == 1)*/
91+ }/*for (ib = 0; ib < nb; ib++)*/
92+ glDisableClientState(GL_NORMAL_ARRAY);
93+ glDisableClientState(GL_COLOR_ARRAY);
94+ /*
95+ Arrow for 3D value
96+ */
97+ if (blackback == 1) glColor4f(white[0], white[1], white[2], white[3]);
98+ else glColor4f(black[0], black[1], black[2], black[3]);
99+ glNormal3f(0.0f, 0.0f, 1.0f);
100+ if (color_scale == 3) {
101+ for (ib = 0; ib < nb; ib++) {
102+ if (draw_band[ib] == 1) {
103+ glVertexPointer(3, GL_FLOAT, 0, arw_rot[ib]);
104+ glDrawArrays(GL_LINES, 0, ntri[ib] * 3 * 2);
105+ }
106+ }
107+ }
108+ /*
109+ Nodeline
110+ */
111+ if (nodeline == 1) {
112+ /*
113+ First, rotate k-vector
114+ */
115+#pragma omp parallel default(none) \
116+ shared(nb,draw_band,nnl,rot,trans,kvnl,kvnl_rot) \
117+ private(ib)
118+ {
119+ int i, j, itri;
120+
121+ for (ib = 0; ib < nb; ib++) {
122+ /**/
123+ if (draw_band[ib] == 1) {
124+#pragma omp for nowait
125+ for (itri = 0; itri < nnl[ib]; ++itri) {
126+ for (i = 0; i < 2; ++i) {
127+ /**/
128+ for (j = 0; j < 3; ++j)
129+ kvnl_rot[ib][j+3*i+6*itri]
130+ = rot[j][0] * kvnl[ib][itri][i][0]
131+ + rot[j][1] * kvnl[ib][itri][i][1]
132+ + rot[j][2] * kvnl[ib][itri][i][2]
133+ + trans[j];
134+ kvnl_rot[ib][2 + 3 * i + 6 * itri] += 0.001f;
135+ }/*for (i = 0; i < 2; ++i)*/
136+ }/*for (itri = 0; itri < nnl[ib]; ++itri)*/
137+ }/*if (draw_band[ib] == 1)*/
138+ }/*for (ib = 0; ib < nb; ib++)*/
139+ }/*End of parallel region*/
140+ /*
141+ Second, draw each lines
142+ */
143+ glColor4f(black[0], black[1], black[2], black[3]);
144+ glNormal3f(0.0f, 0.0f, 1.0f);
145+ for (ib = 0; ib < nb; ib++) {
146+ if (draw_band[ib] == 1) {
147+ glVertexPointer(3, GL_FLOAT, 0, kvnl_rot[ib]);
148+ glDrawArrays(GL_LINES, 0, 2 * nnl[ib]);
149+ }/*if (draw_band[ib] == 1)*/
150+ }/* for (ib = 0; ib < nb; ib++)*/
151+ }/*if (nodeline == 1)*/
152+ /*
153+ Equator
154+ */
155+ if (lequator == 1) {
156+ /*
157+ First, rotate k-vector
158+ */
159+#pragma omp parallel default(none) \
160+ shared(nb,draw_band,nequator,rot,trans,kveq,kveq_rot) \
161+ private(ib)
162+ {
163+ int i, j, itri;
164+
165+ for (ib = 0; ib < nb; ib++) {
166+ /**/
167+ if (draw_band[ib] == 1) {
168+#pragma omp for nowait
169+ for (itri = 0; itri < nequator[ib]; ++itri) {
170+ for (i = 0; i < 2; ++i) {
171+ /**/
172+ for (j = 0; j < 3; ++j)
173+ kveq_rot[ib][j + 3 * i + 6 * itri]
174+ = rot[j][0] * kveq[ib][itri][i][0]
175+ + rot[j][1] * kveq[ib][itri][i][1]
176+ + rot[j][2] * kveq[ib][itri][i][2]
177+ + trans[j];
178+ kveq_rot[ib][2 + 3 * i + 6 * itri] += 0.001f;
179+ }/*for (i = 0; i < 2; ++i)*/
180+ }/*for (itri = 0; itri < nequator[ib]; ++itri)*/
181+ }/*if (draw_band[ib] == 1)*/
182+ }/*for (ib = 0; ib < nb; ib++)*/
183+ }/*End of parallel region*/
184+ /*
185+ Second, draw each lines
186+ */
187+ glColor4f(black[0], black[1], black[2], black[3]);
188+ glNormal3f(0.0f, 0.0f, 1.0f);
189+ for (ib = 0; ib < nb; ib++) {
190+ if (draw_band[ib] == 1) {
191+ glVertexPointer(3, GL_FLOAT, 0, kveq_rot[ib]);
192+ glDrawArrays(GL_LINES, 0, 2 * nequator[ib]);
193+ }/*if (draw_band[ib] == 1)*/
194+ }/* for (ib = 0; ib < nb; ib++)*/
195+ }/*if (nodeline == 1)*/
196+}/*void draw_fermi*/
197+/**
198+ @brief Draw lines of BZ boundaries
199+*/
200+static void draw_bz_lines() {
201+ int ibzl, i, j;
202+ GLfloat bzl2[3], bvec2[3][3], linecolor[4], secvec2[3];
203+ GLfloat vertices[300];
204+ /*
205+ Line color is oposit of BG color
206+ */
207+ if (blackback == 1)
208+ for (i = 0; i<4; i++) linecolor[i] = white[i];
209+ else
210+ for (i = 0; i<4; i++) linecolor[i] = black[i];
211+ /**/
212+ glLineWidth(linewidth);
213+ /*
214+ First Brillouin zone mode
215+ */
216+ if (fbz == 1) {
217+ for (ibzl = 0; ibzl < nbzl; ++ibzl) {
218+ for (i = 0; i< 2; ++i) {
219+ for (j = 0; j < 3; ++j)
220+ bzl2[j] = rot[j][0] * bzl[ibzl][i][0]
221+ + rot[j][1] * bzl[ibzl][i][1]
222+ + rot[j][2] * bzl[ibzl][i][2]
223+ + trans[j];
224+ for (j = 0; j < 3; j++) vertices[j + 3 * i] = bzl2[j];
225+ }/*for (i = 0; i< 2; ++i)*/
226+ glColor4f(linecolor[0], linecolor[1], linecolor[2], linecolor[3]);
227+ glNormal3f(0.0f, 0.0f, 1.0f);
228+ glVertexPointer(3, GL_FLOAT, 0, vertices);
229+ glDrawArrays(GL_LINES, 0, 2);
230+ }/*for (ibzl = 0; ibzl < nbzl; ++ibzl)*/
231+ }/*if (fbz == 1)*/
232+ else {
233+ /*
234+ Premitive BZ mode
235+ */
236+ for (i = 0; i < 3; ++i) {
237+ for (j = 0; j < 3; ++j) {
238+ bvec2[i][j] = rot[j][0] * bvec[i][0]
239+ + rot[j][1] * bvec[i][1]
240+ + rot[j][2] * bvec[i][2];
241+ }/*for (j = 0; j < 3; ++j)*/
242+ }/*for (i = 0; i < 3; ++i)*/
243+ for (i = 0; i<3; ++i) vertices[i] = trans[i];
244+ for (i = 0; i<3; ++i) vertices[i+3] = trans[i] + bvec2[0][i];
245+ for (i = 0; i<3; ++i) vertices[i+3*2] = trans[i] + bvec2[0][i] + bvec2[1][i];
246+ for (i = 0; i<3; ++i) vertices[i+3*3] = trans[i] + bvec2[0][i] + bvec2[1][i] + bvec2[2][i];
247+ for (i = 0; i<3; ++i) vertices[i+3*4] = trans[i] + bvec2[1][i] + bvec2[2][i];
248+ for (i = 0; i<3; ++i) vertices[i+3*5] = trans[i] + bvec2[1][i];
249+ for (i = 0; i<3; ++i) vertices[i+3*6] = trans[i];
250+ for (i = 0; i<3; ++i) vertices[i+3*7] = trans[i] + bvec2[2][i];
251+ for (i = 0; i<3; ++i) vertices[i+3*8] = trans[i] + bvec2[0][i] + bvec2[2][i];
252+ for (i = 0; i<3; ++i) vertices[i+3*9] = trans[i] + bvec2[0][i] + bvec2[1][i] + bvec2[2][i];
253+ for (i = 0; i<3; ++i) vertices[i+3*10] = trans[i] + bvec2[0][i] + bvec2[2][i];
254+ for (i = 0; i<3; ++i) vertices[i+3*11] = trans[i] + bvec2[0][i];
255+ for (i = 0; i<3; ++i) vertices[i+3*12] = trans[i];
256+ for (i = 0; i<3; ++i) vertices[i+3*13] = trans[i] + bvec2[2][i];
257+ for (i = 0; i<3; ++i) vertices[i+3*14] = trans[i] + bvec2[1][i] + bvec2[2][i];
258+ for (i = 0; i<3; ++i) vertices[i+3*15] = trans[i] + bvec2[1][i];
259+ for (i = 0; i<3; ++i) vertices[i+3*16] = trans[i] + bvec2[0][i] + bvec2[1][i];
260+ glColor4f(linecolor[0], linecolor[1], linecolor[2], linecolor[3]);
261+ glNormal3f(0.0f, 0.0f, 1.0f);
262+ glVertexPointer(3, GL_FLOAT, 0, vertices);
263+ glDrawArrays(GL_LINE_STRIP, 0, 17);
264+ }/*if (fbz != 1)*/
265+ /*
266+ Section for the 2D Fermi line
267+ */
268+ if (lsection == 1 && fbz == 1) {
269+ for (j = 0; j < 3; ++j)
270+ secvec2[j] = rot[j][0] * secvec[0]
271+ + rot[j][1] * secvec[1]
272+ + rot[j][2] * secvec[2];
273+ for (ibzl = 0; ibzl < nbzl2d; ++ibzl) {
274+ for (j = 0; j < 3; ++j)
275+ bzl2[j] = rot[j][0] * bzl2d[ibzl][0]
276+ + rot[j][1] * bzl2d[ibzl][1]
277+ + rot[j][2] * bzl2d[ibzl][2]
278+ + trans[j];
279+ for (j = 0; j < 3; j++)vertices[j + 3 * ibzl] = bzl2[j];
280+ }/*for (ibzl = 0; ibzl < nbzl2d; ++ibzl)*/
281+ glColor4f(gray[0], gray[1], gray[2], gray[3]);
282+ glNormal3f(secvec2[0], secvec2[1], secvec2[2]);
283+ glVertexPointer(3, GL_FLOAT, 0, vertices);
284+ glDrawArrays(GL_TRIANGLE_FAN, 0, nbzl2d);
285+ }/*if (lsection == 1)*/
286+}/*draw bz_lines */
287+/**
288+ @brief Draw color-bar or colr-circle (periodic)
289+ as a color scale
290+*/
291+static void draw_colorbar()
292+{
293+ int i, j, k;
294+ GLfloat mat2, vertices[366], colors[488];
295+ /**/
296+ glEnableClientState(GL_COLOR_ARRAY);
297+ if (color_scale == 1 || color_scale == 4) {
298+ for (i = 0; i < 5; i++) {
299+ for (j = 0; j < 2; j++) {
300+ vertices[0 + j * 3 + i * 6] = -1.0f + 0.5f*(GLfloat)i;
301+ vertices[1 + j * 3 + i * 6] = -1.0f - 0.1f*(GLfloat)j;
302+ vertices[2 + j * 3 + i * 6] = 0.0f;
303+ if ( i == 0) for (k = 0; k < 4; k++) colors[k + 4 * j + 8 * i] = blue[k];
304+ else if (i == 1) for (k = 0; k < 4; k++) colors[k + 4 * j + 8 * i] = cyan[k];
305+ else if (i == 2) for (k = 0; k < 4; k++) colors[k + 4 * j + 8 * i] = green[k];
306+ else if (i == 3) for (k = 0; k < 4; k++) colors[k + 4 * j + 8 * i] = yellow[k];
307+ else if (i == 4) for (k = 0; k < 4; k++) colors[k + 4 * j + 8 * i] = red[k];
308+ }
309+ }/*for (i = 0; i < 10; i++)*/
310+ glNormal3f(0.0f, 0.0f, 1.0f);
311+ glVertexPointer(3, GL_FLOAT, 0, vertices);
312+ glColorPointer(4, GL_FLOAT, 0, colors);
313+ glDrawArrays(GL_TRIANGLE_STRIP, 0, 10);
314+ }/*if (color_scale == 1 || color_scale == 4)*/
315+ else if (color_scale == 2) {
316+ /*
317+ Periodic color scale
318+ */
319+ vertices[0] = 0.0f;
320+ vertices[1] = -1.0f;
321+ vertices[2] = 0.0f;
322+ if (blackback == 1)
323+ for (j = 0; j < 4; j++) colors[j] = wgray[j];
324+ else
325+ for (j = 0; j < 4; j++) colors[j] = bgray[j];
326+ /**/
327+ for (i = 0; i <= 60; i++) {
328+ /**/
329+ mat2 = (GLfloat)i / 60.0f * 6.0f;
330+ /**/
331+ if (mat2 <= 1.0) {
332+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = yellow[j] * mat2 + red[j] * (1.0f - mat2);
333+ }
334+ else if (mat2 <= 2.0) {
335+ mat2 = mat2 - 1.0f;
336+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = green[j] * mat2 + yellow[j] * (1.0f - mat2);
337+ }
338+ else if (mat2 <= 3.0) {
339+ mat2 = mat2 - 2.0f;
340+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = cyan[j] * mat2 + green[j] * (1.0f - mat2);
341+ }
342+ else if (mat2 <= 4.0) {
343+ mat2 = mat2 - 3.0f;
344+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = blue[j] * mat2 + cyan[j] * (1.0f - mat2);
345+ }
346+ else if (mat2 <= 5.0) {
347+ mat2 = mat2 - 4.0f;
348+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = magenta[j] * mat2 + blue[j] * (1.0f - mat2);
349+ }
350+ else {
351+ mat2 = mat2 - 5.0f;
352+ for (j = 0; j<4; ++j) colors[j + 4 * (i + 1)] = red[j] * mat2 + magenta[j] * (1.0f - mat2);
353+ }
354+ /**/
355+ vertices[0 + 3 * (i + 1)] = 0.2f * cosf((GLfloat)i / 60.0f * 6.283185307f);
356+ vertices[1 + 3 * (i + 1)] = 0.2f * sinf((GLfloat)i / 60.0f * 6.283185307f) - 1.0f;
357+ vertices[2 + 3 * (i + 1)] = 0.0f;
358+ }/*for (i = 0; i <= 60; i++)*/
359+ glNormal3f(0.0f, 0.0f, 1.0f);
360+ glVertexPointer(3, GL_FLOAT, 0, vertices);
361+ glColorPointer(4, GL_FLOAT, 0, colors);
362+ glDrawArrays(GL_TRIANGLE_FAN, 0, 62);
363+ }/*else if (color_scale == 2)*/
364+ else if (color_scale == 6 || color_scale == 7) {
365+ for (i = 0; i < 2; i++) {
366+ for (j = 0; j < 2; j++) {
367+ vertices[0 + j * 3 + i * 6] = -1.0f + 2.0f*(GLfloat)i;
368+ vertices[1 + j * 3 + i * 6] = -1.0f - 0.1f*(GLfloat)j;
369+ vertices[2 + j * 3 + i * 6] = 0.0f;
370+ if (i == 0) for (k = 0; k < 4; k++) colors[k + 4 * j + 8 * i] = bgray[k];
371+ else if (i == 1) for (k = 0; k < 4; k++) colors[k + 4 * j + 8 * i] = wgray[k];
372+ }
373+ }/*for (i = 0; i < 10; i++)*/
374+ glNormal3f(0.0f, 0.0f, 1.0f);
375+ glVertexPointer(3, GL_FLOAT, 0, vertices);
376+ glColorPointer(4, GL_FLOAT, 0, colors);
377+ glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
378+ }/*if (color_scale == 6 || color_scale == 7)*/
379+ glDisableClientState(GL_COLOR_ARRAY);
380+}/*void draw_colorbar*/
381+/**
382+ @brief Draw eye-points for the stereogram
383+*/
384+static void draw_circles(
385+ GLfloat dx2d //!< [in] Translation used for the section-mode
386+) {
387+ int i;
388+ GLfloat r, vertices[66];
389+ /**/
390+ r = 0.05f;
391+ /**/
392+ vertices[0] = 0.7f - dx2d;
393+ vertices[1] = scl;
394+ vertices[2] = 0.0f;
395+ for (i = 0; i <= 20; i++) {
396+ vertices[0 + (i + 1) * 3] = r * cosf((GLfloat)i / 20.0f * 6.283185307f) + 0.7f - dx2d;
397+ vertices[1 + (i + 1) * 3] = r * sinf((GLfloat)i / 20.0f * 6.283185307f) + scl;
398+ vertices[2 + (i + 1) * 3] = 0.0f;
399+ }/*for (i = 0; i <= 20; i++)*/
400+ glNormal3f(0.0f, 0.0f, 1.0f);
401+ if (blackback == 1) glColor4f(white[0], white[1], white[2], white[3]);
402+ else glColor4f(black[0], black[1], black[2], black[3]);
403+ glVertexPointer(3, GL_FLOAT, 0, vertices);
404+ glDrawArrays(GL_TRIANGLE_FAN, 0, 22);
405+ /**/
406+ for (i = 0; i < 22; i++) vertices[3 * i] += -1.4f;
407+ glVertexPointer(3, GL_FLOAT, 0, vertices);
408+ glDrawArrays(GL_TRIANGLE_FAN, 0, 22);
409+}/*void draw_circles*/
410+/**
411+ @brief Draw 2D Fermi lines
412+*/
413+static void draw_fermi_line() {
414+ int i, ib, ibzl;
415+ GLfloat vertices[60];
416+ /*
417+ Draw 2D BZ lines
418+ */
419+ glLineWidth(linewidth);
420+ for (ibzl = 0; ibzl < nbzl2d; ++ibzl) {
421+ for (i = 0; i < 3; i++) vertices[i + 3 * ibzl] = bzl2d_proj[ibzl][i];
422+ }/*for (ibzl = 0; ibzl < nbzl2d; ++ibzl)*/
423+ glNormal3f(0.0f, 0.0f, 1.0f);
424+ if (blackback == 1) glColor4f(white[0], white[1], white[2], white[3]);
425+ else glColor4f(black[0], black[1], black[2], black[3]);
426+ glVertexPointer(3, GL_FLOAT, 0, vertices);
427+ glDrawArrays(GL_LINE_LOOP, 0, nbzl2d);
428+ /*
429+ Draw Fermi lines
430+ */
431+ glLineWidth(linewidth);
432+ glEnableClientState(GL_COLOR_ARRAY);
433+ glNormal3f(0.0f, 0.0f, 1.0f);
434+ for (ib = 0; ib < nb; ib++) {
435+ if (draw_band[ib] == 1) {
436+ glVertexPointer(3, GL_FLOAT, 0, kv2d[ib]);
437+ glColorPointer(4, GL_FLOAT, 0, clr2d[ib]);
438+ glDrawArrays(GL_LINES, 0, 2 * n2d[ib]);
439+ }/*if (draw_band[ib] == 1)*/
440+ }/* for (ib = 0; ib < nb; ib++)*/
441+ glDisableClientState(GL_COLOR_ARRAY);
442+}/*void draw_fermi_line*/
443+/**
444+ @brief Glut Display function
445+ called by glutDisplayFunc
446+*/
447+void draw()
448+{
449+ GLfloat pos[] = { 1.0f, 1.0f, 1.0f, 0.0f };
450+ GLfloat amb[] = { 0.2f, 0.2f, 0.2f, 0.0f };
451+ GLfloat dx, dx2d, theta, posz, phi;
452+ GLfloat pos1[4], pos2[4];
453+
454+ if (draw_band == NULL) return;
455+
456+ if (lstereo == 2) {
457+ /*
458+ Parallel eyes
459+ */
460+ theta = 3.1416f / 180.0f * 2.5f;
461+ posz = 5.0f;
462+ dx = 0.7f;
463+ phi = atanf(posz / dx) - theta;
464+ theta = (3.1416f * 0.5f - phi) / 3.1416f * 180.0f;
465+ /**/
466+ pos1[0] = posz * cosf(phi) - dx;
467+ pos1[1] = 0.0;
468+ pos1[2] = posz * sinf(phi);
469+ pos1[3] = 1.0;
470+ /**/
471+ pos2[0] = -posz * cosf(phi) + dx;
472+ pos2[1] = 0.0;
473+ pos2[2] = posz * sinf(phi);
474+ pos2[3] = 1.0;
475+ }/*if (lstereo == 2)*/
476+ else if (lstereo == 3) {
477+ /*
478+ Cross eyes
479+ */
480+ theta = -3.1416f / 180.0f * 2.0f;
481+ posz = 5.0;
482+ dx = -0.7f;
483+ /**/
484+ pos1[0] = posz * sinf(theta) - dx;
485+ pos1[1] = 0.0;
486+ pos1[2] = posz * cosf(theta);
487+ pos1[3] = 1.0;
488+ /**/
489+ pos2[0] = -posz * sinf(theta) + dx;
490+ pos2[1] = 0.0;
491+ pos2[2] = posz * cosf(theta);
492+ pos2[3] = 1.0;
493+ }/*if (lstereo == 3)*/
494+ else {
495+ theta = 0.0;
496+ dx = 0.0;
497+ }/*lstero == 1*/
498+ if (lsection == 1 && fbz == 1) dx2d = 0.7f;
499+ else dx2d = 0.0;
500+ /*
501+ Initialize
502+ */
503+ glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
504+ glLoadIdentity();
505+ /*
506+ Set view point & view line
507+ */
508+ glTranslatef(0.0, 0.0, -5.0);
509+ //gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
510+ /*
511+ Set position of light
512+ */
513+ if (lstereo == 1) {
514+ glLightfv(GL_LIGHT0, GL_POSITION, pos);
515+ glTranslatef(-dx2d, 0.0, 0.0);
516+ /*
517+ Draw color scale
518+ */
519+ if (lcolorbar == 1) draw_colorbar();
520+ }
521+ else {
522+ glLightfv(GL_LIGHT0, GL_POSITION, pos1);
523+ draw_circles(dx2d);
524+ glTranslatef(-dx-dx2d, 0.0, 0.0);
525+ glRotatef(theta, 0.0, 1.0, 0.0);
526+ }
527+ glLightfv(GL_LIGHT1, GL_AMBIENT, amb);
528+ /*
529+ Rotation & Zoom
530+ */
531+ glScalef(scl, scl, scl);
532+ /*
533+ Draw Brillouin zone boundaries
534+ */
535+ draw_bz_lines();
536+ /*
537+ Draw Fermi surfaces
538+ */
539+ draw_fermi();
540+ /*
541+ Draw the second object for stereogram
542+ */
543+ if (lstereo != 1) {
544+ glPushMatrix();
545+ glLoadIdentity();
546+ glTranslatef(0.0, 0.0, -5.0);
547+ //gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
548+ glLightfv(GL_LIGHT0, GL_POSITION, pos2);
549+ /**/
550+ glTranslatef(dx-dx2d, 0.0, 0.0);
551+ glRotatef(-theta, 0.0, 1.0, 0.0);
552+ /**/
553+ glScalef(scl, scl, scl);
554+ draw_bz_lines();
555+ draw_fermi();
556+ /**/
557+ glPopMatrix();
558+ }
559+ /*
560+ Draw 2D Fermi line
561+ */
562+ if (lsection == 1 && fbz == 1) {
563+ glPushMatrix();
564+ glLoadIdentity();
565+ glTranslatef(0.0, 0.0, -5.0);
566+ //gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
567+ glLightfv(GL_LIGHT0, GL_POSITION, pos);
568+ /**/
569+ if (lstereo == 1) glTranslatef(dx2d, 0.0, 0.0);
570+ else glTranslatef(2.0f * dx2d, 0.0, 0.0);
571+ /**/
572+ glScalef(scl, scl, scl);
573+ draw_fermi_line();
574+ /**/
575+ glPopMatrix();
576+ }/*if (lsection == 1)*/
577+}/*void display*/
--- /dev/null
+++ b/Android1.NativeActivity/draw.hpp
@@ -0,0 +1,25 @@
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+
25+void draw();
--- /dev/null
+++ b/Android1.NativeActivity/equator.cpp
@@ -0,0 +1,127 @@
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 Compute nodal lines
26+*/
27+#include <GLES/gl.h>
28+#include <stdio.h>
29+#include <stdlib.h>
30+#include <math.h>
31+#include "basic_math.hpp"
32+#include "variable.hpp"
33+/**
34+@brief Compute equator \f$\{\bf v}_{n k} \cdot {\bf k} = 0\f$
35+
36+ Modify : ::ntri_th, ::nequator, ::kveq, ::kveq_rot
37+
38+ If ::query = 1, this routine only compute the number of
39+ line segmants and malloc variables.
40+*/
41+void equator() {
42+ int ib, itri, i, j, ithread;
43+
44+#pragma omp parallel default(none) \
45+ shared(nb,nequator,nmlp,kveq,kvp,ntri,ntri_th,query,eqvec) \
46+ private(ib,itri,i,j,ithread)
47+ {
48+ int sw[3], nequator0;
49+ GLfloat a[3][3], prod[3];
50+
51+ ithread = get_thread();
52+ for (ib = 0; ib < nb; ib++) {
53+ if(query == 1) nequator0 = 0;
54+ else nequator0 = ntri_th[ib][ithread];
55+#pragma omp for
56+ for (itri = 0; itri < ntri[ib]; ++itri) {
57+
58+ for (i = 0; i < 3; ++i) {
59+ prod[i] = 0.0f;
60+ for (j = 0; j < 3; ++j) prod[i] += eqvec[j] * nmlp[ib][itri][i][j];
61+ }/*for (i = 0; i < 3; ++i)*/
62+
63+ eigsort(3, prod, sw);
64+ for (i = 0; i < 3; ++i) {
65+ for (j = 0; j < 3; ++j) {
66+ a[i][j] = (0.f - prod[sw[j]]) / (prod[sw[i]] - prod[sw[j]]);
67+ }/*for (j = 0; j < 3; ++j)*/
68+ }/*for (i = 0; i < 3; ++i)*/
69+
70+ if ((prod[sw[0]] < 0.0 && 0.0 <= prod[sw[1]])
71+ || (prod[sw[0]] <= 0.0 && 0.0 < prod[sw[1]])) {
72+ if (query == 0) {
73+ for (i = 0; i < 3; ++i) {
74+ kveq[ib][nequator0][0][i] = kvp[ib][itri][sw[0]][i] * a[0][1] + kvp[ib][itri][sw[1]][i] * a[1][0];
75+ kveq[ib][nequator0][1][i] = kvp[ib][itri][sw[0]][i] * a[0][2] + kvp[ib][itri][sw[2]][i] * a[2][0];
76+ }/*for (i = 0; i < 3; ++i)*/
77+ }/*if (query == 0)*/
78+ nequator0 += 1;
79+ }/*else if (prod[sw[0]] < 0.0 && 0.0 <= prod[sw[1]])*/
80+ else if ((prod[sw[1]] < 0.0 && 0.0 <= prod[sw[2]])
81+ || (prod[sw[1]] <= 0.0 && 0.0 < prod[sw[2]])) {
82+ if (query == 0) {
83+ for (i = 0; i < 3; ++i) {
84+ kveq[ib][nequator0][0][i] = kvp[ib][itri][sw[0]][i] * a[0][2] + kvp[ib][itri][sw[2]][i] * a[2][0];
85+ kveq[ib][nequator0][1][i] = kvp[ib][itri][sw[1]][i] * a[1][2] + kvp[ib][itri][sw[2]][i] * a[2][1];
86+ }/*for (i = 0; i < 3; ++i)*/
87+ }/*if (query == 0)*/
88+ nequator0 += 1;
89+ }/*else if (prod[sw[1]] < 0.0 && 0.0 <= prod[sw[2]])*/
90+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
91+ if (query == 1) ntri_th[ib][ithread] = nequator0;
92+ }/*for (ib = 0; ib < nb; ib++)*/
93+ }/* End of parallel region */
94+ /*
95+ Allocation of nodeline
96+ */
97+ if (query == 1) {
98+ /*
99+ Sum node-lines in all threads
100+ */
101+ for (ib = 0; ib < nb; ib++) {
102+ for (ithread = 1; ithread < nthreads; ithread++) {
103+ ntri_th[ib][ithread] += ntri_th[ib][ithread - 1];
104+ }
105+ nequator[ib] = ntri_th[ib][nthreads - 1];
106+ for (ithread = nthreads - 1; ithread > 0; ithread--) {
107+ ntri_th[ib][ithread] = ntri_th[ib][ithread - 1];
108+ }
109+ ntri_th[ib][0] = 0;
110+ }/*for (ib = 0; ib < nb; ib++)*/
111+ /*
112+ Allocation of nodeline
113+ */
114+ kveq = new GLfloat***[nb];
115+ kveq_rot = new GLfloat*[nb];
116+ for (ib = 0; ib < nb; ++ib) {
117+ kveq[ib] = new GLfloat**[nequator[ib]];
118+ kveq_rot[ib] = new GLfloat[6 * nequator[ib]];
119+ for (itri = 0; itri < nequator[ib]; ++itri) {
120+ kveq[ib][itri] = new GLfloat*[2];
121+ for (i = 0; i < 2; ++i) {
122+ kveq[ib][itri][i] = new GLfloat[3];
123+ }/*for (j = 0; j < 2; ++j)*/
124+ }/*for (i = 0; i < nequator[ib]; ++i)*/
125+ }/*for (ib = 0; ib < nb; ++ib)*/
126+ }/*if (query == 1)*/
127+}/*void equator()*/
--- /dev/null
+++ b/Android1.NativeActivity/equator.hpp
@@ -0,0 +1,25 @@
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+
25+void equator();
--- /dev/null
+++ b/Android1.NativeActivity/fermi_patch.cpp
@@ -0,0 +1,594 @@
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 Functions for computing patch of Fermi surface
26+*/
27+#include <GLES/gl.h>
28+#include <stdlib.h>
29+#include <stdio.h>
30+#include <math.h>
31+#include "basic_math.hpp"
32+#include "variable.hpp"
33+/**
34+ @brief Store triangle patch
35+
36+ Modify : ::matp, ::kvp
37+
38+ For the 1st BZ mode, this routine cuts triangle recursivly at the
39+ BZ boundary (Bragg plane).
40+ If ::query == 1, this routine only increment the number of patch.
41+ If ::query == 0, store acutually the corner.
42+
43+ - DO @f${\bf l}@f$ in Bragg vector
44+ - @f${p_i = {\bf l}\cdot{\bf k}}@f$
45+ - Sort : @f$p_0<p_1<p_2@f$
46+ - @f[
47+ a_{i j} \equiv \frac{-p_j}{p_i - p_j}
48+ @f]
49+ - if (@f$|{\bf l}| < p_0@f$)
50+ - This patch is not in the 1st BZ
51+ - if (@f$p_0 < |{\bf l}| < p_1@f$)
52+ - @f${\bf k}'_0 = {\bf k}_0@f$
53+ - @f${\bf k}'_1 = {\bf k}_0 a_{0 1} + {\bf k}_1 a_{1 0}@f$
54+ - @f${\bf k}'_2 = {\bf k}_0 a_{0 2} + {\bf k}_2 a_{2 0}@f$
55+ - if (@f$p_1 < |{\bf l}| < p_2@f$)
56+ - @f${\bf k}'_0 = {\bf k}_0@f$
57+ - @f${\bf k}'_1 = {\bf k}_1@f$
58+ - @f${\bf k}'_2 = {\bf k}_0 a_{0 2} + {\bf k}_2 a_{2 0}@f$
59+ - and
60+ - @f${\bf k}'_0 = {\bf k}_1 a_{1 2} + {\bf k}_2 a_{2 1}@f$
61+ - @f${\bf k}'_1 = {\bf k}_1@f$
62+ - @f${\bf k}'_2 = {\bf k}_0 a_{0 2} + {\bf k}_2 a_{2 0}@f$
63+ - if (@f$p_2 < |{\bf l}| < p_3@f$)
64+ - @f${\bf k}'_0 = {\bf k}_0@f$
65+ - @f${\bf k}'_1 = {\bf k}_1@f$
66+ - @f${\bf k}'_2 = {\bf k}_2@f$
67+ - END DO
68+*/
69+static void triangle(
70+ int ib, //!<[in] The band index
71+ int *ntri0, //!<[inout] Index of triangle patch
72+ int nbr, //!<[in] Bragg plane
73+ GLfloat mat1[3][3], //!<[in] The matrix element
74+ GLfloat kvec1[3][3], //!<[in] @f$k@f$-vector of corners
75+ GLfloat vf1[3][3] //!<[in] @f$v_f@f$-vector of corners
76+)
77+{
78+ int ibr, i, j, sw[3];
79+ GLfloat prod[3], thr, mat2[3][3], kvec2[3][3],
80+ vf2[3][3], a[3][3], bshift, vfave[3], norm[3];
81+ /*
82+ If the area is nearly 0, it is ignored.
83+ */
84+ for (i = 0; i < 3; i++)norm[i] = 0.0f;
85+ for (i = 0; i < 3; i++) {
86+ norm[0] += (kvec1[1][i] - kvec1[2][i])*(kvec1[1][i] - kvec1[2][i]);
87+ norm[1] += (kvec1[2][i] - kvec1[0][i])*(kvec1[2][i] - kvec1[0][i]);
88+ norm[2] += (kvec1[0][i] - kvec1[1][i])*(kvec1[0][i] - kvec1[1][i]);
89+ }
90+ for (i = 0; i < 3; i++) {
91+ if (norm[i] < 1.0e-10f*brnrm_min) return;
92+ }
93+ /*
94+ For 1st BZ, it is cut at the BZ boundary.
95+ */
96+ if (fbz == 1) {
97+ /**/
98+ for (ibr = 0; ibr < nbragg; ++ibr) {
99+
100+ thr = brnrm[ibr] * 0.001f;
101+ /**/
102+ for (i = 0; i < 3; ++i)
103+ prod[i] = bragg[ibr][0] * kvec1[i][0]
104+ + bragg[ibr][1] * kvec1[i][1]
105+ + bragg[ibr][2] * kvec1[i][2];
106+ eigsort(3, prod, sw);
107+ for (i = 0; i < 3; ++i) {
108+ for (j = 0; j < 3; ++j) {
109+ a[i][j] = (brnrm[ibr] - prod[sw[j]]) / (prod[sw[i]] - prod[sw[j]]);
110+ }/*for (j = 0; j < 3; ++j)*/
111+ }/*for (i = 0; i < 3; ++i)*/
112+ i = (int)(0.5f * ((prod[sw[2]] / brnrm[ibr]) + 1.0f));
113+ bshift = -2.0f *(GLfloat)i;
114+
115+ if (brnrm[ibr] + thr > prod[sw[2]]) continue;
116+
117+ if (brnrm[ibr] < prod[sw[0]]) {
118+ /*
119+ All corners are outside of the Bragg plane
120+ */
121+ for (i = 0; i < 3; ++i) {
122+ for (j = 0; j < 3; ++j) {
123+ kvec2[i][j] = kvec1[sw[i]][j] + bshift * bragg[ibr][j];
124+ mat2[i][j] = mat1[sw[i]][j];
125+ vf2[i][j] = vf1[sw[i]][j];
126+ }
127+ }
128+ triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2);
129+ return;
130+ }
131+ else if (brnrm[ibr] < prod[sw[1]]) {
132+ /*
133+ Single corner (#0) is inside of the Bragg plane
134+ */
135+ for (i = 0; i < 3; ++i) {
136+ kvec2[0][i] = kvec1[sw[0]][i];
137+ kvec2[1][i] = kvec1[sw[0]][i] * a[0][1] + kvec1[sw[1]][i] * a[1][0];
138+ kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
139+
140+ mat2[0][i] = mat1[sw[0]][i];
141+ mat2[1][i] = mat1[sw[0]][i] * a[0][1] + mat1[sw[1]][i] * a[1][0];
142+ mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
143+
144+ vf2[0][i] = vf1[sw[0]][i];
145+ vf2[1][i] = vf1[sw[0]][i] * a[0][1] + vf1[sw[1]][i] * a[1][0];
146+ vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
147+ }/*for (i = 0; i < 3; ++i)*/
148+ triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2);
149+
150+ for (i = 0; i < 3; ++i) {
151+ kvec2[0][i] = kvec1[sw[0]][i] * a[0][1] + kvec1[sw[1]][i] * a[1][0];
152+ kvec2[1][i] = kvec1[sw[1]][i];
153+ kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
154+
155+ mat2[0][i] = mat1[sw[0]][i] * a[0][1] + mat1[sw[1]][i] * a[1][0];
156+ mat2[1][i] = mat1[sw[1]][i];
157+ mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
158+
159+ vf2[0][i] = vf1[sw[0]][i] * a[0][1] + vf1[sw[1]][i] * a[1][0];
160+ vf2[1][i] = vf1[sw[1]][i];
161+ vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
162+ }/*for (i = 0; i < 3; ++i)*/
163+ for (i = 0; i < 3; ++i) for (j = 0; j < 3; ++j)
164+ kvec2[i][j] += bshift * bragg[ibr][j];
165+ triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2);
166+
167+ for (i = 0; i < 3; ++i) {
168+ kvec2[0][i] = kvec1[sw[2]][i];
169+ kvec2[1][i] = kvec1[sw[1]][i];
170+ kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
171+
172+ mat2[0][i] = mat1[sw[2]][i];
173+ mat2[1][i] = mat1[sw[1]][i];
174+ mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
175+
176+ vf2[0][i] = vf1[sw[2]][i];
177+ vf2[1][i] = vf1[sw[1]][i];
178+ vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
179+ }/*for (i = 0; i < 3; ++i)*/
180+ for (i = 0; i < 3; ++i) for (j = 0; j < 3; ++j)
181+ kvec2[i][j] += bshift * bragg[ibr][j];
182+ triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2);
183+ return;
184+ }
185+ else if (brnrm[ibr] < prod[sw[2]]) {
186+ /*
187+ Two corners (#0, #1) are inside of the Bragg plane
188+ */
189+ for (i = 0; i < 3; ++i) {
190+ kvec2[0][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
191+ kvec2[1][i] = kvec1[sw[1]][i] * a[1][2] + kvec1[sw[2]][i] * a[2][1];
192+ kvec2[2][i] = kvec1[sw[2]][i];
193+
194+ mat2[0][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
195+ mat2[1][i] = mat1[sw[1]][i] * a[1][2] + mat1[sw[2]][i] * a[2][1];
196+ mat2[2][i] = mat1[sw[2]][i];
197+
198+ vf2[0][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
199+ vf2[1][i] = vf1[sw[1]][i] * a[1][2] + vf1[sw[2]][i] * a[2][1];
200+ vf2[2][i] = vf1[sw[2]][i];
201+ }/*for (i = 0; i < 3; ++i)*/
202+ for (i = 0; i < 3; ++i) for (j = 0; j < 3; ++j)
203+ kvec2[i][j] += bshift * bragg[ibr][j];
204+ triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2);
205+
206+ for (i = 0; i < 3; ++i) {
207+ kvec2[0][i] = kvec1[sw[0]][i];
208+ kvec2[1][i] = kvec1[sw[1]][i];
209+ kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
210+
211+ mat2[0][i] = mat1[sw[0]][i];
212+ mat2[1][i] = mat1[sw[1]][i];
213+ mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
214+
215+ vf2[0][i] = vf1[sw[0]][i];
216+ vf2[1][i] = vf1[sw[1]][i];
217+ vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
218+ }/*for (i = 0; i < 3; ++i)*/
219+ triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2);
220+
221+ for (i = 0; i < 3; ++i) {
222+ kvec2[0][i] = kvec1[sw[1]][i] * a[1][2] + kvec1[sw[2]][i] * a[2][1];
223+ kvec2[1][i] = kvec1[sw[1]][i];
224+ kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0];
225+
226+ mat2[0][i] = mat1[sw[1]][i] * a[1][2] + mat1[sw[2]][i] * a[2][1];
227+ mat2[1][i] = mat1[sw[1]][i];
228+ mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0];
229+
230+ vf2[0][i] = vf1[sw[1]][i] * a[1][2] + vf1[sw[2]][i] * a[2][1];
231+ vf2[1][i] = vf1[sw[1]][i];
232+ vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0];
233+ }/*for (i = 0; i < 3; ++i)*/
234+ triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2);
235+ return;
236+ }
237+ else {
238+ /*
239+ All corners are inside of the Bragg plane
240+ */
241+ } /* brnrm[ibr] + thr < prod */
242+ } /* for ibr = 1; ibr < nbragg*/
243+ } /* if fbz == 1 */
244+ /*
245+ If it is not query, store.
246+ If query, only count the number of patches.
247+ */
248+ if (query != 1){
249+ normal_vec(kvec1[0], kvec1[1], kvec1[2], norm);
250+ for (i = 0; i < 3; ++i) {
251+ vfave[i] = 0.0f;
252+ for (j = 0; j < 3; ++j) vfave[i] += vf1[j][i];
253+ }
254+ prod[0] = 0.0f;
255+ for (i = 0; i < 3; ++i) prod[0] += vfave[i] * norm[i];
256+
257+ if (prod[0] < 0.0f) {
258+ for (i = 0; i < 3; ++i) {
259+ for (j = 0; j < 3; ++j) {
260+ kvp[ib][*ntri0][i][j] = kvec1[2-i][j];
261+ matp[ib][*ntri0][i][j] = mat1[2 - i][j];
262+ nmlp[ib][*ntri0][i][j] = vf1[2-i][j];
263+ }
264+ }/*for (i = 0; i < 3; ++i)*/
265+ }
266+ else {
267+ for (i = 0; i < 3; ++i) {
268+ for (j = 0; j < 3; ++j) {
269+ kvp[ib][*ntri0][i][j] = kvec1[i][j];
270+ matp[ib][*ntri0][i][j] = mat1[i][j];
271+ nmlp[ib][*ntri0][i][j] = vf1[i][j];
272+ }
273+ }/*for (i = 0; i < 3; ++i)*/
274+ }
275+ }/*if (query != 1)*/
276+ *ntri0 += 1;
277+}/* triangle */
278+/**
279+@brief Cut triangle patch with the tetrahedron method.
280+
281+ - Sort : @f$\varepsilon_0<\varepsilon_1<\varepsilon_2<\varepsilon_3@f$
282+ - @f[
283+ a_{i j} \equiv \frac{-\varepsilon_j}{\varepsilon_i - \varepsilon_j}
284+ @f]
285+ - if (@f$\varepsilon_0 < 0 < \varepsilon_1@f$)
286+ - @f${\bf k}'_0 = {\bf k}_0 a_{0 1} + {\bf k}_1 a_{1 0}@f$
287+ - @f${\bf k}'_1 = {\bf k}_0 a_{0 2} + {\bf k}_2 a_{2 0}@f$
288+ - @f${\bf k}'_2 = {\bf k}_0 a_{0 3} + {\bf k}_3 a_{3 0}@f$
289+ - if (@f$\varepsilon_1 < 0 < \varepsilon_2@f$)
290+ - @f${\bf k}'_0 = {\bf k}_0 a_{0 2} + {\bf k}_2 a_{2 0}@f$
291+ - @f${\bf k}'_1 = {\bf k}_0 a_{0 3} + {\bf k}_3 a_{3 0}@f$
292+ - @f${\bf k}'_2 = {\bf k}_1 a_{1 2} + {\bf k}_2 a_{2 1}@f$
293+ - and
294+ - @f${\bf k}'_0 = {\bf k}_1 a_{1 3} + {\bf k}_3 a_{3 1}@f$
295+ - @f${\bf k}'_1 = {\bf k}_0 a_{0 3} + {\bf k}_3 a_{3 0}@f$
296+ - @f${\bf k}'_2 = {\bf k}_1 a_{1 2} + {\bf k}_2 a_{2 1}@f$
297+ - if (@f$\varepsilon_2 < 0 < \varepsilon_3@f$)
298+ - @f${\bf k}'_0 = {\bf k}_3 a_{3 0} + {\bf k}_0 a_{0 3}@f$
299+ - @f${\bf k}'_1 = {\bf k}_3 a_{3 1} + {\bf k}_1 a_{1 3}@f$
300+ - @f${\bf k}'_2 = {\bf k}_3 a_{3 2} + {\bf k}_2 a_{2 3}@f$
301+*/
302+static void tetrahedron(
303+ int ib, //!< [in] The band index
304+ int *ntri0, //!< [in,out] Counter of trixngle
305+ GLfloat eig1[8], //!< [in] Orbital energies @f$\varepsilon_{n k}@f$
306+ GLfloat mat1[8][3], //!< [in] Matrix elements @f$\Delta_{n k}@f$
307+ GLfloat kvec1[8][3], //!< [in] @f$k@f$-vectors
308+ GLfloat vf1[8][3] //!< [in] @f$v_f@f$-vectors
309+)
310+{
311+ int it, i, j, sw[4];
312+ GLfloat eig2[4], mat2[4][3], kvec2[4][3], vf2[4][3], a[4][4],
313+ kvec3[3][3], mat3[3][3], vf3[3][3], vol;
314+
315+ for (it = 0; it < 6; ++it) {
316+ /*
317+ Define corners of the tetrahedron
318+ */
319+ for (i = 0; i < 4; ++i) {
320+ eig2[i] = eig1[corner[it][i]];
321+ for (j = 0; j < 3; ++j) {
322+ mat2[i][j] = mat1[corner[it][i]][j];
323+ vf2[i][j] = vf1[corner[it][i]][j];
324+ }
325+ /*
326+ Fractional -> Cartecian
327+ */
328+ for (j = 0; j < 3; ++j)
329+ kvec2[i][j] = bvec[0][j] * kvec1[corner[it][i]][0]
330+ + bvec[1][j] * kvec1[corner[it][i]][1]
331+ + bvec[2][j] * kvec1[corner[it][i]][2];
332+ }/*for (i = 0; i < 4; ++i)*/
333+ eigsort(4, eig2, sw);
334+
335+ for (i = 0; i < 4; ++i) {
336+ for (j = 0; j < 4; ++j) {
337+ a[i][j] = (0.0f - eig2[sw[j]]) / (eig2[sw[i]] - eig2[sw[j]]);
338+ }/*for (j = 0; j < 4; ++j)*/
339+ }/*for (i = 0; i < 4; ++i)*/
340+ /*
341+ Draw triangle in each cases
342+ */
343+ if (eig2[sw[0]] <= 0.0 && 0.0 < eig2[sw[1]]) {
344+ for (i = 0; i < 3; ++i) {
345+ kvec3[0][i] = kvec2[sw[0]][i] * a[0][1] + kvec2[sw[1]][i] * a[1][0];
346+ kvec3[1][i] = kvec2[sw[0]][i] * a[0][2] + kvec2[sw[2]][i] * a[2][0];
347+ kvec3[2][i] = kvec2[sw[0]][i] * a[0][3] + kvec2[sw[3]][i] * a[3][0];
348+
349+ mat3[0][i] = mat2[sw[0]][i] * a[0][1] + mat2[sw[1]][i] * a[1][0];
350+ mat3[1][i] = mat2[sw[0]][i] * a[0][2] + mat2[sw[2]][i] * a[2][0];
351+ mat3[2][i] = mat2[sw[0]][i] * a[0][3] + mat2[sw[3]][i] * a[3][0];
352+
353+ vf3[0][i] = vf2[sw[0]][i] * a[0][1] + vf2[sw[1]][i] * a[1][0];
354+ vf3[1][i] = vf2[sw[0]][i] * a[0][2] + vf2[sw[2]][i] * a[2][0];
355+ vf3[2][i] = vf2[sw[0]][i] * a[0][3] + vf2[sw[3]][i] * a[3][0];
356+ }
357+
358+ vol = a[1][0] * a[2][0] * a[3][0];
359+ triangle(ib, ntri0, 0, mat3, kvec3, vf3);
360+ }
361+ else if (eig2[sw[1]] <= 0.0 && 0.0 < eig2[sw[2]]) {
362+ for (i = 0; i < 3; ++i) {
363+ kvec3[0][i] = kvec2[sw[0]][i] * a[0][2] + kvec2[sw[2]][i] * a[2][0];
364+ kvec3[1][i] = kvec2[sw[0]][i] * a[0][3] + kvec2[sw[3]][i] * a[3][0];
365+ kvec3[2][i] = kvec2[sw[1]][i] * a[1][2] + kvec2[sw[2]][i] * a[2][1];
366+
367+ mat3[0][i] = mat2[sw[0]][i] * a[0][2] + mat2[sw[2]][i] * a[2][0];
368+ mat3[1][i] = mat2[sw[0]][i] * a[0][3] + mat2[sw[3]][i] * a[3][0];
369+ mat3[2][i] = mat2[sw[1]][i] * a[1][2] + mat2[sw[2]][i] * a[2][1];
370+
371+ vf3[0][i] = vf2[sw[0]][i] * a[0][2] + vf2[sw[2]][i] * a[2][0];
372+ vf3[1][i] = vf2[sw[0]][i] * a[0][3] + vf2[sw[3]][i] * a[3][0];
373+ vf3[2][i] = vf2[sw[1]][i] * a[1][2] + vf2[sw[2]][i] * a[2][1];
374+ }
375+
376+ vol = a[1][2] * a[2][0] * a[3][0];
377+ triangle(ib, ntri0, 0, mat3, kvec3, vf3);
378+ /**/
379+ for (i = 0; i < 3; ++i) {
380+ kvec3[0][i] = kvec2[sw[1]][i] * a[1][3] + kvec2[sw[3]][i] * a[3][1];
381+ kvec3[1][i] = kvec2[sw[0]][i] * a[0][3] + kvec2[sw[3]][i] * a[3][0];
382+ kvec3[2][i] = kvec2[sw[1]][i] * a[1][2] + kvec2[sw[2]][i] * a[2][1];
383+
384+ mat3[0][i] = mat2[sw[1]][i] * a[1][3] + mat2[sw[3]][i] * a[3][1];
385+ mat3[1][i] = mat2[sw[0]][i] * a[0][3] + mat2[sw[3]][i] * a[3][0];
386+ mat3[2][i] = mat2[sw[1]][i] * a[1][2] + mat2[sw[2]][i] * a[2][1];
387+
388+ vf3[0][i] = vf2[sw[1]][i] * a[1][3] + vf2[sw[3]][i] * a[3][1];
389+ vf3[1][i] = vf2[sw[0]][i] * a[0][3] + vf2[sw[3]][i] * a[3][0];
390+ vf3[2][i] = vf2[sw[1]][i] * a[1][2] + vf2[sw[2]][i] * a[2][1];
391+ }
392+
393+ vol = a[1][3] * a[3][0] * a[2][1];
394+ triangle(ib, ntri0, 0, mat3, kvec3, vf3);
395+ }
396+ else if (eig2[sw[2]] <= 0.0 && 0.0 < eig2[sw[3]]) {
397+ for (i = 0; i < 3; ++i) {
398+ kvec3[0][i] = kvec2[sw[3]][i] * a[3][0] + kvec2[sw[0]][i] * a[0][3];
399+ kvec3[1][i] = kvec2[sw[3]][i] * a[3][1] + kvec2[sw[1]][i] * a[1][3];
400+ kvec3[2][i] = kvec2[sw[3]][i] * a[3][2] + kvec2[sw[2]][i] * a[2][3];
401+
402+ mat3[0][i] = mat2[sw[3]][i] * a[3][0] + mat2[sw[0]][i] * a[0][3];
403+ mat3[1][i] = mat2[sw[3]][i] * a[3][1] + mat2[sw[1]][i] * a[1][3];
404+ mat3[2][i] = mat2[sw[3]][i] * a[3][2] + mat2[sw[2]][i] * a[2][3];
405+
406+ vf3[0][i] = vf2[sw[3]][i] * a[3][0] + vf2[sw[0]][i] * a[0][3];
407+ vf3[1][i] = vf2[sw[3]][i] * a[3][1] + vf2[sw[1]][i] * a[1][3];
408+ vf3[2][i] = vf2[sw[3]][i] * a[3][2] + vf2[sw[2]][i] * a[2][3];
409+ }
410+
411+ vol = a[0][3] * a[1][3] * a[2][3];
412+ triangle(ib, ntri0, 0, mat3, kvec3, vf3);
413+ }
414+ else {
415+ }
416+ }/*for (it = 0; it < 6; ++it)*/
417+}/* tetrahedron */
418+/**
419+ @brief Compute patches for Fermi surfaces
420+
421+ Modify : ::ntri, ::ntri_th, nmlp, ::matp, ::kvp, ::clr, ::nmlp_rot, ::kvp_rot
422+ If ::query = 1, this routine compute the number of patches and malloc these variables.
423+*/
424+void fermi_patch()
425+{
426+ int ib, i0, i1, j0, start[3], last[3];
427+ int ithread;
428+ /**/
429+ if (fbz == 1) {
430+ for (i0 = 0; i0 < 3; ++i0) {
431+ start[i0] = ng[i0] / 2 - ng[i0];
432+ last[i0] = ng[i0] / 2;
433+ }
434+ }
435+ else {
436+ for (i0 = 0; i0 < 3; ++i0) {
437+ start[i0] = 0;
438+ last[i0] = ng[i0];
439+ }
440+ }
441+ /**/
442+#pragma omp parallel default(none) \
443+ shared(nb,ntri,ntri_th,start,last,ng,ng0,eig,vf,EF,mat,shiftk,query) \
444+ private(ib,j0,i0,i1,ithread)
445+ {
446+ int ntri0, i, j, i2, j1, j2, ii0, ii1, ii2;
447+ GLfloat kvec1[8][3], mat1[8][3], eig1[8], vf1[8][3];
448+
449+ ithread = get_thread();
450+ for (ib = 0; ib < nb; ++ib) {
451+
452+ if(query == 1) ntri0 = 0;
453+ else ntri0 = ntri_th[ib][ithread];
454+
455+#pragma omp for nowait schedule(static, 1)
456+ for (j0 = start[0]; j0 < last[0]; ++j0) {
457+ for (j1 = start[1]; j1 < last[1]; ++j1) {
458+ for (j2 = start[2]; j2 < last[2]; ++j2) {
459+ /**/
460+ i0 = j0;
461+ i1 = j1;
462+ i2 = j2;
463+ ii0 = j0 + 1;
464+ ii1 = j1 + 1;
465+ ii2 = j2 + 1;
466+ /**/
467+ kvec1[0][0] = (GLfloat)i0 / (GLfloat)ng[0];
468+ kvec1[1][0] = (GLfloat)i0 / (GLfloat)ng[0];
469+ kvec1[2][0] = (GLfloat)i0 / (GLfloat)ng[0];
470+ kvec1[3][0] = (GLfloat)i0 / (GLfloat)ng[0];
471+ kvec1[4][0] = (GLfloat)ii0 / (GLfloat)ng[0];
472+ kvec1[5][0] = (GLfloat)ii0 / (GLfloat)ng[0];
473+ kvec1[6][0] = (GLfloat)ii0 / (GLfloat)ng[0];
474+ kvec1[7][0] = (GLfloat)ii0 / (GLfloat)ng[0];
475+ /**/
476+ kvec1[0][1] = (GLfloat)i1 / (GLfloat)ng[1];
477+ kvec1[1][1] = (GLfloat)i1 / (GLfloat)ng[1];
478+ kvec1[2][1] = (GLfloat)ii1 / (GLfloat)ng[1];
479+ kvec1[3][1] = (GLfloat)ii1 / (GLfloat)ng[1];
480+ kvec1[4][1] = (GLfloat)i1 / (GLfloat)ng[1];
481+ kvec1[5][1] = (GLfloat)i1 / (GLfloat)ng[1];
482+ kvec1[6][1] = (GLfloat)ii1 / (GLfloat)ng[1];
483+ kvec1[7][1] = (GLfloat)ii1 / (GLfloat)ng[1];
484+ /**/
485+ kvec1[0][2] = (GLfloat)i2 / (GLfloat)ng[2];
486+ kvec1[1][2] = (GLfloat)ii2 / (GLfloat)ng[2];
487+ kvec1[2][2] = (GLfloat)i2 / (GLfloat)ng[2];
488+ kvec1[3][2] = (GLfloat)ii2 / (GLfloat)ng[2];
489+ kvec1[4][2] = (GLfloat)i2 / (GLfloat)ng[2];
490+ kvec1[5][2] = (GLfloat)ii2 / (GLfloat)ng[2];
491+ kvec1[6][2] = (GLfloat)i2 / (GLfloat)ng[2];
492+ kvec1[7][2] = (GLfloat)ii2 / (GLfloat)ng[2];
493+ /**/
494+ for (i = 0; i < 8; i++)
495+ for (j = 0; j < 3; j++)
496+ kvec1[i][j] = kvec1[i][j] + (GLfloat)shiftk[j] / (GLfloat)(2 * ng0[j]);
497+ /**/
498+ i0 = modulo(i0, ng[0]);
499+ i1 = modulo(i1, ng[1]);
500+ i2 = modulo(i2, ng[2]);
501+ ii0 = modulo(ii0, ng[0]);
502+ ii1 = modulo(ii1, ng[1]);
503+ ii2 = modulo(ii2, ng[2]);
504+ /**/
505+ eig1[0] = eig[ib][i0][i1][i2] - EF;
506+ eig1[1] = eig[ib][i0][i1][ii2] - EF;
507+ eig1[2] = eig[ib][i0][ii1][i2] - EF;
508+ eig1[3] = eig[ib][i0][ii1][ii2] - EF;
509+ eig1[4] = eig[ib][ii0][i1][i2] - EF;
510+ eig1[5] = eig[ib][ii0][i1][ii2] - EF;
511+ eig1[6] = eig[ib][ii0][ii1][i2] - EF;
512+ eig1[7] = eig[ib][ii0][ii1][ii2] - EF;
513+ /**/
514+ for (j = 0; j < 3; j++) {
515+ mat1[0][j] = mat[ib][i0][i1][i2][j];
516+ mat1[1][j] = mat[ib][i0][i1][ii2][j];
517+ mat1[2][j] = mat[ib][i0][ii1][i2][j];
518+ mat1[3][j] = mat[ib][i0][ii1][ii2][j];
519+ mat1[4][j] = mat[ib][ii0][i1][i2][j];
520+ mat1[5][j] = mat[ib][ii0][i1][ii2][j];
521+ mat1[6][j] = mat[ib][ii0][ii1][i2][j];
522+ mat1[7][j] = mat[ib][ii0][ii1][ii2][j];
523+ /**/
524+ vf1[0][j] = vf[ib][i0][i1][i2][j];
525+ vf1[1][j] = vf[ib][i0][i1][ii2][j];
526+ vf1[2][j] = vf[ib][i0][ii1][i2][j];
527+ vf1[3][j] = vf[ib][i0][ii1][ii2][j];
528+ vf1[4][j] = vf[ib][ii0][i1][i2][j];
529+ vf1[5][j] = vf[ib][ii0][i1][ii2][j];
530+ vf1[6][j] = vf[ib][ii0][ii1][i2][j];
531+ vf1[7][j] = vf[ib][ii0][ii1][ii2][j];
532+ }/*for (j = 0; j < 3; j++)*/
533+ /**/
534+ tetrahedron(ib, &ntri0, eig1, mat1, kvec1, vf1);
535+ }/*for (j0 = start[0]; j0 < ng[0]; ++j0)*/
536+ }/*for (j1 = start[1]; j1 < ng[1]; ++j1)*/
537+ }/*for (j0 = start[0]; j0 < ng[0]; ++j0)*/
538+ if (query == 1) ntri_th[ib][ithread] = ntri0;
539+ }/*for (ib = 0; ib < nb; ++ib)*/
540+ } /* End of parallel region */
541+ /*
542+ If this run is workspace query, malloc arrays
543+ */
544+ if (query == 1) {
545+ /*
546+ Sum patches in all threads
547+ */
548+ for (ib = 0; ib < nb; ib++) {
549+ for (ithread = 1; ithread < nthreads; ithread++) {
550+ ntri_th[ib][ithread] += ntri_th[ib][ithread - 1];
551+ }
552+ ntri[ib] = ntri_th[ib][nthreads - 1];
553+ for (ithread = nthreads - 1; ithread > 0; ithread--) {
554+ ntri_th[ib][ithread] = ntri_th[ib][ithread - 1];
555+ }
556+ ntri_th[ib][0] = 0;
557+ }
558+ /*
559+ Allocation of triangler patches
560+ */
561+ matp = new GLfloat***[nb];
562+ clr = new GLfloat*[nb];
563+ kvp = new GLfloat***[nb];
564+ nmlp = new GLfloat***[nb];
565+ arw = new GLfloat * ***[nb];
566+ kvp_rot = new GLfloat*[nb];
567+ nmlp_rot = new GLfloat*[nb];
568+ arw_rot = new GLfloat * [nb];
569+ for (ib = 0; ib < nb; ++ib) {
570+ matp[ib] = new GLfloat**[ntri[ib]];
571+ clr[ib] = new GLfloat[12 * ntri[ib]];
572+ kvp[ib] = new GLfloat**[ntri[ib]];
573+ nmlp[ib] = new GLfloat**[ntri[ib]];
574+ arw[ib] = new GLfloat ***[ntri[ib]];
575+ kvp_rot[ib] = new GLfloat[9 * ntri[ib]];
576+ nmlp_rot[ib] = new GLfloat[9 * ntri[ib]];
577+ arw_rot[ib] = new GLfloat[18 * ntri[ib]];
578+ for (i0 = 0; i0 < ntri[ib]; ++i0) {
579+ matp[ib][i0] = new GLfloat*[3];
580+ kvp[ib][i0] = new GLfloat*[3];
581+ nmlp[ib][i0] = new GLfloat*[3];
582+ arw[ib][i0] = new GLfloat ** [3];
583+ for (i1 = 0; i1 < 3; ++i1) {
584+ matp[ib][i0][i1] = new GLfloat[3];
585+ kvp[ib][i0][i1] = new GLfloat[3];
586+ nmlp[ib][i0][i1] = new GLfloat[3];
587+ arw[ib][i0][i1] = new GLfloat*[2];
588+ for (j0 = 0; j0 < 2; ++j0)
589+ arw[ib][i0][i1][j0] = new GLfloat[3];
590+ }/*for (i1 = 0; i1 < 3; ++i1)*/
591+ }/*for (i0 = 0; i0 < ntri[ib]; ++i0)*/
592+ }/*for (ib = 0; ib < nb; ++ib)*/
593+ }/*if (query == 1)*/
594+} /* fermi_patch */
--- /dev/null
+++ b/Android1.NativeActivity/fermi_patch.hpp
@@ -0,0 +1,25 @@
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+
25+void fermi_patch();
--- /dev/null
+++ b/Android1.NativeActivity/free_patch.cpp
@@ -0,0 +1,529 @@
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 Refresh patch
26+*/
27+#include <stdlib.h>
28+#include <math.h>
29+#include <stdio.h>
30+#include "basic_math.hpp"
31+#include "variable.hpp"
32+
33+/**
34+ @brief Free variables for patch before new patch is computed
35+
36+ Free : ::nmlp, ::matp, ::clr, ::kvp, ::nmlp_rot, ::kvp_rot,
37+ ::kvnl, ::kvnl_rot, ::kv2d, ::clr2d
38+*/
39+void free_patch() {
40+ int ib, i0, i1, i2;
41+ /*
42+ Fermi patch
43+ */
44+ if (refresh_patch == 1) {
45+ for (ib = 0; ib < nb; ++ib) {
46+ for (i0 = 0; i0 < ntri[ib]; ++i0) {
47+ for (i1 = 0; i1 < 3; ++i1) {
48+ for (i2 = 0; i2 < 2; ++i2)
49+ delete[] arw[ib][i0][i1][i2];
50+ delete[] nmlp[ib][i0][i1];
51+ delete[] matp[ib][i0][i1];
52+ delete[] kvp[ib][i0][i1];
53+ delete[] arw[ib][i0][i1];
54+ }
55+ delete[] nmlp[ib][i0];
56+ delete[] matp[ib][i0];
57+ delete[] kvp[ib][i0];
58+ delete[] arw[ib][i0];
59+ }
60+ delete[] nmlp[ib];
61+ delete[] matp[ib];
62+ delete[] clr[ib];
63+ delete[] kvp[ib];
64+ delete[] arw[ib];
65+ delete[] nmlp_rot[ib];
66+ delete[] kvp_rot[ib];
67+ delete[] arw_rot[ib];
68+ }
69+ delete[] nmlp;
70+ delete[] matp;
71+ delete[] clr;
72+ delete[] kvp;
73+ delete[] arw;
74+ delete[] nmlp_rot;
75+ delete[] kvp_rot;
76+ delete[] arw_rot;
77+ }/*if (refresh_patch == 1)*/
78+ /*
79+ Nodal line
80+ */
81+ if (refresh_nodeline == 1) {
82+ for (ib = 0; ib < nb; ++ib) {
83+ for (i0 = 0; i0 < nnl[ib]; ++i0) {
84+ for (i1 = 0; i1 < 2; ++i1) {
85+ delete[] kvnl[ib][i0][i1];
86+ }/*for (i1 = 0; i1 < 2; ++i1)*/
87+ delete[] kvnl[ib][i0];
88+ }/*for (i0 = 0; i0 < nnl[ib]; ++i0)*/
89+ delete[] kvnl[ib];
90+ delete[] kvnl_rot[ib];
91+ }/*for (ib = 0; ib < nb; ++ib)*/
92+ delete[] kvnl;
93+ delete[] kvnl_rot;
94+ }/*if (refresh_nodeline == 1)*/
95+ /*
96+ 2D Fermi line
97+ */
98+ if (refresh_section == 1) {
99+ for (ib = 0; ib < nb; ++ib) {
100+ delete[] kv2d[ib];
101+ delete[] clr2d[ib];
102+ }/*for (ib = 0; ib < nb; ++ib)*/
103+ delete[] kv2d;
104+ delete[] clr2d;
105+ }/*if (refresh_section == 1)*/
106+ /*
107+ equator
108+ */
109+ if (refresh_equator == 1) {
110+ for (ib = 0; ib < nb; ++ib) {
111+ for (i0 = 0; i0 < nequator[ib]; ++i0) {
112+ for (i1 = 0; i1 < 2; ++i1) {
113+ delete[] kveq[ib][i0][i1];
114+ }/*for (i1 = 0; i1 < 2; ++i1)*/
115+ delete[] kveq[ib][i0];
116+ }/*for (i0 = 0; i0 < nequator[ib]; ++i0)*/
117+ delete[] kveq[ib];
118+ delete[] kveq_rot[ib];
119+ }/*for (ib = 0; ib < nb; ++ib)*/
120+ delete[] kveq;
121+ delete[] kveq_rot;
122+ }/*if (refresh_equator == 1)*/
123+}/*void free_patch()*/
124+/**
125+ @brief Compute Max. & Min. of matrix elements.
126+ Compute color of each patch
127+
128+ Modify : ::clr
129+*/
130+void max_and_min()
131+{
132+ int itri, ithread;
133+ GLfloat *max_th, *min_th;
134+
135+ max_th = new GLfloat[nthreads];
136+ min_th = new GLfloat[nthreads];
137+ /*
138+ Search max and min.
139+ */
140+ if (color_scale == 1 || color_scale == 6) {
141+#pragma omp parallel default(none) \
142+shared(nb,ntri,matp,max_th,min_th) private(itri,ithread)
143+ {
144+ int i, ib;
145+
146+ ithread = get_thread();
147+ max_th[ithread] = -1.0e10f;
148+ min_th[ithread] = 1.0e10f;
149+
150+ for (ib = 0; ib < nb; ib++) {
151+#pragma omp for
152+ for (itri = 0; itri < ntri[ib]; ++itri) {
153+ for (i = 0; i < 3; ++i) {
154+ if (matp[ib][itri][i][0] > max_th[ithread]) max_th[ithread] = matp[ib][itri][i][0];
155+ if (matp[ib][itri][i][0] < min_th[ithread]) min_th[ithread] = matp[ib][itri][i][0];
156+ }
157+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
158+ }/*for (ib = 0; ib < nb; ib++)*/
159+ }/*End of parallel region*/
160+ /**/
161+ patch_max = max_th[0];
162+ patch_min = min_th[0];
163+ for (ithread = 1; ithread < nthreads; ithread++) {
164+ if (max_th[ithread] > patch_max) patch_max = max_th[ithread];
165+ if (min_th[ithread] < patch_min) patch_min = min_th[ithread];
166+ }
167+ }/*if (color_scale == 0 || color_scale == 4)*/
168+ else if (color_scale == 2) {
169+#pragma omp parallel default(none) \
170+shared(nb,ntri,matp,max_th,min_th) private(itri,ithread)
171+ {
172+ int i, ib;
173+ GLfloat abs;
174+
175+ ithread = get_thread();
176+ max_th[ithread] = -1.0e10f;
177+ min_th[ithread] = 1.0e10f;
178+
179+ for (ib = 0; ib < nb; ib++) {
180+#pragma omp for
181+ for (itri = 0; itri < ntri[ib]; ++itri) {
182+ for (i = 0; i < 3; ++i) {
183+ abs = sqrtf(matp[ib][itri][i][0] * matp[ib][itri][i][0]
184+ + matp[ib][itri][i][1] * matp[ib][itri][i][1]);
185+ if (abs > max_th[ithread]) max_th[ithread] = abs;
186+ if (abs < min_th[ithread]) min_th[ithread] = abs;
187+ }
188+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
189+ }/*for (ib = 0; ib < nb; ib++)*/
190+ }/*End of parallel region*/
191+ /**/
192+ patch_min = min_th[0];
193+ patch_max = max_th[0];
194+ for (ithread = 1; ithread < nthreads; ithread++) {
195+ if (max_th[ithread] < patch_min) patch_min = max_th[ithread];
196+ if (max_th[ithread] > patch_max) patch_max = max_th[ithread];
197+ }
198+ }/*if (color_scale == 2)*/
199+ else if (color_scale == 3) {
200+#pragma omp parallel default(none) \
201+shared(nb,ntri,matp,min_th,max_th) private(itri,ithread)
202+ {
203+ int i, ib;
204+ GLfloat abs;
205+
206+ ithread = get_thread();
207+ min_th[ithread] = 1.0e10f;
208+ max_th[ithread] = -1.0e10f;
209+
210+ for (ib = 0; ib < nb; ib++) {
211+#pragma omp for
212+ for (itri = 0; itri < ntri[ib]; ++itri) {
213+ for (i = 0; i < 3; ++i) {
214+ abs = sqrtf(matp[ib][itri][i][0] * matp[ib][itri][i][0]
215+ + matp[ib][itri][i][1] * matp[ib][itri][i][1]
216+ + matp[ib][itri][i][2] * matp[ib][itri][i][2]);
217+ if (abs > max_th[ithread]) max_th[ithread] = abs;
218+ if (abs < min_th[ithread]) min_th[ithread] = abs;
219+ }
220+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
221+ }/*for (ib = 0; ib < nb; ib++)*/
222+ }/*End of parallel region*/
223+ /**/
224+ patch_max = max_th[0];
225+ patch_min = min_th[0];
226+ for (ithread = 1; ithread < nthreads; ithread++) {
227+ if (max_th[ithread] > patch_max) patch_max = max_th[ithread];
228+ if (min_th[ithread] < patch_min) patch_min = min_th[ithread];
229+ }
230+ }/*if (color_scale == 3)*/
231+ else if (color_scale == 4 || color_scale == 7) {
232+#pragma omp parallel default(none) \
233+shared(nb,ntri,nmlp,max_th,min_th) private(itri,ithread)
234+ {
235+ int i, j, ib;
236+ GLfloat norm;
237+
238+ ithread = get_thread();
239+ max_th[ithread] = -1.0e10f;
240+ min_th[ithread] = 1.0e10f;
241+
242+ for (ib = 0; ib < nb; ib++) {
243+#pragma omp for
244+ for (itri = 0; itri < ntri[ib]; ++itri) {
245+ for (i = 0; i < 3; ++i) {
246+ norm = 0.0f;
247+ for (j = 0; j < 3; ++j) norm += nmlp[ib][itri][i][j]*nmlp[ib][itri][i][j];
248+ norm = sqrtf(norm);
249+
250+ if (norm > max_th[ithread]) max_th[ithread] = norm;
251+ if (norm < min_th[ithread]) min_th[ithread] = norm;
252+ }
253+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
254+ }/*for (ib = 0; ib < nb; ib++)*/
255+ }/*End of parallel region*/
256+ /**/
257+ patch_max = max_th[0];
258+ patch_min = min_th[0];
259+ for (ithread = 1; ithread < nthreads; ithread++) {
260+ if (max_th[ithread] > patch_max) patch_max = max_th[ithread];
261+ if (min_th[ithread] < patch_min) patch_min = min_th[ithread];
262+ }
263+ }/*if (color_scale == 5 || color_scale == 6)*/
264+
265+ delete[] max_th;
266+ delete[] min_th;
267+}/* max_and_min */
268+ /**
269+ @brief Compute Max. & Min. of matrix elements.
270+ Compute color of each patch
271+
272+ Modify : ::clr
273+ */
274+void paint()
275+{
276+ int itri, j;
277+ GLfloat origin[4];
278+
279+ if (color_scale == 1) {
280+#pragma omp parallel default(none) \
281+shared(nb,ntri,matp,clr,cyan,blue,green,yellow,red,patch_max,patch_min) \
282+private(itri, j)
283+ {
284+ int i, ib;
285+ GLfloat mat2;
286+
287+ for (ib = 0; ib < nb; ib++) {
288+#pragma omp for nowait
289+ for (itri = 0; itri < ntri[ib]; ++itri) {
290+ for (i = 0; i < 3; ++i) {
291+ /**/
292+ mat2 = (matp[ib][itri][i][0] - patch_min) / (patch_max - patch_min);
293+ mat2 = mat2 * 4.0f;
294+ /**/
295+ if (mat2 <= 1.0) {
296+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = cyan[j] * mat2 + blue[j] * (1.0f - mat2);
297+ }
298+ else if (mat2 <= 2.0) {
299+ mat2 = mat2 - 1.0f;
300+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = green[j] * mat2 + cyan[j] * (1.0f - mat2);
301+ }
302+ else if (mat2 <= 3.0) {
303+ mat2 = mat2 - 2.0f;
304+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = yellow[j] * mat2 + green[j] * (1.0f - mat2);
305+ }
306+ else {
307+ mat2 = mat2 - 3.0f;
308+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = red[j] * mat2 + yellow[j] * (1.0f - mat2);
309+ }
310+ }/*for (i = 0; i < 3; ++i)*/
311+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
312+ }/*for (ib = 0; ib < nb; ib++)*/
313+ }/*End of parallel region*/
314+ }/*if (color_scale == 1 || color_scale == 2)*/
315+ else if (color_scale == 2) {
316+
317+ if (blackback == 1) for (j = 0; j < 4; ++j) origin[j] = wgray[j];
318+ else for (j = 0; j < 4; ++j) origin[j] = bgray[j];
319+
320+#pragma omp parallel default(none) \
321+shared(nb,ntri,matp,clr,cyan,blue,green,yellow,red,magenta,bgray,wgray,blackback,patch_max,origin) \
322+private(itri, j)
323+ {
324+ int i, ib;
325+ GLfloat theta, abs, theta2;
326+
327+ for (ib = 0; ib < nb; ib++) {
328+ for (itri = 0; itri < ntri[ib]; ++itri) {
329+ for (i = 0; i < 3; ++i) {
330+ /**/
331+ abs = sqrtf(matp[ib][itri][i][0] * matp[ib][itri][i][0]
332+ + matp[ib][itri][i][1] * matp[ib][itri][i][1]);
333+ if (abs / patch_max < 0.00001) theta = 0.0f;
334+ else if (matp[ib][itri][i][1] > 0.0) theta = acosf(matp[ib][itri][i][0] / abs);
335+ else theta = -acosf(matp[ib][itri][i][0] / abs);
336+ abs /= patch_max;
337+ theta = 3.0f * theta / acosf(-1.0f);
338+ /**/
339+ if (-3.0f <= theta && theta < -2.0f) {
340+ theta2 = theta + 3.0f;
341+ for (j = 0; j < 4; ++j)
342+ clr[ib][j + 4 * i + 12 * itri] = blue[j] * theta2 + cyan[j] * (1.0f - theta2);
343+ }
344+ else if (-2.0f <= theta && theta < -1.0f) {
345+ theta2 = theta + 2.0f;
346+ for (j = 0; j < 4; ++j)
347+ clr[ib][j + 4 * i + 12 * itri] = magenta[j] * theta2 + blue[j] * (1.0f - theta2);
348+ }
349+ else if (-1.0f <= theta && theta < 0.0f) {
350+ theta2 = theta + 1.0f;
351+ for (j = 0; j < 4; ++j)
352+ clr[ib][j + 4 * i + 12 * itri] = red[j] * theta2 + magenta[j] * (1.0f - theta2);
353+ }
354+ else if (0.0f <= theta && theta < 1.0f) {
355+ theta2 = theta;
356+ for (j = 0; j < 4; ++j)
357+ clr[ib][j + 4 * i + 12 * itri] = yellow[j] * theta2 + red[j] * (1.0f - theta2);
358+ }
359+ else if (1.0f <= theta && theta < 2.0f) {
360+ theta2 = theta - 1.0f;
361+ for (j = 0; j < 4; ++j)
362+ clr[ib][j + 4 * i + 12 * itri] = green[j] * theta2 + yellow[j] * (1.0f - theta2);
363+ }
364+ else {
365+ theta2 = theta - 2.0f;
366+ for (j = 0; j < 4; ++j)
367+ clr[ib][j + 4 * i + 12 * itri] = cyan[j] * theta2 + green[j] * (1.0f - theta2);
368+ }
369+ clr[ib][j + 4 * i + 12 * itri] = clr[ib][j + 4 * i + 12 * itri] * abs + origin[j] * (1.0f - abs);
370+
371+ }/*for (i = 0; i < 3; ++i)*/
372+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
373+ }/*for (ib = 0; ib < nb; ib++)*/
374+ }/*End of parallel region*/
375+ }/*if (color_scale == 2)*/
376+ else if (color_scale == 4) {
377+#pragma omp parallel default(none) \
378+shared(nb,ntri,nmlp,clr,cyan,blue,green,yellow,red,patch_max,patch_min) \
379+private(itri, j)
380+ {
381+ int i, ib;
382+ GLfloat mat2;
383+
384+ for (ib = 0; ib < nb; ib++) {
385+ #pragma omp for nowait
386+ for (itri = 0; itri < ntri[ib]; ++itri) {
387+ for (i = 0; i < 3; ++i) {
388+ /**/
389+ mat2 = 0.0f;
390+ for (j = 0; j < 3; ++j) mat2 += nmlp[ib][itri][i][j] * nmlp[ib][itri][i][j];
391+ mat2 = sqrtf(mat2);
392+ mat2 = (mat2 - patch_min) / (patch_max - patch_min);
393+ mat2 = mat2 * 4.0f;
394+ /**/
395+ if (mat2 <= 1.0) {
396+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = cyan[j] * mat2 + blue[j] * (1.0f - mat2);
397+ }
398+ else if (mat2 <= 2.0) {
399+ mat2 = mat2 - 1.0f;
400+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = green[j] * mat2 + cyan[j] * (1.0f - mat2);
401+ }
402+ else if (mat2 <= 3.0) {
403+ mat2 = mat2 - 2.0f;
404+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = yellow[j] * mat2 + green[j] * (1.0f - mat2);
405+ }
406+ else {
407+ mat2 = mat2 - 3.0f;
408+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = red[j] * mat2 + yellow[j] * (1.0f - mat2);
409+ }
410+ }/*for (i = 0; i < 3; ++i)*/
411+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
412+ }/*for (ib = 0; ib < nb; ib++)*/
413+ }/*End of parallel region*/
414+ }/*if (color_scale == 4)*/
415+ else if (color_scale == 3 || color_scale == 5) {
416+#pragma omp parallel default(none) \
417+shared(nb,ntri,matp,clr,cyan,blue,green,yellow,red,color_scale,kvp,arw,patch_max) \
418+private(itri, j)
419+ {
420+ int i, ib;
421+ GLfloat mat2;
422+
423+ for (ib = 0; ib < nb; ib++) {
424+ /**/
425+ if (nb == 1) mat2 = 0.5f;
426+ else mat2 = 1.0f / (GLfloat)(nb - 1) * (GLfloat)ib;
427+ mat2 *= 4.0f;
428+ /**/
429+ if (mat2 <= 1.0) {
430+#pragma omp for nowait
431+ for (itri = 0; itri < ntri[ib]; ++itri) {
432+ for (i = 0; i < 3; ++i) {
433+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = cyan[j] * mat2 + blue[j] * (1.0f - mat2);
434+ }
435+ }
436+ }
437+ else if (mat2 <= 2.0) {
438+ mat2 = mat2 - 1.0f;
439+#pragma omp for nowait
440+ for (itri = 0; itri < ntri[ib]; ++itri) {
441+ for (i = 0; i < 3; ++i) {
442+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = green[j] * mat2 + cyan[j] * (1.0f - mat2);
443+ }
444+ }
445+ }
446+ else if (mat2 <= 3.0) {
447+ mat2 = mat2 - 2.0f;
448+#pragma omp for nowait
449+ for (itri = 0; itri < ntri[ib]; ++itri) {
450+ for (i = 0; i < 3; ++i) {
451+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = yellow[j] * mat2 + green[j] * (1.0f - mat2);
452+ }
453+ }
454+ }
455+ else {
456+ mat2 = mat2 - 3.0f;
457+#pragma omp for nowait
458+ for (itri = 0; itri < ntri[ib]; ++itri) {
459+ for (i = 0; i < 3; ++i) {
460+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = red[j] * mat2 + yellow[j] * (1.0f - mat2);
461+ }
462+ }
463+ }
464+ }/*for (ib = 0; ib < nb; ib++*/
465+
466+ if (color_scale == 3) {
467+ for (ib = 0; ib < nb; ib++) {
468+#pragma omp for nowait
469+ for (itri = 0; itri < ntri[ib]; ++itri) {
470+ for (i = 0; i < 3; ++i) {
471+ for (j = 0; j < 3; ++j) {
472+ arw[ib][itri][i][0][j] = kvp[ib][itri][i][j]
473+ - 0.5f * matp[ib][itri][i][j] / patch_max;
474+ arw[ib][itri][i][1][j] = kvp[ib][itri][i][j]
475+ + 0.5f * matp[ib][itri][i][j] / patch_max;
476+ }/*for (j = 0; j < 3; ++j)*/
477+ }/*for (i = 0; i < 3; ++i)*/
478+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
479+ }/*for (ib = 0; ib < nb; ib++)*/
480+ }/*if (color_scale == 3)*/
481+ }/*End of parallel region*/
482+ }/*if (color_scale == 5)*/
483+ else if (color_scale == 6) {
484+#pragma omp parallel default(none) \
485+shared(nb,ntri,matp,clr,bgray,wgray,patch_max,patch_min) \
486+private(itri, j)
487+ {
488+ int i, ib;
489+ GLfloat mat2;
490+
491+ for (ib = 0; ib < nb; ib++) {
492+#pragma omp for nowait
493+ for (itri = 0; itri < ntri[ib]; ++itri) {
494+ for (i = 0; i < 3; ++i) {
495+ /**/
496+ mat2 = (matp[ib][itri][i][0] - patch_min) / (patch_max - patch_min);
497+ /**/
498+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = wgray[j] * mat2 + bgray[j] * (1.0f - mat2);
499+ }/*for (i = 0; i < 3; ++i)*/
500+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
501+ }/*for (ib = 0; ib < nb; ib++)*/
502+ }/*End of parallel region*/
503+ }/*if (color_scale == 6)*/
504+ else if (color_scale == 7) {
505+#pragma omp parallel default(none) \
506+shared(nb,ntri,nmlp,clr,bgray,wgray,patch_max,patch_min) \
507+private(itri, j)
508+ {
509+ int i, ib;
510+ GLfloat mat2;
511+
512+ for (ib = 0; ib < nb; ib++) {
513+#pragma omp for nowait
514+ for (itri = 0; itri < ntri[ib]; ++itri) {
515+ for (i = 0; i < 3; ++i) {
516+ /**/
517+ mat2 = 0.0f;
518+ for (j = 0; j < 3; ++j) mat2 += nmlp[ib][itri][i][j] * nmlp[ib][itri][i][j];
519+ mat2 = sqrtf(mat2);
520+ mat2 = (mat2 - patch_min) / (patch_max - patch_min);
521+ /**/
522+ for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = wgray[j] * mat2 + bgray[j] * (1.0f - mat2);
523+ }/*for (i = 0; i < 3; ++i)*/
524+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
525+ }/*for (ib = 0; ib < nb; ib++)*/
526+ }/*End of parallel region*/
527+ }/*if (color_scale == 7)*/
528+
529+}/* paint */
--- /dev/null
+++ b/Android1.NativeActivity/free_patch.hpp
@@ -0,0 +1,28 @@
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+
25+void free_patch();
26+void max_and_min();
27+void paint();
28+
--- /dev/null
+++ b/Android1.NativeActivity/initialize.cpp
@@ -0,0 +1,257 @@
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 Functions that initilize variables.
26+*/
27+#include <stdio.h>
28+#include "variable.hpp"
29+
30+/**
31+ @brief Specify corners of tetrahedron
32+
33+ Modify : ::corner
34+*/
35+void init_corner() {
36+ int i, j;
37+ int corner1[16][6][4] = {
38+ /*
39+ [0] min = 0-7
40+ */
41+ { { 0, 1, 3, 7 },
42+ { 0, 1, 5, 7 },
43+ { 0, 2, 3, 7 },
44+ { 0, 2, 6, 7 },
45+ { 0, 4, 5, 7 },
46+ { 0, 4, 6, 7 } },
47+ /*
48+ [1] min = 1-6
49+ */
50+ { { 1, 0, 2, 6 },
51+ { 1, 0, 4, 6 },
52+ { 1, 3, 2, 6 },
53+ { 1, 3, 7, 6 },
54+ { 1, 5, 4, 6 },
55+ { 1, 5, 7, 6 } },
56+ /*
57+ [2] min = 2-5
58+ */
59+ { { 2, 0, 1, 5 },
60+ { 2, 0, 4, 5 },
61+ { 2, 3, 1, 5 },
62+ { 2, 3, 7, 5 },
63+ { 2, 6, 4, 5 },
64+ { 2, 6, 7, 5 } },
65+ /*
66+ [3] min = 3-4
67+ */
68+ { { 4, 0, 1, 3 },
69+ { 4, 0, 2, 3 },
70+ { 4, 5, 1, 3 },
71+ { 4, 5, 7, 3 },
72+ { 4, 6, 2, 3 },
73+ { 4, 6, 7, 3 } },
74+ /*
75+ [4] min = 0-7, max = 3-4
76+ */
77+ { { 0, 4, 5, 6 },
78+ { 1, 2, 3, 7 },
79+ { 0, 7, 2, 6 },
80+ { 0, 7, 1, 5 },
81+ { 0, 7, 1, 2 },
82+ { 0, 7, 6, 5 } },
83+ /*
84+ [5] min = 1-6, max = 3-4
85+ */
86+ { { 0, 4, 5, 6 },
87+ { 1, 2, 3, 7 },
88+ { 1, 6, 5, 7 },
89+ { 1, 6, 7, 2 },
90+ { 1, 6, 2, 0 },
91+ { 1, 6, 0, 5 } },
92+ /*
93+ [6] min = 2-5, max = 3-4
94+ */
95+ { { 0, 4, 5, 6 },
96+ { 1, 2, 3, 7 },
97+ { 2, 5, 7, 6 },
98+ { 2, 5, 6, 0 },
99+ { 2, 5, 0, 1 },
100+ { 2, 5, 1, 7 } },
101+ /*
102+ [7] min = 3-4, max = 0-7
103+ */
104+ { { 0, 1, 2, 4 },
105+ { 7, 3, 5, 6 },
106+ { 3, 4, 1, 5 },
107+ { 3, 4, 5, 6 },
108+ { 3, 4, 6, 2 },
109+ { 3, 4, 2, 1 } },
110+ /*
111+ [8] min = 2-5, max = 0-7
112+ */
113+ { { 0, 1, 2, 4 },
114+ { 7, 3, 5, 6 },
115+ { 2, 5, 1, 3 },
116+ { 2, 5, 3, 6 },
117+ { 2, 5, 6, 4 },
118+ { 2, 5, 4, 1 } },
119+ /*
120+ [9] min = 1-6, max = 0-7
121+ */
122+ { { 0, 1, 2, 4 },
123+ { 7, 3, 5, 6 },
124+ { 1, 6, 2, 3 },
125+ { 1, 6, 3, 5 },
126+ { 1, 6, 5, 4 },
127+ { 1, 6, 4, 2 } },
128+ /*
129+ [10] min = 0-7, max = 1-6
130+ */
131+ { { 1, 0, 3, 5 },
132+ { 6, 2, 4, 7 },
133+ { 0, 7, 2, 3 },
134+ { 0, 7, 3, 5 },
135+ { 0, 7, 5, 4 },
136+ { 0, 7, 4, 2 } },
137+ /*
138+ [11] min = 3-4, max = 1-6
139+ */
140+ { { 1, 0, 3, 5 },
141+ { 6, 2, 4, 7 },
142+ { 3, 4, 0, 2 },
143+ { 3, 4, 2, 7 },
144+ { 3, 4, 7, 5 },
145+ { 3, 4, 5, 0 } },
146+ /*
147+ [12] min = 2-5, max = 1-6
148+ */
149+ { { 1, 0, 3, 5 },
150+ { 6, 2, 4, 7 },
151+ { 2, 5, 0, 3 },
152+ { 2, 5, 3, 7 },
153+ { 2, 5, 7, 4 },
154+ { 2, 5, 4, 0 } },
155+ /*
156+ [13] min = 0-7, max = 2-5
157+ */
158+ { { 2, 0, 3, 6 },
159+ { 5, 1, 4, 7 },
160+ { 0, 7, 1, 3 },
161+ { 0, 7, 3, 6 },
162+ { 0, 7, 6, 4 },
163+ { 0, 7, 4, 1 } },
164+ /*
165+ [14] min = 1-6, max = 2-5
166+ */
167+ { { 2, 0, 3, 6 },
168+ { 5, 1, 4, 7 },
169+ { 1, 6, 0, 3 },
170+ { 1, 6, 3, 7 },
171+ { 1, 6, 7, 4 },
172+ { 1, 6, 4, 0 } },
173+ /*
174+ [15] min = 3-4, max = 2-5
175+ */
176+ { { 2, 0, 3, 6 },
177+ { 5, 1, 4, 7 },
178+ { 3, 4, 0, 1 },
179+ { 3, 4, 1, 7 },
180+ { 3, 4, 7, 6 },
181+ { 3, 4, 6, 0 } },
182+ };
183+ /**/
184+ for (i = 0; i < 6; ++i) {
185+ for (j = 0; j < 4; ++j) {
186+ corner[i][j] = corner1[itet][i][j];
187+ }
188+ }
189+}
190+/**
191+ @brief Compute Bragg vetor
192+
193+ Modify : ::bragg, ::brnrm
194+*/
195+void bragg_vector() {
196+ int i0, i1, i2, i, ibr;
197+ /**/
198+ ibr = 0;
199+ /**/
200+ for (i0 = -1; i0 <= 1; ++i0) {
201+ for (i1 = -1; i1 <= 1; ++i1) {
202+ for (i2 = -1; i2 <= 1; ++i2) {
203+ /*
204+ Excepte Gamma points
205+ */
206+ if (i0 == 0 && i1 == 0 && i2 == 0) continue;
207+ /*
208+ Fractional -> Cartecian
209+ */
210+ for (i = 0; i < 3; ++i)
211+ bragg[ibr][i] = ((GLfloat)i0 * bvec[0][i]
212+ + (GLfloat)i1 * bvec[1][i]
213+ + (GLfloat)i2 * bvec[2][i]) * 0.5f;
214+ /*
215+ And its norm
216+ */
217+ brnrm[ibr] = bragg[ibr][0] * bragg[ibr][0]
218+ + bragg[ibr][1] * bragg[ibr][1]
219+ + bragg[ibr][2] * bragg[ibr][2];
220+ /**/
221+ ibr = ibr + 1;
222+ }/*for (i2 = -1; i2 <= 1; ++i2)*/
223+ }/*for (i1 = -1; i1 <= 1; ++i1)*/
224+ }/*for (i0 = -1; i0 <= 1; ++i0)*/
225+ /*
226+ Search min. of brnrm
227+ */
228+ brnrm_min = brnrm[0];
229+ for (ibr = 1; ibr < 26; ibr++) {
230+ if (brnrm_min > brnrm[ibr]) brnrm_min = brnrm[ibr];
231+ }
232+}/* bragg_vector */
233+/**
234+ @brief Print max and minimum @f$\varepsilon_{n k}, \Delta_{n k}@f$
235+ in the whole Brillouine zone
236+*/
237+void max_and_min_bz() {
238+ int ib, i0, i1, i2;
239+ GLfloat eigmin, eigmax, matmin, matmax;
240+ /**/
241+ for (ib = 0; ib < nb; ib++) {
242+ eigmax = eig0[0][0][0][0];
243+ eigmin = eig0[0][0][0][0];
244+ matmax = mat0[0][0][0][0][0];
245+ matmin = mat0[0][0][0][0][0];
246+ for (i0 = 0; i0 < ng0[0]; ++i0) {
247+ for (i1 = 0; i1 < ng0[1]; ++i1) {
248+ for (i2 = 0; i2 < ng0[2]; ++i2) {
249+ if (eig0[ib][i0][i1][i2] > eigmax) eigmax = eig0[ib][i0][i1][i2];
250+ if (eig0[ib][i0][i1][i2] < eigmin) eigmin = eig0[ib][i0][i1][i2];
251+ if (mat0[ib][i0][i1][i2][0] > matmax) matmax = mat0[ib][i0][i1][i2][0];
252+ if (mat0[ib][i0][i1][i2][0] < matmin) matmin = mat0[ib][i0][i1][i2][0];
253+ }/*for (i2 = 0; i2 < ng0[2]; ++i2)*/
254+ }/*for (i1 = 0; i1 < ng0[1]; ++i1)*/
255+ }/*for (i0 = 0; i0 < ng0[0]; ++i0)*/
256+ }/*for (ib = 0; ib < nb; ib++)*/
257+}/* max_and_min_bz */
--- /dev/null
+++ b/Android1.NativeActivity/initialize.hpp
@@ -0,0 +1,27 @@
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+
25+void init_corner();
26+void bragg_vector();
27+void max_and_min_bz();
--- /dev/null
+++ b/Android1.NativeActivity/kumo.cpp
@@ -0,0 +1,217 @@
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 Compute @f$\varepsilon_{n k}, \Delta_{n k}@f$ on
26+denser @f$k@f$-grid with French-curve (Kumo) interpolation
27+*/
28+#include <GLES/gl.h>
29+#include <stdlib.h>
30+#include <stdio.h>
31+#include "basic_math.hpp"
32+#include "variable.hpp"
33+/**
34+ @brief Compute coefficient for the French-curve (Kumo) interpolation
35+ @f[
36+ A^{\rm intp} = \sum_{i = 1}^4 C_i A_i^{\rm orig}
37+ @f]
38+*/
39+static void kumo_coef(
40+ int j, //!< [in] Interpolated grid index
41+ GLfloat *coef //!< [out] Coefficient of interpolation @f$C_i@f$
42+) {
43+ GLfloat x, mx;
44+ x = (GLfloat)j / (GLfloat)interpol;
45+ mx = 1.0f - x;
46+ coef[0] = -0.5f * x * mx * mx;
47+ coef[1] = mx * (mx*mx + 3.0f* x*mx + 0.5f* x* x);
48+ coef[2] = x * ( x* x + 3.0f*mx* x + 0.5f*mx*mx);
49+ coef[3] = -0.5f * x * x * mx;
50+}
51+/**
52+ @brief Interpolation of energy and matrix
53+ with the French-curve (Kumo) interpolation.
54+
55+ Modify : ::eig, ::mat
56+*/
57+void interpol_energy() {
58+ int ib, i0, i1, i2, ii;
59+ /*
60+ Reallocate
61+ */
62+ for (ib = 0; ib < nb; ib++) {
63+ for (i0 = 0; i0 < ng[0]; i0++) {
64+ for (i1 = 0; i1 < ng[1]; i1++) {
65+ for (i2 = 0; i2 < ng[2]; i2++) {
66+ delete[] vf[ib][i0][i1][i2];
67+ delete[] mat[ib][i0][i1][i2];
68+ }
69+ delete[] eig[ib][i0][i1];
70+ delete[] mat[ib][i0][i1];
71+ delete[] vf[ib][i0][i1];
72+ }/*for (i1 = 0; i1 < ng[1]; i1++)*/
73+ delete[] eig[ib][i0];
74+ delete[] mat[ib][i0];
75+ delete[] vf[ib][i0];
76+ }/*for (i0 = 0; i0 < ng[0]; i0++)*/
77+ delete[] eig[ib];
78+ delete[] mat[ib];
79+ delete[] vf[ib];
80+ }/*for (ib = 0; ib < nb; ib++)*/
81+ for (ii = 0; ii < 3; ii++)ng[ii] = ng0[ii] * interpol;
82+ /**/
83+ for (ib = 0; ib < nb; ib++) {
84+ eig[ib] = new GLfloat**[ng[0]];
85+ mat[ib] = new GLfloat***[ng[0]];
86+ vf[ib] = new GLfloat***[ng[0]];
87+ for (i0 = 0; i0 < ng[0]; i0++) {
88+ eig[ib][i0] = new GLfloat*[ng[1]];
89+ mat[ib][i0] = new GLfloat**[ng[1]];
90+ vf[ib][i0] = new GLfloat**[ng[1]];
91+ for (i1 = 0; i1 < ng[1]; i1++) {
92+ eig[ib][i0][i1] = new GLfloat[ng[2]];
93+ mat[ib][i0][i1] = new GLfloat*[ng[2]];
94+ vf[ib][i0][i1] = new GLfloat*[ng[2]];
95+ for (i2 = 0; i2 < ng[2]; i2++) {
96+ mat[ib][i0][i1][i2] = new GLfloat[3];
97+ vf[ib][i0][i1][i2] = new GLfloat[3];
98+ }
99+ }/*for (i1 = 0; i1 < ng[1]; i1++)*/
100+ }/*for (i0 = 0; i0 < ng[0]; i0++)*/
101+ }/*for (ib = 0; ib < nb; ib++)*/
102+ /*
103+ 3rd order - three dimensional Kumo interpolation
104+ */
105+#pragma omp parallel default(none) \
106+ shared(nb,ng0,ng,eig,eig0,mat,mat0,interpol) \
107+ private (ib,i0,i1,i2,ii)
108+ {
109+ int j0, j1, j2, jj;
110+ GLfloat coef[4],
111+ mat1[4][4][4][3], mat2[4][4][3], mat3[4][3],
112+ eig1[4][4][4], eig2[4][4], eig3[4];
113+
114+ for (ib = 0; ib < nb; ib++) {
115+# pragma omp for nowait
116+ for (i0 = 0; i0 < ng0[0]; i0++) {
117+ //if (ith == 1) continue;
118+ for (i1 = 0; i1 < ng0[1]; i1++) {
119+ for (i2 = 0; i2 < ng0[2]; i2++) {
120+ for (j0 = 0; j0 < 4; j0++) {
121+ for (j1 = 0; j1 < 4; j1++) {
122+ for (j2 = 0; j2 < 4; j2++) {
123+ eig1[j0][j1][j2] = eig0[ib][modulo(i0 + j0 - 1, ng0[0])]
124+ [modulo(i1 + j1 - 1, ng0[1])]
125+ [modulo(i2 + j2 - 1, ng0[2])];
126+ for (jj = 0; jj < 3; jj++) {
127+ mat1[j0][j1][j2][jj] = mat0[ib][modulo(i0 + j0 - 1, ng0[0])]
128+ [modulo(i1 + j1 - 1, ng0[1])]
129+ [modulo(i2 + j2 - 1, ng0[2])][jj];
130+ }
131+ }/*for (j2 = 0; j2 < 4; j2++)*/
132+ }/*for (j1 = 0; j1 < 4; j1++)*/
133+ }/*for (i2 = 0; i2 < ng0[2]; i2++)*/
134+ for (j0 = 0; j0 < interpol; j0++) {
135+ kumo_coef(j0, &coef[0]);
136+ for (j1 = 0; j1 < 4; j1++) {
137+ for (j2 = 0; j2 < 4; j2++) {
138+ eig2[j1][j2] = 0.0;
139+ for (jj = 0; jj < 3; jj++) mat2[j1][j2][jj] = 0.0;
140+ for (ii = 0; ii < 4; ii++) {
141+ eig2[j1][j2] += coef[ii] * eig1[ii][j1][j2];
142+ for (jj = 0; jj < 3; jj++)
143+ mat2[j1][j2][jj] += coef[ii] * mat1[ii][j1][j2][jj];
144+ }/*for (ii = 0; ii < 4; ii++)*/
145+ }/*for (j2 = 0; j2 < 4; j2++)*/
146+ }/*for (j1 = 0; j1 < 4; j1++)*/
147+ for (j1 = 0; j1 < interpol; j1++) {
148+ kumo_coef(j1, &coef[0]);
149+ for (j2 = 0; j2 < 4; j2++) {
150+ eig3[j2] = 0.0;
151+ for (jj = 0; jj < 3; jj++) mat3[j2][jj] = 0.0;
152+ for (ii = 0; ii < 4; ii++) {
153+ eig3[j2] += coef[ii] * eig2[ii][j2];
154+ for (jj = 0; jj < 3; jj++)
155+ mat3[j2][jj] += coef[ii] * mat2[ii][j2][jj];
156+ }/*for (ii = 0; ii < 4; ii++)*/
157+ }/*for (j2 = 0; j2 < 4; j2++)*/
158+ for (j2 = 0; j2 < interpol; j2++) {
159+ kumo_coef(j2, &coef[0]);
160+ eig[ib][i0*interpol + j0]
161+ [i1*interpol + j1]
162+ [i2*interpol + j2] = 0.0;
163+ for (jj = 0; jj < 3; jj++)
164+ mat[ib][i0*interpol + j0]
165+ [i1*interpol + j1]
166+ [i2*interpol + j2][jj] = 0.0;
167+ for (ii = 0; ii < 4; ii++) {
168+ eig[ib][i0*interpol + j0]
169+ [i1*interpol + j1]
170+ [i2*interpol + j2] += coef[ii] * eig3[ii];
171+ for (jj = 0; jj < 3; jj++)
172+ mat[ib][i0*interpol + j0]
173+ [i1*interpol + j1]
174+ [i2*interpol + j2][jj] += coef[ii] * mat3[ii][jj];
175+ }/*for (ii = 0; ii < 4; ii++)*/
176+ }/*for (j2 = 0; j2 < interpol; j2++)*/
177+ }/*for (j1 = 0; j1 < interpol; j1++)*/
178+ }/*for (j0 = 0; j0 < interpol; j0++)*/
179+ }/*for (i2 = 0; i2 < ng0[2]; i2++)*/
180+ }/*for (i1 = 0; i1 < ng0[1]; i1++)*/
181+ }/*for (i0 = 0; i0 < ng0[0]; i0++)*/
182+ }/*for (ib = 0; ib < nb; ib++)*/
183+ }/*End of parallel region*/
184+ /*
185+ Fermi velocity
186+ */
187+#pragma omp parallel default(none) \
188+ shared(nb,ng,eig,vf,avec) \
189+ private (ib,i0,i1,i2,ii)
190+ {
191+ int i0p, i0m, i1p, i1m, i2p, i2m;
192+ GLfloat de[3];
193+
194+ for (ib = 0; ib < nb; ib++) {
195+ for (i0 = 0; i0 < ng[0]; i0++) {
196+ i0p = modulo(i0 + 1, ng[0]);
197+ i0m = modulo(i0 - 1, ng[0]);
198+ for (i1 = 0; i1 < ng[1]; i1++) {
199+ i1p= modulo(i1 + 1, ng[1]);
200+ i1m = modulo(i1 - 1, ng[1]);
201+ for (i2 = 0; i2 < ng[2]; i2++) {
202+ i2p = modulo(i2 + 1, ng[2]);
203+ i2m = modulo(i2 - 1, ng[2]);
204+
205+ de[0] = eig[ib][i0p][i1][i2] - eig[ib][i0m][i1][i2];
206+ de[1] = eig[ib][i0][i1p][i2] - eig[ib][i0][i1m][i2];
207+ de[2] = eig[ib][i0][i1][i2p] - eig[ib][i0][i1][i2m];
208+ for (ii = 0; ii < 3; ii++)de[ii] *= 0.5f * (GLfloat)ng[ii];
209+ for (ii = 0; ii < 3; ii++) vf[ib][i0][i1][i2][ii] =
210+ avec[0][ii] * de[0] + avec[1][ii] * de[1] + avec[2][ii] * de[2];
211+
212+ }/*for (i2 = 0; i2 < ng[2]; i2++)*/
213+ }/*for (i1 = 0; i1 < ng[1]; i1++)*/
214+ }/*for (i0 = 0; i0 < ng[0]; i0++)*/
215+ }/*for (ib = 0; ib < nb; ib++)*/
216+ }/*End of parallel region*/
217+}/*void interpol_energy() */
--- /dev/null
+++ b/Android1.NativeActivity/kumo.hpp
@@ -0,0 +1,25 @@
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+
25+void interpol_energy();
--- a/Android1.NativeActivity/main.cpp
+++ b/Android1.NativeActivity/main.cpp
@@ -14,6 +14,13 @@
1414 * 本ライセンスを参照してください。
1515 *
1616 */
17+#include "read_file.hpp"
18+#include "kumo.hpp"
19+#include "initialize.hpp"
20+#include "bz_lines.hpp"
21+#include "section.hpp"
22+#include "menu.hpp"
23+#include "draw.hpp"
1724 #include <EGL/egl.h>
1825 #include <GLES/gl.h>
1926
@@ -26,9 +33,141 @@
2633
2734 #define LOGI(...) ((void)__android_log_print(ANDROID_LOG_INFO, "AndroidProject1.NativeActivity", __VA_ARGS__))
2835 #define LOGW(...) ((void)__android_log_print(ANDROID_LOG_WARN, "AndroidProject1.NativeActivity", __VA_ARGS__))
29-GLfloat rot[3][3] = { { 1.0, 0.0, 0.0 },
30- { 0.0, 1.0, 0.0 },
36+/*
37+ Input variables
38+*/
39+int ng0[3]; //!< @f$k@f$-point grid in the input file
40+int shiftk[3]; //!< Wherether @f$k@f$-grid is shifted or not
41+int nb; //!< The number of Bands
42+GLfloat avec[3][3]; //!< Direct lattice vector
43+GLfloat bvec[3][3]; //!< Reciprocal lattice vector
44+GLfloat**** eig0; //!< Eigenvalues @f$\varepsilon_{n k}@f$[::nb][::ng0[0]][::ng0[1]][::ng0[2]]
45+GLfloat***** mat0; //!< Matrix element [::nb][::ng0[0]][::ng0[1]][::ng0[2]][3]
46+/*
47+ Interpolation
48+*/
49+int ng[3]; //!< @b Interpolated @f$k@f$-grids
50+GLfloat**** eig; //!< Eigenvalues @f$\varepsilon_{n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]]
51+GLfloat***** mat; //!< Matrix element @f$\delta_{n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]][3]
52+GLfloat***** vf; //!< Matrix element @f$\{\bf v}_{{\rm f} n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]][3]
53+int interpol = 1; //!< Ratio of interpolation
54+/*
55+ Switch for some modes
56+*/
57+int blackback = 1; //!< Switch for black background
58+int color_scale = 1; //!< Switch for full color scale mode
59+int fbz = 1; //!< Switch for 1st Brillouin zone mode
60+int nodeline = 0; //!< Switch for node lines
61+int lcolorbar = 1; //!< Switch for colorbar
62+int lstereo = 1; //!< Switch for the stereogram
63+int lmouse = 1; //!< Switch for the mouse function
64+int lsection = 0; //!< Switch for the 2D Fermi lines
65+int lequator = 0; //!< Switch for equator
66+/*
67+ Variables for Brillouin zone boundaries
68+*/
69+int nbzl; //!< The number of Lines of 1st Brillouin zone
70+GLfloat bzl[676][2][3]; //!< Lines of 1st BZ [nbzl(max:26*26=676)][2][3]
71+GLfloat bragg[26][3]; //!< Bragg plane vectors
72+GLfloat brnrm[26]; //!< Norms of Bragg plane vectors
73+GLfloat brnrm_min; //!< Minimum scale of the reciplocal space
74+int nbragg; //!< Number of Bragg plane og 1st BZ
75+/*
76+ Variables for patchs
77+*/
78+int* ntri; //!< The number of triangle patch [::nb]
79+int** ntri_th; //!< The number of triangle patch in each thread [::nb]
80+int* draw_band; //!< Switch for drawn bands [::nb]
81+GLfloat**** nmlp; //!< Normal vector of patchs [::nb][::ntri][3][3]
82+GLfloat**** kvp; //!< @f$k@f$-vectors of points [::nb][::ntri][3][3]
83+GLfloat***** arw;
84+GLfloat** nmlp_rot; //!< Normal vector of patchs [::nb][::ntri*3*3]
85+GLfloat** kvp_rot; //!< @f$k@f$-vectors of points [::nb][::ntri*3*3]
86+GLfloat** arw_rot;
87+GLfloat**** matp; //!< Matrix elements of points [::nb][::ntri][3][3]
88+GLfloat** clr; //!< Colors of points [::nb][::ntri*3*4]
89+int itet = 0; //!< Counter for tetrahedron
90+GLfloat side = 1.0; //!< Which side is lighted
91+GLfloat patch_max; //!< Max value across patch
92+GLfloat patch_min; //!< Max value across patch
93+/*
94+ Variables for nodeline
95+*/
96+int* nnl; //!< The number of nodeline
97+GLfloat**** kvnl; //!< @f$k@f$-vector of nodeline [::nb][::nnl][2][3]
98+GLfloat** kvnl_rot; //!< @f$k@f$-vector of nodeline [::nb][::nnl*2*3]
99+/*
100+ 2D Fermi line
101+*/
102+GLfloat secvec[3]; //!< @f$k@f$-vector to define section
103+GLfloat secvec_fr[3]; //!< @f$k@f$-vector to define section
104+GLfloat secscale; //!< 0.0 (across @f$\Gamma@f$) or 1.0
105+GLfloat axis2d[2][3]; //!< @f$k@f$-vector to define section
106+int* n2d; //!< Number of line segment
107+GLfloat** kv2d; //!< @f$k@f$-vector for 2D plot [::nb][::n2d*2*3]
108+GLfloat** clr2d; //!< Matrix element for 2D plot [::nb][::n2d*2*4]
109+int nbzl2d; //!< The number of Lines of 1st Brillouin zone
110+GLfloat bzl2d[26][3]; //!< Lines of 1st BZ [::nbzl2d (max:26)][3]
111+GLfloat bzl2d_proj[26][3]; //!< Lines of 1st BZ [::nbzl2d (max:26)][3], projected into 2D plane
112+/*
113+ Equator
114+*/
115+GLfloat eqvec[3]; //!< @f$k@f$-vector for equator
116+GLfloat eqvec_fr[3]; //!< @f$k@f$-vector for equator
117+int* nequator; //!< The number of equator
118+GLfloat**** kveq; //!< @f$k@f$-vector of equator [::nb][::nequator][2][3]
119+GLfloat** kveq_rot; //!< @f$k@f$-vector of equator [::nb][::nequator*2*3]
120+/*
121+ Variables for mouse & cursorkey
122+*/
123+GLfloat sx; //!< Scale of mouse movement
124+GLfloat sy; //!< Scale of mouse movement
125+int cx; //!< Starting point of drug
126+int cy; //!< Starting point of drug
127+GLfloat scl = 1.0; //!< Initial scale
128+GLfloat trans[3] = { 0.0, 0.0, 0.0 }; //!< Translation
129+GLfloat rot[3][3] = { { 1.0, 0.0, 0.0 },
130+ { 0.0, 1.0, 0.0 },
31131 { 0.0, 0.0, 1.0 } }; //!< Rotation matrix
132+GLfloat thetax = 0.0; //!< Rotation angle
133+GLfloat thetay = 0.0; //!< Rotation angle
134+GLfloat thetaz = 0.0; //!< Rotation angle
135+GLfloat linewidth = 3.0; //!< BZ/nodal-line/Fermiline width
136+/*
137+ Colors
138+*/
139+GLfloat black[4] = { 0.0, 0.0, 0.0, 1.0 }; //!< Black color code
140+GLfloat gray[4] = { 0.5f, 0.5f, 0.5f, 1.0 }; //!< Gray color code
141+GLfloat wgray[4] = { 0.9f, 0.9f, 0.9f, 1.0 }; //!< Gray color code
142+GLfloat bgray[4] = { 0.1f, 0.1f, 0.1f, 1.0 }; //!< Gray color code
143+GLfloat white[4] = { 1.0, 1.0, 1.0, 1.0 }; //!< White color code
144+GLfloat cyan[4] = { 0.0, 1.0, 1.0, 1.0 }; //!< Cyan color code
145+GLfloat magenta[4] = { 1.0, 0.0, 1.0, 1.0 }; //!< Magenta color code
146+GLfloat yellow[4] = { 1.0, 1.0, 0.0, 1.0 }; //!< Yellow color code
147+GLfloat red[4] = { 1.0, 0.0, 0.0, 1.0 }; //!< Red color code
148+GLfloat green[4] = { 0.0, 1.0, 0.0, 1.0 }; //!< Green color code
149+GLfloat blue[4] = { 0.0, 0.0, 1.0, 1.0 }; //!< Blue color code
150+/*
151+ Others
152+*/
153+int query; //!< Query switch
154+int corner[6][4]; //!< Corners of tetrahedron
155+GLfloat EF = 0.0; //!< Fermi energy
156+enum
157+{
158+ MOUSE_SCROLL_UP = 3, //!< Mouse wheel up
159+ MOUSE_SCROLL_DOWN = 4 //!< Mouse wheel down
160+};
161+int nthreads;//!< Number of OpenMP threads
162+int refresh_interpol = 0;
163+int refresh_patch = 1;
164+int refresh_color = 1;
165+int refresh_nodeline = 1;
166+int refresh_equator = 1;
167+int refresh_section = 1;
168+int skip_minmax = 0;
169+int windowx = 1100;
170+int windowy = 850;
32171 /**
33172 * 保存状態のデータです。
34173 */
@@ -110,7 +249,7 @@ static int engine_init_display(struct engine* engine) {
110249 // GL の状態を初期化します。
111250 glDisable(GL_DITHER);
112251 glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_FASTEST);
113- glClearColor(1.0f, 0.41f, 0.71f, 1.0f);
252+ glClearColor(0.0f, 0.0f, 0.0f, 1.0f);
114253 //glEnable(GL_CULL_FACE);
115254 glShadeModel(GL_SMOOTH);
116255 glEnable(GL_DEPTH_TEST);
@@ -135,121 +274,53 @@ static int engine_init_display(struct engine* engine) {
135274 * ディスプレイ内の現在のフレームのみ。
136275 */
137276 static void engine_draw_frame(struct engine* engine) {
138- int ntri = 8, itri, i, j;
139- GLfloat kvp_rot[72];
140- GLfloat kvp[8][3][3] = {
141- {{ 1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}},
142- {{ 1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, -1.0}},
143- {{ 1.0, 0.0, 0.0}, {0.0, -1.0, 0.0}, {0.0, 0.0, 1.0}},
144- {{ 1.0, 0.0, 0.0}, {0.0, -1.0, 0.0}, {0.0, 0.0, -1.0}},
145- {{-1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}},
146- {{-1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, -1.0}},
147- {{-1.0, 0.0, 0.0}, {0.0, -1.0, 0.0}, {0.0, 0.0, 1.0}},
148- {{-1.0, 0.0, 0.0}, {0.0, -1.0, 0.0}, {0.0, 0.0, -1.0}}
149- };
150- GLfloat nmlp_rot[72];
151- GLfloat nmlp[8][3][3] = {
152- {{ 1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}},
153- {{ 1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, -1.0}},
154- {{ 1.0, 0.0, 0.0}, {0.0, -1.0, 0.0}, {0.0, 0.0, 1.0}},
155- {{ 1.0, 0.0, 0.0}, {0.0, -1.0, 0.0}, {0.0, 0.0, -1.0}},
156- {{-1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}},
157- {{-1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, -1.0}},
158- {{-1.0, 0.0, 0.0}, {0.0, -1.0, 0.0}, {0.0, 0.0, 1.0}},
159- {{-1.0, 0.0, 0.0}, {0.0, -1.0, 0.0}, {0.0, 0.0, -1.0}}
160- };
161- GLfloat clr[8][3][4] = {
162- {{1.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0}},
163- {{1.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0}, {1.0, 1.0, 0.0, 0.0}},
164- {{1.0, 0.0, 0.0, 0.0}, {1.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 1.0, 0.0}},
165- {{1.0, 0.0, 0.0, 0.0}, {1.0, 0.0, 1.0, 0.0}, {1.0, 1.0, 0.0, 0.0}},
166- {{0.0, 1.0, 1.0, 0.0}, {0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0}},
167- {{0.0, 1.0, 1.0, 0.0}, {0.0, 1.0, 0.0, 0.0}, {1.0, 1.0, 0.0, 0.0}},
168- {{0.0, 1.0, 1.0, 0.0}, {1.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 1.0, 0.0}},
169- {{0.0, 1.0, 1.0, 0.0}, {1.0, 0.0, 1.0, 0.0}, {1.0, 1.0, 0.0, 0.0}}
170- };
171- GLfloat pos[] = { 1.0f, 1.0f, 1.0f, 0.0f };
172- GLfloat amb[] = { 0.2f, 0.2f, 0.2f, 0.0f };
173- GLfloat dx, dy, a, rot0[3][3], rot1[3][3], ax, ay, sx = 1.0/engine->width, sy = 1.0/engine->height;
277+ int i, j;
278+ GLfloat dx, dy, a, rot0[3][3], rot1[3][3], ax, ay, sx = 1.0/engine->width, sy = 1.0/engine->height;
174279
175- /*
176- Translation of mousepointer from starting point
177- */
178- dx = (engine->state.x - engine->state.x0) * sx;
179- dy = (engine->state.y - engine->state.y0) * sy;
180- /*
181- Distanse from starting point
182- */
183- a = sqrtf(dx * dx + dy * dy);
184- /**/
185- if (a != 0.0) {
186- /*
187- Compute rotational matrix from translation of mousepointer
188- */
189- ax = -dy;
190- ay = dx;
191- /**/
192- a = a * 10.0f;
193- /**/
194- rot0[0][0] = (ax * ax + ay * ay * cosf(a)) / (ax * ax + ay * ay);
195- rot0[0][1] = ax * ay * (cosf(a) - 1.0f) / (ax * ax + ay * ay);
196- rot0[0][2] = ay * sinf(a) / sqrtf(ax * ax + ay * ay);
197- rot0[1][0] = ax * ay * (cosf(a) - 1.0f) / (ax * ax + ay * ay);
198- rot0[1][1] = (ax * ax * cosf(a) + ay * ay) / (ax * ax + ay * ay);
199- rot0[1][2] = ax * sinf(a) / sqrtf(ax * ax + ay * ay);
200- rot0[2][0] = -ay * sinf(a) / sqrtf(ax * ax + ay * ay);
201- rot0[2][1] = -ax * sinf(a) / sqrtf(ax * ax + ay * ay);
202- rot0[2][2] = cosf(a);
203- /**/
204- for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) rot1[i][j] = rot[i][j];
205- /**/
206- for (i = 0; i < 3; i++) {
207- for (j = 0; j < 3; j++) {
208- rot[i][j] = rot0[i][0] * rot1[0][j]
209- + rot0[i][1] * rot1[1][j]
210- + rot0[i][2] * rot1[2][j];
211- }
212- }
280+ /*
281+ Translation of mousepointer from starting point
282+ */
283+ dx = (engine->state.x - engine->state.x0) * sx;
284+ dy = (engine->state.y - engine->state.y0) * sy;
285+ /*
286+ Distanse from starting point
287+ */
288+ a = sqrtf(dx * dx + dy * dy);
289+ /**/
290+ if (a != 0.0) {
291+ /*
292+ Compute rotational matrix from translation of mousepointer
293+ */
294+ ax = -dy;
295+ ay = dx;
296+ /**/
297+ a = a * 10.0f;
298+ /**/
299+ rot0[0][0] = (ax * ax + ay * ay * cosf(a)) / (ax * ax + ay * ay);
300+ rot0[0][1] = ax * ay * (cosf(a) - 1.0f) / (ax * ax + ay * ay);
301+ rot0[0][2] = ay * sinf(a) / sqrtf(ax * ax + ay * ay);
302+ rot0[1][0] = ax * ay * (cosf(a) - 1.0f) / (ax * ax + ay * ay);
303+ rot0[1][1] = (ax * ax * cosf(a) + ay * ay) / (ax * ax + ay * ay);
304+ rot0[1][2] = ax * sinf(a) / sqrtf(ax * ax + ay * ay);
305+ rot0[2][0] = -ay * sinf(a) / sqrtf(ax * ax + ay * ay);
306+ rot0[2][1] = -ax * sinf(a) / sqrtf(ax * ax + ay * ay);
307+ rot0[2][2] = cosf(a);
308+ /**/
309+ for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) rot1[i][j] = rot[i][j];
310+ /**/
311+ for (i = 0; i < 3; i++) {
312+ for (j = 0; j < 3; j++) {
313+ rot[i][j] = rot0[i][0] * rot1[0][j]
314+ + rot0[i][1] * rot1[1][j]
315+ + rot0[i][2] * rot1[2][j];
316+ }
317+ }
213318 }
214-
215- for (itri = 0; itri < ntri; ++itri) {
216- for (i = 0; i < 3; ++i) {
217- for (j = 0; j < 3; ++j) {
218- kvp_rot[j + 3 * i + 9 * itri]
219- = rot[j][0] * kvp[itri][i][0]
220- + rot[j][1] * kvp[itri][i][1]
221- + rot[j][2] * kvp[itri][i][2];
222- nmlp_rot[j + 3 * i + 9 * itri]
223- = rot[j][0] * nmlp[itri][i][0]
224- + rot[j][1] * nmlp[itri][i][1]
225- + rot[j][2] * nmlp[itri][i][2];
226- }
227- }/*for (i = 0; i < 3; ++i)*/
228- }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
229-
230319 if (engine->display == NULL) {
231320 // ディスプレイがありません。
232321 return;
233322 }
234-
235- glClearColor(0.0, 0.0, 0.0, 1);
236- glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
237-
238- glLoadIdentity();
239- glTranslatef(0.0, 0.0, -5.0);
240- glLightfv(GL_LIGHT0, GL_POSITION, pos);
241- glLightfv(GL_LIGHT1, GL_AMBIENT, amb);
242- glScalef(2.0, 2.0, 2.0);
243-
244- glEnableClientState(GL_NORMAL_ARRAY);
245- glEnableClientState(GL_COLOR_ARRAY);
246- glVertexPointer(3, GL_FLOAT, 0, kvp_rot);
247- glNormalPointer(GL_FLOAT, 0, nmlp_rot);
248- glColorPointer(4, GL_FLOAT, 0, clr);
249- glDrawArrays(GL_TRIANGLES, 0, ntri * 3);
250- glDisableClientState(GL_NORMAL_ARRAY);
251- glDisableClientState(GL_COLOR_ARRAY);
252-
323+ draw();
253324 eglSwapBuffers(engine->display, engine->surface);
254325 }
255326 /**
@@ -330,6 +401,7 @@ static void engine_handle_cmd(struct android_app* app, int32_t cmd) {
330401 void android_main(struct android_app* state) {
331402 struct engine engine;
332403
404+ nthreads = 1;
333405 memset(&engine, 0, sizeof(engine));
334406 state->userData = &engine;
335407 state->onAppCmd = engine_handle_cmd;
@@ -341,8 +413,19 @@ void android_main(struct android_app* state) {
341413 engine.state = *(struct saved_state*)state->savedState;
342414 }
343415
344- // ループはスタッフによる開始を待っています。
345-
416+ color_scale = read_file();
417+ interpol_energy();
418+ init_corner();
419+ bragg_vector();
420+ /*
421+ Brillouin zone
422+ */
423+ bz_lines();
424+ calc_2dbz();
425+ /**/
426+ max_and_min_bz();
427+ /**/
428+ compute_patch_segment();
346429 while (1) {
347430 // 保留中のすべてのイベントを読み取ります。
348431 int ident;
--- /dev/null
+++ b/Android1.NativeActivity/menu.cpp
@@ -0,0 +1,79 @@
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 Create & modify right-click menu. And operate their function.
26+*/
27+#include "free_patch.hpp"
28+#include "fermi_patch.hpp"
29+#include "calc_nodeline.hpp"
30+#include "kumo.hpp"
31+#include "section.hpp"
32+#include "equator.hpp"
33+#include "draw.hpp"
34+#include "variable.hpp"
35+
36+void compute_patch_segment() {
37+ if (refresh_interpol == 1){
38+ interpol_energy();
39+ refresh_patch = 1;
40+ refresh_interpol = 0;
41+ }
42+ if (refresh_patch == 1) {
43+ query = 1; fermi_patch();
44+ query = 0; fermi_patch();
45+ refresh_color = 1;
46+ refresh_section = 1;
47+ refresh_equator = 1;
48+ refresh_nodeline = 1;
49+ refresh_patch = 0;
50+ }
51+ if (refresh_color == 1) {
52+ if (skip_minmax == 1) skip_minmax = 0;
53+ else max_and_min();
54+ paint();
55+ refresh_section = 1;
56+ refresh_color = 0;
57+ }
58+ if (refresh_nodeline == 1) {
59+ query = 1; calc_nodeline();
60+ query = 0; calc_nodeline();
61+ refresh_nodeline = 0;
62+ }
63+ if (refresh_section == 1) {
64+ calc_2dbz();
65+ query = 1; calc_section();
66+ query = 0; calc_section();
67+ refresh_section = 0;
68+ }
69+ if (refresh_equator == 1) {
70+ query = 1; equator();
71+ query = 0; equator();
72+ refresh_equator = 0;
73+ }
74+}
75+
76+void refresh_patch_segment() {
77+ free_patch();
78+ compute_patch_segment();
79+}
--- /dev/null
+++ b/Android1.NativeActivity/menu.hpp
@@ -0,0 +1,30 @@
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+#if ! defined(FERMISURFER_MENU_H)
25+#define FERMISURFER_MENU_H
26+
27+void compute_patch_segment();
28+void refresh_patch_segment();
29+
30+#endif // FERMISURFER_MENU_H
--- /dev/null
+++ b/Android1.NativeActivity/read_file.cpp
@@ -0,0 +1,196 @@
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 Read .frmsf file
26+*/
27+#include <stdlib.h>
28+#include <stdio.h>
29+#include <math.h>
30+#include <ctype.h>
31+#include <string.h>
32+#include <GLES/gl.h>
33+
34+#include "basic_math.hpp"
35+#include "variable.hpp"
36+/*
37+ Allocation of Kohn-Sham energies $ matrix elements
38+*/
39+static void allocate_griddata()
40+{
41+ int i, j, ib, i0, i1, i2;
42+
43+ for (i = 0; i < 3; i++) ng[i] = ng0[i];
44+
45+ ntri = new int[nb];
46+ ntri_th = new int* [nb];
47+ for (ib = 0; ib < nb; ib++) ntri_th[ib] = new int[nthreads];
48+ nnl = new int[nb];
49+ n2d = new int[nb];
50+ nequator = new int[nb];
51+ draw_band = new int[nb];
52+ for (ib = 0; ib < nb; ib++) draw_band[ib] = 1;
53+
54+ scl /= sqrtf(bvec[0][0] * bvec[0][0] + bvec[0][1] * bvec[0][1] + bvec[0][2] * bvec[0][2]);
55+ linewidth /= scl;
56+ /*
57+ Direct lattice vector
58+ */
59+ for (i = 0; i < 3; ++i) {
60+ for (j = 0; j < 3; ++j) avec[i][j] = 0.0f;
61+ avec[i][i] = 1.0f;
62+ solve3(bvec, avec[i]);
63+ }/*for (i = 0; i < 3; ++i)*/
64+ for (i = 0; i < 3; ++i) {
65+ secvec[i] = bvec[2][i];
66+ eqvec[i] = bvec[2][i];
67+ eqvec_fr[i] = 0.0;
68+ secvec_fr[i] = 0.0;
69+ }
70+ eqvec_fr[2] = 1.0;
71+ secvec_fr[2] = 1.0;
72+ secscale = 0.0;
73+
74+ eig0 = new GLfloat * **[nb];
75+ eig = new GLfloat * **[nb];
76+ mat0 = new GLfloat * ***[nb];
77+ mat = new GLfloat * ***[nb];
78+ vf = new GLfloat * ***[nb];
79+ for (ib = 0; ib < nb; ib++) {
80+ eig0[ib] = new GLfloat * *[ng0[0]];
81+ eig[ib] = new GLfloat * *[ng0[0]];
82+ mat0[ib] = new GLfloat * **[ng0[0]];
83+ mat[ib] = new GLfloat * **[ng0[0]];
84+ vf[ib] = new GLfloat * **[ng0[0]];
85+ for (i0 = 0; i0 < ng0[0]; i0++) {
86+ eig0[ib][i0] = new GLfloat * [ng0[1]];
87+ eig[ib][i0] = new GLfloat * [ng0[1]];
88+ mat0[ib][i0] = new GLfloat * *[ng0[1]];
89+ mat[ib][i0] = new GLfloat * *[ng0[1]];
90+ vf[ib][i0] = new GLfloat * *[ng0[1]];
91+ for (i1 = 0; i1 < ng0[1]; i1++) {
92+ eig0[ib][i0][i1] = new GLfloat[ng0[2]];
93+ eig[ib][i0][i1] = new GLfloat[ng0[2]];
94+ mat0[ib][i0][i1] = new GLfloat * [ng0[2]];
95+ mat[ib][i0][i1] = new GLfloat * [ng0[2]];
96+ vf[ib][i0][i1] = new GLfloat * [ng0[2]];
97+ for (i2 = 0; i2 < ng0[2]; ++i2) {
98+ mat0[ib][i0][i1][i2] = new GLfloat[3];
99+ mat[ib][i0][i1][i2] = new GLfloat[3];
100+ vf[ib][i0][i1][i2] = new GLfloat[3];
101+ }
102+ }
103+ }
104+ }
105+}
106+/**
107+ @brief Input from Fermi surface file
108+*/
109+int read_file()
110+{
111+ int ib, i, i0, i1, i2, ii0, ii1, ii2, ierr, iaxis;
112+ FILE *fp;
113+ char* ctemp1;
114+ char ctemp2[256];
115+ int lshift; //!< Switch for shifted Brillouin zone
116+ /*
117+ Open input file.
118+ */
119+ fp = fopen("/sdcard/Download/vfermi.frmsf", "r");
120+ /*
121+ k-point grid
122+ */
123+ ctemp1 = fgets(ctemp2, 256, fp);
124+ ierr = sscanf(ctemp2, "%d%d%d", &ng0[0], &ng0[1], &ng0[2]);
125+ /*
126+ Shift of k-point grid
127+ */
128+ ierr = fscanf(fp, "%d", &lshift);
129+ if (lshift == 0) {
130+ for (i = 0; i < 3; i++) shiftk[i] = (ng0[i] + 1) % 2;
131+ }
132+ else if (lshift == 1) {
133+ for (i = 0; i < 3; i++) shiftk[i] = 0;
134+ }
135+ else if (lshift == 2) {
136+ for (i = 0; i < 3; i++) shiftk[i] = 1;
137+ }
138+ else {
139+ exit(0);
140+ }
141+ /*
142+ # of bands
143+ */
144+ ierr = fscanf(fp, "%d", &nb);
145+ /*
146+ Reciplocal lattice vectors
147+ */
148+ for (i = 0; i < 3; ++i) {
149+ ierr = fscanf(fp, "%e%e%e", &bvec[i][0], &bvec[i][1], &bvec[i][2]);
150+ }/*for (i = 0; i < 3; ++i)*/
151+ allocate_griddata();
152+ /*
153+ Kohn-Sham energies
154+ */
155+ for (ib = 0; ib < nb; ++ib) {
156+ for (i0 = 0; i0 < ng0[0]; ++i0) {
157+ if (lshift != 0) ii0 = i0;
158+ else ii0 = modulo(i0 + (ng0[0] + 1) / 2, ng0[0]);
159+ for (i1 = 0; i1 < ng0[1]; ++i1) {
160+ if (lshift != 0) ii1 = i1;
161+ else ii1 = modulo(i1 + (ng0[1] + 1) / 2, ng0[1]);
162+ for (i2 = 0; i2 < ng0[2]; ++i2) {
163+ if (lshift != 0) ii2 = i2;
164+ else ii2 = modulo(i2 + (ng0[2] + 1) / 2, ng0[2]);
165+ ierr = fscanf(fp, "%e", &eig0[ib][ii0][ii1][ii2]);
166+ }
167+ }
168+ }
169+ }
170+ /*
171+ Matrix elements
172+ */
173+ for (iaxis = 0; iaxis < 3; iaxis++) {
174+ for (ib = 0; ib < nb; ++ib) {
175+ for (i0 = 0; i0 < ng0[0]; ++i0) {
176+ if (lshift != 0) ii0 = i0;
177+ else ii0 = modulo(i0 + (ng0[0] + 1) / 2, ng0[0]);
178+ for (i1 = 0; i1 < ng0[1]; ++i1) {
179+ if (lshift != 0) ii1 = i1;
180+ else ii1 = modulo(i1 + (ng0[1] + 1) / 2, ng0[1]);
181+ for (i2 = 0; i2 < ng0[2]; ++i2) {
182+ if (lshift != 0) ii2 = i2;
183+ else ii2 = modulo(i2 + (ng0[2] + 1) / 2, ng0[2]);
184+ ierr = fscanf(fp, "%e", &mat0[ib][ii0][ii1][ii2][iaxis]);
185+ if (ierr == EOF) {
186+ fclose(fp);
187+ return iaxis;
188+ }
189+ }/*for (i2 = 0; i2 < ng0[2]; ++i2)*/
190+ }/*for (i1 = 0; i1 < ng0[1]; ++i1)*/
191+ }/*for (i0 = 0; i0 < ng0[0]; ++i0)*/
192+ }/*for (ib = 0; ib < nb; ++ib)*/
193+ }
194+ fclose(fp);
195+ return 3;
196+} /* read_file */
--- /dev/null
+++ b/Android1.NativeActivity/read_file.hpp
@@ -0,0 +1,24 @@
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+int read_file();
--- /dev/null
+++ b/Android1.NativeActivity/section.cpp
@@ -0,0 +1,366 @@
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+/**
25+@file
26+@brief Functions for the 2D Fermi line
27+*/
28+#include <GLES/gl.h>
29+#include <stdio.h>
30+#include <stdlib.h>
31+#include <math.h>
32+#include "basic_math.hpp"
33+#include "variable.hpp"
34+/**
35+ @brief Project 3D \f$k\f$-vector into 2D plane.
36+
37+ Modify: ::kv2d
38+*/
39+static void proj_2d(
40+ GLfloat vec[3] //!< [in,out] Line ends to be projected
41+)
42+{
43+ int ii, kk;
44+ GLfloat vec0[3];
45+
46+ for (kk = 0; kk < 2; kk++) {
47+ vec0[kk] = 0.0;
48+ for (ii = 0; ii < 3; ii++)
49+ vec0[kk] += axis2d[kk][ii] * vec[ii];
50+ }/*for (kk = 0; kk < 2; kk++)*/
51+ vec0[2] = 0.0;
52+ for (kk = 0; kk < 3; kk++)vec[kk] = vec0[kk];
53+}/*proj_2d*/
54+/**
55+ @brief Set Projection axis for 2D plane
56+
57+ Modify : ::axis2d
58+*/
59+static void set2daxis() {
60+ int ii, jj;
61+ GLfloat snorm, norm, thr = 0.001f;
62+
63+ snorm = 0.0;
64+ for (ii = 0; ii < 3; ii++) snorm += secvec[ii] * secvec[ii];
65+ /*
66+ Define the first axis
67+ */
68+ for (ii = 0; ii < 3; ii++) {
69+ for (jj = 0; jj < 3; jj++) axis2d[0][jj] = 0.0;
70+ axis2d[0][ii] = 1.0;
71+ for (jj = 0; jj < 3; jj++) axis2d[0][jj] += - secvec[ii] * secvec[jj] / snorm;
72+
73+ norm = 0.0;
74+ for (jj = 0; jj < 3; jj++) norm += axis2d[0][jj] * axis2d[0][jj];
75+ norm = sqrtf(norm);
76+ if (norm > thr) {
77+ for (jj = 0; jj < 3; jj++) axis2d[0][jj] /= norm;
78+ break;
79+ }/*if (norm > 0.thr)*/
80+ }/*for (ii = 0; ii < 3; ii++)*/
81+ /*
82+ Define the second axis with outor product
83+ */
84+ axis2d[1][0] = secvec[1] * axis2d[0][2] - secvec[2] * axis2d[0][1];
85+ axis2d[1][1] = secvec[2] * axis2d[0][0] - secvec[0] * axis2d[0][2];
86+ axis2d[1][2] = secvec[0] * axis2d[0][1] - secvec[1] * axis2d[0][0];
87+ norm = 0.0;
88+ for (jj = 0; jj < 3; jj++) norm += axis2d[1][jj] * axis2d[1][jj];
89+ norm = sqrtf(norm);
90+ for (jj = 0; jj < 3; jj++) axis2d[1][jj] /= norm;
91+}/*static void set_2daxis*/
92+/**
93+ @brief Judge wheser this line is the edge of 1st BZ (or the premitive BZ)
94+*/
95+int bragg_vert2d(
96+ int jbr, //!< [in] Index of a Bragg plane
97+ int nbr, //!< [in]
98+ GLfloat vert[3], //!< [inout] start point of line
99+ GLfloat vert2[3] //!< [in] end point of line
100+)
101+{
102+ int kbr, i, lbr, nbr0;
103+ GLfloat bmat[3][3], rhs[3], prod, thr, det;
104+ /**/
105+ nbr0 = nbr;
106+ /**/
107+ for (kbr = nbr0; kbr < nbragg; ++kbr) {
108+ /**/
109+ /**/
110+ for (i = 0; i<3; ++i) bmat[0][i] = secvec[i];
111+ for (i = 0; i<3; ++i) bmat[1][i] = bragg[jbr][i];
112+ for (i = 0; i<3; ++i) bmat[2][i] = bragg[kbr][i];
113+ /**/
114+ rhs[0] = 0.0;
115+ for (i = 0; i < 3; ++i)rhs[0] += secvec[i] * secvec[i];
116+ rhs[1] = brnrm[jbr];
117+ rhs[2] = brnrm[kbr];
118+ thr = sqrtf(rhs[0] * rhs[1] * rhs[2]) * 0.001f;
119+ rhs[0] *= secscale;
120+ /*
121+ if Bragg planes do not cross, roop next kbr
122+ */
123+ det = solve3(bmat, rhs);
124+ if (fabsf(det) < thr) continue;
125+ /*
126+ if vert0 = vert1, roop next kbr
127+ */
128+ prod = (vert2[0] - rhs[0]) * (vert2[0] - rhs[0])
129+ + (vert2[1] - rhs[1]) * (vert2[1] - rhs[1])
130+ + (vert2[2] - rhs[2]) * (vert2[2] - rhs[2]);
131+ if (prod < thr) continue;
132+ /*
133+ is corner really in 1st BZ ?
134+ */
135+ i = 0;
136+ for (lbr = 0; lbr < nbragg; ++lbr) {
137+ prod = bragg[lbr][0] * rhs[0]
138+ + bragg[lbr][1] * rhs[1]
139+ + bragg[lbr][2] * rhs[2];
140+ /**/
141+ if (prod > brnrm[lbr] + thr) {
142+ i = 1;
143+ break;
144+ }
145+ }/*for (lbr = 0; lbr < nbragg; ++lbr)*/
146+ if (i == 1) {
147+ }
148+ else {
149+ for (i = 0; i<3; ++i) vert[i] = rhs[i];
150+ return kbr + 1;
151+ }
152+ }/*for (kbr = nbr0; kbr < nbragg; ++kbr)*/
153+ /*
154+ this line is not a BZ boundary
155+ */
156+ return 0;
157+}/* bragg_vert2d */
158+/**
159+ @brief Compute boundary of 2D BZ
160+
161+ Modify : ::nbzl2d, ::bzl2d_proj
162+*/
163+void calc_2dbz() {
164+ int jbr, nbr, i, j, lvert, ibzl;
165+ GLfloat vert[2][3], vec[26][2][3], prod, thr;
166+
167+ if (fbz != 1)return;
168+ /*
169+ Set Projection axis for 2D plane
170+ */
171+ set2daxis();
172+
173+ nbzl2d = 0;
174+
175+ for (jbr = 0; jbr < nbragg; ++jbr) {
176+ /**/
177+ for (i = 0; i < 3; ++i) vert[1][i] = 0.0;
178+ nbr = 0;
179+ lvert = bragg_vert2d(jbr, nbr, vert[0], vert[1]);
180+ if (lvert == 0) continue;
181+ nbr = lvert;
182+ /**/
183+ lvert = bragg_vert2d(jbr, nbr, vert[1], vert[0]);
184+ if (lvert == 0) continue;
185+ /**/
186+ for (i = 0; i < 2; ++i) for (j = 0; j < 3; ++j) vec[nbzl2d][i][j] = vert[i][j];
187+ nbzl2d = nbzl2d + 1;
188+ }/*for (jbr = 0; jbr < nbragg; ++jbr)*/
189+ /*
190+ Order bz lines
191+ */
192+ for (i = 0; i < 3; i++) bzl2d[0][i] = vec[0][0][i];
193+ for (i = 0; i < 3; i++) bzl2d[1][i] = vec[0][1][i];
194+ for (ibzl = 0; ibzl < nbzl2d; ibzl++) {
195+
196+ thr = 0.0f;
197+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) thr += bzl2d[j][i] * bzl2d[j][i];
198+ thr *= 0.001f;
199+
200+ prod = 0.0;
201+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++)
202+ prod += (bzl2d[j][i] - vec[ibzl][j][i]) * (bzl2d[j][i] - vec[ibzl][j][i]);
203+ if (prod < thr)
204+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) vec[ibzl][j][i] = 0.0;
205+
206+ prod = 0.0;
207+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++)
208+ prod += (bzl2d[1 - j][i] - vec[ibzl][j][i]) * (bzl2d[1 - j][i] - vec[ibzl][j][i]);
209+ if (prod < thr)
210+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) vec[ibzl][j][i] = 0.0;
211+ }/*for (ibzl = 1; ibzl < nbzl2d; ibzl++)*/
212+
213+ for (jbr = 1; jbr < nbzl2d - 1; jbr++) {
214+
215+ thr = 0.0f;
216+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) thr += bzl2d[jbr][i] * bzl2d[jbr][i];
217+ thr *= 0.001f;
218+
219+ prod = 0.0;
220+ for (ibzl = 0; ibzl < nbzl2d; ibzl++) for (i = 0; i < 3; i++)
221+ prod += vec[ibzl][0][i] * vec[ibzl][0][i];
222+ if (prod < thr) {
223+ nbzl2d = jbr + 1;
224+ break;
225+ }
226+
227+ for (ibzl = 1; ibzl < nbzl2d; ibzl++) {
228+ prod = (bzl2d[jbr][0] - vec[ibzl][0][0]) * (bzl2d[jbr][0] - vec[ibzl][0][0])
229+ + (bzl2d[jbr][1] - vec[ibzl][0][1]) * (bzl2d[jbr][1] - vec[ibzl][0][1])
230+ + (bzl2d[jbr][2] - vec[ibzl][0][2]) * (bzl2d[jbr][2] - vec[ibzl][0][2]);
231+ if (prod < thr) {
232+ for (i = 0; i < 3; i++) bzl2d[jbr + 1][i] = vec[ibzl][1][i];
233+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) vec[ibzl][j][i] = 0.0;
234+ }
235+
236+ prod = (bzl2d[jbr][0] - vec[ibzl][1][0]) * (bzl2d[jbr][0] - vec[ibzl][1][0])
237+ + (bzl2d[jbr][1] - vec[ibzl][1][1]) * (bzl2d[jbr][1] - vec[ibzl][1][1])
238+ + (bzl2d[jbr][2] - vec[ibzl][1][2]) * (bzl2d[jbr][2] - vec[ibzl][1][2]);
239+ if (prod < thr) {
240+ for (i = 0; i < 3; i++) bzl2d[jbr + 1][i] = vec[ibzl][0][i];
241+ for (j = 0; j < 2; j++) for (i = 0; i < 3; i++) vec[ibzl][j][i] = 0.0;
242+ }
243+ }/*for (ibzl = 1; ibzl < nbzl2d; ibzl++)*/
244+ }/*for (ibzl = 1; ibzl < nbzl; ibzl++)*/
245+ /*
246+ Project into 2D plane
247+ */
248+ for (ibzl = 0; ibzl < nbzl2d; ibzl++) {
249+ for (i = 0; i < 3; i++) bzl2d_proj[ibzl][i] = bzl2d[ibzl][i];
250+ proj_2d(bzl2d_proj[ibzl]);
251+ }/*for (ibzl = 0; ibzl < nbzl2d; ibzl++)*/
252+}/*calc_2dbz()*/
253+/**
254+ @brief Compute Fermi-line
255+
256+ Modify : ::ntri_th, ::n2d, ::clr2d, ::kv2d
257+
258+ If ::query = 1, this routine only compute the number of line segment and
259+ malloc variables.
260+
261+ If ::query = 0, actually compute lines and store them.
262+*/
263+void calc_section() {
264+ int i, ib, itri, ithread;
265+
266+ if (fbz != 1)return;
267+
268+#pragma omp parallel default(none) \
269+ shared(nb,n2d,clr,clr2d,kvp,kv2d,ntri,ntri_th,secvec,secscale,query) \
270+ private(ib,itri,i,ithread)
271+ {
272+ int j, n2d0, sw[3];
273+ GLfloat norm[3], a[3][3];
274+
275+ ithread = get_thread();
276+ for (ib = 0; ib < nb; ib++) {
277+ if (query == 1) n2d0 = 0;
278+ else n2d0 = ntri_th[ib][ithread];
279+#pragma omp for
280+ for (itri = 0; itri < ntri[ib]; ++itri) {
281+ /**/
282+ for (i = 0; i < 3; i++) {
283+ norm[i] = 0.0;
284+ for (j = 0; j < 3; j++) norm[i] += secvec[j] * (kvp[ib][itri][i][j] - secscale * secvec[j]);
285+ }
286+ eigsort(3, norm, sw);
287+ for (i = 0; i < 3; ++i) {
288+ for (j = 0; j < 3; ++j) {
289+ a[i][j] = (0.0f - norm[sw[j]]) / (norm[sw[i]] - norm[sw[j]]);
290+ }/*for (j = 0; j < 3; ++j)*/
291+ }/*for (i = 0; i < 3; ++i)*/
292+ /**/
293+ if ((norm[sw[0]] < 0.0 && 0.0 <= norm[sw[1]]) || (norm[sw[0]] <= 0.0 && 0.0 < norm[sw[1]])) {
294+ if (query == 0) {
295+ for (i = 0; i < 3; ++i) {
296+ kv2d[ib][i + 0 * 3 + 6 * n2d0]
297+ = kvp[ib][itri][sw[1]][i] * a[1][0] + kvp[ib][itri][sw[0]][i] * a[0][1];
298+ kv2d[ib][i + 1 * 3 + 6 * n2d0]
299+ = kvp[ib][itri][sw[2]][i] * a[2][0] + kvp[ib][itri][sw[0]][i] * a[0][2];
300+ }/*for (i = 0; i < 3; ++i)*/
301+ for (i = 0; i < 4; ++i) {
302+ clr2d[ib][i + 0 * 4 + 8 * n2d0]
303+ = clr[ib][i + 4 * sw[1] + 12 * itri] * a[1][0]
304+ + clr[ib][i + 4 * sw[0] + 12 * itri] * a[0][1];
305+ clr2d[ib][i + 1 * 4 + 8 * n2d0]
306+ = clr[ib][i + 4 * sw[2] + 12 * itri] * a[2][0]
307+ + clr[ib][i + 4 * sw[0] + 12 * itri] * a[0][2];
308+ }/*for (i = 0; i < 4; ++i)*/
309+ proj_2d(&kv2d[ib][0 * 3 + 6 * n2d0]);
310+ proj_2d(&kv2d[ib][1 * 3 + 6 * n2d0]);
311+ }/*if (query == 0)*/
312+ n2d0 += 1;
313+ }/*else if (norm[sw[0]] < 0.0 && 0.0 <= norm[sw[1]])*/
314+ else if ((norm[sw[1]] < 0.0 && 0.0 <= norm[sw[2]]) || (norm[sw[1]] <= 0.0 && 0.0 < norm[sw[2]])) {
315+ if (query == 0) {
316+ for (i = 0; i < 3; ++i) {
317+ kv2d[ib][i + 0 * 3 + 6 * n2d0]
318+ = kvp[ib][itri][sw[2]][i] * a[2][0] + kvp[ib][itri][sw[0]][i] * a[0][2];
319+ kv2d[ib][i + 1 * 3 + 6 * n2d0]
320+ = kvp[ib][itri][sw[2]][i] * a[2][1] + kvp[ib][itri][sw[1]][i] * a[1][2];
321+ }/*for (i = 0; i < 3; ++i)*/
322+ for (i = 0; i < 4; ++i) {
323+ clr2d[ib][i + 0 * 4 + 8 * n2d0]
324+ = clr[ib][i + 4 * sw[2] + 12 * itri] * a[2][0]
325+ + clr[ib][i + 4 * sw[0] + 12 * itri] * a[0][2];
326+ clr2d[ib][i + 1 * 4 + 8 * n2d0]
327+ = clr[ib][i + 4 * sw[2] + 12 * itri] * a[2][1]
328+ + clr[ib][i + 4 * sw[1] + 12 * itri] * a[1][2];
329+ }/*for (i = 0; i < 4; ++i)*/
330+ proj_2d(&kv2d[ib][0 * 3 + 6 * n2d0]);
331+ proj_2d(&kv2d[ib][1 * 3 + 6 * n2d0]);
332+ }/*if (query == 0)*/
333+ n2d0 += 1;
334+ }/*else if (norm[sw[1]] < 0.0 && 0.0 <= norm[sw[2]])*/
335+ }/*for (itri = 0; itri < ntri[ib]; ++itri)*/
336+ if(query == 1) ntri_th[ib][ithread] = n2d0;
337+ }/*for (ib = 0; ib < nb; ib++)*/
338+ }/* End of parallel region */
339+ /*
340+ Allocation of Fermi-lines
341+ */
342+ if (query == 1) {
343+ /*
344+ Sum Fermi-lines in all threads
345+ */
346+ for (ib = 0; ib < nb; ib++) {
347+ for (ithread = 1; ithread < nthreads; ithread++) {
348+ ntri_th[ib][ithread] += ntri_th[ib][ithread - 1];
349+ }
350+ n2d[ib] = ntri_th[ib][nthreads - 1];
351+ for (ithread = nthreads - 1; ithread > 0; ithread--) {
352+ ntri_th[ib][ithread] = ntri_th[ib][ithread - 1];
353+ }
354+ ntri_th[ib][0] = 0;
355+ }
356+ /*
357+ Allocation of Fermi-lines
358+ */
359+ kv2d = new GLfloat*[nb];
360+ clr2d = new GLfloat*[nb];
361+ for (ib = 0; ib < nb; ++ib) {
362+ kv2d[ib] = new GLfloat[6 * n2d[ib]];
363+ clr2d[ib] = new GLfloat[8 * n2d[ib]];
364+ }/*for (ib = 0; ib < nb; ++ib)*/
365+ }/*if (query == 0)*/
366+}/*void calc_nodeline()*/
--- /dev/null
+++ b/Android1.NativeActivity/section.hpp
@@ -0,0 +1,26 @@
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+
25+void calc_2dbz();
26+void calc_section();
--- /dev/null
+++ b/Android1.NativeActivity/variable.hpp
@@ -0,0 +1,172 @@
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+/**
25+@file
26+@brief Global variables
27+*/
28+#if ! defined(FERMISURFER_VARIABLE_H)
29+#define FERMISURFER_VARIABLE_H
30+
31+#include <GLES/gl.h>
32+
33+#include "menu.hpp"
34+
35+/*
36+ Input variables
37+*/
38+extern int ng0[3]; //!< @f$k@f$-point grid in the input file
39+extern int shiftk[3]; //!< Wherether @f$k@f$-grid is shifted or not
40+extern int nb; //!< The number of Bands
41+extern GLfloat avec[3][3]; //!< Direct lattice vector
42+extern GLfloat bvec[3][3]; //!< Reciprocal lattice vector
43+extern GLfloat ****eig0; //!< Eigenvalues @f$\varepsilon_{n k}@f$[::nb][::ng0[0]][::ng0[1]][::ng0[2]]
44+extern GLfloat *****mat0; //!< Matrix element [::nb][::ng0[0]][::ng0[1]][::ng0[2]][3]
45+/*
46+ Interpolation
47+*/
48+extern int ng[3]; //!< @b Interpolated @f$k@f$-grids
49+extern GLfloat ****eig; //!< Eigenvalues @f$\varepsilon_{n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]]
50+extern GLfloat *****mat; //!< Matrix element @f$\delta_{n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]][3]
51+extern GLfloat *****vf; //!< Matrix element @f$\{\bf v}_{{\rm f} n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]][3]
52+extern int interpol; //!< Ratio of interpolation
53+/*
54+ Switch for some modes
55+*/
56+extern int blackback; //!< Switch for black background
57+extern int color_scale; //!< Switch for full color scale mode
58+extern int fbz; //!< Switch for 1st Brillouin zone mode
59+extern int nodeline; //!< Switch for node lines
60+extern int lcolorbar; //!< Switch for colorbar
61+extern int lstereo; //!< Switch for the stereogram
62+extern int lmouse; //!< Switch for the mouse function
63+extern int lsection; //!< Switch for the 2D Fermi lines
64+extern int lequator; //!< Switch for equator
65+/*
66+ Variables for Brillouin zone boundaries
67+*/
68+extern int nbzl; //!< The number of Lines of 1st Brillouin zone
69+extern GLfloat bzl[676][2][3]; //!< Lines of 1st BZ [nbzl(max:26*26=676)][2][3]
70+extern GLfloat bragg[26][3]; //!< Bragg plane vectors
71+extern GLfloat brnrm[26]; //!< Norms of Bragg plane vectors
72+extern GLfloat brnrm_min; //!< Minimum scale of the reciplocal space
73+extern int nbragg; //!< Number of Bragg plane og 1st BZ
74+/*
75+ Variables for patchs
76+*/
77+extern int *ntri; //!< The number of triangle patch [::nb]
78+extern int **ntri_th; //!< The number of triangle patch in each thread [::nb]
79+extern int *draw_band; //!< Switch for drawn bands [::nb]
80+extern GLfloat ****nmlp; //!< Normal vector of patchs [::nb][::ntri][3][3]
81+extern GLfloat***** arw;
82+extern GLfloat ****kvp; //!< @f$k@f$-vectors of points [::nb][::ntri][3][3]
83+extern GLfloat **nmlp_rot; //!< Normal vector of patchs [::nb][::ntri*3*3]
84+extern GLfloat **kvp_rot; //!< @f$k@f$-vectors of points [::nb][::ntri*3*3]
85+extern GLfloat **arw_rot;
86+extern GLfloat ****matp; //!< Matrix elements of points [::nb][::ntri][3][3]
87+extern GLfloat **clr; //!< Colors of points [::nb][::ntri*3*4]
88+extern int itet; //!< Counter for tetrahedron
89+extern GLfloat side; //!< Which side is lighted
90+extern GLfloat patch_max; //!< Max value across patch
91+extern GLfloat patch_min; //!< Max value across patch
92+/*
93+ Variables for nodeline
94+*/
95+extern int *nnl; //!< The number of nodeline
96+extern GLfloat ****kvnl; //!< @f$k@f$-vector of nodeline [::nb][::nnl][2][3]
97+extern GLfloat **kvnl_rot; //!< @f$k@f$-vector of nodeline [::nb][::nnl*2*3]
98+/*
99+ 2D Fermi line
100+*/
101+extern GLfloat secvec[3]; //!< @f$k@f$-vector to define section
102+extern GLfloat secvec_fr[3]; //!< @f$k@f$-vector to define section
103+extern GLfloat secscale; //!< 0.0 (across @f$\Gamma@f$) or 1.0
104+extern GLfloat axis2d[2][3]; //!< @f$k@f$-vector to define section
105+extern int *n2d; //!< Number of line segment
106+extern GLfloat **kv2d; //!< @f$k@f$-vector for 2D plot [::nb][::n2d*2*3]
107+extern GLfloat **clr2d; //!< Matrix element for 2D plot [::nb][::n2d*2*4]
108+extern int nbzl2d; //!< The number of Lines of 1st Brillouin zone
109+extern GLfloat bzl2d[26][3]; //!< Lines of 1st BZ [::nbzl2d (max:26)][3]
110+extern GLfloat bzl2d_proj[26][3]; //!< Lines of 1st BZ [::nbzl2d (max:26)][3], projected into 2D plane
111+/*
112+ Equator
113+*/
114+extern GLfloat eqvec[3]; //!< @f$k@f$-vector for equator
115+extern GLfloat eqvec_fr[3]; //!< @f$k@f$-vector for equator
116+extern int *nequator; //!< The number of equator
117+extern GLfloat ****kveq; //!< @f$k@f$-vector of equator [::nb][::nequator][2][3]
118+extern GLfloat **kveq_rot; //!< @f$k@f$-vector of equator [::nb][::nequator*2*3]
119+/*
120+ Variables for mouse & cursorkey
121+*/
122+extern GLfloat sx; //!< Scale of mouse movement
123+extern GLfloat sy; //!< Scale of mouse movement
124+extern int cx; //!< Starting point of drug
125+extern int cy; //!< Starting point of drug
126+extern GLfloat scl; //!< Initial scale
127+extern GLfloat trans[3]; //!< Translation
128+extern GLfloat rot[3][3]; //!< Rotation matrix
129+extern GLfloat thetax; //!< Rotation angle
130+extern GLfloat thetay; //!< Rotation angle
131+extern GLfloat thetaz; //!< Rotation angle
132+extern GLfloat linewidth; //!< BZ/nodal-line/Fermiline width
133+/*
134+ Colors
135+*/
136+extern GLfloat black[4]; //!< Black color code
137+extern GLfloat gray[4]; //!< Gray color code
138+extern GLfloat wgray[4]; //!< Gray color code
139+extern GLfloat bgray[4]; //!< Gray color code
140+extern GLfloat white[4]; //!< White color code
141+extern GLfloat cyan[4]; //!< Cyan color code
142+extern GLfloat magenta[4]; //!< Magenta color code
143+extern GLfloat yellow[4]; //!< Yellow color code
144+extern GLfloat red[4]; //!< Red color code
145+extern GLfloat green[4]; //!< Green color code
146+extern GLfloat blue[4]; //!< Blue color code
147+/*
148+ Others
149+*/
150+extern int query; //!< Query switch
151+extern int corner[6][4]; //!< Corners of tetrahedron
152+extern GLfloat EF; //!< Fermi energy
153+enum
154+{
155+ MOUSE_SCROLL_UP = 3, //!< Mouse wheel up
156+ MOUSE_SCROLL_DOWN = 4 //!< Mouse wheel down
157+};
158+extern int nthreads;//!< Number of OpenMP threads
159+/*
160+Batch mode
161+*/
162+extern int lbatch;
163+
164+extern int refresh_interpol;
165+extern int refresh_patch;
166+extern int refresh_color;
167+extern int refresh_nodeline;
168+extern int refresh_equator;
169+extern int refresh_section;
170+extern int skip_minmax;
171+
172+#endif
Show on old repository browser