1! Copyright Eyvaz Isaev, PWSCF group, 2009-2010
2!
3! Department of Physics, Chemistry and Biology (IFM)
4! Linkoping Univrsity, Sweden
5!
6! Theoretical Physics Department
7! Moscow State Institute of Steel and Alloys, Russia
8!
9! Program prepares a part of gnuplot file
10! 1-defines axis labels
11! 2-defines axis
12!
13!	implicit real*8(a-h,o-z)
14	implicit none
15
16!        real*8  kpoint(3,1000)
17        integer maxlin, idiv, integs(100)
18        integer i, j, k, nbnd, kpnt
19        real*8 ax(3),by(3),cz(3)
20	real*8 dl,kx,ky,kz
21	real*8 E_min, E_max, THz, e_min_low
22        character*10 dummy
23        character*1 dummy1
24        character*11 dummy3
25        character*3 dummy2
26	character*65 plotchar
27
28	real*8, allocatable :: e(:)
29	real*8, allocatable :: kpoint(:,:)
30!
31	read(5,'(12x,i4,6X,i5,1X)') nbnd, kpnt
32	allocate (e(nbnd))
33	allocate (kpoint(3,kpnt))
34
35	do i=1,kpnt
36	read(5,*)
37	read(5,*) (e(j),j=1,nbnd)
38
39        if(i.eq.1) then
40        E_min=e(1)
41        E_max=e(1)
42        endif
43
44	do k=1,nbnd
45        if(e(k).le.E_min) E_min=e(k)
46        if(e(k).ge.E_max) E_max=e(k)
47        enddo
48
49	enddo
50
51	e_min=dint(e_min*1.2)
52	e_max=dint(e_max*1.2)
53        e_min_low=e_min-33.0/2
54
55	open(12,file='Freq_plot_unit')
56	open(15,file='plot.GNU.tmp')
57
58	read(12,*) THz
59
60	write(15,*)
61	write(15,'("THz=", f8.4)') THz
62	write(15,*)
63	write(15,'("E_min=",f8.1)') E_min/Thz
64	write(15,'("E_max=",f8.1)') E_max/Thz
65	write(15,'("Label_position=",f8.1)') E_min_low/Thz
66	write(15,*)
67
68        close(12)
69
70! Writing labels
71!
72	idiv=1
73	maxlin=0
74
75	open(9,file='kpts')
76	open(10,file='K_points')
77
78        read(10,*) dummy
79        print*, dummy
80
81        read(9,*)
82	read(9,*)
83
84	read(10,*) ax(1),by(1),cz(1)
85	read(10,*) ax(2),by(2),cz(2)
86	read(10,*) ax(3),by(3),cz(3)
87
88	read(10,*)
8911	read(10,*,end=99) integs(idiv), (kpoint(j,idiv),j=1,3), dummy1
90	print*, dummy1
91	if(dummy1.eq.'G') dummy3="{/Symbol G}"
92	dummy2="//dummy1//"
93!	write(6,'(i6,2x,3f10.6)')  integs(idiv), (kpoint(j,idiv),j=1,3)
94
95	read(9,*) dl, kx,ky,kz
96
97	dl=dl-0.03
98	if(dummy1.eq.'G') then
99	write(15,'("set label  ", 1h", a, 1h", "  at  ", f6.2, " , ", "Label_position")') dummy3, dl
100	else
101	write(15,'("set label  ", 1h",  a, 1h", "  at  ", f6.2, " , ", "Label_position")') dummy1, dl
102	endif
103
104	idiv=idiv+1
105	maxlin=maxlin+1
106
107	goto 11
10899	continue
109
110	write(15,*)
111        close(9)
112
113        open(9,file='kpts')
114	open(15,file='plot.GNU.tmp')
115
116	read(9,'(1x,I5)') maxlin
117	read(9,*)
118
119	do i=1,maxlin
120	read(9,*) dl,kx,ky,kz
121!
122        if(i.ne.1.and.i.ne.maxlin) then
123	 write(15,'("set arrow nohead from  ",f8.5, ", E_min", "   to  ", f8.5,", E_max", "  lw 3")') dl,  dl
124        endif
125!
126	enddo
127
128	close(9)
129	close(10)
130	close(14)
131
132	write(15,*)
133
134        plotchar="plot [:] [E_min:E_max] 'Frequency' u 1:($2/THz) w lines lt 1 lw 3"
135        write(15,'(a)') plotchar
136	write(15,*)
137
138!	deallocate (e(nbnd))
139!	deallocate (kpoint(3,kpnt))
140
141	stop
142	end
143
144