OSDN > Finden Software > Multimedia > Graphics > fermisurfer > SCM > git > fermisurfer > Commit

fermisurfer
Fork
fermisurfer
android

R/O
HTTP
SSH
HTTPS

fermisurfer: Commit

fermisurfer Git

Commit MetaInfo

Revision	35f669eb33fc9b9858c3d07b80ff7aacbb094fc8 (tree)
Zeit	2017-04-14 15:24:58
Autor	mitsuaki1987 <kawamitsuaki@gmai...>
Commiter	mitsuaki1987

Log Message

Separate routines in each file
Add interpolation function

Ändern Zusammenfassung

modified: .gitignore (diff)
modified: doc/fermisurfer.doxygen (diff)
modified: src/Makefile (diff)
add: src/basic_math.c (diff)
add: src/basic_math.h (diff)
add: src/bz_lines.c (diff)
add: src/bz_lines.h (diff)
add: src/calc_nodeline.c (diff)
add: src/calc_nodeline.h (diff)
add: src/draw.c (diff)
add: src/draw.h (diff)
add: src/fermi_patch.c (diff)
add: src/fermi_patch.h (diff)
modified: src/fermisurfer.c (diff)
add: src/free_patch.c (diff)
add: src/free_patch.h (diff)
add: src/initialize.c (diff)
add: src/initialize.h (diff)
add: src/kumo.c (diff)
add: src/kumo.h (diff)
add: src/menu.c (diff)
add: src/menu.h (diff)
add: src/operation.c (diff)
add: src/operation.h (diff)
add: src/read_file.c (diff)
add: src/read_file.h (diff)
add: src/variable.h (diff)

Diff

--- a/.gitignore

+++ b/.gitignore

		@@ -28,3 +28,5 @@ src/bxsf2frmsf
28	28	src/resource.h
29	29	doc/ja/_build/
30	30	doc/en/_build/
	31	+*.o
	32	+doc/html/
		\ No newline at end of file

--- a/doc/fermisurfer.doxygen

+++ b/doc/fermisurfer.doxygen

		@@ -38,13 +38,13 @@ PROJECT_NAME = fermisurfer
38	38	# could be handy for archiving the generated documentation or if some version
39	39	# control system is used.
40	40
41		-PROJECT_NUMBER =
	41	+PROJECT_NUMBER = 1.8
42	42
43	43	# Using the PROJECT_BRIEF tag one can provide an optional one line description
44	44	# for a project that appears at the top of each page and should give viewer a
45	45	# quick idea about the purpose of the project. Keep the description short.
46	46
47		-PROJECT_BRIEF =
	47	+PROJECT_BRIEF = fermisurfer
48	48
49	49	# With the PROJECT_LOGO tag one can specify a logo or an icon that is included
50	50	# in the documentation. The maximum height of the logo should not exceed 55

		@@ -58,7 +58,7 @@ PROJECT_LOGO =
58	58	# entered, it will be relative to the location where doxygen was started. If
59	59	# left blank the current directory will be used.
60	60
61		-OUTPUT_DIRECTORY = C:/Users/kawamuura/xubuntu/docment/public_html/fermisurfer
	61	+OUTPUT_DIRECTORY = .
62	62
63	63	# If the CREATE_SUBDIRS tag is set to YES then doxygen will create 4096 sub-
64	64	# directories (in 2 levels) under the output directory of each output format and

		@@ -768,7 +768,7 @@ WARN_LOGFILE =
768	768	# spaces.
769	769	# Note: If this tag is empty the current directory is searched.
770	770
771		-INPUT = C:/Users/kawamuura/xubuntu/programs/fermisurfer/src
	771	+INPUT = ../src
772	772
773	773	# This tag can be used to specify the character encoding of the source files
774	774	# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses

--- a/src/Makefile

+++ b/src/Makefile

		@@ -1,16 +1,49 @@
1	1
2	2	CC = gcc
3	3
4		-CFLAGS=-O3 -lglut -lGLU -lGL -lm -g
	4	+CFLAGS=-O3 -lglut -lGLU -lGL -lm -g -fopenmp
5	5	# for mac
6	6	#CFLAGS=-D MAC -O3 -framework OpenGL -framework GLUT -lm
7	7
	8	+OBJS= \
	9	+basic_math.o \
	10	+bz_lines.o \
	11	+calc_nodeline.o \
	12	+draw.o \
	13	+fermisurfer.o \
	14	+fermi_patch.o \
	15	+free_patch.o \
	16	+initialize.o \
	17	+kumo.o \
	18	+menu.o \
	19	+operation.o \
	20	+read_file.o
	21	+
	22	+HEADERS= \
	23	+basic_math.h \
	24	+bz_lines.h \
	25	+calc_nodeline.h \
	26	+draw.h \
	27	+fermi_patch.h \
	28	+free_patch.h \
	29	+initialize.h \
	30	+kumo.h \
	31	+menu.h \
	32	+operation.h \
	33	+read_file.h \
	34	+variable.h
	35	+
8	36	all:fermisurfer bxsf2frmsf
9	37
10		-fermisurfer:fermisurfer.c
11		- $(CC) $< $(CFLAGS) -o $@
	38	+SUFFIXES: .o .c
	39	+
	40	+.c.o:
	41	+ $(CC) $(CFLAGS) -c $<
	42	+
	43	+fermisurfer:$(OBJS)
	44	+ $(CC) $(OBJS) $(CFLAGS) -o $@
12	45
13		-bxsf2frmsf:bxsf2frmsf.c
	46	+bxsf2frmsf:bxsf2frmsf.o
14	47	$(CC) $< $(CFLAGS) -o $@
15	48
16	49	install:

		@@ -20,5 +53,17 @@ uninstall:
20	53	sudo rm /usr/local/bin/fermisurfer
21	54
22	55	clean:
23		- rm -rf fermisurfer bxsf2frmsf
	56	+ rm -rf *.o fermisurfer bxsf2frmsf
24	57
	58	+basic_math.o:$(HEADERS)
	59	+bz_lines.o:$(HEADERS)
	60	+calc_nodeline.o:$(HEADERS)
	61	+draw.o:$(HEADERS)
	62	+fermisurfer.o:$(HEADERS)
	63	+fermi_patch.o:$(HEADERS)
	64	+free_patch.o:$(HEADERS)
	65	+initialize.o:$(HEADERS)
	66	+kumo.o:$(HEADERS)
	67	+menu.o:$(HEADERS)
	68	+operation.o:$(HEADERS)
	69	+read_file.o:$(HEADERS)

--- /dev/null

+++ b/src/basic_math.c

		@@ -0,0 +1,124 @@
	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	+#include <math.h>
	26	+#if defined(MAC)
	27	+#include <GLUT/glut.h>
	28	+#else
	29	+#include <GL/glut.h>
	30	+#endif
	31	+
	32	+/**
	33	+*Work as Modulo function of fortran
	34	+*/
	35	+int modulo(int i, int n) {
	36	+ int j;
	37	+ j = (i + 100 * n) % n;
	38	+ return j;
	39	+}/modulo(int i, int n)/
	40	+ /**
	41	+ * Solve linear system
	42	+ */
	43	+GLfloat solve3(
	44	+ GLfloat a[3][3] /*< [in] Matix/,
	45	+ GLfloat b[3] /*< [inout] Right hand side vector/)
	46	+{
	47	+ int i;
	48	+ GLfloat det, c[3];
	49	+ /**/
	50	+ det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
	51	+ + a[0][1] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
	52	+ + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
	53	+ /**/
	54	+ c[0] = b[0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
	55	+ + b[1] * (a[0][2] * a[2][1] - a[0][1] * a[2][2])
	56	+ + b[2] * (a[0][1] * a[1][2] - a[0][2] * a[1][1]);
	57	+ /**/
	58	+ c[1] = b[0] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
	59	+ + b[1] * (a[0][0] * a[2][2] - a[0][2] * a[2][0])
	60	+ + b[2] * (a[0][2] * a[1][0] - a[0][0] * a[1][2]);
	61	+ /**/
	62	+ c[2] = b[0] * (a[1][0] * a[2][1] - a[1][1] * a[2][0])
	63	+ + b[1] * (a[0][1] * a[2][0] - a[0][0] * a[2][1])
	64	+ + b[2] * (a[0][0] * a[1][1] - a[0][1] * a[1][0]);
	65	+ /**/
	66	+ for (i = 0; i<3; ++i) b[i] = c[i] / det;
	67	+ return det;
	68	+ /**/
	69	+}
	70	+/**
	71	+* Sort eigenvalues
	72	+*/
	73	+void eigsort(
	74	+ int n /*< [in] the number of components/,
	75	+ GLfloat* eig2 /*< [inout] the orbital energy/,
	76	+ GLfloat* mat2 /*< [inout] the matrix element/,
	77	+ GLfloat kvec2[][3] /*< [inout] k-vectors of corners/)
	78	+{
	79	+ int i, j, k;
	80	+ GLfloat tmp;
	81	+ /**/
	82	+ for (i = 0; i < n - 1; ++i) {
	83	+ for (j = i + 1; j < n; ++j) {
	84	+ if (eig2[j] < eig2[i]) {
	85	+ tmp = eig2[j];
	86	+ eig2[j] = eig2[i];
	87	+ eig2[i] = tmp;
	88	+ /**/
	89	+ tmp = mat2[j];
	90	+ mat2[j] = mat2[i];
	91	+ mat2[i] = tmp;
	92	+ /**/
	93	+ for (k = 0; k < 3; ++k) {
	94	+ tmp = kvec2[j][k];
	95	+ kvec2[j][k] = kvec2[i][k];
	96	+ kvec2[i][k] = tmp;
	97	+ }
	98	+ }
	99	+ }
	100	+ }
	101	+} /* eigsort */
	102	+ /**
	103	+ * Calculate normal vector
	104	+ */
	105	+void normal_vec(
	106	+ GLfloat in1[3] /*< [in] Corner 1/,
	107	+ GLfloat in2[3] /*< [in] Corner 2/,
	108	+ GLfloat in3[3] /*< [in] Corner 3/,
	109	+ GLfloat out[3] /*< [out] The normal vector/)
	110	+{
	111	+ int i;
	112	+ GLfloat norm;
	113	+ out[0] = in1[1] * in2[2] - in1[2] * in2[1]
	114	+ + in2[1] * in3[2] - in2[2] * in3[1]
	115	+ + in3[1] * in1[2] - in3[2] * in1[1];
	116	+ out[1] = in1[2] * in2[0] - in1[0] * in2[2]
	117	+ + in2[2] * in3[0] - in2[0] * in3[2]
	118	+ + in3[2] * in1[0] - in3[0] * in1[2];
	119	+ out[2] = in1[0] * in2[1] - in1[1] * in2[0]
	120	+ + in2[0] * in3[1] - in2[1] * in3[0]
	121	+ + in3[0] * in1[1] - in3[1] * in1[0];
	122	+ norm = sqrtf(out[0] * out[0] + out[1] * out[1] + out[2] * out[2]);
	123	+ for (i = 0; i<3; i++) out[i] = out[i] / norm;
	124	+} /* normal_vec */

--- /dev/null

+++ b/src/basic_math.h

		@@ -0,0 +1,44 @@
	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	+#if defined(MAC)
	26	+#include <GLUT/glut.h>
	27	+#else
	28	+#include <GL/glut.h>
	29	+#endif
	30	+
	31	+int modulo(int i, int n);
	32	+GLfloat solve3(
	33	+ GLfloat a[3][3] /*< [in] Matix/,
	34	+ GLfloat b[3] /*< [inout] Right hand side vector/);
	35	+void eigsort(
	36	+ int n /*< [in] the number of components/,
	37	+ GLfloat* eig2 /*< [inout] the orbital energy/,
	38	+ GLfloat* mat2 /*< [inout] the matrix element/,
	39	+ GLfloat kvec2[][3] /*< [inout] k-vectors of corners/);
	40	+void normal_vec(
	41	+ GLfloat in1[3] /*< [in] Corner 1/,
	42	+ GLfloat in2[3] /*< [in] Corner 2/,
	43	+ GLfloat in3[3] /*< [in] Corner 3/,
	44	+ GLfloat out[3] /*< [out] The normal vector/);
		\ No newline at end of file

--- /dev/null

+++ b/src/bz_lines.c

		@@ -0,0 +1,146 @@
	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	+#include <stdlib.h>
	26	+#include <stdio.h>
	27	+#include <math.h>
	28	+#include "variable.h"
	29	+#include "basic_math.h"
	30	+
	31	+#if defined(MAC)
	32	+#include <GLUT/glut.h>
	33	+#else
	34	+#include <GL/glut.h>
	35	+#endif
	36	+
	37	+/**
	38	+* Judge wheser this line is the edge of 1st BZ
	39	+*/
	40	+int bragg_vert(
	41	+ int ibr /*< [in] Index of a Bragg plane/,
	42	+ int jbr /*< [in] Index of a Bragg plane/,
	43	+ int nbr /*< [in] /,
	44	+ GLfloat vert[3] /*< [in] start point of line/,
	45	+ GLfloat vert2[3] /*< [in] end point of line/)
	46	+{
	47	+ int kbr, i, lbr, nbr0;
	48	+ GLfloat bmat[3][3], rhs[3], prod, thr = 0.0001, det;
	49	+ /**/
	50	+ nbr0 = nbr;
	51	+ /**/
	52	+ for (kbr = nbr0; kbr < 26; ++kbr) {
	53	+ /**/
	54	+ for (i = 0; i<3; ++i) bmat[0][i] = bragg[ibr][i];
	55	+ for (i = 0; i<3; ++i) bmat[1][i] = bragg[jbr][i];
	56	+ for (i = 0; i<3; ++i) bmat[2][i] = bragg[kbr][i];
	57	+ /**/
	58	+ rhs[0] = brnrm[ibr];
	59	+ rhs[1] = brnrm[jbr];
	60	+ rhs[2] = brnrm[kbr];
	61	+ /*
	62	+ if Bragg planes do not cross, roop next kbr
	63	+ */
	64	+ det = solve3(bmat, rhs);
	65	+ if (fabsf(det) < thr) continue;
	66	+ /*
	67	+ if vert0 = vert1, roop next kbr
	68	+ */
	69	+ prod = (vert2[0] - rhs[0]) * (vert2[0] - rhs[0])
	70	+ + (vert2[1] - rhs[1]) * (vert2[1] - rhs[1])
	71	+ + (vert2[2] - rhs[2]) * (vert2[2] - rhs[2]);
	72	+ if (prod < thr) continue;
	73	+ /*
	74	+ is corner really in 1st BZ ?
	75	+ */
	76	+ i = 0;
	77	+ for (lbr = 0; lbr < 26; ++lbr) {
	78	+ prod = bragg[lbr][0] * rhs[0]
	79	+ + bragg[lbr][1] * rhs[1]
	80	+ + bragg[lbr][2] * rhs[2];
	81	+ /**/
	82	+ if (prod > brnrm[lbr] + thr) {
	83	+ i = 1;
	84	+ break;
	85	+ }
	86	+ }
	87	+ if (i == 1) {
	88	+ }
	89	+ else {
	90	+ for (i = 0; i<3; ++i) vert[i] = rhs[i];
	91	+ return kbr + 1;
	92	+ }
	93	+ }
	94	+ /*
	95	+ this line is not a BZ boundary
	96	+ */
	97	+ return 0;
	98	+ /**/
	99	+}/* bragg_vert */
	100	+ /**
	101	+ * Compute Brillouin zone boundariy lines
	102	+ */
	103	+void bz_lines() {
	104	+ /**/
	105	+ int ibr, jbr, nbr, ibzl, i, j, lvert;
	106	+ GLfloat vert[2][3];
	107	+ /**/
	108	+ ibzl = 0;
	109	+ /**/
	110	+ for (ibr = 0; ibr < 26; ++ibr) {
	111	+ for (jbr = 0; jbr < 26; ++jbr) {
	112	+ /**/
	113	+ for (i = 0; i<3; ++i) vert[1][i] = 0.0;
	114	+ nbr = 0;
	115	+ lvert = bragg_vert(ibr, jbr, nbr, vert[0], vert[1]);
	116	+ if (lvert == 0) continue;
	117	+ nbr = lvert;
	118	+ /**/
	119	+ lvert = bragg_vert(ibr, jbr, nbr, vert[1], vert[0]);
	120	+ if (lvert == 0) continue;
	121	+ /**/
	122	+ if (query == 1) {
	123	+ nbzl = nbzl + 1;
	124	+ }
	125	+ else {
	126	+ for (i = 0; i<2; ++i) for (j = 0; j<3; ++j) bzl[ibzl][i][j] = vert[i][j];
	127	+ ibzl = ibzl + 1;
	128	+ }
	129	+ /**/
	130	+ }
	131	+ }
	132	+ /**/
	133	+ if (query == 1) {
	134	+ printf("# of lines for BZ : %d \n", nbzl);
	135	+ /**/
	136	+ bzl = (GLfloat**)malloc(nbzl sizeof(GLfloat*));
	137	+ for (ibzl = 0; ibzl < nbzl; ++ibzl) {
	138	+ bzl[ibzl] = (GLfloat*)malloc(2 sizeof(GLfloat*));
	139	+ for (i = 0; i < 2; ++i) {
	140	+ bzl[ibzl][i] = (GLfloat)malloc(3 sizeof(GLfloat));
	141	+ }
	142	+ }
	143	+ /**/
	144	+ }
	145	+ /**/
	146	+} /* bz_lines */

--- /dev/null

+++ b/src/bz_lines.h

		@@ -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/src/calc_nodeline.c

		@@ -0,0 +1,158 @@
	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	+#include <stdio.h>
	26	+#include <math.h>
	27	+#include "variable.h"
	28	+
	29	+#if defined(MAC)
	30	+#include <GLUT/glut.h>
	31	+#else
	32	+#include <GL/glut.h>
	33	+#endif
	34	+
	35	+/**
	36	+* Node line
	37	+*/
	38	+void calc_nodeline() {
	39	+ int ib, itri, i, j;
	40	+ GLfloat mprod[2];
	41	+ /*
	42	+ Query
	43	+ */
	44	+#pragma omp parallel default(none) \
	45	+ shared(nb,nnl,matp,ntri) \
	46	+ private(ib,itri,mprod)
	47	+ {
	48	+#pragma omp for
	49	+ for (ib = 0; ib < nb; ib++) {
	50	+ nnl[ib] = 0;
	51	+ for (itri = 0; itri < ntri[ib]; ++itri) {
	52	+ /**/
	53	+ mprod[0] = matp[ib][itri][0] * matp[ib][itri][1];
	54	+ mprod[1] = matp[ib][itri][1] * matp[ib][itri][2];
	55	+ /**/
	56	+ if (fabsf(matp[ib][itri][0]) < 0.00001 && fabsf(matp[ib][itri][1]) < 0.00001) {
	57	+ nnl[ib] = nnl[ib] + 1;
	58	+ }
	59	+ else if (fabsf(matp[ib][itri][0]) < 0.00001 && fabsf(matp[ib][itri][2]) < 0.00001) {
	60	+ nnl[ib] = nnl[ib] + 1;
	61	+ }
	62	+ else if (fabsf(matp[ib][itri][1]) < 0.00001 && fabsf(matp[ib][itri][2]) < 0.00001) {
	63	+ nnl[ib] = nnl[ib] + 1;
	64	+ }
	65	+ else if (mprod[0] < 0.0) {
	66	+ if (mprod[1] < 0.0) {
	67	+ nnl[ib] = nnl[ib] + 1;
	68	+ }
	69	+ else {
	70	+ nnl[ib] = nnl[ib] + 1;
	71	+ }
	72	+ }
	73	+ else if (mprod[1] < 0.0) {
	74	+ nnl[ib] = nnl[ib] + 1;
	75	+ }
	76	+ }
	77	+ }
	78	+ } /* End of parallel region */
	79	+ /**/
	80	+ printf("band # of nodeline \n");
	81	+ for (ib = 0; ib < nb; ib++) {
	82	+ printf("%d %d \n", ib + 1, nnl[ib]);
	83	+ }
	84	+ printf("\n");
	85	+ /*
	86	+ Allocation of nodeline
	87	+ */
	88	+ kvnl = (GLfloat***)malloc(nb sizeof(GLfloat***));
	89	+ for (ib = 0; ib < nb; ++ib) {
	90	+ kvnl[ib] = (GLfloat**)malloc(nnl[ib] sizeof(GLfloat**));
	91	+ for (i = 0; i < nnl[ib]; ++i) {
	92	+ kvnl[ib][i] = (GLfloat*)malloc(2 sizeof(GLfloat*));
	93	+ for (j = 0; j < 2; ++j) {
	94	+ kvnl[ib][i][j] = (GLfloat)malloc(3 sizeof(GLfloat));
	95	+ }
	96	+ }
	97	+ }
	98	+ /**/
	99	+#pragma omp parallel default(none) \
	100	+ shared(nb,nnl,matp,kvnl,kvp,ntri) \
	101	+ private(ib,itri,mprod,i)
	102	+ {
	103	+#pragma omp for
	104	+ for (ib = 0; ib < nb; ib++) {
	105	+ nnl[ib] = 0;
	106	+ for (itri = 0; itri < ntri[ib]; ++itri) {
	107	+ /**/
	108	+ mprod[0] = matp[ib][itri][0] * matp[ib][itri][1];
	109	+ mprod[1] = matp[ib][itri][1] * matp[ib][itri][2];
	110	+ /**/
	111	+ if (fabsf(matp[ib][itri][0]) < 0.00001 && fabsf(matp[ib][itri][1]) < 0.00001) {
	112	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][0][i] = kvp[ib][itri][0][i];
	113	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][1][i] = kvp[ib][itri][1][i];
	114	+ nnl[ib] = nnl[ib] + 1;
	115	+ }
	116	+ else if (fabsf(matp[ib][itri][0]) < 0.00001 && fabsf(matp[ib][itri][2]) < 0.00001) {
	117	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][0][i] = kvp[ib][itri][0][i];
	118	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][1][i] = kvp[ib][itri][2][i];
	119	+ nnl[ib] = nnl[ib] + 1;
	120	+ }
	121	+ else if (fabsf(matp[ib][itri][1]) < 0.00001 && fabsf(matp[ib][itri][2]) < 0.00001) {
	122	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][0][i] = kvp[ib][itri][1][i];
	123	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][1][i] = kvp[ib][itri][2][i];
	124	+ nnl[ib] = nnl[ib] + 1;
	125	+ }
	126	+ else if (mprod[0] < 0.0) {
	127	+ if (mprod[1] < 0.0) {
	128	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][0][i] = 1.0 / (matp[ib][itri][0] - matp[ib][itri][1])
	129	+ * ((0.0 - matp[ib][itri][1]) * kvp[ib][itri][0][i]
	130	+ + (matp[ib][itri][0] - 0.0) * kvp[ib][itri][1][i]);
	131	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][1][i] = 1.0 / (matp[ib][itri][1] - matp[ib][itri][2])
	132	+ * ((0.0 - matp[ib][itri][2]) * kvp[ib][itri][1][i]
	133	+ + (matp[ib][itri][1] - 0.0) * kvp[ib][itri][2][i]);
	134	+ nnl[ib] = nnl[ib] + 1;
	135	+ }
	136	+ else {
	137	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][0][i] = 1.0 / (matp[ib][itri][0] - matp[ib][itri][1])
	138	+ * ((0.0 - matp[ib][itri][1]) * kvp[ib][itri][0][i]
	139	+ + (matp[ib][itri][0] - 0.0) * kvp[ib][itri][1][i]);
	140	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][1][i] = 1.0 / (matp[ib][itri][0] - matp[ib][itri][2])
	141	+ * ((0.0 - matp[ib][itri][2]) * kvp[ib][itri][0][i]
	142	+ + (matp[ib][itri][0] - 0.0) * kvp[ib][itri][2][i]);
	143	+ nnl[ib] = nnl[ib] + 1;
	144	+ }
	145	+ }
	146	+ else if (mprod[1] < 0.0) {
	147	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][0][i] = 1.0 / (matp[ib][itri][0] - matp[ib][itri][2])
	148	+ * ((0.0 - matp[ib][itri][2]) * kvp[ib][itri][0][i]
	149	+ + (matp[ib][itri][0] - 0.0) * kvp[ib][itri][2][i]);
	150	+ for (i = 0; i<3; ++i)kvnl[ib][nnl[ib]][1][i] = 1.0 / (matp[ib][itri][1] - matp[ib][itri][2])
	151	+ * ((0.0 - matp[ib][itri][2]) * kvp[ib][itri][1][i]
	152	+ + (matp[ib][itri][1] - 0.0) * kvp[ib][itri][2][i]);
	153	+ nnl[ib] = nnl[ib] + 1;
	154	+ }
	155	+ }
	156	+ }
	157	+ } /* End of parallel region */
	158	+}

