OSDN -- Open-Source-Software-Entwicklung und Downloads
Übersetzungsstatus von Deutsch translation status for de
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

fermisurfer: Commit

fermisurfer Git


Commit MetaInfo

Revision2f09127057cc816c7b429797591f17d4baa4f45c (tree)
Zeit2015-04-15 16:33:44
AutorMitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

TEST

Ändern Zusammenfassung

Diff

--- /dev/null
+++ b/fermisurfer.c
@@ -0,0 +1,2068 @@
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 <stdlib.h>
25+#include <stdio.h>
26+#include <math.h>
27+#include <GL/glut.h>
28+#include <omp.h>
29+/*
30+ Input variables
31+*/
32+int ng[3]; /* BZ grids */
33+int lshift; /* Switch for shifted Brillouin zone */
34+int nb; /* # of Bands */
35+GLfloat bvec[3][3]; /* Resiplocal lattice vector */
36+GLfloat ****eig; /* Eigenvalues [nb][ng[0]][ng[1]][ng[2]] */
37+GLfloat ****mat; /* Matrix element [nb][ng[0]][ng[1]][ng[2]] */
38+/*
39+ Switch for some modes
40+*/
41+int blackback = 1; /* Switch for black background */
42+int fcscl = 1; /* Switch for full color scale mode */
43+int fbz = 1; /* Switch for 1st Brillouin zone mode */
44+int nodeline = 0; /* Switch for node lines */
45+int lcolorbar = 1; /* Switch for colorbar */
46+int lstereo = 1; /* Switch for the stereogram */
47+/*
48+ Variables for Brillouin zone boundaries
49+*/
50+int nbzl; /* # of Lines of 1st Brillouin zone */
51+GLfloat ***bzl; /* Lines of 1st BZ [nbzl][2][3] */
52+GLfloat bragg[26][3]; /* Bragg plane vectors */
53+GLfloat brnrm[26]; /* Norms of Bragg plane vectors */
54+/*
55+ Variables for patchs
56+*/
57+int *ntri; /* # of triangle patch [nb] */
58+int *draw_band; /* Switch for drawn bands [nb] */
59+GLfloat ***nmlp; /* Normal vector of patchs [nb][ntri][3] */
60+GLfloat ****kvp; /* K-vectors of points [nb][ntri][3][3] */
61+GLfloat ***matp; /* Matrix elements of points [nb][ntri][3] */
62+GLfloat ****clr; /* Colors of points [nb][ntri][3][4] */
63+int itet = 0;
64+/*
65+ Variables for nodeline
66+*/
67+int *nnl; /* # of nodeline */
68+GLfloat ****kvnl; /* K-vector of nodeline [nb][nnl][2][3] */
69+/*
70+ Variables for mouce & cursorkey
71+*/
72+GLfloat sx, sy; /* Scale of mouce movement */
73+int cx, cy; /* Starting point of drug */
74+GLfloat scl = 1.0; /* Initial scale */
75+GLfloat trans[3] = {0.0, 0.0, 0.0}; /* Translation */
76+GLfloat rot[3][3] = {{1.0, 0.0, 0.0}, /* Rotation matrix */
77+ {0.0, 1.0, 0.0},
78+ {0.0, 0.0, 1.0}};
79+/*
80+ Colors
81+ */
82+GLfloat black[] = {0.0, 0.0, 0.0, 1.0};
83+GLfloat white[] = {1.0, 1.0, 1.0, 1.0};
84+GLfloat cyan[] = {0.0, 1.0, 1.0, 1.0};
85+GLfloat magenta[] = {1.0, 0.0, 1.0, 1.0};
86+GLfloat yellow[] = {1.0, 1.0, 0.0, 1.0};
87+GLfloat red[] = {1.0, 0.0, 0.0, 1.0};
88+GLfloat green[] = {0.0, 1.0, 0.0, 1.0};
89+GLfloat blue[] = {0.0, 0.0, 1.0, 1.0};
90+/*
91+ Others
92+*/
93+int query; /* Query switch */
94+int corner[6][4]; /* Corners of tetrahedron */
95+GLfloat def = 0.0; /* Shift of Fermi energy */
96+enum
97+ {
98+ MOUSE_LEFT_BUTTON = 0,
99+ MOUSE_MIDDLE_BUTTON = 1,
100+ MOUSE_RIGHT_BUTTON = 2,
101+ MOUSE_SCROLL_UP = 3,
102+ MOUSE_SCROLL_DOWN = 4
103+ };
104+/*
105+ Input from Fermi surface file
106+*/
107+void read_file(char *fname)
108+{
109+ int ib, i, i1, i2, i3, ierr;
110+ FILE *fp;
111+ /*
112+ Open input file.
113+ */
114+ if ((fp = fopen(fname, "r")) == NULL) {
115+ printf("file open error!!\n");
116+ exit(EXIT_FAILURE);
117+ }
118+ printf("\n##### Brillouin zone informations ##### \n\n");
119+ /*
120+ k-point grid
121+ */
122+ ierr = fscanf(fp, "%d%d%d", &ng[0], &ng[1], &ng[2]);
123+ if(ierr == 0) printf("error ! reading ng");
124+ printf("k point grid : %d %d %d \n", ng[0], ng[1], ng[2]);
125+ /*
126+ Shift of k-point grid
127+ */
128+ ierr = fscanf(fp, "%d", &lshift);
129+ if(ierr == 0) printf("error ! reading lshift");
130+ if(lshift == 0){
131+ printf("k point grid is not shifted \n");
132+ }
133+ else if(lshift == 1){
134+ printf("k point grid is shifted on 0.5, 0.5, 0.5 . \n");
135+ }
136+ else{
137+ exit(0);
138+ }
139+ /*
140+ # of bands
141+ */
142+ ierr = fscanf(fp, "%d", &nb);
143+ if(ierr == 0) printf("error ! reading nb");
144+ printf("# of bands : %d \n", nb);
145+ ntri = (int*)malloc(nb * sizeof(int));
146+ nnl = (int*)malloc(nb * sizeof(int));
147+ draw_band = (int*)malloc(nb * sizeof(int));
148+ for (ib = 0; ib < nb; ib++) draw_band[ib] = 1;
149+ /*
150+ Reciplocal lattice vectors
151+ */
152+ for (i = 0; i < 3; ++i) {
153+ ierr = fscanf(fp, "%e%e%e", &bvec[i][0], &bvec[i][1], &bvec[i][2]);
154+ if(ierr == 0) printf("error ! reading bvec");
155+ printf("bvec %d : %f %f %f \n", i + 1, bvec[i][0], bvec[i][1], bvec[i][2]);
156+ }
157+ /*
158+ Allocation of Kohn-Sham energies $ matrix elements
159+ */
160+ eig = (GLfloat****)malloc(nb * sizeof(GLfloat***));
161+ mat = (GLfloat****)malloc(nb * sizeof(GLfloat***));
162+ for (ib = 0; ib < nb; ib++){
163+ eig[ib] = (GLfloat***)malloc(ng[0] * sizeof(GLfloat**));
164+ mat[ib] = (GLfloat***)malloc(ng[0] * sizeof(GLfloat**));
165+ for (i1 = 0; i1 < ng[0]; i1++){
166+ eig[ib][i1] = (GLfloat**)malloc(ng[1] * sizeof(GLfloat*));
167+ mat[ib][i1] = (GLfloat**)malloc(ng[1] * sizeof(GLfloat*));
168+ for (i2 = 0; i2 < ng[1]; i2++){
169+ eig[ib][i1][i2] = (GLfloat*)malloc(ng[2] * sizeof(GLfloat));
170+ mat[ib][i1][i2] = (GLfloat*)malloc(ng[2] * sizeof(GLfloat));
171+ }
172+ }
173+ }
174+ /*
175+ Kohn-Sham energies
176+ */
177+ for (ib = 0; ib < nb; ++ib) {
178+ for (i1 = 0; i1 < ng[0]; ++i1) {
179+ for (i2 = 0; i2 < ng[1]; ++i2) {
180+ for (i3 = 0; i3 < ng[2]; ++i3) {
181+ ierr = fscanf(fp, "%e", &eig[ib][i1][i2][i3]);
182+ }
183+ }
184+ }
185+ }
186+ /*
187+ Matrix elements
188+ */
189+ for (ib = 0; ib < nb; ++ib) {
190+ for (i1 = 0; i1 < ng[0]; ++i1) {
191+ for (i2 = 0; i2 < ng[1]; ++i2) {
192+ for (i3 = 0; i3 < ng[2]; ++i3) {
193+ ierr = fscanf(fp, "%e", &mat[ib][i1][i2][i3]);
194+ }
195+ }
196+ }
197+ }
198+ fclose(fp);
199+ /**/
200+} /* read_file */
201+/*
202+ Initialize corners of tetrahedron
203+*/
204+void init_corner(){
205+ int i, j;
206+ int corner1[16][6][4] = {
207+ /*
208+ [0] min = 0-7
209+ */
210+ {{0, 1, 3, 7},
211+ {0, 1, 5, 7},
212+ {0, 2, 3, 7},
213+ {0, 2, 6, 7},
214+ {0, 4, 5, 7},
215+ {0, 4, 6, 7}},
216+ /*
217+ [1] min = 1-6
218+ */
219+ {{1, 0, 2, 6},
220+ {1, 0, 4, 6},
221+ {1, 3, 2, 6},
222+ {1, 3, 7, 6},
223+ {1, 5, 4, 6},
224+ {1, 5, 7, 6}},
225+ /*
226+ [2] min = 2-5
227+ */
228+ {{2, 0, 1, 5},
229+ {2, 0, 4, 5},
230+ {2, 3, 1, 5},
231+ {2, 3, 7, 5},
232+ {2, 6, 4, 5},
233+ {2, 6, 7, 5}},
234+ /*
235+ [3] min = 3-4
236+ */
237+ {{4, 0, 1, 3},
238+ {4, 0, 2, 3},
239+ {4, 5, 1, 3},
240+ {4, 5, 7, 3},
241+ {4, 6, 2, 3},
242+ {4, 6, 7, 3}},
243+ /*
244+ [4] min = 0-7, max = 3-4
245+ */
246+ {{0, 4, 5, 6},
247+ {1, 2, 3, 7},
248+ {0, 7, 2, 6},
249+ {0, 7, 1, 5},
250+ {0, 7, 1, 2},
251+ {0, 7, 6, 5}},
252+ /*
253+ [5] min = 1-6, max = 3-4
254+ */
255+ {{0, 4, 5, 6},
256+ {1, 2, 3, 7},
257+ {1, 6, 5, 7},
258+ {1, 6, 7, 2},
259+ {1, 6, 2, 0},
260+ {1, 6, 0, 5}},
261+ /*
262+ [6] min = 2-5, max = 3-4
263+ */
264+ {{0, 4, 5, 6},
265+ {1, 2, 3, 7},
266+ {2, 5, 7, 6},
267+ {2, 5, 6, 0},
268+ {2, 5, 0, 1},
269+ {2, 5, 1, 7}},
270+ /*
271+ [7] min = 3-4, max = 0-7
272+ */
273+ {{0, 1, 2, 4},
274+ {7, 3, 5, 6},
275+ {3, 4, 1, 5},
276+ {3, 4, 5, 6},
277+ {3, 4, 6, 2},
278+ {3, 4, 2, 1}},
279+ /*
280+ [8] min = 2-5, max = 0-7
281+ */
282+ {{0, 1, 2, 4},
283+ {7, 3, 5, 6},
284+ {2, 5, 1, 3},
285+ {2, 5, 3, 6},
286+ {2, 5, 6, 4},
287+ {2, 5, 4, 1}},
288+ /*
289+ [9] min = 1-6, max = 0-7
290+ */
291+ {{0, 1, 2, 4},
292+ {7, 3, 5, 6},
293+ {1, 6, 2, 3},
294+ {1, 6, 3, 5},
295+ {1, 6, 5, 4},
296+ {1, 6, 4, 2}},
297+ /*
298+ [10] min = 0-7, max = 1-6
299+ */
300+ {{1, 0, 3, 5},
301+ {6, 2, 4, 7},
302+ {0, 7, 2, 3},
303+ {0, 7, 3, 5},
304+ {0, 7, 5, 4},
305+ {0, 7, 4, 2}},
306+ /*
307+ [11] min = 3-4, max = 1-6
308+ */
309+ {{1, 0, 3, 5},
310+ {6, 2, 4, 7},
311+ {3, 4, 0, 2},
312+ {3, 4, 2, 7},
313+ {3, 4, 7, 5},
314+ {3, 4, 5, 0}},
315+ /*
316+ [12] min = 2-5, max = 1-6
317+ */
318+ {{1, 0, 3, 5},
319+ {6, 2, 4, 7},
320+ {2, 5, 0, 3},
321+ {2, 5, 3, 7},
322+ {2, 5, 7, 4},
323+ {2, 5, 4, 0}},
324+ /*
325+ [13] min = 0-7, max = 2-5
326+ */
327+ {{2, 0, 3, 6},
328+ {5, 1, 4, 7},
329+ {0, 7, 1, 3},
330+ {0, 7, 3, 6},
331+ {0, 7, 6, 4},
332+ {0, 7, 4, 1}},
333+ /*
334+ [14] min = 1-6, max = 2-5
335+ */
336+ {{2, 0, 3, 6},
337+ {5, 1, 4, 7},
338+ {1, 6, 0, 3},
339+ {1, 6, 3, 7},
340+ {1, 6, 7, 4},
341+ {1, 6, 4, 0}},
342+ /*
343+ [15] min = 3-4, max = 2-5
344+ */
345+ {{2, 0, 3, 6},
346+ {5, 1, 4, 7},
347+ {3, 4, 0, 1},
348+ {3, 4, 1, 7},
349+ {3, 4, 7, 6},
350+ {3, 4, 6, 0}},
351+ };
352+ /**/
353+ for (i = 0; i < 6; ++i){
354+ for (j = 0; j < 4; ++j){
355+ corner[i][j] = corner1[itet][i][j];
356+ }
357+ }
358+}
359+/*
360+ Compute Bragg vetor
361+*/
362+void bragg_vector(){
363+ int i1, i2, i3, i, ibr;
364+ /**/
365+ ibr = 0;
366+ /**/
367+ for(i1 = -1; i1 <= 1; ++i1){
368+ for(i2 = -1; i2 <= 1; ++i2){
369+ for(i3 = -1; i3 <= 1; ++i3){
370+ /**/
371+ if(i1 == 0 && i2 == 0 && i3 ==0){
372+ }
373+ else {
374+ for(i = 0; i < 3; ++i) bragg[ibr][i] = (( GLfloat)i1 * bvec[0][i]
375+ + (GLfloat)i2 * bvec[1][i]
376+ + (GLfloat)i3 * bvec[2][i]) * 0.5;
377+ /**/
378+ brnrm[ibr] = bragg[ibr][0] * bragg[ibr][0]
379+ + bragg[ibr][1] * bragg[ibr][1]
380+ + bragg[ibr][2] * bragg[ibr][2];
381+ /**/
382+ ibr = ibr + 1;
383+ }
384+ }
385+ }
386+ }
387+} /* bragg_vector */
388+/*
389+ Solve linear system
390+*/
391+GLfloat solve3(GLfloat a[3][3], GLfloat b[3]){
392+ int i;
393+ GLfloat det, c[3];
394+ /**/
395+ det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
396+ + a[0][1] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
397+ + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
398+ /**/
399+ c[0] = b[0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
400+ + b[1] * (a[0][2] * a[2][1] - a[0][1] * a[2][2])
401+ + b[2] * (a[0][1] * a[1][2] - a[0][2] * a[1][1]);
402+ /**/
403+ c[1] = b[0] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
404+ + b[1] * (a[0][0] * a[2][2] - a[0][2] * a[2][0])
405+ + b[2] * (a[0][2] * a[1][0] - a[0][0] * a[1][2]);
406+ /**/
407+ c[2] = b[0] * (a[1][0] * a[2][1] - a[1][1] * a[2][0])
408+ + b[1] * (a[0][1] * a[2][0] - a[0][0] * a[2][1])
409+ + b[2] * (a[0][0] * a[1][1] - a[0][1] * a[1][0]);
410+ /**/
411+ for(i=0;i<3;++i) b[i] = c[i] / det;
412+ return det;
413+ /**/
414+}
415+/*
416+ Corner of Bragg plane
417+*/
418+int bragg_vert(int ibr, int jbr, int nbr, GLfloat vert[3], GLfloat vert2[3]){
419+ int kbr, i, lbr, nbr0;
420+ GLfloat bmat[3][3], rhs[3], prod, thr = 0.0001, det;
421+ /**/
422+ nbr0 = nbr;
423+ /**/
424+ for(kbr = nbr0; kbr < 26; ++kbr){
425+ /**/
426+ for(i=0;i<3;++i) bmat[0][i] = bragg[ibr][i];
427+ for(i=0;i<3;++i) bmat[1][i] = bragg[jbr][i];
428+ for(i=0;i<3;++i) bmat[2][i] = bragg[kbr][i];
429+ /**/
430+ rhs[0] = brnrm[ibr];
431+ rhs[1] = brnrm[jbr];
432+ rhs[2] = brnrm[kbr];
433+ /*
434+ if Bragg planes do not cross, roop next kbr
435+ */
436+ det = solve3(bmat, rhs);
437+ if(fabsf(det) < thr) continue;
438+ /*
439+ if vert0 = vert1, roop next kbr
440+ */
441+ prod = (vert2[0] - rhs[0]) * (vert2[0] - rhs[0])
442+ + (vert2[1] - rhs[1]) * (vert2[1] - rhs[1])
443+ + (vert2[2] - rhs[2]) * (vert2[2] - rhs[2]);
444+ if(prod < thr) continue;
445+ /*
446+ is corner really in 1st BZ ?
447+ */
448+ i = 0;
449+ for(lbr = 0; lbr < 26; ++lbr){
450+ prod = bragg[lbr][0] * rhs[0]
451+ + bragg[lbr][1] * rhs[1]
452+ + bragg[lbr][2] * rhs[2];
453+ /**/
454+ if(prod > brnrm[lbr] + thr){
455+ i = 1;
456+ break;
457+ }
458+ }
459+ if(i == 1) {
460+ }
461+ else{
462+ for(i=0;i<3;++i) vert[i] = rhs[i];
463+ return kbr + 1;
464+ }
465+ }
466+ /*
467+ this line is not a BZ boundary
468+ */
469+ return 0;
470+ /**/
471+}/* bragg_vert */
472+/*
473+ Compute Brillouin zone boundariy lines
474+*/
475+void bz_lines(){
476+ /**/
477+ int ibr, jbr, nbr, ibzl, i, j, lvert;
478+ GLfloat vert[2][3];
479+ /**/
480+ ibzl = 0;
481+ /**/
482+ for(ibr = 0; ibr < 26; ++ibr){
483+ for(jbr = 0; jbr < 26; ++jbr){
484+ /**/
485+ for(i=0;i<3;++i) vert[1][i] = 0.0;
486+ nbr = 0;
487+ lvert = bragg_vert(ibr, jbr, nbr, vert[0], vert[1]);
488+ if(lvert == 0) continue;
489+ nbr = lvert;
490+ /**/
491+ lvert = bragg_vert(ibr, jbr, nbr, vert[1], vert[0]);
492+ if(lvert == 0) continue;
493+ /**/
494+ if(query == 1){
495+ nbzl = nbzl + 1;
496+ }
497+ else{
498+ for(i=0;i<2;++i) for(j=0;j<3;++j) bzl[ibzl][i][j] = vert[i][j];
499+ ibzl = ibzl + 1;
500+ }
501+ /**/
502+ }
503+ }
504+ /**/
505+ if(query == 1){
506+ printf("# of lines for BZ : %d \n", nbzl);
507+ /**/
508+ bzl = (GLfloat***)malloc(nbzl * sizeof(GLfloat*));
509+ for(ibzl = 0; ibzl < nbzl; ++ibzl){
510+ bzl[ibzl] = (GLfloat**)malloc(2 * sizeof(GLfloat*));
511+ for(i = 0; i < 2; ++i){
512+ bzl[ibzl][i] = (GLfloat*)malloc(3 * sizeof(GLfloat));
513+ }
514+ }
515+ /**/
516+ }
517+ /**/
518+} /* bz_lines */
519+/*
520+ Max and Minimum in Brillouine zone
521+*/
522+void max_and_min_bz(){
523+ int ib, i1, i2, i3;
524+ GLfloat eigmin, eigmax, matmin, matmax;
525+ /**/
526+ printf("\n##### Max. and Min. of each bands ##### \n\n");
527+ printf("Band Eig_Min. Eig_Max Mat_Min Mat_Max \n");
528+ for(ib = 0; ib < nb; ib++){
529+ eigmax = - 100000000.0000;
530+ eigmin = 100000000.0000;
531+ matmax = - 100000000.0000;
532+ matmin = 100000000.0000;
533+ for (i1 = 0; i1 < ng[0]; ++i1) {
534+ for (i2 = 0; i2 < ng[1]; ++i2) {
535+ for (i3 = 0; i3 < ng[2]; ++i3) {
536+ if(eig[ib][i1][i2][i3] > eigmax) eigmax = eig[ib][i1][i2][i3];
537+ if(eig[ib][i1][i2][i3] < eigmin) eigmin = eig[ib][i1][i2][i3];
538+ if(mat[ib][i1][i2][i3] > matmax) matmax = mat[ib][i1][i2][i3];
539+ if(mat[ib][i1][i2][i3] < matmin) matmin = mat[ib][i1][i2][i3];
540+ }
541+ }
542+ }
543+ printf("%d %f %f %f %f \n", ib + 1, eigmin, eigmax, matmin, matmax);
544+ }
545+ /**/
546+}/* max_and_min_bz */
547+/*
548+ Sort eigenvalues
549+*/
550+void eigsort(int n, GLfloat* eig2, GLfloat* mat2, GLfloat kvec2[][3]){
551+ int i, j, k;
552+ GLfloat tmp;
553+ /**/
554+ for (i = 0; i < n - 1; ++i){
555+ for (j = i + 1; j < n; ++j){
556+ if(eig2[j] < eig2[i]) {
557+ tmp = eig2[j];
558+ eig2[j] = eig2[i];
559+ eig2[i] = tmp;
560+ /**/
561+ tmp = mat2[j];
562+ mat2[j] = mat2[i];
563+ mat2[i] = tmp;
564+ /**/
565+ for (k = 0; k < 3; ++k){
566+ tmp = kvec2[j][k];
567+ kvec2[j][k] = kvec2[i][k];
568+ kvec2[i][k] = tmp;
569+ }
570+ }
571+ }
572+ }
573+} /* eigsort */
574+/*
575+ Calculate normal vector
576+*/
577+void normal_vec(GLfloat in1[3], GLfloat in2[3], GLfloat in3[3], GLfloat out[3]){
578+ int i;
579+ GLfloat norm;
580+ out[0] = in1[1] * in2[2] - in1[2] * in2[1]
581+ + in2[1] * in3[2] - in2[2] * in3[1]
582+ + in3[1] * in1[2] - in3[2] * in1[1];
583+ out[1] = in1[2] * in2[0] - in1[0] * in2[2]
584+ + in2[2] * in3[0] - in2[0] * in3[2]
585+ + in3[2] * in1[0] - in3[0] * in1[2];
586+ out[2] = in1[0] * in2[1] - in1[1] * in2[0]
587+ + in2[0] * in3[1] - in2[1] * in3[0]
588+ + in3[0] * in1[1] - in3[1] * in1[0];
589+ norm = sqrtf(out[0]*out[0] + out[1]*out[1] + out[2]*out[2]);
590+ for(i=0;i<3;i++) out[i] = out[i] / norm;
591+} /* normal_vec */
592+/*
593+ Store triangle patch
594+*/
595+void triangle(int ib, int nbr, GLfloat mat1[3], GLfloat kvec1[3][3]){
596+ /**/
597+ int ibr, i, j;
598+ GLfloat prod[3], thr = 0.0000, mat2[3], kvec2[3][3];
599+ /**/
600+ if(fbz == 1){
601+ /**/
602+ for(ibr = nbr; ibr < 26; ++ibr){
603+ /**/
604+ for(i = 0; i < 3; ++i){
605+ prod[i] = bragg[ibr][0] * kvec1[i][0]
606+ + bragg[ibr][1] * kvec1[i][1]
607+ + bragg[ibr][2] * kvec1[i][2];
608+ }
609+ eigsort(3, prod, mat1, kvec1);
610+ /**/
611+ if(brnrm[ibr] + thr < prod[0]) {
612+ return;
613+ }
614+ else if(brnrm[ibr] + thr < prod[1]){
615+ mat2[0] = mat1[0];
616+ mat2[1] = mat1[0] * (brnrm[ibr] - prod[1]) / (prod[0] - prod[1])
617+ + mat1[1] * (brnrm[ibr] - prod[0]) / (prod[1] - prod[0]);
618+ mat2[2] = mat1[0] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
619+ + mat1[2] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
620+ for(i=0;i<3;++i) kvec2[0][i] = kvec1[0][i];
621+ for(i=0;i<3;++i) kvec2[1][i] = kvec1[0][i] * (brnrm[ibr] - prod[1]) / (prod[0] - prod[1])
622+ + kvec1[1][i] * (brnrm[ibr] - prod[0]) / (prod[1] - prod[0]);
623+ for(i=0;i<3;++i) kvec2[2][i] = kvec1[0][i] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
624+ + kvec1[2][i] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
625+ triangle(ib, ibr + 1, mat2, kvec2);
626+ return;
627+ }
628+ else if(brnrm[ibr] + thr < prod[2]){
629+ mat2[0] = mat1[0];
630+ mat2[1] = mat1[1];
631+ mat2[2] = mat1[0] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
632+ + mat1[2] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
633+ for(i=0;i<3;++i) kvec2[0][i] = kvec1[0][i];
634+ for(i=0;i<3;++i) kvec2[1][i] = kvec1[1][i];
635+ for(i=0;i<3;++i) kvec2[2][i] = kvec1[0][i] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
636+ + kvec1[2][i] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
637+ triangle(ib, ibr + 1, mat2, kvec2);
638+ /**/
639+ mat2[0] = mat1[1] * (brnrm[ibr] - prod[2]) / (prod[1] - prod[2])
640+ + mat1[2] * (brnrm[ibr] - prod[1]) / (prod[2] - prod[1]);
641+ mat2[1] = mat1[1];
642+ mat2[2] = mat1[0] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
643+ + mat1[2] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
644+ for(i=0;i<3;++i) kvec2[0][i] = kvec1[1][i] * (brnrm[ibr] - prod[2]) / (prod[1] - prod[2])
645+ + kvec1[2][i] * (brnrm[ibr] - prod[1]) / (prod[2] - prod[1]);
646+ for(i=0;i<3;++i) kvec2[1][i] = kvec1[1][i];
647+ for(i=0;i<3;++i) kvec2[2][i] = kvec1[0][i] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
648+ + kvec1[2][i] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
649+ triangle(ib, ibr + 1, mat2, kvec2);
650+ return;
651+ }
652+ else{
653+ } /* brnrm[ibr] + thr < prod */
654+ } /* for ibr = 1; ibr < 26*/
655+ } /* if fbz == 1 */
656+ /**/
657+ if(query == 1){
658+ ntri[ib] = ntri[ib] + 1;
659+ }
660+ else{
661+ normal_vec(kvec1[0], kvec1[1], kvec1[2], nmlp[ib][ntri[ib]]);
662+ for(i = 0; i < 3; ++i){
663+ matp[ib][ntri[ib]][i] = mat1[i];
664+ for(j = 0; j < 3; ++j){
665+ kvp[ib][ntri[ib]][i][j] = kvec1[i][j];
666+ }
667+ }
668+ ntri[ib] = ntri[ib] + 1;
669+ }
670+ /**/
671+}/* triangle */
672+/*
673+ Tetrahedrron method
674+*/
675+void tetrahedron(int ib, GLfloat eig1[8], GLfloat mat1[8], GLfloat kvec1[8][3]){
676+ /**/
677+ int it, i, j;
678+ GLfloat eig2[4], mat2[4], kvec2[4][3], a[4][4], kvec3[3][3], mat3[3];
679+ /**/
680+ for (it = 0; it < 6; ++it) {
681+ /*
682+ Define corners of the tetrahedron
683+ */
684+ for (i = 0; i < 4; ++i) {
685+ eig2[ i] = eig1[corner[it][i]];
686+ mat2[ i] = mat1[corner[it][i]];
687+ /**/
688+ for(j = 0; j < 3; ++j) kvec2[i][j] = bvec[0][j] * kvec1[corner[it][i]][0]
689+ + bvec[1][j] * kvec1[corner[it][i]][1]
690+ + bvec[2][j] * kvec1[corner[it][i]][2];
691+ }
692+ /*
693+ Sort of eig1
694+ */
695+ eigsort(4, eig2, mat2, kvec2);
696+ for (i = 0; i < 4; ++i){
697+ for (j = 0; j < 4; ++j){
698+ a[i][j] = (0.00 - eig2[j]) / (eig2[i] - eig2[j]);
699+ }
700+ }
701+ /*
702+ Draw triangle in each cases
703+ */
704+ if(eig2[0] <= 0.00 && 0.00 < eig2[1]) {
705+ for(i = 0; i < 3; ++i) kvec3[0][i] = kvec2[0][i] * a[0][1] + kvec2[1][i] * a[1][0];
706+ for(i = 0; i < 3; ++i) kvec3[1][i] = kvec2[0][i] * a[0][2] + kvec2[2][i] * a[2][0];
707+ for(i = 0; i < 3; ++i) kvec3[2][i] = kvec2[0][i] * a[0][3] + kvec2[3][i] * a[3][0];
708+ mat3[0] = mat2[0] * a[0][1] + mat2[1] * a[1][0];
709+ mat3[1] = mat2[0] * a[0][2] + mat2[2] * a[2][0];
710+ mat3[2] = mat2[0] * a[0][3] + mat2[3] * a[3][0];
711+ triangle(ib, 0, mat3, kvec3);
712+ }
713+ else if(eig2[1] <= 0.00 && 0.00 < eig2[2]){
714+ for(i = 0; i < 3; ++i) kvec3[0][i] = kvec2[0][i] * a[0][2] + kvec2[2][i] * a[2][0];
715+ for(i = 0; i < 3; ++i) kvec3[1][i] = kvec2[0][i] * a[0][3] + kvec2[3][i] * a[3][0];
716+ for(i = 0; i < 3; ++i) kvec3[2][i] = kvec2[1][i] * a[1][2] + kvec2[2][i] * a[2][1];
717+ mat3[0] = mat2[0] * a[0][2] + mat2[2] * a[2][0];
718+ mat3[1] = mat2[0] * a[0][3] + mat2[3] * a[3][0];
719+ mat3[2] = mat2[1] * a[1][2] + mat2[2] * a[2][1];
720+ triangle(ib, 0, mat3, kvec3);
721+ /**/
722+ for(i = 0; i < 3; ++i) kvec3[0][i] = kvec2[1][i] * a[1][3] + kvec2[3][i] * a[3][1];
723+ for(i = 0; i < 3; ++i) kvec3[1][i] = kvec2[0][i] * a[0][3] + kvec2[3][i] * a[3][0];
724+ for(i = 0; i < 3; ++i) kvec3[2][i] = kvec2[1][i] * a[1][2] + kvec2[2][i] * a[2][1];
725+ mat3[0] = mat2[1] * a[1][3] + mat2[3] * a[3][1];
726+ mat3[1] = mat2[0] * a[0][3] + mat2[3] * a[3][0];
727+ mat3[2] = mat2[1] * a[1][2] + mat2[2] * a[2][1];
728+ triangle(ib, 0, mat3, kvec3);
729+ }
730+ else if(eig2[2] <= 0.00 && 0.00 < eig2[3]){
731+ for(i = 0; i < 3; ++i) kvec3[0][i] = kvec2[3][i] * a[3][0] + kvec2[0][i] * a[0][3];
732+ for(i = 0; i < 3; ++i) kvec3[1][i] = kvec2[3][i] * a[3][1] + kvec2[1][i] * a[1][3];
733+ for(i = 0; i < 3; ++i) kvec3[2][i] = kvec2[3][i] * a[3][2] + kvec2[2][i] * a[2][3];
734+ mat3[0] = mat2[3] * a[3][0] + mat2[0] * a[0][3];
735+ mat3[1] = mat2[3] * a[3][1] + mat2[1] * a[1][3];
736+ mat3[2] = mat2[3] * a[3][2] + mat2[2] * a[2][3];
737+ triangle(ib, 0, mat3, kvec3);
738+ }
739+ else{
740+ }
741+ }
742+}/* tetrahedron */
743+/*
744+ Patches for FSs
745+*/
746+void fermi_patch()
747+{
748+ int ib, i1, i2, i3, ii1, ii2, ii3, j1, j2, j3, start[3], i, j;
749+ GLfloat kvec1[8][3], eig1[8], mat1[8];
750+ /**/
751+ if(fbz == 1){
752+ for(i1 = 0; i1 < 3;++i1) start[i1] = - ng[i1];
753+ }
754+ else{
755+ for(i1 = 0; i1 < 3;++i1) start[i1] = 0;
756+ }
757+ /**/
758+#pragma omp parallel default(none) \
759+ shared(nb,ntri,start,ng,eig,mat,lshift) \
760+ private(ib,j1,j2,j3,i1,i2,i3,ii1,ii2,ii3,kvec1,eig1,mat1,i,j)
761+ {
762+#pragma omp for nowait
763+ for (ib = 0; ib < nb; ++ib) {
764+ ntri[ib] = 0;
765+ for (j1 = start[0]; j1 < ng[0]; ++j1) {
766+ for (j2 = start[1]; j2 < ng[1]; ++j2) {
767+ for (j3 = start[2]; j3 < ng[2]; ++j3) {
768+ /**/
769+ i1 = j1;
770+ i2 = j2;
771+ i3 = j3;
772+ ii1 = j1 + 1;
773+ ii2 = j2 + 1;
774+ ii3 = j3 + 1;
775+ /**/
776+ kvec1[0][0] = (GLfloat)i1 / (GLfloat)ng[0];
777+ kvec1[1][0] = (GLfloat)i1 / (GLfloat)ng[0];
778+ kvec1[2][0] = (GLfloat)i1 / (GLfloat)ng[0];
779+ kvec1[3][0] = (GLfloat)i1 / (GLfloat)ng[0];
780+ kvec1[4][0] = (GLfloat)ii1 / (GLfloat)ng[0];
781+ kvec1[5][0] = (GLfloat)ii1 / (GLfloat)ng[0];
782+ kvec1[6][0] = (GLfloat)ii1 / (GLfloat)ng[0];
783+ kvec1[7][0] = (GLfloat)ii1 / (GLfloat)ng[0];
784+ /**/
785+ kvec1[0][1] = (GLfloat)i2 / (GLfloat)ng[1];
786+ kvec1[1][1] = (GLfloat)i2 / (GLfloat)ng[1];
787+ kvec1[2][1] = (GLfloat)ii2 / (GLfloat)ng[1];
788+ kvec1[3][1] = (GLfloat)ii2 / (GLfloat)ng[1];
789+ kvec1[4][1] = (GLfloat)i2 / (GLfloat)ng[1];
790+ kvec1[5][1] = (GLfloat)i2 / (GLfloat)ng[1];
791+ kvec1[6][1] = (GLfloat)ii2 / (GLfloat)ng[1];
792+ kvec1[7][1] = (GLfloat)ii2 / (GLfloat)ng[1];
793+ /**/
794+ kvec1[0][2] = (GLfloat)i3 / (GLfloat)ng[2];
795+ kvec1[1][2] = (GLfloat)ii3 / (GLfloat)ng[2];
796+ kvec1[2][2] = (GLfloat)i3 / (GLfloat)ng[2];
797+ kvec1[3][2] = (GLfloat)ii3 / (GLfloat)ng[2];
798+ kvec1[4][2] = (GLfloat)i3 / (GLfloat)ng[2];
799+ kvec1[5][2] = (GLfloat)ii3 / (GLfloat)ng[2];
800+ kvec1[6][2] = (GLfloat)i3 / (GLfloat)ng[2];
801+ kvec1[7][2] = (GLfloat)ii3 / (GLfloat)ng[2];
802+ /**/
803+ if(lshift == 1){
804+ for(i=0;i<8;i++)for(j=0;j<3;j++) kvec1[i][j] = kvec1[i][j] + 0.50 / (GLfloat)ng[j];
805+ }
806+ /**/
807+ if(i1 < 0) i1 = i1 + ng[0];
808+ if(i2 < 0) i2 = i2 + ng[1];
809+ if(i3 < 0) i3 = i3 + ng[2];
810+ if(ii1 < 0) ii1 = ii1 + ng[0];
811+ if(ii2 < 0) ii2 = ii2 + ng[1];
812+ if(ii3 < 0) ii3 = ii3 + ng[2];
813+ /**/
814+ if(ii1 >= ng[0]) ii1 = 0;
815+ if(ii2 >= ng[1]) ii2 = 0;
816+ if(ii3 >= ng[2]) ii3 = 0;
817+ /**/
818+ eig1[0] = eig[ib][ i1][ i2][ i3];
819+ eig1[1] = eig[ib][ i1][ i2][ii3];
820+ eig1[2] = eig[ib][ i1][ii2][ i3];
821+ eig1[3] = eig[ib][ i1][ii2][ii3];
822+ eig1[4] = eig[ib][ii1][ i2][ i3];
823+ eig1[5] = eig[ib][ii1][ i2][ii3];
824+ eig1[6] = eig[ib][ii1][ii2][ i3];
825+ eig1[7] = eig[ib][ii1][ii2][ii3];
826+ /**/
827+ mat1[0] = mat[ib][ i1][ i2][ i3];
828+ mat1[1] = mat[ib][ i1][ i2][ii3];
829+ mat1[2] = mat[ib][ i1][ii2][ i3];
830+ mat1[3] = mat[ib][ i1][ii2][ii3];
831+ mat1[4] = mat[ib][ii1][ i2][ i3];
832+ mat1[5] = mat[ib][ii1][ i2][ii3];
833+ mat1[6] = mat[ib][ii1][ii2][ i3];
834+ mat1[7] = mat[ib][ii1][ii2][ii3];
835+ /**/
836+ tetrahedron(ib, eig1, mat1, kvec1);
837+ }
838+ }
839+ }
840+ }
841+ } /** End of parallel region **/
842+ /**/
843+ if(query == 1){
844+ printf("band # of patchs \n");
845+ for(ib =0; ib < nb; ib++){
846+ printf("%d %d \n", ib + 1, ntri[ib]);
847+ }
848+ printf("\n");
849+ /*
850+ Allocation of triangler patches
851+ */
852+ nmlp = (GLfloat***)malloc(nb * sizeof(GLfloat**));
853+ matp = (GLfloat***)malloc(nb * sizeof(GLfloat**));
854+ clr = (GLfloat****)malloc(nb * sizeof(GLfloat***));
855+ kvp = (GLfloat****)malloc(nb * sizeof(GLfloat***));
856+ for (ib = 0; ib < nb; ++ib) {
857+ nmlp[ib] = (GLfloat**)malloc(ntri[ib] * sizeof(GLfloat*));
858+ matp[ib] = (GLfloat**)malloc(ntri[ib] * sizeof(GLfloat*));
859+ clr[ib] = (GLfloat***)malloc(ntri[ib] * sizeof(GLfloat**));
860+ kvp[ib] = (GLfloat***)malloc(ntri[ib] * sizeof(GLfloat**));
861+ for (i1 = 0; i1 < ntri[ib]; ++i1){
862+ nmlp[ib][i1] = (GLfloat*)malloc(3 * sizeof(GLfloat));
863+ matp[ib][i1] = (GLfloat*)malloc(3 * sizeof(GLfloat));
864+ clr[ib][i1] = (GLfloat**)malloc(3 * sizeof(GLfloat*));
865+ kvp[ib][i1] = (GLfloat**)malloc(3 * sizeof(GLfloat*));
866+ for (i2 = 0; i2 < 3; ++i2){
867+ kvp[ib][i1][i2] = (GLfloat*)malloc(3 * sizeof(GLfloat));
868+ clr[ib][i1][i2] = (GLfloat*)malloc(4 * sizeof(GLfloat));
869+ }
870+ }
871+ }
872+ /**/
873+ }
874+ /**/
875+} /* fermi_patch */
876+/*
877+ Max. & Min. of matrix elements.
878+*/
879+void max_and_min(){
880+ int ib, itri, i, j, ierr;
881+ GLfloat matmax, matmin, mat2;
882+ /**/
883+ matmax = - 100000000.0000;
884+ matmin = 100000000.0000;
885+ /**/
886+ for(ib = 0; ib < nb; ib++){
887+ for(itri = 0; itri < ntri[ib]; ++itri){
888+ for(i = 0; i < 3; ++i){
889+ if(matp[ib][itri][i] > matmax) matmax = matp[ib][itri][i];
890+ if(matp[ib][itri][i] < matmin) matmin = matp[ib][itri][i];
891+ }
892+ }
893+ }
894+ /**/
895+ printf("Max. value : %f \n", matmax);
896+ printf("Min. value : %f \n \n", matmin);
897+ /**/
898+ if(fcscl == 2){
899+ printf("Set min. value : ");
900+ ierr = scanf("%f", &matmin);
901+ if(ierr == 0) printf("error ! reading min");
902+ printf("Set max. value : ");
903+ ierr = scanf("%f", &matmax);
904+ if(ierr == 0) printf("error ! reading max");
905+ }
906+ /**/
907+ if(fcscl == 1 || fcscl == 2){
908+ for(ib = 0; ib < nb; ib++){
909+ for(itri = 0; itri < ntri[ib]; ++itri){
910+ for (i = 0; i < 3; ++i){
911+ /**/
912+ mat2 = (matp[ib][itri][i] - matmin) / (matmax - matmin);
913+ mat2 = mat2 * 4.0;
914+ /**/
915+ if(mat2 <= 1.0) {
916+ for(j=0;j<4;++j) clr[ib][itri][i][j] = cyan[j] * mat2 + blue[j] * (1.0 - mat2);
917+ }
918+ else if(mat2 <= 2.0){
919+ mat2 = mat2 - 1.0;
920+ for(j=0;j<4;++j) clr[ib][itri][i][j] = green[j] * mat2 + cyan[j] * (1.0 - mat2);
921+ }
922+ else if(mat2 <= 3.0){
923+ mat2 = mat2 - 2.0;
924+ for(j=0;j<4;++j) clr[ib][itri][i][j] = yellow[j] * mat2 + green[j] * (1.0 - mat2);
925+ }
926+ else{
927+ mat2 = mat2 - 3.0;
928+ for(j=0;j<4;++j) clr[ib][itri][i][j] = red[j] * mat2 + yellow[j] * (1.0 - mat2);
929+ }
930+ }
931+ }
932+ }
933+ }
934+ else if(fcscl == 4){
935+ for(ib = 0; ib < nb; ib++){
936+ for(itri = 0; itri < ntri[ib]; ++itri){
937+ for (i = 0; i < 3; ++i){
938+ /**/
939+ mat2 = matp[ib][itri][i] / 6.283185307;
940+ mat2 = mat2 - floorf(mat2);
941+ mat2 = mat2 * 6.0;
942+ /**/
943+ if(mat2 <= 1.0) {
944+ for(j=0;j<4;++j) clr[ib][itri][i][j] = yellow[j] * mat2 + red[j] * (1.0 - mat2);
945+ }
946+ else if(mat2 <= 2.0){
947+ mat2 = mat2 - 1.0;
948+ for(j=0;j<4;++j) clr[ib][itri][i][j] = green[j] * mat2 + yellow[j] * (1.0 - mat2);
949+ }
950+ else if(mat2 <= 3.0){
951+ mat2 = mat2 - 2.0;
952+ for(j=0;j<4;++j) clr[ib][itri][i][j] = cyan[j] * mat2 + green[j] * (1.0 - mat2);
953+ }
954+ else if(mat2 <= 4.0){
955+ mat2 = mat2 - 3.0;
956+ for(j=0;j<4;++j) clr[ib][itri][i][j] = blue[j] * mat2 + cyan[j] * (1.0 - mat2);
957+ }
958+ else if(mat2 <= 5.0){
959+ mat2 = mat2 - 4.0;
960+ for(j=0;j<4;++j) clr[ib][itri][i][j] = magenta[j] * mat2 + blue[j] * (1.0 - mat2);
961+ }
962+ else{
963+ mat2 = mat2 - 5.0;
964+ for(j=0;j<4;++j) clr[ib][itri][i][j] = red[j] * mat2 + magenta[j] * (1.0 - mat2);
965+ }
966+ }
967+ }
968+ }
969+ }
970+ else{
971+ for(ib = 0; ib < nb; ib++){
972+ /**/
973+ mat2 = 1.0 / (GLfloat)(nb - 1) * (GLfloat)ib;
974+ mat2 = mat2 * 4.0;
975+ /**/
976+ if(mat2 <= 1.0) {
977+ for(itri = 0; itri < ntri[ib]; ++itri){
978+ for (i = 0; i < 3; ++i){
979+ for(j=0;j<4;++j) clr[ib][itri][i][j] = cyan[j] * mat2 + blue[j] * (1.0 - mat2);
980+ }
981+ }
982+ }
983+ else if(mat2 <= 2.0){
984+ mat2 = mat2 - 1.0;
985+ for(itri = 0; itri < ntri[ib]; ++itri){
986+ for (i = 0; i < 3; ++i){
987+ for(j=0;j<4;++j) clr[ib][itri][i][j] = green[j] * mat2 + cyan[j] * (1.0 - mat2);
988+ }
989+ }
990+ }
991+ else if(mat2 <= 3.0){
992+ mat2 = mat2 - 2.0;
993+ for(itri = 0; itri < ntri[ib]; ++itri){
994+ for (i = 0; i < 3; ++i){
995+ for(j=0;j<4;++j) clr[ib][itri][i][j] = yellow[j] * mat2 + green[j] * (1.0 - mat2);
996+ }
997+ }
998+ }
999+ else{
1000+ mat2 = mat2 - 3.0;
1001+ for(itri = 0; itri < ntri[ib]; ++itri){
1002+ for (i = 0; i < 3; ++i){
1003+ for(j=0;j<4;++j) clr[ib][itri][i][j] = red[j] * mat2 + yellow[j] * (1.0 - mat2);
1004+ }
1005+ }
1006+ }
1007+ }
1008+ }
1009+ /**/
1010+} /* max_and_min */
1011+/*
1012+ Node line
1013+*/
1014+void calc_nodeline(){
1015+ int ib, itri, i, j;
1016+ GLfloat mprod[2];
1017+ /*
1018+ Query
1019+ */
1020+#pragma omp parallel default(none) \
1021+ shared(nb,nnl,matp,ntri) \
1022+ private(ib,itri,mprod)
1023+ {
1024+#pragma omp for
1025+ for(ib = 0; ib < nb; ib++){
1026+ nnl[ib] = 0;
1027+ for(itri = 0; itri < ntri[ib]; ++itri){
1028+ /**/
1029+ mprod[0] = matp[ib][itri][0] * matp[ib][itri][1];
1030+ mprod[1] = matp[ib][itri][1] * matp[ib][itri][2];
1031+ /**/
1032+ if( fabsf(matp[ib][itri][0]) < 0.00001 && fabsf(matp[ib][itri][1]) < 0.00001){
1033+ nnl[ib] = nnl[ib] + 1;
1034+ }
1035+ else if(fabsf(matp[ib][itri][0]) < 0.00001 && fabsf(matp[ib][itri][2]) < 0.00001){
1036+ nnl[ib] = nnl[ib] + 1;
1037+ }
1038+ else if(fabsf(matp[ib][itri][1]) < 0.00001 && fabsf(matp[ib][itri][2]) < 0.00001){
1039+ nnl[ib] = nnl[ib] + 1;
1040+ }
1041+ else if(mprod[0] < 0.0){
1042+ if(mprod[1] < 0.0){
1043+ nnl[ib] = nnl[ib] + 1;
1044+ }
1045+ else{
1046+ nnl[ib] = nnl[ib] + 1;
1047+ }
1048+ }
1049+ else if(mprod[1] < 0.0){
1050+ nnl[ib] = nnl[ib] + 1;
1051+ }
1052+ }
1053+ }
1054+ } /** End of parallel region **/
1055+ /**/
1056+ printf("band # of nodeline \n");
1057+ for(ib =0; ib < nb; ib++){
1058+ printf("%d %d \n", ib + 1, nnl[ib]);
1059+ }
1060+ printf("\n");
1061+ /*
1062+ Allocation of nodeline
1063+ */
1064+ kvnl = (GLfloat****)malloc(nb * sizeof(GLfloat***));
1065+ for (ib = 0; ib < nb; ++ib) {
1066+ kvnl[ib] = (GLfloat***)malloc(nnl[ib] * sizeof(GLfloat**));
1067+ for (i = 0; i < nnl[ib]; ++i){
1068+ kvnl[ib][i] = (GLfloat**)malloc(2 * sizeof(GLfloat*));
1069+ for (j = 0; j < 2; ++j){
1070+ kvnl[ib][i][j] = (GLfloat*)malloc(3 * sizeof(GLfloat));
1071+ }
1072+ }
1073+ }
1074+ /**/
1075+#pragma omp parallel default(none) \
1076+ shared(nb,nnl,matp,kvnl,kvp,ntri) \
1077+ private(ib,itri,mprod,i)
1078+ {
1079+#pragma omp for
1080+ for(ib = 0; ib < nb; ib++){
1081+ nnl[ib] = 0;
1082+ for(itri = 0; itri < ntri[ib]; ++itri){
1083+ /**/
1084+ mprod[0] = matp[ib][itri][0] * matp[ib][itri][1];
1085+ mprod[1] = matp[ib][itri][1] * matp[ib][itri][2];
1086+ /**/
1087+ if( fabsf(matp[ib][itri][0]) < 0.00001 && fabsf(matp[ib][itri][1]) < 0.00001){
1088+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][0][i] = kvp[ib][itri][0][i];
1089+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][1][i] = kvp[ib][itri][1][i];
1090+ nnl[ib] = nnl[ib] + 1;
1091+ }
1092+ else if(fabsf(matp[ib][itri][0]) < 0.00001 && fabsf(matp[ib][itri][2]) < 0.00001){
1093+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][0][i] = kvp[ib][itri][0][i];
1094+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][1][i] = kvp[ib][itri][2][i];
1095+ nnl[ib] = nnl[ib] + 1;
1096+ }
1097+ else if(fabsf(matp[ib][itri][1]) < 0.00001 && fabsf(matp[ib][itri][2]) < 0.00001){
1098+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][0][i] = kvp[ib][itri][1][i];
1099+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][1][i] = kvp[ib][itri][2][i];
1100+ nnl[ib] = nnl[ib] + 1;
1101+ }
1102+ else if(mprod[0] < 0.0){
1103+ if(mprod[1] < 0.0){
1104+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][0][i] = 1.0 / (matp[ib][itri][0] - matp[ib][itri][1])
1105+ * (( 0.0 - matp[ib][itri][1]) * kvp[ib][itri][0][i]
1106+ + (matp[ib][itri][0] - 0.0) * kvp[ib][itri][1][i]);
1107+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][1][i] = 1.0 / (matp[ib][itri][1] - matp[ib][itri][2])
1108+ * (( 0.0 - matp[ib][itri][2]) * kvp[ib][itri][1][i]
1109+ + (matp[ib][itri][1] - 0.0) * kvp[ib][itri][2][i]);
1110+ nnl[ib] = nnl[ib] + 1;
1111+ }
1112+ else{
1113+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][0][i] = 1.0 / (matp[ib][itri][0] - matp[ib][itri][1])
1114+ * (( 0.0 - matp[ib][itri][1]) * kvp[ib][itri][0][i]
1115+ + (matp[ib][itri][0] - 0.0) * kvp[ib][itri][1][i]);
1116+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][1][i] = 1.0 / (matp[ib][itri][0] - matp[ib][itri][2])
1117+ * (( 0.0 - matp[ib][itri][2]) * kvp[ib][itri][0][i]
1118+ + (matp[ib][itri][0] - 0.0) * kvp[ib][itri][2][i]);
1119+ nnl[ib] = nnl[ib] + 1;
1120+ }
1121+ }
1122+ else if(mprod[1] < 0.0){
1123+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][0][i] = 1.0 / (matp[ib][itri][0] - matp[ib][itri][2])
1124+ * (( 0.0 - matp[ib][itri][2]) * kvp[ib][itri][0][i]
1125+ + (matp[ib][itri][0] - 0.0) * kvp[ib][itri][2][i]);
1126+ for(i=0;i<3;++i)kvnl[ib][nnl[ib]][1][i] = 1.0 / (matp[ib][itri][1] - matp[ib][itri][2])
1127+ * (( 0.0 - matp[ib][itri][2]) * kvp[ib][itri][1][i]
1128+ + (matp[ib][itri][1] - 0.0) * kvp[ib][itri][2][i]);
1129+ nnl[ib] = nnl[ib] + 1;
1130+ }
1131+ }
1132+ }
1133+ } /** End of parallel region **/
1134+}
1135+/*
1136+ Draw Fermi surfaces
1137+*/
1138+void draw_fermi(){
1139+ /**/
1140+ int i, j, ib, itri;
1141+ GLfloat vect[3];
1142+ /*
1143+ Triagle patchs
1144+ */
1145+ for(ib = 0;ib < nb; ib++){
1146+ /**/
1147+ if(draw_band[ib] == 1){
1148+#pragma omp parallel default(none) \
1149+ shared(ib,ntri,rot,trans,nmlp,kvp,clr) \
1150+ private(itri,i,j,vect)
1151+ glBegin(GL_TRIANGLES);
1152+ {
1153+#pragma omp for
1154+ for(itri = 0; itri < ntri[ib]; ++itri){
1155+ /**/
1156+ for(j = 0; j < 3; ++j) vect[j] = rot[j][0] * nmlp[ib][itri][0]
1157+ + rot[j][1] * nmlp[ib][itri][1]
1158+ + rot[j][2] * nmlp[ib][itri][2];
1159+ glNormal3fv(vect);
1160+ /**/
1161+ for (i = 0; i < 3; ++i){
1162+ /**/
1163+ for(j = 0; j < 3; ++j) vect[j] = rot[j][0] * kvp[ib][itri][i][0]
1164+ + rot[j][1] * kvp[ib][itri][i][1]
1165+ + rot[j][2] * kvp[ib][itri][i][2]
1166+ + trans[j];
1167+ /**/
1168+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, clr[ib][itri][i]);
1169+ glVertex3fv(vect);
1170+ }
1171+ }
1172+ glEnd();
1173+ } /** End of parallel region **/
1174+ }
1175+ /**/
1176+ }
1177+ /*
1178+ Nodeline
1179+ */
1180+ if(nodeline == 1){
1181+ glLineWidth(10.0);
1182+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black);
1183+ glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
1184+ glBegin(GL_LINES);
1185+ for(ib = 0;ib < nb; ib++){
1186+ /**/
1187+ if(draw_band[ib] == 1){
1188+ for(itri = 0; itri < nnl[ib]; ++itri){
1189+ for (i = 0; i < 2; ++i){
1190+ /**/
1191+ for(j = 0; j < 3; ++j) vect[j] = rot[j][0] * kvnl[ib][itri][i][0]
1192+ + rot[j][1] * kvnl[ib][itri][i][1]
1193+ + rot[j][2] * kvnl[ib][itri][i][2]
1194+ + trans[j] + 0.001;
1195+ /**/
1196+ glVertex3fv(vect);
1197+ }
1198+ }
1199+ }
1200+ /**/
1201+ }
1202+ glEnd();
1203+ }
1204+ /**/
1205+} /* draw_ferm */
1206+/*
1207+ Draw lines of BZ boundaries
1208+*/
1209+void draw_bz_lines(){
1210+ /**/
1211+ int ibzl, i, j;
1212+ GLfloat bzl2[3], bvec2[3][3], linecolor[4];
1213+ /**/
1214+ if(blackback == 1){
1215+ for(i=0;i<4;i++) linecolor[i] = white[i];
1216+ }
1217+ else{
1218+ for(i=0;i<4;i++) linecolor[i] = black[i];
1219+ }
1220+ /**/
1221+ glLineWidth(2.0);
1222+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, linecolor);
1223+ glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, linecolor);
1224+ /**/
1225+ if(fbz == 1){
1226+ /**/
1227+ glBegin(GL_LINES);
1228+ for(ibzl = 0; ibzl < nbzl; ++ibzl){
1229+ for(i = 0; i< 2; ++i){
1230+ for(j = 0; j < 3; ++j) bzl2[j] = rot[j][0] * bzl[ibzl][i][0]
1231+ + rot[j][1] * bzl[ibzl][i][1]
1232+ + rot[j][2] * bzl[ibzl][i][2]
1233+ + trans[j];
1234+ glVertex3fv(bzl2);
1235+ }
1236+ }
1237+ glEnd();
1238+ /**/
1239+ }
1240+ else{
1241+ /**/
1242+ for(i = 0; i < 3; ++i){
1243+ for(j = 0; j < 3; ++j){
1244+ bvec2[i][j] = rot[j][0] * bvec[i][0]
1245+ + rot[j][1] * bvec[i][1]
1246+ + rot[j][2] * bvec[i][2];
1247+ }
1248+ }
1249+ /**/
1250+ glBegin(GL_LINE_STRIP);
1251+ for(i=0;i<3;++i) bzl2[i] = trans[i] ; glVertex3fv(bzl2);
1252+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[0][i] ; glVertex3fv(bzl2);
1253+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[1][i] ; glVertex3fv(bzl2);
1254+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[1][i] + bvec2[2][i]; glVertex3fv(bzl2);
1255+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[1][i] + bvec2[2][i]; glVertex3fv(bzl2);
1256+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[1][i] ; glVertex3fv(bzl2);
1257+ for(i=0;i<3;++i) bzl2[i] = trans[i] ; glVertex3fv(bzl2);
1258+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[2][i]; glVertex3fv(bzl2);
1259+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[2][i]; glVertex3fv(bzl2);
1260+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[1][i] + bvec2[2][i]; glVertex3fv(bzl2);
1261+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[2][i]; glVertex3fv(bzl2);
1262+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[0][i] ; glVertex3fv(bzl2);
1263+ for(i=0;i<3;++i) bzl2[i] = trans[i] ; glVertex3fv(bzl2);
1264+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[2][i]; glVertex3fv(bzl2);
1265+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[1][i] + bvec2[2][i]; glVertex3fv(bzl2);
1266+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[1][i] ; glVertex3fv(bzl2);
1267+ for(i=0;i<3;++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[1][i] ; glVertex3fv(bzl2);
1268+ glEnd();
1269+ /**/
1270+ }
1271+ /**/
1272+ glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
1273+ /**/
1274+} /* draw bz_lines */
1275+/*
1276+ Draw color scale
1277+*/
1278+void draw_colorbar()
1279+{
1280+ int i, j;
1281+ GLfloat mat2, barcolor[4];
1282+ GLfloat colorbar[13][3] = {
1283+ { 0.0, 0.0, 1.0 },
1284+ { -1.0, -1.0, 0.0 },
1285+ { -1.0, -1.0 - 0.1, 0.0 },
1286+ { -0.5, -1.0, 0.0 },
1287+ { -0.5, -1.0 - 0.1, 0.0 },
1288+ { 0.0, -1.0, 0.0 },
1289+ { 0.0, -1.0 - 0.1, 0.0 },
1290+ { 0.5, -1.0, 0.0 },
1291+ { 0.5, -1.0 - 0.1, 0.0 },
1292+ { 1.0, -1.0, 0.0 },
1293+ { 1.0, -1.0 - 0.1, 0.0 },
1294+ { 0.0, -1.0, 0.00001 },
1295+ { 0.0, -1.0 - 0.1, 0.00001 }
1296+ };
1297+ /**/
1298+ if(fcscl == 1 || fcscl == 2){
1299+ glBegin(GL_TRIANGLE_STRIP);
1300+ glNormal3fv(colorbar[0]);
1301+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, blue);
1302+ glVertex3fv(colorbar[1]);
1303+ glVertex3fv(colorbar[2]);
1304+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, cyan);
1305+ glVertex3fv(colorbar[3]);
1306+ glVertex3fv(colorbar[4]);
1307+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, green);
1308+ glVertex3fv(colorbar[5]);
1309+ glVertex3fv(colorbar[6]);
1310+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, yellow);
1311+ glVertex3fv(colorbar[7]);
1312+ glVertex3fv(colorbar[8]);
1313+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
1314+ glVertex3fv(colorbar[9]);
1315+ glVertex3fv(colorbar[10]);
1316+ glEnd();
1317+ /**/
1318+ /* if(nodeline == 1){ */
1319+ /* glLineWidth(2.0); */
1320+ /* glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black); */
1321+ /* glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black); */
1322+ /* glBegin(GL_LINES); */
1323+ /* glVertex3fv(colorbar[11]); */
1324+ /* glVertex3fv(colorbar[12]); */
1325+ /* glEnd(); */
1326+ /* } */
1327+ }
1328+ else if(fcscl == 4){
1329+ for(i = 0; i <= 60; i ++){
1330+ /**/
1331+ mat2 = (GLfloat)i / 60.0 * 6.0;
1332+ /**/
1333+ if(mat2 <= 1.0) {
1334+ for(j=0;j<4;++j) barcolor[j] = yellow[j] * mat2 + red[j] * (1.0 - mat2);
1335+ }
1336+ else if(mat2 <= 2.0){
1337+ mat2 = mat2 - 1.0;
1338+ for(j=0;j<4;++j) barcolor[j] = green[j] * mat2 + yellow[j] * (1.0 - mat2);
1339+ }
1340+ else if(mat2 <= 3.0){
1341+ mat2 = mat2 - 2.0;
1342+ for(j=0;j<4;++j) barcolor[j] = cyan[j] * mat2 + green[j] * (1.0 - mat2);
1343+ }
1344+ else if(mat2 <= 4.0){
1345+ mat2 = mat2 - 3.0;
1346+ for(j=0;j<4;++j) barcolor[j] = blue[j] * mat2 + cyan[j] * (1.0 - mat2);
1347+ }
1348+ else if(mat2 <= 5.0){
1349+ mat2 = mat2 - 4.0;
1350+ for(j=0;j<4;++j) barcolor[j] = magenta[j] * mat2 + blue[j] * (1.0 - mat2);
1351+ }
1352+ else{
1353+ mat2 = mat2 - 5.0;
1354+ for(j=0;j<4;++j) barcolor[j] = red[j] * mat2 + magenta[j] * (1.0 - mat2);
1355+ }
1356+ /**/
1357+ glBegin(GL_TRIANGLES);
1358+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, barcolor);
1359+ glNormal3fv(colorbar[0]);
1360+ glVertex3f(0.15 * cos((GLfloat)(i+1) / 60.0 * 6.283185307),
1361+ 0.15 * sin((GLfloat)(i+1) / 60.0 * 6.283185307) - 1.0, 0.0);
1362+ glVertex3f(0.15 * cos((GLfloat)i / 60.0 * 6.283185307),
1363+ 0.15 * sin((GLfloat)i / 60.0 * 6.283185307) - 1.0, 0.0);
1364+ glVertex3f(0.2 * cos((GLfloat)i / 60.0 * 6.283185307),
1365+ 0.2 * sin((GLfloat)i / 60.0 * 6.283185307) - 1.0, 0.0);
1366+ glVertex3f(0.2 * cos((GLfloat)i / 60.0 * 6.283185307),
1367+ 0.2 * sin((GLfloat)i / 60.0 * 6.283185307) - 1.0, 0.0);
1368+ glVertex3f(0.2 * cos((GLfloat)(i+1) / 60.0 * 6.283185307),
1369+ 0.2 * sin((GLfloat)(i+1) / 60.0 * 6.283185307) - 1.0, 0.0);
1370+ glVertex3f(0.15 * cos((GLfloat)(i+1) / 60.0 * 6.283185307),
1371+ 0.15 * sin((GLfloat)(i+1) / 60.0 * 6.283185307) - 1.0, 0.0);
1372+ glEnd();
1373+ }
1374+ }
1375+ else{
1376+ }
1377+} /* draw_colorbar */
1378+/*
1379+Draw points for the stereogram
1380+*/
1381+void draw_circles(){
1382+ int i;
1383+ GLfloat r;
1384+ /**/
1385+ r = 0.05;
1386+ /**/
1387+ if(blackback == 1){
1388+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, white);
1389+ glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, white);
1390+ }
1391+ else{
1392+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black);
1393+ glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
1394+ }
1395+ /**/
1396+ glBegin(GL_TRIANGLE_FAN);
1397+ glNormal3f(0.0, 0.0, 1.0);
1398+ glVertex3f(0.7, scl, 0.0);
1399+ for(i = 0; i <= 20; i ++){
1400+ glVertex3f(r * cos((GLfloat)i / 20.0 * 6.283185307) + 0.7,
1401+ r * sin((GLfloat)i / 20.0 * 6.283185307) + scl, 0.0);
1402+ }
1403+ glEnd();
1404+ /**/
1405+ glBegin(GL_TRIANGLE_FAN);
1406+ glNormal3f(0.0, 0.0, 1.0);
1407+ glVertex3f(-0.7, scl, 0.0);
1408+ for(i = 0; i <= 20; i ++){
1409+ glVertex3f(r * cos((GLfloat)i / 20.0 * 6.283185307) - 0.7,
1410+ r * sin((GLfloat)i / 20.0 * 6.283185307) + scl, 0.0);
1411+ }
1412+ glEnd();
1413+}
1414+/*
1415+ Glut Display function
1416+*/
1417+void display()
1418+{
1419+ GLfloat pos[] = { 1.0, 1.0, 1.0, 0.0 };
1420+ GLfloat amb[] = {0.2, 0.2, 0.2, 0.0};
1421+ double dx, theta, posz, phi;
1422+ GLfloat pos1[4], pos2[4];
1423+ /**/
1424+ if(lstereo == 2){
1425+ theta = 3.1416 /180.0 * 2.50;
1426+ posz = 5.0;
1427+ dx = 0.7;
1428+ phi = atan(posz/dx) -theta;
1429+ theta = (3.1416 * 0.50 - phi) / 3.1416 * 180.0;
1430+ /**/
1431+ pos1[0] = posz * cos(phi) - dx;
1432+ pos1[1] = 0.0;
1433+ pos1[2] = posz * sin(phi);
1434+ pos1[3] = 1.0;
1435+ /**/
1436+ pos2[0] = - posz * cos(phi) + dx;
1437+ pos2[1] = 0.0;
1438+ pos2[2] = posz * sin(phi);
1439+ pos2[3] = 1.0;
1440+ }
1441+ else if(lstereo == 3){
1442+ theta = -3.1416 /180.0 * 2.0;
1443+ posz = 5.0;
1444+ dx = - 0.7;
1445+ /**/
1446+ pos1[0] = posz * sin(theta) - dx;
1447+ pos1[1] = 0.0;
1448+ pos1[2] = posz * cos(theta);
1449+ pos1[3] = 1.0;
1450+ /**/
1451+ pos2[0] = - posz * sin(theta) + dx;
1452+ pos2[1] = 0.0;
1453+ pos2[2] = posz * cos(theta);
1454+ pos2[3] = 1.0;
1455+ }
1456+ else{
1457+ theta = 0.0;
1458+ dx = 0.0;
1459+ }
1460+ /*
1461+ Initialize
1462+ */
1463+ glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
1464+ glLoadIdentity();
1465+ /*
1466+ Set view point & view line
1467+ */
1468+ gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
1469+ /*
1470+ Set position of light
1471+ */
1472+ if(lstereo == 1){
1473+ glLightfv(GL_LIGHT0, GL_POSITION, pos);
1474+ /*
1475+ Draw color scale
1476+ */
1477+ if(lcolorbar == 1) draw_colorbar();
1478+ }
1479+ else{
1480+ glLightfv(GL_LIGHT0, GL_POSITION, pos1);
1481+ draw_circles();
1482+ glTranslated(-dx, 0.0, 0.0);
1483+ glRotated(theta, 0.0, 1.0, 0.0);
1484+ }
1485+ glLightfv(GL_LIGHT1, GL_AMBIENT, amb);
1486+ /*
1487+ Rotation & Zoom
1488+ */
1489+ glScaled(scl, scl, scl);
1490+ /*
1491+ Draw Brillouin zone boundaries
1492+ */
1493+ draw_bz_lines();
1494+ /*
1495+ Draw Fermi surfaces
1496+ */
1497+ draw_fermi();
1498+ /*
1499+ Draw the second
1500+ */
1501+ if(lstereo != 1){
1502+ glPushMatrix();
1503+ glLoadIdentity();
1504+ gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
1505+ glLightfv(GL_LIGHT0, GL_POSITION, pos2);
1506+ /**/
1507+ glTranslated(dx, 0.0, 0.0);
1508+ glRotated(-theta, 0.0, 1.0, 0.0);
1509+ /**/
1510+ glScaled(scl, scl, scl);
1511+ draw_bz_lines();
1512+ draw_fermi();
1513+ /**/
1514+ glPopMatrix();
1515+ }
1516+ /**/
1517+ glutSwapBuffers();
1518+} /* display */
1519+/*
1520+ Window resization
1521+*/
1522+void resize(int w, int h)
1523+{
1524+ /*
1525+ Scale of translation of mousepointer
1526+ */
1527+ sx = 1.0 / (GLfloat)w;
1528+ sy = 1.0 / (GLfloat)h;
1529+ /**/
1530+ glViewport(0, 0, w, h);
1531+ /**/
1532+ glMatrixMode(GL_PROJECTION);
1533+ /**/
1534+ glLoadIdentity();
1535+ gluPerspective(30.0, (GLfloat)w / (GLfloat)h, 1.0, 100.0);
1536+ /**/
1537+ glMatrixMode(GL_MODELVIEW);
1538+} /* end resize */
1539+/*
1540+ Idling
1541+*/
1542+void idle(void)
1543+{
1544+ glutPostRedisplay();
1545+} /* idle */
1546+/*
1547+ Glut mouse function
1548+*/
1549+void mouse(int button, int state, int x, int y)
1550+{
1551+ switch (button) {
1552+ /*
1553+ Drag
1554+ */
1555+ case GLUT_LEFT_BUTTON:
1556+ switch (state) {
1557+ case GLUT_DOWN:
1558+ /* Start of animation */
1559+ glutIdleFunc(idle);
1560+ /* Record drag start point */
1561+ cx = x;
1562+ cy = y;
1563+ break;
1564+ case GLUT_UP:
1565+ /* End of animation */
1566+ glutIdleFunc(0);
1567+ break;
1568+ default:
1569+ break;
1570+ }
1571+ break;
1572+ /*
1573+ Zoom up
1574+ */
1575+ case MOUSE_SCROLL_UP:
1576+ switch (state) {
1577+ case GLUT_DOWN:
1578+ scl = scl * 1.1;
1579+ glutPostRedisplay();
1580+ break;
1581+ case GLUT_UP:
1582+ break;
1583+ default:
1584+ break;
1585+ }
1586+ break;
1587+ /*
1588+ Zoom down
1589+ */
1590+ case MOUSE_SCROLL_DOWN:
1591+ switch (state) {
1592+ case GLUT_DOWN:
1593+ scl = scl * 0.9;
1594+ glutPostRedisplay();
1595+ break;
1596+ case GLUT_UP:
1597+ break;
1598+ default:
1599+ break;
1600+ }
1601+ break;
1602+ /*
1603+ No pushing
1604+ */
1605+ default:
1606+ break;
1607+ }
1608+} /* end mouse */
1609+/*
1610+ Glut motion function
1611+*/
1612+void motion(int x, int y)
1613+{
1614+ int i, j;
1615+ GLfloat dx, dy, a, rot0[3][3], rot1[3][3], ax, ay;
1616+ /*
1617+ Translation of mousepointer from starting point
1618+ */
1619+ dx = (x - cx) * sx;
1620+ dy = (y - cy) * sy;
1621+ /**/
1622+ /*
1623+ Distanse from starting point
1624+ */
1625+ a = sqrt(dx * dx + dy * dy);
1626+ /**/
1627+ if (a != 0.0) {
1628+ /*
1629+ Compute rotational matrix from translation of mousepointer
1630+ */
1631+ ax = - dy;
1632+ ay = dx;
1633+ /**/
1634+ a = a * 10.0;
1635+ /**/
1636+ rot0[0][0] = (ax * ax + ay * ay * cos(a)) / (ax * ax + ay * ay);
1637+ rot0[0][1] = ax * ay * (cos(a) - 1.0) / (ax * ax + ay * ay);
1638+ rot0[0][2] = ay * sin(a) / sqrtf(ax * ax + ay * ay);
1639+ rot0[1][0] = ax * ay * (cos(a) - 1.0) / (ax * ax + ay * ay);
1640+ rot0[1][1] = (ax * ax * cos(a) + ay * ay) / (ax * ax + ay * ay);
1641+ rot0[1][2] = ax * sin(a) / sqrtf(ax * ax + ay * ay);
1642+ rot0[2][0] = - ay * sin(a) / sqrtf(ax * ax + ay * ay);
1643+ rot0[2][1] = - ax * sin(a) / sqrtf(ax * ax + ay * ay);
1644+ rot0[2][2] = cos(a);
1645+ /**/
1646+ for(i = 0; i < 3; i++) for(j = 0; j < 3; j++) rot1[i][j] = rot[i][j];
1647+ /**/
1648+ for(i = 0; i < 3; i++){
1649+ for(j = 0; j < 3; j++){
1650+ rot[i][j] = rot0[i][0] * rot1[0][j]
1651+ + rot0[i][1] * rot1[1][j]
1652+ + rot0[i][2] * rot1[2][j];
1653+ }
1654+ }
1655+ }
1656+ cx = x;
1657+ cy = y;
1658+ /**/
1659+} /* motion */
1660+/*
1661+ Glut keyboard function
1662+*/
1663+void keyboard(unsigned char key, int x, int y)
1664+{
1665+ switch (key) {
1666+ }
1667+} /* keyboard */
1668+/*
1669+ Glut special key function
1670+ */
1671+void special_key(int key, int x, int y)
1672+{
1673+ switch (key) {
1674+ case GLUT_KEY_UP:
1675+ trans[1] = trans[1] + 0.1;
1676+ glutPostRedisplay();
1677+ break;
1678+ case GLUT_KEY_DOWN:
1679+ trans[1] = trans[1] - 0.1;
1680+ glutPostRedisplay();
1681+ break;
1682+ case GLUT_KEY_RIGHT:
1683+ /**/
1684+ trans[0] = trans[0] + 0.1;
1685+ glutPostRedisplay();
1686+ break;
1687+ /**/
1688+ case GLUT_KEY_LEFT:
1689+ /**/
1690+ trans[0] = trans[0] - 0.1;
1691+ glutPostRedisplay();
1692+ break;
1693+ /**/
1694+ }
1695+} /* special_key */
1696+/*
1697+ Main menu
1698+*/
1699+void main_menu(int value){
1700+ /**/
1701+ int ib, i1, i2, i3, ierr;
1702+ GLfloat emin, emax;
1703+ /**/
1704+ if(value == 2){
1705+ /*
1706+ Shift Fermi energy
1707+ */
1708+ emin = 100000.0;
1709+ emax = -100000.0;
1710+ for (ib = 0; ib < nb; ++ib) {
1711+ for (i1 = 0; i1 < ng[0]; ++i1) {
1712+ for (i2 = 0; i2 < ng[1]; ++i2) {
1713+ for (i3 = 0; i3 < ng[2]; ++i3) {
1714+ eig[ib][i1][i2][i3] = eig[ib][i1][i2][i3] + def;
1715+ if(emin > eig[ib][i1][i2][i3]) emin = eig[ib][i1][i2][i3];
1716+ if(emax < eig[ib][i1][i2][i3]) emax = eig[ib][i1][i2][i3];
1717+ }
1718+ }
1719+ }
1720+ }
1721+ printf("Min Max E_F \n");
1722+ printf("%f %f %f \n",emin, emax, def);
1723+ printf("New Fermi energy : ");
1724+ ierr = scanf("%f", &def);
1725+ if(ierr == 0) printf("error ! reading ef");
1726+ for (ib = 0; ib < nb; ++ib) {
1727+ for (i1 = 0; i1 < ng[0]; ++i1) {
1728+ for (i2 = 0; i2 < ng[1]; ++i2) {
1729+ for (i3 = 0; i3 < ng[2]; ++i3) {
1730+ eig[ib][i1][i2][i3] = eig[ib][i1][i2][i3] - def;
1731+ }
1732+ }
1733+ }
1734+ }
1735+ query = 1;
1736+ fermi_patch();
1737+ query = 0;
1738+ fermi_patch();
1739+ calc_nodeline();
1740+ /**/
1741+ max_and_min();
1742+ /**/
1743+ glutPostRedisplay();
1744+ }
1745+ else if(value == 9){
1746+ /*
1747+ Exit
1748+ */
1749+ printf("\nExit. \n\n");
1750+ free(eig);
1751+ free(mat);
1752+ free(ntri);
1753+ free(draw_band);
1754+ free(nmlp);
1755+ free(kvp);
1756+ free(matp);
1757+ free(clr);
1758+ free(bzl);
1759+ free(nnl);
1760+ free(kvnl);
1761+ exit(0);
1762+ }
1763+}
1764+/*
1765+ On / Off band
1766+*/
1767+void menu_band(int value){
1768+ /**/
1769+ if(draw_band[value] == 0){
1770+ printf("band # %d : On \n", value + 1);
1771+ draw_band[value] = 1;
1772+ }
1773+ else{
1774+ printf("band # %d : Off \n", value + 1);
1775+ draw_band[value] = 0;
1776+ }
1777+ glutPostRedisplay();
1778+ /**/
1779+} /* menu_band */
1780+/*
1781+ Change background color
1782+*/
1783+void menu_bgcolor(int value){
1784+ if(value == 1 && blackback != 1){
1785+ printf("Background color becomes BLACK. \n");
1786+ glClearColor(0.0, 0.0, 0.0, 0.0);
1787+ blackback = 1;
1788+ glutPostRedisplay();
1789+ }
1790+ else if(value == 2 && blackback != 0){
1791+ printf("Background color becomes WHITE. \n");
1792+ glClearColor(1.0, 1.0, 1.0, 0.0);
1793+ blackback = 0;
1794+ }
1795+/**/
1796+}/* bgcolor change*/
1797+/*
1798+ Change color scale mode
1799+*/
1800+void menu_colorscale(int value){
1801+ /**/
1802+ if(value == 1 && fcscl != 1){
1803+ fcscl = 1;
1804+ printf("##### Automatic color scale mode ##### \n");
1805+ max_and_min();
1806+ glutPostRedisplay();
1807+ }
1808+ else if(value == 2 && fcscl != 2){
1809+ fcscl = 2;
1810+ printf("##### Manual color scale mode ##### \n");
1811+ max_and_min();
1812+ glutPostRedisplay();
1813+ }
1814+ else if(value == 3 && fcscl != 3){
1815+ fcscl = 3;
1816+ printf("##### Unicolor scale mode ##### \n");
1817+ max_and_min();
1818+ glutPostRedisplay();
1819+ }
1820+ else if(value == 4 && fcscl != 4){
1821+ fcscl = 4;
1822+ printf("##### Periodic color scale mode ##### \n");
1823+ max_and_min();
1824+ glutPostRedisplay();
1825+ }
1826+ /**/
1827+} /* menu_colorscale */
1828+/*
1829+ Change Brillouin zone
1830+*/
1831+void menu_bzmode(int value){
1832+ if(value == 1 && fbz != 1){
1833+ fbz = 1;
1834+ printf("##### First Brillouin zone mode ##### \n\n");
1835+ free(kvp);
1836+ free(nmlp);
1837+ free(clr);
1838+ free(matp);
1839+ free(kvnl);
1840+ /**/
1841+ query = 1;
1842+ fermi_patch();
1843+ query = 0;
1844+ fermi_patch();
1845+ calc_nodeline();
1846+ /**/
1847+ max_and_min();
1848+ /**/
1849+ glutPostRedisplay();
1850+ }
1851+ else if(value == 2 && fbz != -1){
1852+ fbz = -1;
1853+ printf("##### Primitive Brillouin zone mode ##### \n\n");
1854+ free(kvp);
1855+ free(nmlp);
1856+ free(clr);
1857+ free(matp);
1858+ free(kvnl);
1859+ /**/
1860+ query = 1;
1861+ fermi_patch();
1862+ query = 0;
1863+ fermi_patch();
1864+ calc_nodeline();
1865+ /**/
1866+ max_and_min();
1867+ /**/
1868+ glutPostRedisplay();
1869+ }
1870+} /* menu_bzmode */
1871+/*
1872+ On/Off Node line
1873+*/
1874+void menu_nodeline(int value){
1875+ if(value == 1 && nodeline != 1){
1876+ printf("Nodeline : on \n\n");
1877+ nodeline = 1;
1878+ glutPostRedisplay();
1879+ }
1880+ else if(value == 2 && nodeline != 0){
1881+ printf("Nodeline : off \n\n");
1882+ nodeline = 0;
1883+ glutPostRedisplay();
1884+ }
1885+ /**/
1886+} /* menu_nodeline */
1887+/*
1888+ Tern stereogram
1889+*/
1890+void menu_stereo(int value){
1891+ if(value == 1 && lstereo != 1){
1892+ printf("Stereo : Off. \n\n");
1893+ lstereo = 1;
1894+ glutPostRedisplay();
1895+ }
1896+ if(value == 2 && lstereo != 2){
1897+ printf("Stereo : On. Parallel eyes \n\n");
1898+ lstereo = 2;
1899+ glutPostRedisplay();
1900+ }
1901+ if(value == 3 && lstereo != 3){
1902+ printf("Stereo : On. Crossed eyes \n\n");
1903+ lstereo = 3;
1904+ glutPostRedisplay();
1905+ }
1906+} /* menu_stereo */
1907+/*
1908+ On/Off Colorbar
1909+*/
1910+void menu_colorbar(int value){
1911+ if(value == 1 && lcolorbar != 1){
1912+ printf("Color bar : on \n\n");
1913+ lcolorbar = 1;
1914+ glutPostRedisplay();
1915+ }
1916+ if(value == 2 && lcolorbar != 0){
1917+ printf("Color bar : off \n\n");
1918+ lcolorbar = 0;
1919+ glutPostRedisplay();
1920+ }
1921+} /* menu_colorbar */
1922+/*
1923+ Change tetrahedron
1924+*/
1925+void menu_tetra(int value){
1926+ /**/
1927+ if(value != itet){
1928+ printf("Tetra patern %d \n", value + 1);
1929+ itet = value;
1930+ free(kvp);
1931+ free(nmlp);
1932+ free(clr);
1933+ free(matp);
1934+ free(kvnl);
1935+ init_corner();
1936+ query = 1;
1937+ fermi_patch();
1938+ query = 0;
1939+ fermi_patch();
1940+ calc_nodeline();
1941+ max_and_min();
1942+ glutPostRedisplay();
1943+ }
1944+} /* menu_tetra */
1945+/*
1946+ Glut init function
1947+*/
1948+void init(void)
1949+{
1950+ int ib;
1951+ int ibandmenu, ibgmenu, icsmenu, ibzmenu, inlmenu,
1952+ icbmenu, itetmenu, istereo;
1953+ char ibstr[20]={0};
1954+ /**/
1955+ glClearColor(0.0, 0.0, 0.0, 0.0);
1956+ glEnable(GL_DEPTH_TEST);
1957+ glEnable(GL_LIGHTING);
1958+ glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,GL_TRUE );
1959+ glEnable(GL_LIGHT0);
1960+ glEnable(GL_LIGHT1);
1961+ glEnable(GL_NORMALIZE);
1962+ /* Band menu */
1963+ ibandmenu = glutCreateMenu(menu_band);
1964+ for(ib =0;ib < nb; ib++){
1965+ sprintf(ibstr,"band # %d",ib + 1);
1966+ glutAddMenuEntry(ibstr, ib);
1967+ }
1968+ /* Background color */
1969+ ibgmenu = glutCreateMenu(menu_bgcolor);
1970+ glutAddMenuEntry("Black",1);
1971+ glutAddMenuEntry("White",2);
1972+ /* Color scale mode */
1973+ icsmenu = glutCreateMenu(menu_colorscale);
1974+ glutAddMenuEntry("Auto",1);
1975+ glutAddMenuEntry("Manual",2);
1976+ glutAddMenuEntry("Unicolor",3);
1977+ glutAddMenuEntry("Periodic",4);
1978+ /* Brillouin zone */
1979+ ibzmenu = glutCreateMenu(menu_bzmode);
1980+ glutAddMenuEntry("First Brillouin zone",1);
1981+ glutAddMenuEntry("Primitive Brillouin zone",2);
1982+ /* Nodeline on/off */
1983+ inlmenu = glutCreateMenu(menu_nodeline);
1984+ glutAddMenuEntry("On",1);
1985+ glutAddMenuEntry("Off",2);
1986+ /* Colorbar on/off */
1987+ icbmenu = glutCreateMenu(menu_colorbar);
1988+ glutAddMenuEntry("On",1);
1989+ glutAddMenuEntry("Off",2);
1990+ /* Stereogram */
1991+ istereo = glutCreateMenu(menu_stereo);
1992+ glutAddMenuEntry("None",1);
1993+ glutAddMenuEntry("Parallel",2);
1994+ glutAddMenuEntry("Cross",3);
1995+ /* Tetrahedron */
1996+ itetmenu = glutCreateMenu(menu_tetra);
1997+ for(ib =0;ib < 16; ib++){
1998+ sprintf(ibstr,"tetra # %d",ib + 1);
1999+ glutAddMenuEntry(ibstr,ib);
2000+ }
2001+ /*
2002+ Main menu
2003+ */
2004+ glutCreateMenu(main_menu);
2005+ glutAddSubMenu("Band",ibandmenu);
2006+ glutAddMenuEntry("Shift Fermi energy",2);
2007+ glutAddSubMenu("Background color",ibgmenu);
2008+ glutAddSubMenu("Color scale mode",icsmenu);
2009+ glutAddSubMenu("Brillouin zone",ibzmenu);
2010+ glutAddSubMenu("Node lines",inlmenu);
2011+ glutAddSubMenu("Color bar On/Off",icbmenu);
2012+ glutAddSubMenu("Stereogram",istereo);
2013+ glutAddSubMenu("Tetrahedron",itetmenu);
2014+ glutAddMenuEntry("Exit",9);
2015+ glutAttachMenu(GLUT_RIGHT_BUTTON);
2016+} /* init */
2017+/*
2018+ Main routine
2019+*/
2020+int main(int argc, char *argv[])
2021+{
2022+ /**/
2023+ read_file(argv[1]);
2024+ init_corner();
2025+ bragg_vector();
2026+ /**/
2027+ query = 1;
2028+ bz_lines();
2029+ query = 0;
2030+ bz_lines();
2031+ /**/
2032+ max_and_min_bz();
2033+ /**/
2034+#pragma omp parallel
2035+ printf("Threads : %d \n", omp_get_num_threads());
2036+ /**/
2037+ printf("\n##### First Brillouin zone mode ##### \n\n");
2038+ query = 1;
2039+ fermi_patch();
2040+ query = 0;
2041+ fermi_patch();
2042+ calc_nodeline();
2043+ /**/
2044+ printf("\n##### Full color scale ##### \n\n");
2045+ max_and_min();
2046+ /*
2047+ Description
2048+ */
2049+ printf("\n##### How to handle ##### \n\n");
2050+ printf(" mouse drag : Rotate objects \n");
2051+ printf(" mousewheel : Resize objects \n");
2052+ printf(" cursorkey : Move objects \n");
2053+ printf(" mouse right button : Menu \n");
2054+ printf("\n");
2055+ /**/
2056+ glutInit(&argc, argv);
2057+ glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH);
2058+ glutCreateWindow(argv[0]);
2059+ glutDisplayFunc(display);
2060+ glutReshapeFunc(resize);
2061+ glutMouseFunc(mouse);
2062+ glutMotionFunc(motion);
2063+ glutKeyboardFunc(keyboard);
2064+ glutSpecialFunc(special_key);
2065+ init();
2066+ glutMainLoop();
2067+ return 0;
2068+} /* main */
Show on old repository browser