--- /dev/null

+++ b/src/calc_nodeline.h

		@@ -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/src/draw.c

		@@ -0,0 +1,415 @@
	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	+#include <math.h>
	26	+#include "variable.h"
	27	+
	28	+#if defined(MAC)
	29	+#include <GLUT/glut.h>
	30	+#else
	31	+#include <GL/glut.h>
	32	+#endif
	33	+
	34	+/**
	35	+* Draw Fermi surfaces
	36	+*/
	37	+void draw_fermi() {
	38	+ /**/
	39	+ int i, j, ib, itri;
	40	+ GLfloat vect[3];
	41	+ /*
	42	+ Triagle patchs
	43	+ */
	44	+ for (ib = 0; ib < nb; ib++) {
	45	+ /**/
	46	+ if (draw_band[ib] == 1) {
	47	+ glBegin(GL_TRIANGLES);
	48	+ for (itri = 0; itri < ntri[ib]; ++itri) {
	49	+ /**/
	50	+ for (j = 0; j < 3; ++j)
	51	+ vect[j] = rot[j][0] * nmlp[ib][itri][0]
	52	+ + rot[j][1] * nmlp[ib][itri][1]
	53	+ + rot[j][2] * nmlp[ib][itri][2];
	54	+ glNormal3fv(vect);
	55	+ /**/
	56	+ for (i = 0; i < 3; ++i) {
	57	+ /**/
	58	+ for (j = 0; j < 3; ++j)
	59	+ vect[j] = rot[j][0] * kvp[ib][itri][i][0]
	60	+ + rot[j][1] * kvp[ib][itri][i][1]
	61	+ + rot[j][2] * kvp[ib][itri][i][2]
	62	+ + trans[j];
	63	+ /**/
	64	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, clr[ib][itri][i]);
	65	+ glVertex3fv(vect);
	66	+ }
	67	+ }
	68	+ glEnd();
	69	+ }
	70	+ /**/
	71	+ }
	72	+ /*
	73	+ Nodeline
	74	+ */
	75	+ if (nodeline == 1) {
	76	+ glLineWidth(10.0);
	77	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black);
	78	+ glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
	79	+ glBegin(GL_LINES);
	80	+ for (ib = 0; ib < nb; ib++) {
	81	+ /**/
	82	+ if (draw_band[ib] == 1) {
	83	+ for (itri = 0; itri < nnl[ib]; ++itri) {
	84	+ for (i = 0; i < 2; ++i) {
	85	+ /**/
	86	+ for (j = 0; j < 3; ++j)
	87	+ vect[j] = rot[j][0] * kvnl[ib][itri][i][0]
	88	+ + rot[j][1] * kvnl[ib][itri][i][1]
	89	+ + rot[j][2] * kvnl[ib][itri][i][2]
	90	+ + trans[j] + 0.001;
	91	+ /**/
	92	+ glVertex3fv(vect);
	93	+ }
	94	+ }
	95	+ }
	96	+ /**/
	97	+ }
	98	+ glEnd();
	99	+ }
	100	+ /**/
	101	+} /* draw_ferm */
	102	+ /**
	103	+ * Draw lines of BZ boundaries
	104	+ */
	105	+void draw_bz_lines() {
	106	+ /**/
	107	+ int ibzl, i, j;
	108	+ GLfloat bzl2[3], bvec2[3][3], linecolor[4];
	109	+ /**/
	110	+ if (blackback == 1) {
	111	+ for (i = 0; i<4; i++) linecolor[i] = white[i];
	112	+ }
	113	+ else {
	114	+ for (i = 0; i<4; i++) linecolor[i] = black[i];
	115	+ }
	116	+ /**/
	117	+ glLineWidth(2.0);
	118	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, linecolor);
	119	+ glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, linecolor);
	120	+ /**/
	121	+ if (fbz == 1) {
	122	+ /**/
	123	+ glBegin(GL_LINES);
	124	+ for (ibzl = 0; ibzl < nbzl; ++ibzl) {
	125	+ for (i = 0; i< 2; ++i) {
	126	+ for (j = 0; j < 3; ++j)
	127	+ bzl2[j] = rot[j][0] * bzl[ibzl][i][0]
	128	+ + rot[j][1] * bzl[ibzl][i][1]
	129	+ + rot[j][2] * bzl[ibzl][i][2]
	130	+ + trans[j];
	131	+ glVertex3fv(bzl2);
	132	+ }
	133	+ }
	134	+ glEnd();
	135	+ /**/
	136	+ }
	137	+ else {
	138	+ /**/
	139	+ for (i = 0; i < 3; ++i) {
	140	+ for (j = 0; j < 3; ++j) {
	141	+ bvec2[i][j] = rot[j][0] * bvec[i][0]
	142	+ + rot[j][1] * bvec[i][1]
	143	+ + rot[j][2] * bvec[i][2];
	144	+ }
	145	+ }
	146	+ /**/
	147	+ glBegin(GL_LINE_STRIP);
	148	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i]; glVertex3fv(bzl2);
	149	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[0][i]; glVertex3fv(bzl2);
	150	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[1][i]; glVertex3fv(bzl2);
	151	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[1][i] + bvec2[2][i]; glVertex3fv(bzl2);
	152	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[1][i] + bvec2[2][i]; glVertex3fv(bzl2);
	153	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[1][i]; glVertex3fv(bzl2);
	154	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i]; glVertex3fv(bzl2);
	155	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[2][i]; glVertex3fv(bzl2);
	156	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[2][i]; glVertex3fv(bzl2);
	157	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[1][i] + bvec2[2][i]; glVertex3fv(bzl2);
	158	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[2][i]; glVertex3fv(bzl2);
	159	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[0][i]; glVertex3fv(bzl2);
	160	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i]; glVertex3fv(bzl2);
	161	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[2][i]; glVertex3fv(bzl2);
	162	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[1][i] + bvec2[2][i]; glVertex3fv(bzl2);
	163	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[1][i]; glVertex3fv(bzl2);
	164	+ for (i = 0; i<3; ++i) bzl2[i] = trans[i] + bvec2[0][i] + bvec2[1][i]; glVertex3fv(bzl2);
	165	+ glEnd();
	166	+ /**/
	167	+ }
	168	+ /**/
	169	+ glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
	170	+ /**/
	171	+} /* draw bz_lines */
	172	+ /**
	173	+ * Draw color scale
	174	+ */
	175	+void draw_colorbar()
	176	+{
	177	+ int i, j;
	178	+ GLfloat mat2, barcolor[4];
	179	+ GLfloat colorbar[13][3] = {
	180	+ { 0.0, 0.0, 1.0 },
	181	+ { -1.0, -1.0, 0.0 },
	182	+ { -1.0, -1.0 - 0.1, 0.0 },
	183	+ { -0.5, -1.0, 0.0 },
	184	+ { -0.5, -1.0 - 0.1, 0.0 },
	185	+ { 0.0, -1.0, 0.0 },
	186	+ { 0.0, -1.0 - 0.1, 0.0 },
	187	+ { 0.5, -1.0, 0.0 },
	188	+ { 0.5, -1.0 - 0.1, 0.0 },
	189	+ { 1.0, -1.0, 0.0 },
	190	+ { 1.0, -1.0 - 0.1, 0.0 },
	191	+ { 0.0, -1.0, 0.00001 },
	192	+ { 0.0, -1.0 - 0.1, 0.00001 }
	193	+ };
	194	+ /**/
	195	+ if (fcscl == 1 \|\| fcscl == 2) {
	196	+ glBegin(GL_TRIANGLE_STRIP);
	197	+ glNormal3fv(colorbar[0]);
	198	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, blue);
	199	+ glVertex3fv(colorbar[1]);
	200	+ glVertex3fv(colorbar[2]);
	201	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, cyan);
	202	+ glVertex3fv(colorbar[3]);
	203	+ glVertex3fv(colorbar[4]);
	204	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, green);
	205	+ glVertex3fv(colorbar[5]);
	206	+ glVertex3fv(colorbar[6]);
	207	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, yellow);
	208	+ glVertex3fv(colorbar[7]);
	209	+ glVertex3fv(colorbar[8]);
	210	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
	211	+ glVertex3fv(colorbar[9]);
	212	+ glVertex3fv(colorbar[10]);
	213	+ glEnd();
	214	+ /**/
	215	+ /* if(nodeline == 1){ */
	216	+ /* glLineWidth(2.0); */
	217	+ /* glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black); */
	218	+ /* glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black); */
	219	+ /* glBegin(GL_LINES); */
	220	+ /* glVertex3fv(colorbar[11]); */
	221	+ /* glVertex3fv(colorbar[12]); */
	222	+ /* glEnd(); */
	223	+ /* } */
	224	+ }
	225	+ else if (fcscl == 4) {
	226	+ for (i = 0; i <= 60; i++) {
	227	+ /**/
	228	+ mat2 = (GLfloat)i / 60.0 * 6.0;
	229	+ /**/
	230	+ if (mat2 <= 1.0) {
	231	+ for (j = 0; j<4; ++j) barcolor[j] = yellow[j] * mat2 + red[j] * (1.0 - mat2);
	232	+ }
	233	+ else if (mat2 <= 2.0) {
	234	+ mat2 = mat2 - 1.0;
	235	+ for (j = 0; j<4; ++j) barcolor[j] = green[j] * mat2 + yellow[j] * (1.0 - mat2);
	236	+ }
	237	+ else if (mat2 <= 3.0) {
	238	+ mat2 = mat2 - 2.0;
	239	+ for (j = 0; j<4; ++j) barcolor[j] = cyan[j] * mat2 + green[j] * (1.0 - mat2);
	240	+ }
	241	+ else if (mat2 <= 4.0) {
	242	+ mat2 = mat2 - 3.0;
	243	+ for (j = 0; j<4; ++j) barcolor[j] = blue[j] * mat2 + cyan[j] * (1.0 - mat2);
	244	+ }
	245	+ else if (mat2 <= 5.0) {
	246	+ mat2 = mat2 - 4.0;
	247	+ for (j = 0; j<4; ++j) barcolor[j] = magenta[j] * mat2 + blue[j] * (1.0 - mat2);
	248	+ }
	249	+ else {
	250	+ mat2 = mat2 - 5.0;
	251	+ for (j = 0; j<4; ++j) barcolor[j] = red[j] * mat2 + magenta[j] * (1.0 - mat2);
	252	+ }
	253	+ /**/
	254	+ glBegin(GL_TRIANGLES);
	255	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, barcolor);
	256	+ glNormal3fv(colorbar[0]);
	257	+ glVertex3f(0.15 * cos((GLfloat)(i + 1) / 60.0 * 6.283185307),
	258	+ 0.15 * sin((GLfloat)(i + 1) / 60.0 * 6.283185307) - 1.0, 0.0);
	259	+ glVertex3f(0.15 * cos((GLfloat)i / 60.0 * 6.283185307),
	260	+ 0.15 * sin((GLfloat)i / 60.0 * 6.283185307) - 1.0, 0.0);
	261	+ glVertex3f(0.2 * cos((GLfloat)i / 60.0 * 6.283185307),
	262	+ 0.2 * sin((GLfloat)i / 60.0 * 6.283185307) - 1.0, 0.0);
	263	+ glVertex3f(0.2 * cos((GLfloat)i / 60.0 * 6.283185307),
	264	+ 0.2 * sin((GLfloat)i / 60.0 * 6.283185307) - 1.0, 0.0);
	265	+ glVertex3f(0.2 * cos((GLfloat)(i + 1) / 60.0 * 6.283185307),
	266	+ 0.2 * sin((GLfloat)(i + 1) / 60.0 * 6.283185307) - 1.0, 0.0);
	267	+ glVertex3f(0.15 * cos((GLfloat)(i + 1) / 60.0 * 6.283185307),
	268	+ 0.15 * sin((GLfloat)(i + 1) / 60.0 * 6.283185307) - 1.0, 0.0);
	269	+ glEnd();
	270	+ }
	271	+ }
	272	+ else {
	273	+ }
	274	+} /* draw_colorbar */
	275	+ /**
	276	+ * Draw points for the stereogram
	277	+ */
	278	+void draw_circles() {
	279	+ int i;
	280	+ GLfloat r;
	281	+ /**/
	282	+ r = 0.05;
	283	+ /**/
	284	+ if (blackback == 1) {
	285	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, white);
	286	+ glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, white);
	287	+ }
	288	+ else {
	289	+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black);
	290	+ glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
	291	+ }
	292	+ /**/
	293	+ glBegin(GL_TRIANGLE_FAN);
	294	+ glNormal3f(0.0, 0.0, 1.0);
	295	+ glVertex3f(0.7, scl, 0.0);
	296	+ for (i = 0; i <= 20; i++) {
	297	+ glVertex3f(r * cos((GLfloat)i / 20.0 * 6.283185307) + 0.7,
	298	+ r * sin((GLfloat)i / 20.0 * 6.283185307) + scl, 0.0);
	299	+ }
	300	+ glEnd();
	301	+ /**/
	302	+ glBegin(GL_TRIANGLE_FAN);
	303	+ glNormal3f(0.0, 0.0, 1.0);
	304	+ glVertex3f(-0.7, scl, 0.0);
	305	+ for (i = 0; i <= 20; i++) {
	306	+ glVertex3f(r * cos((GLfloat)i / 20.0 * 6.283185307) - 0.7,
	307	+ r * sin((GLfloat)i / 20.0 * 6.283185307) + scl, 0.0);
	308	+ }
	309	+ glEnd();
	310	+}
	311	+/*!
	312	+Glut Display function
	313	+*/
	314	+void display()
	315	+{
	316	+ GLfloat pos[] = { 1.0, 1.0, 1.0, 0.0 };
	317	+ GLfloat amb[] = { 0.2, 0.2, 0.2, 0.0 };
	318	+ double dx, theta, posz, phi;
	319	+ GLfloat pos1[4], pos2[4];
	320	+ /**/
	321	+ if (lstereo == 2) {
	322	+ theta = 3.1416 / 180.0 * 2.50;
	323	+ posz = 5.0;
	324	+ dx = 0.7;
	325	+ phi = atan(posz / dx) - theta;
	326	+ theta = (3.1416 * 0.50 - phi) / 3.1416 * 180.0;
	327	+ /**/
	328	+ pos1[0] = posz * cos(phi) - dx;
	329	+ pos1[1] = 0.0;
	330	+ pos1[2] = posz * sin(phi);
	331	+ pos1[3] = 1.0;
	332	+ /**/
	333	+ pos2[0] = -posz * cos(phi) + dx;
	334	+ pos2[1] = 0.0;
	335	+ pos2[2] = posz * sin(phi);
	336	+ pos2[3] = 1.0;
	337	+ }
	338	+ else if (lstereo == 3) {
	339	+ theta = -3.1416 / 180.0 * 2.0;
	340	+ posz = 5.0;
	341	+ dx = -0.7;
	342	+ /**/
	343	+ pos1[0] = posz * sin(theta) - dx;
	344	+ pos1[1] = 0.0;
	345	+ pos1[2] = posz * cos(theta);
	346	+ pos1[3] = 1.0;
	347	+ /**/
	348	+ pos2[0] = -posz * sin(theta) + dx;
	349	+ pos2[1] = 0.0;
	350	+ pos2[2] = posz * cos(theta);
	351	+ pos2[3] = 1.0;
	352	+ }
	353	+ else {
	354	+ theta = 0.0;
	355	+ dx = 0.0;
	356	+ }
	357	+ /*
	358	+ Initialize
	359	+ */
	360	+ glClear(GL_COLOR_BUFFER_BIT \| GL_DEPTH_BUFFER_BIT);
	361	+ glLoadIdentity();
	362	+ /*
	363	+ Set view point & view line
	364	+ */
	365	+ gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
	366	+ /*
	367	+ Set position of light
	368	+ */
	369	+ if (lstereo == 1) {
	370	+ glLightfv(GL_LIGHT0, GL_POSITION, pos);
	371	+ /*
	372	+ Draw color scale
	373	+ */
	374	+ if (lcolorbar == 1) draw_colorbar();
	375	+ }
	376	+ else {
	377	+ glLightfv(GL_LIGHT0, GL_POSITION, pos1);
	378	+ draw_circles();
	379	+ glTranslated(-dx, 0.0, 0.0);
	380	+ glRotated(theta, 0.0, 1.0, 0.0);
	381	+ }
	382	+ glLightfv(GL_LIGHT1, GL_AMBIENT, amb);
	383	+ /*
	384	+ Rotation & Zoom
	385	+ */
	386	+ glScaled(scl, scl, scl);
	387	+ /*
	388	+ Draw Brillouin zone boundaries
	389	+ */
	390	+ draw_bz_lines();
	391	+ /*
	392	+ Draw Fermi surfaces
	393	+ */
	394	+ draw_fermi();
	395	+ /*
	396	+ Draw the second
	397	+ */
	398	+ if (lstereo != 1) {
	399	+ glPushMatrix();
	400	+ glLoadIdentity();
	401	+ gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
	402	+ glLightfv(GL_LIGHT0, GL_POSITION, pos2);
	403	+ /**/
	404	+ glTranslated(dx, 0.0, 0.0);
	405	+ glRotated(-theta, 0.0, 1.0, 0.0);
	406	+ /**/
	407	+ glScaled(scl, scl, scl);
	408	+ draw_bz_lines();
	409	+ draw_fermi();
	410	+ /**/
	411	+ glPopMatrix();
	412	+ }
	413	+ /**/
	414	+ glutSwapBuffers();
	415	+} /* display */

--- /dev/null

+++ b/src/draw.h

		@@ -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 display();

--- /dev/null

+++ b/src/fermi_patch.c

		@@ -0,0 +1,331 @@
	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	+#include <stdlib.h>
	26	+#include <stdio.h>
	27	+#include "variable.h"
	28	+#include "basic_math.h"
	29	+
	30	+#if defined(MAC)
	31	+#include <GLUT/glut.h>
	32	+#else
	33	+#include <GL/glut.h>
	34	+#endif
	35	+
	36	+/**
	37	+* Store triangle patch
	38	+*/
	39	+void triangle(
	40	+ int ib /*<[in] The band index/,
	41	+ int nbr /*<[in] Bragg plane/,
	42	+ GLfloat mat1[3] /*<[in] The matrix element/,
	43	+ GLfloat kvec1[3][3] /*<[in] k-vector of corners/)
	44	+{
	45	+ /**/
	46	+ int ibr, i, j;
	47	+ GLfloat prod[3], thr = 0.0000, mat2[3], kvec2[3][3];
	48	+ /**/
	49	+ if (fbz == 1) {
	50	+ /**/
	51	+ for (ibr = nbr; ibr < 26; ++ibr) {
	52	+ /**/
	53	+ for (i = 0; i < 3; ++i) {
	54	+ prod[i] = bragg[ibr][0] * kvec1[i][0]
	55	+ + bragg[ibr][1] * kvec1[i][1]
	56	+ + bragg[ibr][2] * kvec1[i][2];
	57	+ }
	58	+ eigsort(3, prod, mat1, kvec1);
	59	+ /**/
	60	+ if (brnrm[ibr] + thr < prod[0]) {
	61	+ return;
	62	+ }
	63	+ else if (brnrm[ibr] + thr < prod[1]) {
	64	+ mat2[0] = mat1[0];
	65	+ mat2[1] = mat1[0] * (brnrm[ibr] - prod[1]) / (prod[0] - prod[1])
	66	+ + mat1[1] * (brnrm[ibr] - prod[0]) / (prod[1] - prod[0]);
	67	+ mat2[2] = mat1[0] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
	68	+ + mat1[2] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
	69	+ for (i = 0; i<3; ++i) kvec2[0][i] = kvec1[0][i];
	70	+ for (i = 0; i<3; ++i)
	71	+ kvec2[1][i] = kvec1[0][i] * (brnrm[ibr] - prod[1]) / (prod[0] - prod[1])
	72	+ + kvec1[1][i] * (brnrm[ibr] - prod[0]) / (prod[1] - prod[0]);
	73	+ for (i = 0; i<3; ++i)
	74	+ kvec2[2][i] = kvec1[0][i] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
	75	+ + kvec1[2][i] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
	76	+ triangle(ib, ibr + 1, mat2, kvec2);
	77	+ return;
	78	+ }
	79	+ else if (brnrm[ibr] + thr < prod[2]) {
	80	+ mat2[0] = mat1[0];
	81	+ mat2[1] = mat1[1];
	82	+ mat2[2] = mat1[0] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
	83	+ + mat1[2] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
	84	+ for (i = 0; i<3; ++i) kvec2[0][i] = kvec1[0][i];
	85	+ for (i = 0; i<3; ++i) kvec2[1][i] = kvec1[1][i];
	86	+ for (i = 0; i<3; ++i)
	87	+ kvec2[2][i] = kvec1[0][i] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
	88	+ + kvec1[2][i] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
	89	+ triangle(ib, ibr + 1, mat2, kvec2);
	90	+ /**/
	91	+ mat2[0] = mat1[1] * (brnrm[ibr] - prod[2]) / (prod[1] - prod[2])
	92	+ + mat1[2] * (brnrm[ibr] - prod[1]) / (prod[2] - prod[1]);
	93	+ mat2[1] = mat1[1];
	94	+ mat2[2] = mat1[0] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
	95	+ + mat1[2] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
	96	+ for (i = 0; i<3; ++i)
	97	+ kvec2[0][i] = kvec1[1][i] * (brnrm[ibr] - prod[2]) / (prod[1] - prod[2])
	98	+ + kvec1[2][i] * (brnrm[ibr] - prod[1]) / (prod[2] - prod[1]);
	99	+ for (i = 0; i<3; ++i) kvec2[1][i] = kvec1[1][i];
	100	+ for (i = 0; i<3; ++i)
	101	+ kvec2[2][i] = kvec1[0][i] * (brnrm[ibr] - prod[2]) / (prod[0] - prod[2])
	102	+ + kvec1[2][i] * (brnrm[ibr] - prod[0]) / (prod[2] - prod[0]);
	103	+ triangle(ib, ibr + 1, mat2, kvec2);
	104	+ return;
	105	+ }
	106	+ else {
	107	+ } /* brnrm[ibr] + thr < prod */
	108	+ } /* for ibr = 1; ibr < 26*/
	109	+ } /* if fbz == 1 */
	110	+ /**/
	111	+ if (query == 1) {
	112	+ ntri[ib] = ntri[ib] + 1;
	113	+ }
	114	+ else {
	115	+ normal_vec(kvec1[0], kvec1[1], kvec1[2], nmlp[ib][ntri[ib]]);
	116	+ for (i = 0; i < 3; ++i) {
	117	+ matp[ib][ntri[ib]][i] = mat1[i];
	118	+ for (j = 0; j < 3; ++j) {
	119	+ kvp[ib][ntri[ib]][i][j] = kvec1[i][j];
	120	+ }
	121	+ }
	122	+ ntri[ib] = ntri[ib] + 1;
	123	+ }
	124	+ /**/
	125	+}/* triangle */
	126	+ /**
	127	+ * Tetrahedrron method
	128	+ */
	129	+void tetrahedron(
	130	+ int ib /*< [in] The band index/,
	131	+ GLfloat eig1[8] /*< [in] Orbital energies/,
	132	+ GLfloat mat1[8] /*< [in] Matrix elements/,
	133	+ GLfloat kvec1[8][3] /*< [in] k vectors/)
	134	+{
	135	+ /**/
	136	+ int it, i, j;
	137	+ GLfloat eig2[4], mat2[4], kvec2[4][3], a[4][4], kvec3[3][3], mat3[3];
	138	+ /**/
	139	+ for (it = 0; it < 6; ++it) {
	140	+ /*
	141	+ Define corners of the tetrahedron
	142	+ */
	143	+ for (i = 0; i < 4; ++i) {
	144	+ eig2[i] = eig1[corner[it][i]];
	145	+ mat2[i] = mat1[corner[it][i]];
	146	+ /**/
	147	+ for (j = 0; j < 3; ++j)
	148	+ kvec2[i][j] = bvec[0][j] * kvec1[corner[it][i]][0]
	149	+ + bvec[1][j] * kvec1[corner[it][i]][1]
	150	+ + bvec[2][j] * kvec1[corner[it][i]][2];
	151	+ }
	152	+ /*
	153	+ Sort of eig1
	154	+ */
	155	+ eigsort(4, eig2, mat2, kvec2);
	156	+ for (i = 0; i < 4; ++i) {
	157	+ for (j = 0; j < 4; ++j) {
	158	+ a[i][j] = (0.00 - eig2[j]) / (eig2[i] - eig2[j]);
	159	+ }
	160	+ }
	161	+ /*
	162	+ Draw triangle in each cases
	163	+ */
	164	+ if (eig2[0] <= 0.00 && 0.00 < eig2[1]) {
	165	+ for (i = 0; i < 3; ++i) kvec3[0][i] = kvec2[0][i] * a[0][1] + kvec2[1][i] * a[1][0];
	166	+ for (i = 0; i < 3; ++i) kvec3[1][i] = kvec2[0][i] * a[0][2] + kvec2[2][i] * a[2][0];
	167	+ for (i = 0; i < 3; ++i) kvec3[2][i] = kvec2[0][i] * a[0][3] + kvec2[3][i] * a[3][0];
	168	+ mat3[0] = mat2[0] * a[0][1] + mat2[1] * a[1][0];
	169	+ mat3[1] = mat2[0] * a[0][2] + mat2[2] * a[2][0];
	170	+ mat3[2] = mat2[0] * a[0][3] + mat2[3] * a[3][0];
	171	+ triangle(ib, 0, mat3, kvec3);
	172	+ }
	173	+ else if (eig2[1] <= 0.00 && 0.00 < eig2[2]) {
	174	+ for (i = 0; i < 3; ++i) kvec3[0][i] = kvec2[0][i] * a[0][2] + kvec2[2][i] * a[2][0];
	175	+ for (i = 0; i < 3; ++i) kvec3[1][i] = kvec2[0][i] * a[0][3] + kvec2[3][i] * a[3][0];
	176	+ for (i = 0; i < 3; ++i) kvec3[2][i] = kvec2[1][i] * a[1][2] + kvec2[2][i] * a[2][1];
	177	+ mat3[0] = mat2[0] * a[0][2] + mat2[2] * a[2][0];
	178	+ mat3[1] = mat2[0] * a[0][3] + mat2[3] * a[3][0];
	179	+ mat3[2] = mat2[1] * a[1][2] + mat2[2] * a[2][1];
	180	+ triangle(ib, 0, mat3, kvec3);
	181	+ /**/
	182	+ for (i = 0; i < 3; ++i) kvec3[0][i] = kvec2[1][i] * a[1][3] + kvec2[3][i] * a[3][1];
	183	+ for (i = 0; i < 3; ++i) kvec3[1][i] = kvec2[0][i] * a[0][3] + kvec2[3][i] * a[3][0];
	184	+ for (i = 0; i < 3; ++i) kvec3[2][i] = kvec2[1][i] * a[1][2] + kvec2[2][i] * a[2][1];
	185	+ mat3[0] = mat2[1] * a[1][3] + mat2[3] * a[3][1];
	186	+ mat3[1] = mat2[0] * a[0][3] + mat2[3] * a[3][0];
	187	+ mat3[2] = mat2[1] * a[1][2] + mat2[2] * a[2][1];
	188	+ triangle(ib, 0, mat3, kvec3);
	189	+ }
	190	+ else if (eig2[2] <= 0.00 && 0.00 < eig2[3]) {
	191	+ for (i = 0; i < 3; ++i) kvec3[0][i] = kvec2[3][i] * a[3][0] + kvec2[0][i] * a[0][3];
	192	+ for (i = 0; i < 3; ++i) kvec3[1][i] = kvec2[3][i] * a[3][1] + kvec2[1][i] * a[1][3];
	193	+ for (i = 0; i < 3; ++i) kvec3[2][i] = kvec2[3][i] * a[3][2] + kvec2[2][i] * a[2][3];
	194	+ mat3[0] = mat2[3] * a[3][0] + mat2[0] * a[0][3];
	195	+ mat3[1] = mat2[3] * a[3][1] + mat2[1] * a[1][3];
	196	+ mat3[2] = mat2[3] * a[3][2] + mat2[2] * a[2][3];
	197	+ triangle(ib, 0, mat3, kvec3);
	198	+ }
	199	+ else {
	200	+ }
	201	+ }
	202	+}/* tetrahedron */
	203	+ /**
	204	+ * Patches for FSs
	205	+ */
	206	+void fermi_patch()
	207	+{
	208	+ int ib, i0, i1, i2, ii0, ii1, ii2, j0, j1, j2, start[3], i, j;
	209	+ GLfloat kvec1[8][3], eig1[8], mat1[8];
	210	+ /**/
	211	+ if (fbz == 1) {
	212	+ for (i0 = 0; i0 < 3; ++i0) start[i0] = - ng[i0];
	213	+ }
	214	+ else {
	215	+ for (i0 = 0; i0 < 3; ++i0) start[i0] = 0;
	216	+ }
	217	+ /**/
	218	+#pragma omp parallel default(none) \
	219	+ shared(nb,ntri,start,ng,ng0,eig,EF,mat,shiftk) \
	220	+ private(ib,j0,j1,j2,i0,i1,i2,ii0,ii1,ii2,kvec1,eig1,mat1,i,j)
	221	+ {
	222	+#pragma omp for nowait
	223	+ for (ib = 0; ib < nb; ++ib) {
	224	+ ntri[ib] = 0;
	225	+ for (j0 = start[0]; j0 < ng[0]; ++j0) {
	226	+ for (j1 = start[1]; j1 < ng[1]; ++j1) {
	227	+ for (j2 = start[2]; j2 < ng[2]; ++j2) {
	228	+ /**/
	229	+ i0 = j0;
	230	+ i1 = j1;
	231	+ i2 = j2;
	232	+ ii0 = j0 + 1;
	233	+ ii1 = j1 + 1;
	234	+ ii2 = j2 + 1;
	235	+ /**/
	236	+ kvec1[0][0] = (GLfloat)i0 / (GLfloat)ng[0];
	237	+ kvec1[1][0] = (GLfloat)i0 / (GLfloat)ng[0];
	238	+ kvec1[2][0] = (GLfloat)i0 / (GLfloat)ng[0];
	239	+ kvec1[3][0] = (GLfloat)i0 / (GLfloat)ng[0];
	240	+ kvec1[4][0] = (GLfloat)ii0 / (GLfloat)ng[0];
	241	+ kvec1[5][0] = (GLfloat)ii0 / (GLfloat)ng[0];
	242	+ kvec1[6][0] = (GLfloat)ii0 / (GLfloat)ng[0];
	243	+ kvec1[7][0] = (GLfloat)ii0 / (GLfloat)ng[0];
	244	+ /**/
	245	+ kvec1[0][1] = (GLfloat)i1 / (GLfloat)ng[1];
	246	+ kvec1[1][1] = (GLfloat)i1 / (GLfloat)ng[1];
	247	+ kvec1[2][1] = (GLfloat)ii1 / (GLfloat)ng[1];
	248	+ kvec1[3][1] = (GLfloat)ii1 / (GLfloat)ng[1];
	249	+ kvec1[4][1] = (GLfloat)i1 / (GLfloat)ng[1];
	250	+ kvec1[5][1] = (GLfloat)i1 / (GLfloat)ng[1];
	251	+ kvec1[6][1] = (GLfloat)ii1 / (GLfloat)ng[1];
	252	+ kvec1[7][1] = (GLfloat)ii1 / (GLfloat)ng[1];
	253	+ /**/
	254	+ kvec1[0][2] = (GLfloat)i2 / (GLfloat)ng[2];
	255	+ kvec1[1][2] = (GLfloat)ii2 / (GLfloat)ng[2];
	256	+ kvec1[2][2] = (GLfloat)i2 / (GLfloat)ng[2];
	257	+ kvec1[3][2] = (GLfloat)ii2 / (GLfloat)ng[2];
	258	+ kvec1[4][2] = (GLfloat)i2 / (GLfloat)ng[2];
	259	+ kvec1[5][2] = (GLfloat)ii2 / (GLfloat)ng[2];
	260	+ kvec1[6][2] = (GLfloat)i2 / (GLfloat)ng[2];
	261	+ kvec1[7][2] = (GLfloat)ii2 / (GLfloat)ng[2];
	262	+ /**/
	263	+ for (i = 0; i < 8; i++)
	264	+ for (j = 0; j < 3; j++)
	265	+ kvec1[i][j] = kvec1[i][j] + (double)shiftk[j] / (GLfloat)(2 * ng0[j]);
	266	+ /**/
	267	+ i0 = modulo(i0, ng[0]);
	268	+ i1 = modulo(i1, ng[1]);
	269	+ i2 = modulo(i2, ng[2]);
	270	+ ii0 = modulo(ii0, ng[0]);
	271	+ ii1 = modulo(ii1, ng[1]);
	272	+ ii2 = modulo(ii2, ng[2]);
	273	+ /**/
	274	+ eig1[0] = eig[ib][i0][i1][i2] - EF;
	275	+ eig1[1] = eig[ib][i0][i1][ii2] - EF;
	276	+ eig1[2] = eig[ib][i0][ii1][i2] - EF;
	277	+ eig1[3] = eig[ib][i0][ii1][ii2] - EF;
	278	+ eig1[4] = eig[ib][ii0][i1][i2] - EF;
	279	+ eig1[5] = eig[ib][ii0][i1][ii2] - EF;
	280	+ eig1[6] = eig[ib][ii0][ii1][i2] - EF;
	281	+ eig1[7] = eig[ib][ii0][ii1][ii2] - EF;
	282	+ /**/
	283	+ mat1[0] = mat[ib][i0][i1][i2];
	284	+ mat1[1] = mat[ib][i0][i1][ii2];
	285	+ mat1[2] = mat[ib][i0][ii1][i2];
	286	+ mat1[3] = mat[ib][i0][ii1][ii2];
	287	+ mat1[4] = mat[ib][ii0][i1][i2];
	288	+ mat1[5] = mat[ib][ii0][i1][ii2];
	289	+ mat1[6] = mat[ib][ii0][ii1][i2];
	290	+ mat1[7] = mat[ib][ii0][ii1][ii2];
	291	+ /**/
	292	+ tetrahedron(ib, eig1, mat1, kvec1);
	293	+ }
	294	+ }
	295	+ }
	296	+ }
	297	+ } /* End of parallel region */
	298	+ /**/
	299	+ if (query == 1) {
	300	+ printf("band # of patchs \n");
	301	+ for (ib = 0; ib < nb; ib++) {
	302	+ printf("%d %d \n", ib + 1, ntri[ib]);
	303	+ }
	304	+ printf("\n");
	305	+ /*
	306	+ Allocation of triangler patches
	307	+ */
	308	+ nmlp = (GLfloat**)malloc(nb sizeof(GLfloat**));
	309	+ matp = (GLfloat**)malloc(nb sizeof(GLfloat**));
	310	+ clr = (GLfloat***)malloc(nb sizeof(GLfloat***));
	311	+ kvp = (GLfloat***)malloc(nb sizeof(GLfloat***));
	312	+ for (ib = 0; ib < nb; ++ib) {
	313	+ nmlp[ib] = (GLfloat*)malloc(ntri[ib] sizeof(GLfloat*));
	314	+ matp[ib] = (GLfloat*)malloc(ntri[ib] sizeof(GLfloat*));
	315	+ clr[ib] = (GLfloat**)malloc(ntri[ib] sizeof(GLfloat**));
	316	+ kvp[ib] = (GLfloat**)malloc(ntri[ib] sizeof(GLfloat**));
	317	+ for (i0 = 0; i0 < ntri[ib]; ++i0) {
	318	+ nmlp[ib][i0] = (GLfloat)malloc(3 sizeof(GLfloat));
	319	+ matp[ib][i0] = (GLfloat)malloc(3 sizeof(GLfloat));
	320	+ clr[ib][i0] = (GLfloat*)malloc(3 sizeof(GLfloat*));
	321	+ kvp[ib][i0] = (GLfloat*)malloc(3 sizeof(GLfloat*));
	322	+ for (i1 = 0; i1 < 3; ++i1) {
	323	+ kvp[ib][i0][i1] = (GLfloat)malloc(3 sizeof(GLfloat));
	324	+ clr[ib][i0][i1] = (GLfloat)malloc(4 sizeof(GLfloat));
	325	+ }
	326	+ }
	327	+ }
	328	+ /**/
	329	+ }
	330	+ /**/
	331	+} /* fermi_patch */

--- /dev/null

+++ b/src/fermi_patch.h

		@@ -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();

--- a/src/fermisurfer.c

+++ b/src/fermisurfer.c

		@@ -36,7 +36,16 @@
36	36	*/
37	37	#include <stdlib.h>
38	38	#include <stdio.h>
39		-#include <math.h>
	39	+#include "variable.h"
	40	+#include "read_file.h"
	41	+#include "menu.h"
	42	+#include "operation.h"
	43	+#include "initialize.h"
	44	+#include "draw.h"
	45	+#include "fermi_patch.h"
	46	+#include "calc_nodeline.h"
	47	+#include "bz_lines.h"
	48	+#include "free_patch.h"
40	49
41	50	#if defined(MAC)
42	51	#include <GLUT/glut.h>

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

		@@ -2394,7 +85,7 @@ int main(
2394	85	getchar();
2395	86	exit(-1);
2396	87	}
2397		-
	88	+ initialize_val();
2398	89	/**/
2399	90	read_file(argv[1]);
2400	91	init_corner();

		@@ -2407,10 +98,12 @@ int main(
2407	98	/**/
2408	99	max_and_min_bz();
2409	100	/**/
2410		-#pragma omp parallel
2411	101	#if defined(_OPENMP)
2412		- printf("Threads : %d \n", omp_get_num_threads());
	102	+#pragma omp parallel
	103	+#pragma omp master
	104	+ nthreads = omp_get_num_threads();
2413	105	#endif
	106	+ printf("Number of threads : %d \n", nthreads);
2414	107	/**/
2415	108	printf("\n##### First Brillouin zone mode ##### \n\n");
2416	109	query = 1;

--- /dev/null

+++ b/src/free_patch.c

		@@ -0,0 +1,201 @@
	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	+#include <stdlib.h>
	26	+#include <math.h>
	27	+#include <stdio.h>
	28	+#include "variable.h"
	29	+
	30	+/**
	31	+* Free variables for patch
	32	+*/
	33	+void free_patch() {
	34	+ int ib, i0, i1;
	35	+
	36	+ for (ib = 0; ib < nb; ++ib) {
	37	+ for (i0 = 0; i0 < ntri[ib]; ++i0) {
	38	+ for (i1 = 0; i1 < 3; ++i1) {
	39	+ free(kvp[ib][i0][i1]);
	40	+ free(clr[ib][i0][i1]);
	41	+ }
	42	+ free(nmlp[ib][i0]);
	43	+ free(matp[ib][i0]);
	44	+ free(clr[ib][i0]);
	45	+ free(kvp[ib][i0]);
	46	+ }
	47	+ free(nmlp[ib]);
	48	+ free(matp[ib]);
	49	+ free(clr[ib]);
	50	+ free(kvp[ib]);
	51	+ }
	52	+ free(nmlp);
	53	+ free(matp);
	54	+ free(clr);
	55	+ free(kvp);
	56	+
	57	+ for (ib = 0; ib < nb; ++ib) {
	58	+ for (i0 = 0; i0 < nnl[ib]; ++i0) {
	59	+ for (i1 = 0; i1 < 2; ++i1) {
	60	+ free(kvnl[ib][i0][i1]);
	61	+ }
	62	+ free(kvnl[ib][i0]);
	63	+ }
	64	+ free(kvnl[ib]);
	65	+ }
	66	+ free(kvnl);
	67	+}
	68	+/**
	69	+* Max. & Min. of matrix elements.
	70	+*/
	71	+void max_and_min() {
	72	+ int ib, itri, i, j, ierr;
	73	+ GLfloat matmax, matmin, mat2;
	74	+ /**/
	75	+ matmax = -100000000.0000;
	76	+ matmin = 100000000.0000;
	77	+ /**/
	78	+ for (ib = 0; ib < nb; ib++) {
	79	+ for (itri = 0; itri < ntri[ib]; ++itri) {
	80	+ for (i = 0; i < 3; ++i) {
	81	+ if (matp[ib][itri][i] > matmax) matmax = matp[ib][itri][i];
	82	+ if (matp[ib][itri][i] < matmin) matmin = matp[ib][itri][i];
	83	+ }
	84	+ }
	85	+ }
	86	+ /**/
	87	+ printf("Max. value : %f \n", matmax);
	88	+ printf("Min. value : %f \n \n", matmin);
	89	+ /**/
	90	+ if (fcscl == 2) {
	91	+ printf("Set min. value : ");
	92	+ ierr = scanf("%f", &matmin);
	93	+ if (ierr == 0) printf("error ! reading min");
	94	+ printf("Set max. value : ");
	95	+ ierr = scanf("%f", &matmax);
	96	+ if (ierr == 0) printf("error ! reading max");
	97	+ }
	98	+ /**/
	99	+ if (fcscl == 1 \|\| fcscl == 2) {
	100	+ for (ib = 0; ib < nb; ib++) {
	101	+ for (itri = 0; itri < ntri[ib]; ++itri) {
	102	+ for (i = 0; i < 3; ++i) {
	103	+ /**/
	104	+ mat2 = (matp[ib][itri][i] - matmin) / (matmax - matmin);
	105	+ mat2 = mat2 * 4.0;
	106	+ /**/
	107	+ if (mat2 <= 1.0) {
	108	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = cyan[j] * mat2 + blue[j] * (1.0 - mat2);
	109	+ }
	110	+ else if (mat2 <= 2.0) {
	111	+ mat2 = mat2 - 1.0;
	112	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = green[j] * mat2 + cyan[j] * (1.0 - mat2);
	113	+ }
	114	+ else if (mat2 <= 3.0) {
	115	+ mat2 = mat2 - 2.0;
	116	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = yellow[j] * mat2 + green[j] * (1.0 - mat2);
	117	+ }
	118	+ else {
	119	+ mat2 = mat2 - 3.0;
	120	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = red[j] * mat2 + yellow[j] * (1.0 - mat2);
	121	+ }
	122	+ }
	123	+ }
	124	+ }
	125	+ }
	126	+ else if (fcscl == 4) {
	127	+ for (ib = 0; ib < nb; ib++) {
	128	+ for (itri = 0; itri < ntri[ib]; ++itri) {
	129	+ for (i = 0; i < 3; ++i) {
	130	+ /**/
	131	+ mat2 = matp[ib][itri][i] / 6.283185307;
	132	+ mat2 = mat2 - floorf(mat2);
	133	+ mat2 = mat2 * 6.0;
	134	+ /**/
	135	+ if (mat2 <= 1.0) {
	136	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = yellow[j] * mat2 + red[j] * (1.0 - mat2);
	137	+ }
	138	+ else if (mat2 <= 2.0) {
	139	+ mat2 = mat2 - 1.0;
	140	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = green[j] * mat2 + yellow[j] * (1.0 - mat2);
	141	+ }
	142	+ else if (mat2 <= 3.0) {
	143	+ mat2 = mat2 - 2.0;
	144	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = cyan[j] * mat2 + green[j] * (1.0 - mat2);
	145	+ }
	146	+ else if (mat2 <= 4.0) {
	147	+ mat2 = mat2 - 3.0;
	148	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = blue[j] * mat2 + cyan[j] * (1.0 - mat2);
	149	+ }
	150	+ else if (mat2 <= 5.0) {
	151	+ mat2 = mat2 - 4.0;
	152	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = magenta[j] * mat2 + blue[j] * (1.0 - mat2);
	153	+ }
	154	+ else {
	155	+ mat2 = mat2 - 5.0;
	156	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = red[j] * mat2 + magenta[j] * (1.0 - mat2);
	157	+ }
	158	+ }
	159	+ }
	160	+ }
	161	+ }
	162	+ else {
	163	+ for (ib = 0; ib < nb; ib++) {
	164	+ /**/
	165	+ mat2 = 1.0 / (GLfloat)(nb - 1) * (GLfloat)ib;
	166	+ mat2 = mat2 * 4.0;
	167	+ /**/
	168	+ if (mat2 <= 1.0) {
	169	+ for (itri = 0; itri < ntri[ib]; ++itri) {
	170	+ for (i = 0; i < 3; ++i) {
	171	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = cyan[j] * mat2 + blue[j] * (1.0 - mat2);
	172	+ }
	173	+ }
	174	+ }
	175	+ else if (mat2 <= 2.0) {
	176	+ mat2 = mat2 - 1.0;
	177	+ for (itri = 0; itri < ntri[ib]; ++itri) {
	178	+ for (i = 0; i < 3; ++i) {
	179	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = green[j] * mat2 + cyan[j] * (1.0 - mat2);
	180	+ }
	181	+ }
	182	+ }
	183	+ else if (mat2 <= 3.0) {
	184	+ mat2 = mat2 - 2.0;
	185	+ for (itri = 0; itri < ntri[ib]; ++itri) {
	186	+ for (i = 0; i < 3; ++i) {
	187	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = yellow[j] * mat2 + green[j] * (1.0 - mat2);
	188	+ }
	189	+ }
	190	+ }
	191	+ else {
	192	+ mat2 = mat2 - 3.0;
	193	+ for (itri = 0; itri < ntri[ib]; ++itri) {
	194	+ for (i = 0; i < 3; ++i) {
	195	+ for (j = 0; j<4; ++j) clr[ib][itri][i][j] = red[j] * mat2 + yellow[j] * (1.0 - mat2);
	196	+ }
	197	+ }
	198	+ }
	199	+ }
	200	+ }
	201	+} /* max_and_min */

--- /dev/null

+++ b/src/free_patch.h

		@@ -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 free_patch();
	26	+void max_and_min();

--- /dev/null

+++ b/src/initialize.c

		@@ -0,0 +1,273 @@
	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	+#include <stdio.h>
	26	+#include "variable.h"
	27	+
	28	+/**
	29	+* Initialize corners of tetrahedron
	30	+*/
	31	+void init_corner() {
	32	+ int i, j;
	33	+ int corner1[16][6][4] = {
	34	+ /*
	35	+ [0] min = 0-7
	36	+ */
	37	+ { { 0, 1, 3, 7 },
	38	+ { 0, 1, 5, 7 },
	39	+ { 0, 2, 3, 7 },
	40	+ { 0, 2, 6, 7 },
	41	+ { 0, 4, 5, 7 },
	42	+ { 0, 4, 6, 7 } },
	43	+ /*
	44	+ [1] min = 1-6
	45	+ */
	46	+ { { 1, 0, 2, 6 },
	47	+ { 1, 0, 4, 6 },
	48	+ { 1, 3, 2, 6 },
	49	+ { 1, 3, 7, 6 },
	50	+ { 1, 5, 4, 6 },
	51	+ { 1, 5, 7, 6 } },
	52	+ /*
	53	+ [2] min = 2-5
	54	+ */
	55	+ { { 2, 0, 1, 5 },
	56	+ { 2, 0, 4, 5 },
	57	+ { 2, 3, 1, 5 },
	58	+ { 2, 3, 7, 5 },
	59	+ { 2, 6, 4, 5 },
	60	+ { 2, 6, 7, 5 } },
	61	+ /*
	62	+ [3] min = 3-4
	63	+ */
	64	+ { { 4, 0, 1, 3 },
	65	+ { 4, 0, 2, 3 },
	66	+ { 4, 5, 1, 3 },
	67	+ { 4, 5, 7, 3 },
	68	+ { 4, 6, 2, 3 },
	69	+ { 4, 6, 7, 3 } },
	70	+ /*
	71	+ [4] min = 0-7, max = 3-4
	72	+ */
	73	+ { { 0, 4, 5, 6 },
	74	+ { 1, 2, 3, 7 },
	75	+ { 0, 7, 2, 6 },
	76	+ { 0, 7, 1, 5 },
	77	+ { 0, 7, 1, 2 },
	78	+ { 0, 7, 6, 5 } },
	79	+ /*
	80	+ [5] min = 1-6, max = 3-4
	81	+ */
	82	+ { { 0, 4, 5, 6 },
	83	+ { 1, 2, 3, 7 },
	84	+ { 1, 6, 5, 7 },
	85	+ { 1, 6, 7, 2 },
	86	+ { 1, 6, 2, 0 },
	87	+ { 1, 6, 0, 5 } },
	88	+ /*
	89	+ [6] min = 2-5, max = 3-4
	90	+ */
	91	+ { { 0, 4, 5, 6 },
	92	+ { 1, 2, 3, 7 },
	93	+ { 2, 5, 7, 6 },
	94	+ { 2, 5, 6, 0 },
	95	+ { 2, 5, 0, 1 },
	96	+ { 2, 5, 1, 7 } },
	97	+ /*
	98	+ [7] min = 3-4, max = 0-7
	99	+ */
	100	+ { { 0, 1, 2, 4 },
	101	+ { 7, 3, 5, 6 },
	102	+ { 3, 4, 1, 5 },
	103	+ { 3, 4, 5, 6 },
	104	+ { 3, 4, 6, 2 },
	105	+ { 3, 4, 2, 1 } },
	106	+ /*
	107	+ [8] min = 2-5, max = 0-7
	108	+ */
	109	+ { { 0, 1, 2, 4 },
	110	+ { 7, 3, 5, 6 },
	111	+ { 2, 5, 1, 3 },
	112	+ { 2, 5, 3, 6 },
	113	+ { 2, 5, 6, 4 },
	114	+ { 2, 5, 4, 1 } },
	115	+ /*
	116	+ [9] min = 1-6, max = 0-7
	117	+ */
	118	+ { { 0, 1, 2, 4 },
	119	+ { 7, 3, 5, 6 },
	120	+ { 1, 6, 2, 3 },
	121	+ { 1, 6, 3, 5 },
	122	+ { 1, 6, 5, 4 },
	123	+ { 1, 6, 4, 2 } },
	124	+ /*
	125	+ [10] min = 0-7, max = 1-6
	126	+ */
	127	+ { { 1, 0, 3, 5 },
	128	+ { 6, 2, 4, 7 },
	129	+ { 0, 7, 2, 3 },
	130	+ { 0, 7, 3, 5 },
	131	+ { 0, 7, 5, 4 },
	132	+ { 0, 7, 4, 2 } },
	133	+ /*
	134	+ [11] min = 3-4, max = 1-6
	135	+ */
	136	+ { { 1, 0, 3, 5 },
	137	+ { 6, 2, 4, 7 },
	138	+ { 3, 4, 0, 2 },
	139	+ { 3, 4, 2, 7 },
	140	+ { 3, 4, 7, 5 },
	141	+ { 3, 4, 5, 0 } },
	142	+ /*
	143	+ [12] min = 2-5, max = 1-6
	144	+ */
	145	+ { { 1, 0, 3, 5 },
	146	+ { 6, 2, 4, 7 },
	147	+ { 2, 5, 0, 3 },
	148	+ { 2, 5, 3, 7 },
	149	+ { 2, 5, 7, 4 },
	150	+ { 2, 5, 4, 0 } },
	151	+ /*
	152	+ [13] min = 0-7, max = 2-5
	153	+ */
	154	+ { { 2, 0, 3, 6 },
	155	+ { 5, 1, 4, 7 },
	156	+ { 0, 7, 1, 3 },
	157	+ { 0, 7, 3, 6 },
	158	+ { 0, 7, 6, 4 },
	159	+ { 0, 7, 4, 1 } },
	160	+ /*
	161	+ [14] min = 1-6, max = 2-5
	162	+ */
	163	+ { { 2, 0, 3, 6 },
	164	+ { 5, 1, 4, 7 },
	165	+ { 1, 6, 0, 3 },
	166	+ { 1, 6, 3, 7 },
	167	+ { 1, 6, 7, 4 },
	168	+ { 1, 6, 4, 0 } },
	169	+ /*
	170	+ [15] min = 3-4, max = 2-5
	171	+ */
	172	+ { { 2, 0, 3, 6 },
	173	+ { 5, 1, 4, 7 },
	174	+ { 3, 4, 0, 1 },
	175	+ { 3, 4, 1, 7 },
	176	+ { 3, 4, 7, 6 },
	177	+ { 3, 4, 6, 0 } },
	178	+ };
	179	+ /**/
	180	+ for (i = 0; i < 6; ++i) {
	181	+ for (j = 0; j < 4; ++j) {
	182	+ corner[i][j] = corner1[itet][i][j];
	183	+ }
	184	+ }
	185	+}
	186	+/**
	187	+* Compute Bragg vetor
	188	+*/
	189	+void bragg_vector() {
	190	+ int i0, i1, i2, i, ibr;
	191	+ /**/
	192	+ ibr = 0;
	193	+ /**/
	194	+ for (i0 = -1; i0 <= 1; ++i0) {
	195	+ for (i1 = -1; i1 <= 1; ++i1) {
	196	+ for (i2 = -1; i2 <= 1; ++i2) {
	197	+ /**/
	198	+ if (i0 == 0 && i1 == 0 && i2 == 0) {
	199	+ }
	200	+ else {
	201	+ for (i = 0; i < 3; ++i)
	202	+ bragg[ibr][i] = ((GLfloat)i0 * bvec[0][i]
	203	+ + (GLfloat)i1 * bvec[1][i]
	204	+ + (GLfloat)i2 * bvec[2][i]) * 0.5;
	205	+ /**/
	206	+ brnrm[ibr] = bragg[ibr][0] * bragg[ibr][0]
	207	+ + bragg[ibr][1] * bragg[ibr][1]
	208	+ + bragg[ibr][2] * bragg[ibr][2];
	209	+ /**/
	210	+ ibr = ibr + 1;
	211	+ }
	212	+ }
	213	+ }
	214	+ }
	215	+} /* bragg_vector */
	216	+ /**
	217	+ * Max and Minimum in Brillouine zone
	218	+ */
	219	+void max_and_min_bz() {
	220	+ int ib, i0, i1, i2;
	221	+ GLfloat eigmin, eigmax, matmin, matmax;
	222	+ /**/
	223	+ printf("\n##### Max. and Min. of each bands ##### \n\n");
	224	+ printf("Band Eig_Min. Eig_Max Mat_Min Mat_Max \n");
	225	+ for (ib = 0; ib < nb; ib++) {
	226	+ eigmax = -100000000.0000;
	227	+ eigmin = 100000000.0000;
	228	+ matmax = -100000000.0000;
	229	+ matmin = 100000000.0000;
	230	+ for (i0 = 0; i0 < ng0[0]; ++i0) {
	231	+ for (i1 = 0; i1 < ng0[1]; ++i1) {
	232	+ for (i2 = 0; i2 < ng0[2]; ++i2) {
	233	+ if (eig0[ib][i0][i1][i2] > eigmax) eigmax = eig0[ib][i0][i1][i2];
	234	+ if (eig0[ib][i0][i1][i2] < eigmin) eigmin = eig0[ib][i0][i1][i2];
	235	+ if (mat0[ib][i0][i1][i2] > matmax) matmax = mat0[ib][i0][i1][i2];
	236	+ if (mat0[ib][i0][i1][i2] < matmin) matmin = mat0[ib][i0][i1][i2];
	237	+ }
	238	+ }
	239	+ }
	240	+ printf("%d %f %f %f %f \n", ib + 1, eigmin, eigmax, matmin, matmax);
	241	+ }
	242	+ /**/
	243	+}/* max_and_min_bz */
	244	+
	245	+void initialize_val() {
	246	+ blackback = 1;
	247	+ fcscl = 1;
	248	+ fbz = 1;
	249	+ nodeline = 0;
	250	+ lcolorbar = 1;
	251	+ lstereo = 1;
	252	+ lmouse = 1;
	253	+ itet = 0;
	254	+ scl = 1.0;
	255	+ trans[0] = 0.0; trans[1] = 0.0; trans[2] = 0.0;
	256	+ rot[0][0] = 1.0; rot[0][1] = 0.0; rot[0][2] = 0.0;
	257	+ rot[1][0] = 0.0; rot[1][1] = 1.0; rot[1][2] = 0.0;
	258	+ rot[2][0] = 0.0; rot[2][1] = 0.0; rot[2][2] = 1.0;
	259	+ thetax = 0.0, thetay = 0.0, thetaz = 0.0;
	260	+ /**
	261	+ * Colors
	262	+ */
	263	+ black[0] = 0.0; black[1] = 0.0; black[2] = 0.0; black[3] = 1.0;
	264	+ white[0] = 1.0; white[1] = 1.0; white[2] = 1.0; white[3] = 1.0;
	265	+ cyan[0] = 0.0; cyan[1] = 1.0; cyan[2] = 1.0; cyan[3] = 1.0;
	266	+ magenta[0] = 1.0; magenta[1] = 0.0; magenta[2] = 1.0; magenta[3] = 1.0;
	267	+ yellow[0] = 1.0; yellow[1] = 1.0; yellow[2] = 0.0; yellow[3] = 1.0;
	268	+ red[0] = 1.0; red[1] = 0.0; red[2] = 0.0; red[3] = 1.0;
	269	+ green[0] = 0.0; green[1] = 1.0; green[2] = 0.0; green[3] = 1.0;
	270	+ blue[0] = 0.0; blue[1] = 0.0; blue[2] = 1.0; blue[3] = 1.0;
	271	+ EF = 0.0;
	272	+ interpol = 1;
	273	+}

--- /dev/null

+++ b/src/initialize.h

		@@ -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 init_corner();
	26	+void bragg_vector();
	27	+void max_and_min_bz();
	28	+void initialize_val();

--- /dev/null

+++ b/src/kumo.c

		@@ -0,0 +1,148 @@
	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	+#include <stdlib.h>
	26	+#include "basic_math.h"
	27	+#include "variable.h"
	28	+
	29	+#if defined(MAC)
	30	+#include <GLUT/glut.h>
	31	+#else
	32	+#include <GL/glut.h>
	33	+#endif
	34	+
	35	+/**
	36	+*Compute coefficient for the Kumo interpolation
	37	+*/
	38	+void kumo_coef(int j, GLfloat *coef) {
	39	+ GLfloat x, mx;
	40	+ x = (GLfloat)j / (GLfloat)interpol;
	41	+ mx = 1.0 - x;
	42	+ coef[0] = -0.5 * x * mx * mx;
	43	+ coef[1] = mx * (mxmx + 3.0xmx + 0.5x*x);
	44	+ coef[2] = x * (xx + 3.0mxx + 0.5mx*mx);
	45	+ coef[3] = -0.5 * x * x * mx;
	46	+}
	47	+/**
	48	+* Interpolation of energy and matrix
	49	+*/
	50	+void interpol_energy() {
	51	+ int ib, i0, i1, ii, i2, j0, j1, j2;
	52	+ GLfloat coef[4];
	53	+ GLfloat eig1[4][4][4], mat1[4][4][4], eig2[4][4], mat2[4][4], eig3[4], mat3[4];
	54	+ /*
	55	+ Reallocate
	56	+ */
	57	+ for (ib = 0; ib < nb; ib++) {
	58	+ for (i0 = 0; i0 < ng[0]; i0++) {
	59	+ for (i1 = 0; i1 < ng[1]; i1++) {
	60	+ free(eig[ib][i0][i1]);
	61	+ free(mat[ib][i0][i1]);
	62	+ }
	63	+ free(eig[ib][i0]);
	64	+ free(mat[ib][i0]);
	65	+ }
	66	+ free(eig[ib]);
	67	+ free(mat[ib]);
	68	+ }
	69	+ for (ii = 0; ii < 3; ii++)ng[ii] = ng0[ii] * interpol;
	70	+ /**/
	71	+ for (ib = 0; ib < nb; ib++) {
	72	+ eig[ib] = (GLfloat**)malloc(ng[0] sizeof(GLfloat**));
	73	+ mat[ib] = (GLfloat**)malloc(ng[0] sizeof(GLfloat**));
	74	+ for (i0 = 0; i0 < ng[0]; i0++) {
	75	+ eig[ib][i0] = (GLfloat*)malloc(ng[1] sizeof(GLfloat*));
	76	+ mat[ib][i0] = (GLfloat*)malloc(ng[1] sizeof(GLfloat*));
	77	+ for (i1 = 0; i1 < ng[1]; i1++) {
	78	+ eig[ib][i0][i1] = (GLfloat)malloc(ng[2] sizeof(GLfloat));
	79	+ mat[ib][i0][i1] = (GLfloat)malloc(ng[2] sizeof(GLfloat));
	80	+ }
	81	+ }
	82	+ }
	83	+ /*
	84	+ 3rd order - three dimensional Kumo interpolation
	85	+ */
	86	+ for (ib = 0; ib < nb; ib++) {
	87	+ for (i0 = 0; i0 < ng0[0]; i0++) {
	88	+ for (i1 = 0; i1 < ng0[1]; i1++) {
	89	+ for (i2 = 0; i2 < ng0[2]; i2++) {
	90	+ for (j0 = 0; j0 < 4; j0++) {
	91	+ for (j1 = 0; j1 < 4; j1++) {
	92	+ for (j2 = 0; j2 < 4; j2++) {
	93	+ eig1[j0][j1][j2] = eig0[ib][modulo(i0 + j0 - 1, ng0[0])]
	94	+ [modulo(i1 + j1 - 1, ng0[1])]
	95	+ [modulo(i2 + j2 - 1, ng0[2])];
	96	+ mat1[j0][j1][j2] = mat0[ib][modulo(i0 + j0 - 1, ng0[0])]
	97	+ [modulo(i1 + j1 - 1, ng0[1])]
	98	+ [modulo(i2 + j2 - 1, ng0[2])];
	99	+ }
	100	+ }
	101	+ }
	102	+
	103	+ for (j0 = 0; j0 < interpol; j0++) {
	104	+ kumo_coef(j0, &coef[0]);
	105	+ for (j1 = 0; j1 < 4; j1++) {
	106	+ for (j2 = 0; j2 < 4; j2++) {
	107	+ eig2[j1][j2] = 0.0;
	108	+ mat2[j1][j2] = 0.0;
	109	+ for (ii = 0; ii < 4; ii++) {
	110	+ eig2[j1][j2] += coef[ii] * eig1[ii][j1][j2];
	111	+ mat2[j1][j2] += coef[ii] * mat1[ii][j1][j2];
	112	+ }
	113	+ }
	114	+ }
	115	+ for (j1 = 0; j1 < interpol; j1++) {
	116	+ kumo_coef(j1, &coef[0]);
	117	+ for (j2 = 0; j2 < 4; j2++) {
	118	+ eig3[j2] = 0.0;
	119	+ mat3[j2] = 0.0;
	120	+ for (ii = 0; ii < 4; ii++) {
	121	+ eig3[j2] += coef[ii] * eig2[ii][j2];
	122	+ mat3[j2] += coef[ii] * mat2[ii][j2];
	123	+ }
	124	+ }
	125	+ for (j2 = 0; j2 < interpol; j2++) {
	126	+ kumo_coef(j2, &coef[0]);
	127	+ eig[ib][i0*interpol + j0]
	128	+ [i1*interpol + j1]
	129	+ [i2*interpol + j2] = 0.0;
	130	+ mat[ib][i0*interpol + j0]
	131	+ [i1*interpol + j1]
	132	+ [i2*interpol + j2] = 0.0;
	133	+ for (ii = 0; ii < 4; ii++) {
	134	+ eig[ib][i0*interpol + j0]
	135	+ [i1*interpol + j1]
	136	+ [i2interpol + j2] += coef[ii] eig3[ii];
	137	+ mat[ib][i0*interpol + j0]
	138	+ [i1*interpol + j1]
	139	+ [i2interpol + j2] += coef[ii] mat3[ii];
	140	+ }
	141	+ }
	142	+ }
	143	+ }
	144	+ }
	145	+ }
	146	+ }
	147	+ }
	148	+}/void interpol_energy() /

--- /dev/null

+++ b/src/kumo.h

		@@ -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();

--- /dev/null

+++ b/src/menu.c

		@@ -0,0 +1,618 @@
	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	+#include <stdlib.h>
	26	+#include <stdio.h>
	27	+#include <math.h>
	28	+#include "variable.h"
	29	+#include "free_patch.h"
	30	+#include "fermi_patch.h"
	31	+#include "calc_nodeline.h"
	32	+#include "kumo.h"
	33	+
	34	+#if defined(MAC)
	35	+#include <GLUT/glut.h>
	36	+#else
	37	+#include <GL/glut.h>
	38	+#endif
	39	+
	40	+/**
	41	+ * Main menu
	42	+ */
	43	+void main_menu(int value /*< [in] Selected menu/) {
	44	+ int i0, i1, ib, ibzl, i;
	45	+ /**/
	46	+ if (value == 9) {
	47	+ /*
	48	+ Exit
	49	+ */
	50	+ printf("\nExit. \n\n");
	51	+ free_patch();
	52	+ for (ib = 0; ib < nb; ib++) {
	53	+ for (i0 = 0; i0 < ng[0]; i0++) {
	54	+ for (i1 = 0; i1 < ng[1]; i1++) {
	55	+ free(eig[ib][i0][i1]);
	56	+ free(mat[ib][i0][i1]);
	57	+ }
	58	+ free(eig[ib][i0]);
	59	+ free(mat[ib][i0]);
	60	+ }
	61	+ free(eig[ib]);
	62	+ free(mat[ib]);
	63	+ }
	64	+ free(eig);
	65	+ free(mat);
	66	+ free(ntri);
	67	+ free(draw_band);
	68	+ for (ibzl = 0; ibzl < nbzl; ++ibzl) {
	69	+ for (i = 0; i < 2; ++i) {
	70	+ free(bzl[ibzl][i]);
	71	+ }
	72	+ free(bzl[ibzl]);
	73	+ }
	74	+ free(bzl);
	75	+ free(nnl);
	76	+ exit(0);
	77	+ }
	78	+}
	79	+/*
	80	+ * Shift Fermi energy
	81	+ */
	82	+void menu_shiftEF(int value /*< [in] Selected menu/)
	83	+{
	84	+ int ib, i0, i1, i2, ierr;
	85	+ GLfloat emin, emax;
	86	+
	87	+ if (value == 1) {
	88	+ emin = 100000.0;
	89	+ emax = -100000.0;
	90	+ for (ib = 0; ib < nb; ++ib) {
	91	+ for (i0 = 0; i0 < ng[0]; ++i0) {
	92	+ for (i1 = 0; i1 < ng[1]; ++i1) {
	93	+ for (i2 = 0; i2 < ng[2]; ++i2) {
	94	+ if (emin > eig[ib][i0][i1][i2]) emin = eig[ib][i0][i1][i2];
	95	+ if (emax < eig[ib][i0][i1][i2]) emax = eig[ib][i0][i1][i2];
	96	+ }
	97	+ }
	98	+ }
	99	+ }
	100	+ printf("Min Max E_F \n");
	101	+ printf("%f %f %f \n", emin, emax, EF);
	102	+ printf("New Fermi energy : ");
	103	+ //
	104	+ ierr = scanf("%f", &EF);
	105	+ if (ierr != 1) printf("error ! reading ef");
	106	+ /**/
	107	+ free_patch();
	108	+ query = 1;
	109	+ fermi_patch();
	110	+ query = 0;
	111	+ fermi_patch();
	112	+ calc_nodeline();
	113	+ /**/
	114	+ max_and_min();
	115	+ /**/
	116	+ glutPostRedisplay();
	117	+ }
	118	+}
	119	+/*
	120	+ * Shift Fermi energy
	121	+ */
	122	+void menu_interpol(int value /*< [in] Selected menu/)
	123	+{
	124	+ int ib, i0, i1, i2, ierr;
	125	+ GLfloat emin, emax;
	126	+
	127	+ if (value == 1) {
	128	+ printf("Old interpolation ratio : %d\n", interpol);
	129	+ printf("New interpolation ratio : ");
	130	+ //
	131	+ ierr = scanf("%d", &interpol);
	132	+ if (ierr != 1) printf("error ! reading interpol");
	133	+ /**/
	134	+ interpol_energy();
	135	+ free_patch();
	136	+ query = 1;
	137	+ fermi_patch();
	138	+ query = 0;
	139	+ fermi_patch();
	140	+ calc_nodeline();
	141	+ /**/
	142	+ max_and_min();
	143	+ /**/
	144	+ glutPostRedisplay();
	145	+ }
	146	+}
	147	+/*
	148	+ * Setting of view
	149	+ */
	150	+void menu_view(int value /*< [in] Selected menu/)
	151	+{
	152	+ int ierr;
	153	+
	154	+ if (value == 1) {
	155	+
	156	+ printf(" Current Scale : %f\n", scl);
	157	+ printf(" New Scale : ");
	158	+ ierr = scanf("%f", &scl);
	159	+
	160	+ }
	161	+ else if (value == 2) {
	162	+
	163	+ printf(" Current Position(x y) : %f %f\n", trans[0], trans[1]);
	164	+ printf(" New Position(x y) : ");
	165	+ ierr = scanf("%f %f", &trans[0], &trans[1]);
	166	+
	167	+ }
	168	+ else if (value == 3) {
	169	+
	170	+ /**/
	171	+ thetay = asin(rot[0][2]);
	172	+ if (cosf(thetay) != 0.0) {
	173	+ if (-rot[1][2] / cosf(thetay) >= 0.0) thetax = acosf(rot[2][2] / cosf(thetay));
	174	+ else thetax = 6.283185307 - acosf(rot[2][2] / cosf(thetay));
	175	+ if (-rot[0][1] / cosf(thetay) >= 0.0) thetaz = acosf(rot[0][0] / cosf(thetay));
	176	+ else thetaz = 6.283185307 - acosf(rot[0][0] / cosf(thetay));
	177	+ }
	178	+ else {
	179	+ thetax = 0.0;
	180	+ if (rot[1][0] >= 0.0) thetaz = acosf(rot[1][1]);
	181	+ else thetaz = 6.283185307 - acosf(rot[1][1]);
	182	+ }
	183	+ thetax = 180.0 / 3.14159265 * thetax;
	184	+ thetay = 180.0 / 3.14159265 * thetay;
	185	+ thetaz = 180.0 / 3.14159265 * thetaz;
	186	+ printf(" Current Rotation (theta_x theta_y teta_z) in degree : %f %f %f\n", thetax, thetay, thetaz);
	187	+ printf(" New Rotation (theta_x theta_y teta_z) in degree : ");
	188	+ ierr = scanf("%f %f %f", &thetax, &thetay, &thetaz);
	189	+ thetax = 3.14159265 / 180.0 * thetax;
	190	+ thetay = 3.14159265 / 180.0 * thetay;
	191	+ thetaz = 3.14159265 / 180.0 * thetaz;
	192	+
	193	+ rot[0][0] = cosf(thetay)* cosf(thetaz);
	194	+ rot[0][1] = -cosf(thetay)* sin(thetaz);
	195	+ rot[0][2] = sinf(thetay);
	196	+ rot[1][0] = cosf(thetaz)* sinf(thetax)* sinf(thetay) + cosf(thetax)* sinf(thetaz);
	197	+ rot[1][1] = cosf(thetax) * cosf(thetaz) - sinf(thetax)* sinf(thetay)* sinf(thetaz);
	198	+ rot[1][2] = -cosf(thetay)* sinf(thetax);
	199	+ rot[2][0] = -cosf(thetax)* cosf(thetaz)* sinf(thetay) + sinf(thetax)* sinf(thetaz);
	200	+ rot[2][1] = cosf(thetaz)* sinf(thetax) + cosf(thetax)* sinf(thetay)* sinf(thetaz);
	201	+ rot[2][2] = cos(thetax)* cosf(thetay);
	202	+
	203	+ }
	204	+
	205	+ glutPostRedisplay();
	206	+
	207	+}
	208	+/**
	209	+ * Change mouse function
	210	+ */
	211	+void menu_mouse(int value /*< [in] Selected menu/) {
	212	+ /**/
	213	+ if (value == 1 && lmouse != 1) {
	214	+ lmouse = 1;
	215	+ glutPostRedisplay();
	216	+ }
	217	+ if (value == 2 && lmouse != 2) {
	218	+ lmouse = 2;
	219	+ glutPostRedisplay();
	220	+ }
	221	+ if (value == 3 && lmouse != 3) {
	222	+ lmouse = 3;
	223	+ glutPostRedisplay();
	224	+ }
	225	+ /**/
	226	+} /* menu_band */
	227	+/**
	228	+ * On / Off band
	229	+ */
	230	+void menu_band(int value /*< [in] Selected menu/) {
	231	+ /**/
	232	+ if (draw_band[value] == 0) {
	233	+ draw_band[value] = 1;
	234	+ }
	235	+ else {
	236	+ draw_band[value] = 0;
	237	+ }
	238	+ glutPostRedisplay();
	239	+ /**/
	240	+} /* menu_band */
	241	+/**
	242	+ * Change background color
	243	+ */
	244	+void menu_bgcolor(int value /*<[in] Selected menu/) {
	245	+ if (value == 1 && blackback != 1) {
	246	+ glClearColor(0.0, 0.0, 0.0, 0.0);
	247	+ blackback = 1;
	248	+ glutPostRedisplay();
	249	+ }
	250	+ else if (value == 2 && blackback != 0) {
	251	+ glClearColor(1.0, 1.0, 1.0, 0.0);
	252	+ blackback = 0;
	253	+ glutPostRedisplay();
	254	+ }
	255	+ /**/
	256	+}/* bgcolor change*/
	257	+/**
	258	+ * Change color scale mode
	259	+ */
	260	+void menu_colorscale(int value /*<[in] Selected menu/) {
	261	+ /**/
	262	+ if (value == 1 && fcscl != 1) {
	263	+ fcscl = 1;
	264	+ max_and_min();
	265	+ glutPostRedisplay();
	266	+ }
	267	+ else if (value == 2 && fcscl != 2) {
	268	+ fcscl = 2;
	269	+ max_and_min();
	270	+ glutPostRedisplay();
	271	+ }
	272	+ else if (value == 3 && fcscl != 3) {
	273	+ fcscl = 3;
	274	+ max_and_min();
	275	+ glutPostRedisplay();
	276	+ }
	277	+ else if (value == 4 && fcscl != 4) {
	278	+ fcscl = 4;
	279	+ max_and_min();
	280	+ glutPostRedisplay();
	281	+ }
	282	+ /**/
	283	+} /* menu_colorscale */
	284	+/**
	285	+ * Change Brillouin zone
	286	+ */
	287	+void menu_bzmode(int value /*<[in] Selected menu/) {
	288	+ if (value == 1 && fbz != 1) {
	289	+ fbz = 1;
	290	+ /**/
	291	+ free_patch();
	292	+ query = 1;
	293	+ fermi_patch();
	294	+ query = 0;
	295	+ fermi_patch();
	296	+ calc_nodeline();
	297	+ /**/
	298	+ max_and_min();
	299	+ /**/
	300	+ glutPostRedisplay();
	301	+ }
	302	+ else if (value == 2 && fbz != -1) {
	303	+ fbz = -1;
	304	+ /**/
	305	+ free_patch();
	306	+ query = 1;
	307	+ fermi_patch();
	308	+ query = 0;
	309	+ fermi_patch();
	310	+ calc_nodeline();
	311	+ /**/
	312	+ max_and_min();
	313	+ /**/
	314	+ glutPostRedisplay();
	315	+ }
	316	+} /* menu_bzmode */
	317	+/**
	318	+ * On/Off Node line
	319	+ */
	320	+void menu_nodeline(int value /*<[in] Selected menu/) {
	321	+ if (value == 1 && nodeline != 1) {
	322	+ nodeline = 1;
	323	+ glutPostRedisplay();
	324	+ }
	325	+ else if (value == 2 && nodeline != 0) {
	326	+ nodeline = 0;
	327	+ glutPostRedisplay();
	328	+ }
	329	+ /**/
	330	+} /* menu_nodeline */
	331	+/**
	332	+ * Tern stereogram
	333	+ */
	334	+void menu_stereo(int value /*<[in] Selected menu/) {
	335	+ if (value == 1 && lstereo != 1) {
	336	+ lstereo = 1;
	337	+ glutPostRedisplay();
	338	+ }
	339	+ if (value == 2 && lstereo != 2) {
	340	+ lstereo = 2;
	341	+ glutPostRedisplay();
	342	+ }
	343	+ if (value == 3 && lstereo != 3) {
	344	+ lstereo = 3;
	345	+ glutPostRedisplay();
	346	+ }
	347	+} /* menu_stereo */
	348	+/**
	349	+ * On/Off Colorbar
	350	+ */
	351	+void menu_colorbar(int value /*<[in] Selected menu/) {
	352	+ if (value == 1 && lcolorbar != 1) {
	353	+ lcolorbar = 1;
	354	+ glutPostRedisplay();
	355	+ }
	356	+ if (value == 2 && lcolorbar != 0) {
	357	+ lcolorbar = 0;
	358	+ glutPostRedisplay();
	359	+ }
	360	+} /* menu_colorbar */
	361	+/**
	362	+ * Change tetrahedron
	363	+ */
	364	+void menu_tetra(int value) /*<[in] Selected menu/ {
	365	+ /**/
	366	+ if (value != itet) {
	367	+ printf("Tetra patern %d \n", value + 1);
	368	+ itet = value;
	369	+ init_corner();
	370	+ free_patch();
	371	+ query = 1;
	372	+ fermi_patch();
	373	+ query = 0;
	374	+ fermi_patch();
	375	+ calc_nodeline();
	376	+ max_and_min();
	377	+ glutPostRedisplay();
	378	+ }
	379	+} /* menu_tetra */
	380	+/**
	381	+ * Munu
	382	+ */
	383	+void FS_ModifyMenu(int status)
	384	+{
	385	+ int ib;
	386	+ char menu_str[20] = { 0 };
	387	+
	388	+ if (status == GLUT_MENU_IN_USE) {
	389	+ glutPostRedisplay();
	390	+ }
	391	+ else {
	392	+ /**
	393	+ * Operation with mouse drag
	394	+ */
	395	+ glutSetMenu(imenu_mouse);
	396	+ for (ib = 0; ib < 3; ib++) glutRemoveMenuItem(1);
	397	+ if (lmouse == 1) glutAddMenuEntry("[x] Rotate", 1);
	398	+ else glutAddMenuEntry("[ ] Rotate", 1);
	399	+ if (lmouse == 2) glutAddMenuEntry("[x] Scale", 2);
	400	+ else glutAddMenuEntry("[ ] Scale", 2);
	401	+ if (lmouse == 3) glutAddMenuEntry("[x] Translate", 3);
	402	+ else glutAddMenuEntry("[ ] Translate", 3);
	403	+ /**
	404	+ * Band menu
	405	+ */
	406	+ glutSetMenu(imenu_band);
	407	+ for (ib = 0; ib < nb; ib++) glutRemoveMenuItem(1);
	408	+ for (ib = 0; ib < nb; ib++) {
	409	+ if (draw_band[ib] == 1) sprintf(menu_str, "[x] band # %d", ib + 1);
	410	+ else sprintf(menu_str, "[ ] band # %d", ib + 1);
	411	+ glutAddMenuEntry(menu_str, ib);
	412	+ }
	413	+ /**
	414	+ * Operation with mouse drag
	415	+ */
	416	+ glutSetMenu(imenu_interpol);
	417	+ glutRemoveMenuItem(1);
	418	+ sprintf(menu_str, "Ratio : %d", interpol);
	419	+ glutAddMenuEntry(menu_str, 1);
	420	+ /**
	421	+ * Background color
	422	+ */
	423	+ glutSetMenu(imenu_bgcolor);
	424	+ for (ib = 0; ib < 2; ib++) glutRemoveMenuItem(1);
	425	+ if (blackback == 1) glutAddMenuEntry("[x] Black", 1);
	426	+ else glutAddMenuEntry("[ ] Black", 1);
	427	+ if (blackback == 0) glutAddMenuEntry("[x] White", 2);
	428	+ else glutAddMenuEntry("[ ] White", 2);
	429	+ /**
	430	+ * Color scale mode
	431	+ */
	432	+ glutSetMenu(imenu_colorscale);
	433	+ for (ib = 0; ib < 4; ib++) glutRemoveMenuItem(1);
	434	+ if (fcscl == 1) glutAddMenuEntry("[x] Auto", 1);
	435	+ else glutAddMenuEntry("[ ] Auto", 1);
	436	+ if (fcscl == 2) glutAddMenuEntry("[x] Manual", 2);
	437	+ else glutAddMenuEntry("[ ] Manual", 2);
	438	+ if (fcscl == 3) glutAddMenuEntry("[x] Unicolor", 3);
	439	+ else glutAddMenuEntry("[ ] Unicolor", 3);
	440	+ if (fcscl == 4) glutAddMenuEntry("[x] Periodic", 4);
	441	+ else glutAddMenuEntry("[ ] Periodic", 4);
	442	+ /**
	443	+ * Brillouin zone
	444	+ */
	445	+ glutSetMenu(imenu_bzmode);
	446	+ for (ib = 0; ib < 2; ib++) glutRemoveMenuItem(1);
	447	+ if (fbz == 1) glutAddMenuEntry("[x] First Brillouin zone", 1);
	448	+ else glutAddMenuEntry("[ ] First Brillouin zone", 1);
	449	+ if (fbz == -1) glutAddMenuEntry("[x] Primitive Brillouin zone", 2);
	450	+ else glutAddMenuEntry("[ ] Primitive Brillouin zone", 2);
	451	+ /**
	452	+ * Nodeline on/off
	453	+ */
	454	+ glutSetMenu(imenu_nodeline);
	455	+ for (ib = 0; ib < 2; ib++) glutRemoveMenuItem(1);
	456	+ if (nodeline == 1) glutAddMenuEntry("[x] On", 1);
	457	+ else glutAddMenuEntry("[ ] On", 1);
	458	+ if (nodeline == 0) glutAddMenuEntry("[x] Off", 2);
	459	+ else glutAddMenuEntry("[ ] Off", 2);
	460	+ /**
	461	+ * Colorbar on/off
	462	+ */
	463	+ glutSetMenu(imenu_colorbar);
	464	+ for (ib = 0; ib < 2; ib++) glutRemoveMenuItem(1);
	465	+ if (lcolorbar == 1) glutAddMenuEntry("[x] On", 1);
	466	+ else glutAddMenuEntry("[ ] On", 1);
	467	+ if (lcolorbar == 0) glutAddMenuEntry("[x] Off", 2);
	468	+ else glutAddMenuEntry("[ ] Off", 2);
	469	+ /**
	470	+ * Stereogram
	471	+ */
	472	+ glutSetMenu(imenu_stereo);
	473	+ for (ib = 0; ib < 3; ib++) glutRemoveMenuItem(1);
	474	+ if (lstereo == 1) glutAddMenuEntry("[x] None", 1);
	475	+ else glutAddMenuEntry("[ ] None", 1);
	476	+ if (lstereo == 2) glutAddMenuEntry("[x] Parallel", 2);
	477	+ else glutAddMenuEntry("[ ] Parallel", 2);
	478	+ if (lstereo == 3) glutAddMenuEntry("[x] Cross", 3);
	479	+ else glutAddMenuEntry("[ ] Cross", 3);
	480	+ /**
	481	+ * Tetrahedron
	482	+ */
	483	+ glutSetMenu(imenu_tetra);
	484	+ for (ib = 0; ib < 16; ib++) glutRemoveMenuItem(1);
	485	+ for (ib = 0; ib < 16; ib++) {
	486	+ if (itet == ib) sprintf(menu_str, "[x] tetra # %d", ib + 1);
	487	+ else sprintf(menu_str, "[ ] tetra # %d", ib + 1);
	488	+ glutAddMenuEntry(menu_str, ib);
	489	+ }
	490	+ glutPostRedisplay();
	491	+ }
	492	+}
	493	+/**
	494	+ * Munu
	495	+ */
	496	+void FS_CreateMenu()
	497	+{
	498	+ int ib, imenu_shiftEF, imenu_view;
	499	+ char menu_str[20] = { 0 };
	500	+ /**
	501	+ * Mouse drag works as ...
	502	+ */
	503	+ imenu_mouse = glutCreateMenu(menu_mouse);
	504	+ if (lmouse == 1) glutAddMenuEntry("[x] Rotate", 1);
	505	+ else glutAddMenuEntry("[ ] Rotate", 1);
	506	+ if (lmouse == 2) glutAddMenuEntry("[x] Scale", 2);
	507	+ else glutAddMenuEntry("[ ] Scale", 2);
	508	+ if (lmouse == 3) glutAddMenuEntry("[x] Translate", 3);
	509	+ else glutAddMenuEntry("[ ] Translate", 3);
	510	+ /**
	511	+ * On/Off each band
	512	+ */
	513	+ imenu_band = glutCreateMenu(menu_band);
	514	+ for (ib = 0; ib < nb; ib++) {
	515	+ if (draw_band[ib] == 1) sprintf(menu_str, "[x] band # %d", ib + 1);
	516	+ else sprintf(menu_str, "[ ] band # %d", ib + 1);
	517	+ glutAddMenuEntry(menu_str, ib);
	518	+ }
	519	+ /**
	520	+ * Shift Fermi energy
	521	+ */
	522	+ imenu_shiftEF = glutCreateMenu(menu_shiftEF);
	523	+ glutAddMenuEntry("Shift Fermi energy", 1);
	524	+ /**
	525	+ * Modify interpolation ratio
	526	+ */
	527	+ sprintf(menu_str, "Ratio : %d", interpol);
	528	+ imenu_interpol = glutCreateMenu(menu_interpol);
	529	+ glutAddMenuEntry(menu_str, 1);
	530	+ /**
	531	+ * Set view
	532	+ */
	533	+ imenu_view = glutCreateMenu(menu_view);
	534	+ glutAddMenuEntry("Scale", 1);
	535	+ glutAddMenuEntry("Position", 2);
	536	+ glutAddMenuEntry("Rotation", 3);
	537	+ /**
	538	+ * Background color
	539	+ */
	540	+ imenu_bgcolor = glutCreateMenu(menu_bgcolor);
	541	+ if (blackback == 1) glutAddMenuEntry("[x] Black", 1);
	542	+ else glutAddMenuEntry("[ ] Black", 1);
	543	+ if (blackback == 0) glutAddMenuEntry("[x] White", 2);
	544	+ else glutAddMenuEntry("[ ] White", 2);
	545	+ /**
	546	+ * Color scale mode
	547	+ */
	548	+ imenu_colorscale = glutCreateMenu(menu_colorscale);
	549	+ if (fcscl == 1) glutAddMenuEntry("[x] Auto", 1);
	550	+ else glutAddMenuEntry("[ ] Auto", 1);
	551	+ if (fcscl == 2) glutAddMenuEntry("[x] Manual", 2);
	552	+ else glutAddMenuEntry("[ ] Manual", 2);
	553	+ if (fcscl == 3) glutAddMenuEntry("[x] Unicolor", 3);
	554	+ else glutAddMenuEntry("[ ] Unicolor", 3);
	555	+ if (fcscl == 4) glutAddMenuEntry("[x] Periodic", 4);
	556	+ else glutAddMenuEntry("[ ] Periodic", 4);
	557	+ /**
	558	+ * Brillouin zone
	559	+ */
	560	+ imenu_bzmode = glutCreateMenu(menu_bzmode);
	561	+ if (fbz == 1) glutAddMenuEntry("[x] First Brillouin zone", 1);
	562	+ else glutAddMenuEntry("[ ] First Brillouin zone", 1);
	563	+ if (fbz == -1) glutAddMenuEntry("[x] Primitive Brillouin zone", 2);
	564	+ else glutAddMenuEntry("[ ] Primitive Brillouin zone", 2);
	565	+ /**
	566	+ * Nodeline on/off
	567	+ */
	568	+ imenu_nodeline = glutCreateMenu(menu_nodeline);
	569	+ if (nodeline == 1) glutAddMenuEntry("[x] On", 1);
	570	+ else glutAddMenuEntry("[ ] On", 1);
	571	+ if (nodeline == 0) glutAddMenuEntry("[x] Off", 2);
	572	+ else glutAddMenuEntry("[ ] Off", 2);
	573	+ /**
	574	+ * Colorbar on/off
	575	+ */
	576	+ imenu_colorbar = glutCreateMenu(menu_colorbar);
	577	+ if (lcolorbar == 1) glutAddMenuEntry("[x] On", 1);
	578	+ else glutAddMenuEntry("[ ] On", 1);
	579	+ if (lcolorbar == 0) glutAddMenuEntry("[x] Off", 2);
	580	+ else glutAddMenuEntry("[ ] Off", 2);
	581	+ /**
	582	+ * Stereogram
	583	+ */
	584	+ imenu_stereo = glutCreateMenu(menu_stereo);
	585	+ if (lstereo == 1) glutAddMenuEntry("[x] None", 1);
	586	+ else glutAddMenuEntry("[ ] None", 1);
	587	+ if (lstereo == 2) glutAddMenuEntry("[x] Parallel", 2);
	588	+ else glutAddMenuEntry("[ ] Parallel", 2);
	589	+ if (lstereo == 3) glutAddMenuEntry("[x] Cross", 3);
	590	+ else glutAddMenuEntry("[ ] Cross", 3);
	591	+ /**
	592	+ * Tetrahedron
	593	+ */
	594	+ imenu_tetra = glutCreateMenu(menu_tetra);
	595	+ for (ib = 0; ib < 16; ib++) {
	596	+ if (itet == ib) sprintf(menu_str, "[x] tetra # %d", ib + 1);
	597	+ else sprintf(menu_str, "[ ] tetra # %d", ib + 1);
	598	+ glutAddMenuEntry(menu_str, ib);
	599	+ }
	600	+ /**
	601	+ * Main menu
	602	+ */
	603	+ imenu = glutCreateMenu(main_menu);
	604	+ glutAddSubMenu("Band", imenu_band);
	605	+ glutAddSubMenu("Mouse Drag", imenu_mouse);
	606	+ glutAddSubMenu("Shift Fermi energy", imenu_shiftEF);
	607	+ glutAddSubMenu("Interpolation", imenu_interpol);
	608	+ glutAddSubMenu("Set view", imenu_view);
	609	+ glutAddSubMenu("Background color", imenu_bgcolor);
	610	+ glutAddSubMenu("Color scale mode", imenu_colorscale);
	611	+ glutAddSubMenu("Brillouin zone", imenu_bzmode);
	612	+ glutAddSubMenu("Node lines", imenu_nodeline);
	613	+ glutAddSubMenu("Color bar On/Off", imenu_colorbar);
	614	+ glutAddSubMenu("Stereogram", imenu_stereo);
	615	+ glutAddSubMenu("Tetrahedron", imenu_tetra);
	616	+ glutAddMenuEntry("Exit", 9);
	617	+ glutAttachMenu(GLUT_RIGHT_BUTTON);
	618	+}

--- /dev/null

+++ b/src/menu.h

		@@ -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 FS_CreateMenu();
	26	+void FS_ModifyMenu(int status);

--- /dev/null

+++ b/src/operation.c

		@@ -0,0 +1,234 @@
	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	+#include <math.h>
	26	+#include "variable.h"
	27	+
	28	+#if defined(MAC)
	29	+#include <GLUT/glut.h>
	30	+#else
	31	+#include <GL/glut.h>
	32	+#endif
	33	+
	34	+/**
	35	+ * Window resize
	36	+ */
	37	+void resize(
	38	+ int w /*<[in] Window width/,
	39	+ int h /*<[in] Window height/)
	40	+{
	41	+ /*
	42	+ * Scale of translation of mousepointer
	43	+ */
	44	+ sx = 1.0 / (GLfloat)w;
	45	+ sy = 1.0 / (GLfloat)h;
	46	+ /**/
	47	+ glViewport(0, 0, w, h);
	48	+ /**/
	49	+ glMatrixMode(GL_PROJECTION);
	50	+ /**/
	51	+ glLoadIdentity();
	52	+ gluPerspective(30.0, (GLfloat)w / (GLfloat)h, 1.0, 100.0);
	53	+ /**/
	54	+ glMatrixMode(GL_MODELVIEW);
	55	+} /* end resize */
	56	+/**
	57	+ * Idling
	58	+ */
	59	+void idle(void)
	60	+{
	61	+ glutPostRedisplay();
	62	+} /* idle */
	63	+/**
	64	+ * Glut mouse function
	65	+ */
	66	+void mouse(
	67	+ int button /*< [in] pushed button/,
	68	+ int state /*< [in] down or up or ?/,
	69	+ int x /*< [in] position of mouse cursor/,
	70	+ int y /*< [in] position of mouse cursor/)
	71	+{
	72	+ switch (button) {
	73	+ /*
	74	+ * Drag
	75	+ */
	76	+ case GLUT_LEFT_BUTTON:
	77	+ switch (state) {
	78	+ case GLUT_DOWN:
	79	+ /* Start of animation */
	80	+ glutIdleFunc(idle);
	81	+ /* Record drag start point */
	82	+ cx = x;
	83	+ cy = y;
	84	+ break;
	85	+ case GLUT_UP:
	86	+ /* End of animation */
	87	+ glutIdleFunc(0);
	88	+ break;
	89	+ default:
	90	+ break;
	91	+ }
	92	+ break;
	93	+ /*
	94	+ Zoom up
	95	+ */
	96	+ case MOUSE_SCROLL_UP:
	97	+ switch (state) {
	98	+ case GLUT_DOWN:
	99	+ scl = scl * 1.1;
	100	+ glutPostRedisplay();
	101	+ break;
	102	+ case GLUT_UP:
	103	+ break;
	104	+ default:
	105	+ break;
	106	+ }
	107	+ break;
	108	+ /*
	109	+ Zoom down
	110	+ */
	111	+ case MOUSE_SCROLL_DOWN:
	112	+ switch (state) {
	113	+ case GLUT_DOWN:
	114	+ scl = scl * 0.9;
	115	+ glutPostRedisplay();
	116	+ break;
	117	+ case GLUT_UP:
	118	+ break;
	119	+ default:
	120	+ break;
	121	+ }
	122	+ break;
	123	+ /*
	124	+ No pushing
	125	+ */
	126	+ default:
	127	+ break;
	128	+ }
	129	+} /* end mouse */
	130	+/**
	131	+ * Glut motion function
	132	+ */
	133	+void motion(
	134	+ int x /*< [in] position of cursor/,
	135	+ int y /*< [in] position of cursor/)
	136	+{
	137	+ int i, j;
	138	+ GLfloat dx, dy, a, rot0[3][3], rot1[3][3], ax, ay;
	139	+ /*
	140	+ Translation of mousepointer from starting point
	141	+ */
	142	+ dx = (x - cx) * sx;
	143	+ dy = (y - cy) * sy;
	144	+ /**/
	145	+ /*
	146	+ Distanse from starting point
	147	+ */
	148	+ a = sqrt(dx * dx + dy * dy);
	149	+ /**/
	150	+ if (lmouse == 1) {
	151	+ /**/
	152	+ if (a != 0.0) {
	153	+ /*
	154	+ * Compute rotational matrix from translation of mousepointer
	155	+ */
	156	+ ax = -dy;
	157	+ ay = dx;
	158	+ /**/
	159	+ a = a * 10.0;
	160	+ /**/
	161	+ rot0[0][0] = (ax * ax + ay * ay * cos(a)) / (ax * ax + ay * ay);
	162	+ rot0[0][1] = ax * ay * (cos(a) - 1.0) / (ax * ax + ay * ay);
	163	+ rot0[0][2] = ay * sin(a) / sqrtf(ax * ax + ay * ay);
	164	+ rot0[1][0] = ax * ay * (cos(a) - 1.0) / (ax * ax + ay * ay);
	165	+ rot0[1][1] = (ax * ax * cos(a) + ay * ay) / (ax * ax + ay * ay);
	166	+ rot0[1][2] = ax * sin(a) / sqrtf(ax * ax + ay * ay);
	167	+ rot0[2][0] = -ay * sin(a) / sqrtf(ax * ax + ay * ay);
	168	+ rot0[2][1] = -ax * sin(a) / sqrtf(ax * ax + ay * ay);
	169	+ rot0[2][2] = cos(a);
	170	+ /**/
	171	+ for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) rot1[i][j] = rot[i][j];
	172	+ /**/
	173	+ for (i = 0; i < 3; i++) {
	174	+ for (j = 0; j < 3; j++) {
	175	+ rot[i][j] = rot0[i][0] * rot1[0][j]
	176	+ + rot0[i][1] * rot1[1][j]
	177	+ + rot0[i][2] * rot1[2][j];
	178	+ }
	179	+ }
	180	+ }
	181	+ }
	182	+ else if (lmouse == 2) {
	183	+ scl = scl * exp(-dy);
	184	+ }
	185	+ else {
	186	+ trans[0] = trans[0] + dx;
	187	+ trans[1] = trans[1] - dy;
	188	+ }
	189	+ cx = x;
	190	+ cy = y;
	191	+ /**/
	192	+} /* motion */
	193	+/*
	194	+ * Glut keyboard function
	195	+ */
	196	+void keyboard(
	197	+ unsigned char key /*< [in] Typed key/,
	198	+ int x /*< [in]/,
	199	+ int y /*< [in]/)
	200	+{
	201	+ switch (key) {
	202	+ }
	203	+} /* keyboard */
	204	+/**
	205	+ * Glut special key function
	206	+ */
	207	+void special_key(
	208	+ int key /*< [in] typed special key/,
	209	+ int x /*< [in]/,
	210	+ int y /*< [in]/)
	211	+{
	212	+ switch (key) {
	213	+ case GLUT_KEY_UP:
	214	+ trans[1] = trans[1] + 0.1;
	215	+ glutPostRedisplay();
	216	+ break;
	217	+ case GLUT_KEY_DOWN:
	218	+ trans[1] = trans[1] - 0.1;
	219	+ glutPostRedisplay();
	220	+ break;
	221	+ case GLUT_KEY_RIGHT:
	222	+ /**/
	223	+ trans[0] = trans[0] + 0.1;
	224	+ glutPostRedisplay();
	225	+ break;
	226	+ /**/
	227	+ case GLUT_KEY_LEFT:
	228	+ /**/
	229	+ trans[0] = trans[0] - 0.1;
	230	+ glutPostRedisplay();
	231	+ break;
	232	+ /**/
	233	+ }
	234	+} /* special_key */

--- /dev/null

+++ b/src/operation.h

		@@ -0,0 +1,43 @@
	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 resize(
	26	+ int w /*<[in] Window width/,
	27	+ int h /*<[in] Window height/);
	28	+void mouse(
	29	+ int button /*< [in] pushed button/,
	30	+ int state /*< [in] down or up or ?/,
	31	+ int x /*< [in] position of mouse cursor/,
	32	+ int y /*< [in] position of mouse cursor/);
	33	+void motion(
	34	+ int x /*< [in] position of cursor/,
	35	+ int y /*< [in] position of cursor/);
	36	+void keyboard(
	37	+ unsigned char key /*< [in] Typed key/,
	38	+ int x /*< [in]/,
	39	+ int y /*< [in]/);
	40	+void special_key(
	41	+ int key /*< [in] typed special key/,
	42	+ int x /*< [in]/,
	43	+ int y /*< [in]/);

--- /dev/null

+++ b/src/read_file.c

		@@ -0,0 +1,158 @@
	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	+#include <stdlib.h>
	26	+#include <stdio.h>
	27	+#include "variable.h"
	28	+#include "basic_math.h"
	29	+
	30	+/**
	31	+* Input from Fermi surface file
	32	+*/
	33	+void read_file(char fname/<[in] fname Input file name/)
	34	+{
	35	+ int ib, i, i0, i1, i2, ii0, ii1, ii2, ierr;
	36	+ FILE *fp;
	37	+ /*
	38	+ Open input file.
	39	+ */
	40	+ if ((fp = fopen(fname, "r")) == NULL) {
	41	+ printf("file open error!!\n");
	42	+ printf(" Press any key to exit.\n");
	43	+ getchar();
	44	+ exit(EXIT_FAILURE);
	45	+ }
	46	+ printf("\n##### Brillouin zone informations ##### \n\n");
	47	+ /*
	48	+ k-point grid
	49	+ */
	50	+ ierr = fscanf(fp, "%d%d%d", &ng0[0], &ng0[1], &ng0[2]);
	51	+ if (ierr == 0) printf("error ! reading ng");
	52	+ printf("k point grid : %d %d %d \n", ng0[0], ng0[1], ng0[2]);
	53	+ for (i = 0; i < 3; i++) ng[i] = ng0[i];
	54	+ /*
	55	+ Shift of k-point grid
	56	+ */
	57	+ ierr = fscanf(fp, "%d", &lshift);
	58	+ if (ierr == 0) printf("error ! reading lshift");
	59	+
	60	+ if (lshift == 0) {
	61	+ printf("k point grid is the Monkhorst-Pack grid. \n");
	62	+ for (i = 0; i < 3; i++) shiftk[i] = (ng0[i] + 1) % 2;
	63	+ }
	64	+ else if (lshift == 1) {
	65	+ printf("k point grid starts from Gamma. \n");
	66	+ for (i = 0; i < 3; i++) shiftk[i] = 0;
	67	+ }
	68	+ else if (lshift == 2) {
	69	+ printf("k point grid starts from Gamma + a half grid. \n");
	70	+ for (i = 0; i < 3; i++) shiftk[i] = 1;
	71	+ }
	72	+ else {
	73	+ exit(0);
	74	+ }
	75	+ /*
	76	+ * # of bands
	77	+ */
	78	+ ierr = fscanf(fp, "%d", &nb);
	79	+ if (ierr == 0) printf("error ! reading nb");
	80	+ printf("# of bands : %d \n", nb);
	81	+ ntri = (int)malloc(nb sizeof(int));
	82	+ nnl = (int)malloc(nb sizeof(int));
	83	+ draw_band = (int)malloc(nb sizeof(int));
	84	+ for (ib = 0; ib < nb; ib++) draw_band[ib] = 1;
	85	+ /*
	86	+ Reciplocal lattice vectors
	87	+ */
	88	+ for (i = 0; i < 3; ++i) {
	89	+ ierr = fscanf(fp, "%e%e%e", &bvec[i][0], &bvec[i][1], &bvec[i][2]);
	90	+ if (ierr == 0) printf("error ! reading bvec");
	91	+ printf("bvec %d : %f %f %f \n", i + 1, bvec[i][0], bvec[i][1], bvec[i][2]);
	92	+ }
	93	+ /*
	94	+ Allocation of Kohn-Sham energies $ matrix elements
	95	+ */
	96	+ eig0 = (GLfloat***)malloc(nb sizeof(GLfloat***));
	97	+ mat0 = (GLfloat***)malloc(nb sizeof(GLfloat***));
	98	+ eig = (GLfloat***)malloc(nb sizeof(GLfloat***));
	99	+ mat = (GLfloat***)malloc(nb sizeof(GLfloat***));
	100	+ for (ib = 0; ib < nb; ib++) {
	101	+ eig0[ib] = (GLfloat**)malloc(ng0[0] sizeof(GLfloat**));
	102	+ mat0[ib] = (GLfloat**)malloc(ng0[0] sizeof(GLfloat**));
	103	+ eig[ib] = (GLfloat**)malloc(ng0[0] sizeof(GLfloat**));
	104	+ mat[ib] = (GLfloat**)malloc(ng0[0] sizeof(GLfloat**));
	105	+ for (i0 = 0; i0 < ng0[0]; i0++) {
	106	+ eig0[ib][i0] = (GLfloat*)malloc(ng0[1] sizeof(GLfloat*));
	107	+ mat0[ib][i0] = (GLfloat*)malloc(ng0[1] sizeof(GLfloat*));
	108	+ eig[ib][i0] = (GLfloat*)malloc(ng0[1] sizeof(GLfloat*));
	109	+ mat[ib][i0] = (GLfloat*)malloc(ng0[1] sizeof(GLfloat*));
	110	+ for (i1 = 0; i1 < ng0[1]; i1++) {
	111	+ eig0[ib][i0][i1] = (GLfloat)malloc(ng0[2] sizeof(GLfloat));
	112	+ mat0[ib][i0][i1] = (GLfloat)malloc(ng0[2] sizeof(GLfloat));
	113	+ eig[ib][i0][i1] = (GLfloat)malloc(ng0[2] sizeof(GLfloat));
	114	+ mat[ib][i0][i1] = (GLfloat)malloc(ng0[2] sizeof(GLfloat));
	115	+ }
	116	+ }
	117	+ }
	118	+ /*
	119	+ Kohn-Sham energies
	120	+ */
	121	+ for (ib = 0; ib < nb; ++ib) {
	122	+ for (i0 = 0; i0 < ng0[0]; ++i0) {
	123	+ if (lshift != 0) ii0 = i0;
	124	+ else ii0 = modulo(i0 + (ng0[0] + 1) / 2, ng0[0]);
	125	+ for (i1 = 0; i1 < ng0[1]; ++i1) {
	126	+ if (lshift != 0) ii1 = i1;
	127	+ else ii1 = modulo(i1 + (ng0[1] + 1) / 2, ng0[1]);
	128	+ for (i2 = 0; i2 < ng0[2]; ++i2) {
	129	+ if (lshift != 0) ii2 = i2;
	130	+ else ii2 = modulo(i2 + (ng0[2] + 1) / 2, ng0[2]);
	131	+ ierr = fscanf(fp, "%e", &eig0[ib][ii0][ii1][ii2]);
	132	+ eig[ib][ii0][ii1][ii2] = eig0[ib][ii0][ii1][ii2];
	133	+ }
	134	+ }
	135	+ }
	136	+ }
	137	+ /*
	138	+ Matrix elements
	139	+ */
	140	+ for (ib = 0; ib < nb; ++ib) {
	141	+ for (i0 = 0; i0 < ng0[0]; ++i0) {
	142	+ if (lshift != 0) ii0 = i0;
	143	+ else ii0 = modulo(i0 + (ng0[0] + 1) / 2, ng0[0]);
	144	+ for (i1 = 0; i1 < ng0[1]; ++i1) {
	145	+ if (lshift != 0) ii1 = i1;
	146	+ else ii1 = modulo(i1 + (ng0[1] + 1) / 2, ng0[1]);
	147	+ for (i2 = 0; i2 < ng0[2]; ++i2) {
	148	+ if (lshift != 0) ii2 = i2;
	149	+ else ii2 = modulo(i2 + (ng0[2] + 1) / 2, ng0[2]);
	150	+ ierr = fscanf(fp, "%e", &mat0[ib][ii0][ii1][ii2]);
	151	+ mat[ib][ii0][ii1][ii2] = mat0[ib][ii0][ii1][ii2];
	152	+ }
	153	+ }
	154	+ }
	155	+ }
	156	+ fclose(fp);
	157	+ /**/
	158	+} /* read_file */

--- /dev/null

+++ b/src/read_file.h

		@@ -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 read_file(char fname/<[in] fname Input file name/);

--- /dev/null

+++ b/src/variable.h

		@@ -0,0 +1,119 @@
	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	+#pragma once
	25	+
	26	+#if defined(MAC)
	27	+#include <GLUT/glut.h>
	28	+#else
	29	+#include <GL/glut.h>
	30	+#endif
	31	+
	32	+/**
	33	+* Input variables
	34	+*/
	35	+int ng0[3];
	36	+int lshift; /*< Switch for shifted Brillouin zone /
	37	+int shiftk[3];
	38	+int nb; /*< The number of Bands /
	39	+GLfloat bvec[3][3]; /*< Resiplocal lattice vector /
	40	+GLfloat **eig0; /< Eigenvalues [nb][ng0[0]][ng0[1]][ng0[2]] */
	41	+GLfloat **mat0; /< Matrix element [nb][ng0[0]][ng0[1]][ng0[2]] */
	42	+/*
	43	+ * Interpolation
	44	+ */
	45	+int ng[3]; /*< BZ grids /
	46	+GLfloat **eig; /< Eigenvalues [nb][ng[0]][ng[1]][ng[2]] */
	47	+GLfloat **mat; /< Matrix element [nb][ng[0]][ng[1]][ng[2]] */
	48	+int interpol;
	49	+/**
	50	+ * Switch for some modes
	51	+ */
	52	+int blackback; /*< Switch for black background /
	53	+int fcscl; /*< Switch for full color scale mode /
	54	+int fbz; /*< Switch for 1st Brillouin zone mode /
	55	+int nodeline; /*< Switch for node lines /
	56	+int lcolorbar; /*< Switch for colorbar /
	57	+int lstereo; /*< Switch for the stereogram /
	58	+int lmouse; /*< Switch for the mouse function /
	59	+/*
	60	+ * Menu
	61	+ */
	62	+int imenu_band, imenu_interpol, imenu_bgcolor, imenu_colorscale, imenu_bzmode,
	63	+imenu_nodeline, imenu_colorbar, imenu_tetra, imenu_stereo, imenu_mouse, imenu;
	64	+/**
	65	+* Variables for Brillouin zone boundaries
	66	+*/
	67	+int nbzl; /*< The number of Lines of 1st Brillouin zone /
	68	+GLfloat *bzl; /< Lines of 1st BZ [nbzl][2][3] */
	69	+GLfloat bragg[26][3]; /*< Bragg plane vectors /
	70	+GLfloat brnrm[26]; /*< Norms of Bragg plane vectors /
	71	+/**
	72	+ * Variables for patchs
	73	+ */
	74	+int ntri; /< The number of triangle patch [nb] /
	75	+int draw_band; /< Switch for drawn bands [nb] /
	76	+GLfloat *nmlp; /< Normal vector of patchs [nb][ntri][3] */
	77	+GLfloat **kvp; /< K-vectors of points [nb][ntri][3][3] */
	78	+GLfloat *matp; /< Matrix elements of points [nb][ntri][3] */
	79	+GLfloat **clr; /< Colors of points [nb][ntri][3][4] */
	80	+int itet; /*< Counter for tetrahedron /
	81	+/**
	82	+ * Variables for nodeline
	83	+ */
	84	+int nnl; /< The number of nodeline /
	85	+GLfloat **kvnl; /< K-vector of nodeline [nb][nnl][2][3] */
	86	+/**
	87	+ * Variables for mouse & cursorkey
	88	+ */
	89	+GLfloat sx; /*< Scale of mouse movement /
	90	+GLfloat sy; /*< Scale of mouse movement /
	91	+int cx; /*< Starting point of drug /
	92	+int cy; /*< Starting point of drug /
	93	+GLfloat scl; /*< Initial scale /
	94	+GLfloat trans[3]; /*< Translation /
	95	+GLfloat rot[3][3]; /*< Rotation matrix /
	96	+GLfloat thetax, thetay, thetaz;
	97	+/**
	98	+* Colors
	99	+*/
	100	+GLfloat black[4]; /*< Black color code /
	101	+GLfloat white[4]; /*< White color code /
	102	+GLfloat cyan[4]; /*< Cyan color code /
	103	+GLfloat magenta[4]; /*< Magenta color code /
	104	+GLfloat yellow[4]; /*< Yellow color code /
	105	+GLfloat red[4]; /*< Red color code /
	106	+GLfloat green[4]; /*< Green color code /
	107	+GLfloat blue[4]; /*< Blue color code /
	108	+/*
	109	+ * Others
	110	+ */
	111	+int query; /*< Query switch /
	112	+int corner[6][4]; /*< Corners of tetrahedron /
	113	+GLfloat EF; /*< Fermi energy /
	114	+enum
	115	+{
	116	+ MOUSE_SCROLL_UP = 3, /*< Mouse wheel up /
	117	+ MOUSE_SCROLL_DOWN = 4 /*< Mouse wheel down /
	118	+};
	119	+int nthreads;/*< Number of OpenMP threads /

Show on old repository browser

fermisurfer Forkfermisurferandroid

fermisurfer: Commit

Commit MetaInfo

Log Message

Ändern Zusammenfassung

Diff

fermisurfer
Fork
fermisurfer
android