***************************************************************
*
*     Programm zum zeichen vom Bandstruktur
*
***************************************************************
                                           
      PROGRAM PLOTBAND

      CHARACTER NAME*40,HEAD*40,DTXT1*40,DTXT2*40,LABEL*40
      PARAMETER (NDIR=8,NBANDDIM=499,NNDIM=34,NPRODIM=11)
      DIMENSION STEP(NDIR+1),FE(NDIR*NNDIM,NBANDDIM*NPRODIM),XP(NNDIM)
      DIMENSION FEM(NDIR*NNDIM,2*NBANDDIM)
      CHARACTER*10 TEXT(NDIR+1)
      DIMENSION QBEG(3),QEND(3)
      DIMENSION DSH1(3),DSH2(3),DSH3(3)

 
      WRITE(*,*) 'PROGRAM PLOTBAND: zeichnet Bandstrucktur' 
      WRITE(*,*) 'Plotlabel:'
      READ(*,'(A)') LABEL

C-----Beschriftung einlesen
      OPEN(1,FILE='BANDCTR',STATUS='OLD',ERR=100)
      READ(1,'(A)') NAME
      READ(1,*) ND
      READ(1,*) NN
      DO 20 I=1,ND
        READ(1,*) TEXT(I),STEP(I)
        READ(1,*) QBEG(1),QBEG(2),QBEG(3)
        READ(1,*) QEND(1),QEND(2),QEND(3)
   20 CONTINUE
      STEP(ND+1)=0.
      READ(1,*) TEXT(ND+1)
      EFERMI=0.
      NT=0
      NPROJ=0
      READ(1,*,ERR=25,END=25) EFERMI
      READ(1,*,ERR=25,END=25) NT
      READ(1,*,ERR=25,END=25) ID1,ID2,ID3,RDUM1,RDUM2,NPROJ
   25 CLOSE(1)
      
C-----Art des Plots
      IF (NPROJ.EQ.0) THEN
         METHOD=0
         NPROJ=1
      ELSE
         SIGMA=0
         WRITE(*,*) 'smearing in eV (0 for true curves):'
         READ(*,*) SIGMA
         WRITE(*,*) 'Plot min/max with:'
         WRITE(*,*) '  (1) no interpolation'
         WRITE(*,*) '  (2) quadratic interpolation with'//
     &                   ' periodic boundary'
         WRITE(*,*) '  (3) quadratic interpolation without'//
     &                   ' periodic boundary'
         READ(5,*) METHOD
         IF ((METHOD.LT.1).OR.(METHOD.GT.3)) THEN
            WRITE(*,*) 'Invalid choice !'
            STOP
         ENDIF
         IF ((METHOD.EQ.2).OR.(METHOD.EQ.3)) THEN
            WRITE(*,*) 'EDIFF (for warning, e.g. 0.1):'
            READ(*,*) EDIFF
            WRITE(*,*)
            WRITE(*,*) 'WARNINGS (rel. to E(fermi)):'
            WRITE(*,*) 'what  Point1   Point2   Point3'//
     &                 '   E(old)   E(new)   DE      kPt Band kz'
         ENDIF
      ENDIF

C-----Bandstruktur einlesen 
      NPOINTS=ND*NN
      OPEN(1,FILE='EIGENVAL',STATUS='OLD',ERR=100)  

      READ(1,*) NDUM1,NDUM2,NDUM3,NDUM4
      READ(1,*) RDUM1,RDUM2,RDUM3,RDUM4,RDUM5
      READ(1,*) RDUM6
      READ(1,'(A)') DTXT1
      READ(1,'(A)') DTXT2
      READ(1,*) NEL,NKPT,NBAND
      IF (NT.LE.0) NT=NEL/2+4
      NT=MIN(NT,NBAND)
      IF (NKPT.NE.NPOINTS*NPROJ) THEN
         WRITE(*,*) 'Error on input: ND*NN*NPROJ from file BAND does ',
     &              'not equal NKPT from EIGENVAL'
         WRITE(*,*) ND*NN*NPROJ, NKPT
         STOP
      ENDIF
      EFTEST=-1E30
      DO 30 I=1,NPOINTS
      DO 27 NB=1,NBAND
         FEM(I,NB*2-1)=100.
         FEM(I,NB*2)=-100.
27    CONTINUE  
      DO 30 NP=0,NPROJ-1
C fuer EIGENVAL aus OUTCAR
C         read(1,*)
C ende fuer EIGENVAL aus OUTCAR
         READ(1,*)
         READ(1,*)
         IF (NT.GE.NBAND) NT=NBAND
         
         DO 35 NB=1,NBAND
            READ(1,*) NCOUNT,EIGEN
            IF (NB.LE.NT) THEN
               FE(I,NT*NP+NB)=EIGEN
               FEM(I,2*NB-1)=MIN(FEM(I,2*NB-1),EIGEN)
               FEM(I,2*NB)  =MAX(FEM(I,2*NB),EIGEN)
            ENDIF
	    IF (NB.EQ.(NEL/2)) EFTEST=MAX(EFTEST,EIGEN)
   35    CONTINUE
   30 CONTINUE 
      IF (EFERMI.LT.-1E10) EFERMI=EFTEST
      PRINT*,EFERMI

C calculate quadratic interpolation
      IF(METHOD.GE.2) THEN
       DO I=1,NPOINTS
        DO NB=1,NT
         EMAX=FEM(I,2*NB)
         EMIN=FEM(I,2*NB-1)
         DO NP=0,NPROJ-1

C get the 3 Points E1, E2 and E3
          IF (NP.EQ.0) THEN
           IF(METHOD.EQ.2) THEN
            E1=FE(I,NT*(NPROJ-1)+NB)
           ELSE
            E1=FE(I,NT*(NP+1)+NB)
           ENDIF
          ELSE
           E1=FE(I,NT*(NP-1)+NB)
          ENDIF
          E2=FE(I,NT*NP+NB)
          IF (NP.EQ.NPROJ-1) THEN
           IF(METHOD.EQ.2) THEN
            E3=FE(I,NB)
           ELSE
            E3=E1
           ENDIF
          ELSE
           E3=FE(I,NT*(NP+1)+NB)
          ENDIF

          IF (E1.LE.E2.AND.E3.LE.E2) THEN
           ETEST=E2-(E3-E1)**2/8/(E3+E1-2*E2)
           IF ((ETEST-EMAX).GT.EDIFF) THEN
            WRITE (*,'(A,6F9.5,3I4)') ' MAX',E1-EFERMI,E2-EFERMI,
     &            E3-EFERMI,EMAX-EFERMI,ETEST-EFERMI,ETEST-EMAX,I,NB,NP
           ENDIF
           FEM(I,2*NB)=MAX(FEM(I,2*NB),ETEST)
          ELSE IF (E1.GE.E2.AND.E3.GE.E2) THEN
           ETEST=E2-(E3-E1)**2/8/(E3+E1-2*E2)
           IF ((EMIN-ETEST).GT.EDIFF) THEN
            WRITE (*,'(A,6F9.5,3I4)') ' MIN',E1-EFERMI,E2-EFERMI,
     &            E3-EFERMI,EMIN-EFERMI,ETEST-EFERMI,ETEST-EMIN,I,NB,NP
           ENDIF
           FEM(I,2*NB-1)=MIN(FEM(I,2*NB-1),ETEST)
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDIF

C write BAND or BBAND (for mathematica)
      IF (NPROJ.EQ.1) THEN
        OPEN(3,FILE='BAND',STATUS='UNKNOWN')
        WRITE(3,'(9(I5,F10.4))') ((I,FE(I,N),I=1,NPOINTS),N=1,NT)
        CLOSE(3)
      ELSE
        OPEN(3,FILE='BBAND',STATUS='UNKNOWN')
        WRITE(3,'(9(I5,F10.4))') ((I,FEM(I,N),I=1,NPOINTS),N=1,2*NT)
        CLOSE(3)
      ENDIF

C write BAND.xmgr and BBAND.xmgr (for xmgr)
C if this is a projected bandstructure, then BAND.xmgr is just
C for the style - you can use it if you don't want to draw the
C bandstructure from your surface
      OPEN(3,FILE='BAND.xmgr',STATUS='UNKNOWN')
      WRITE(3,'(3A)') '@ title "',LABEL(1:ITEXT(LABEL,40)),'"'
      WRITE(3,'(A)') '@ type xy'
      IF (NPROJ.EQ.1) THEN
       DO N=1,NT
        X=0
        DO I=1,ND
         DO M=1,NN
           XP(M)=X+(M-1)*STEP(I)/(NN-1)
           WRITE(3,'(2F10.4)') XP(M),FE(M+(I-1)*NN,N)-EFERMI
         ENDDO
         X=X+STEP(I)
         ENDDO
         WRITE(3,'(2F10.4)') X,0
         WRITE(3,'(2F10.4)') 1,0
       ENDDO
      ELSE
C to keep set 0 alive
       WRITE(3,'(2F10.4)') -1,0
      ENDIF
      WRITE(3,'(A)') '&'
      X=0
      DO I=1,ND-1
        X=X+STEP(I)
        WRITE(3,'(A)') '@ type xy'
        WRITE(3,'(2F12.4)') X,-100.
        WRITE(3,'(2F12.4)') X,100.
        WRITE(3,'(A)') '&' 
      ENDDO
C-- search max and min values
      OPEN(15,FILE='BANDSCALE',STATUS='OLD',ERR=789)
      READ(15,*,ERR=789) FMIN,FMAX
      CLOSE(15)
      GOTO 788
  789 CONTINUE
      FMAX=-1000
      FMIN= 1000                                
      DO 50 I=1,NPOINTS
       DO 50 NB=1,NT
        IF (FEM(I,2*NB).GT.FMAX) FMAX=FEM(I,2*NB)
        IF (FEM(I,2*NB-1).LT.FMIN) FMIN=FEM(I,2*NB-1)
   50 CONTINUE  
      FMIN=FMIN-.5-EFERMI
      FMAX=FMAX+.5-EFERMI
  788 CONTINUE
      WRITE(3,'(A)') '@ with g0'
      WRITE(3,'(A,F8.3)') '@ world xmin',0
      WRITE(3,'(A,F8.3)') '@ world xmax',X+STEP(ND)
      WRITE(3,'(A,F8.3)') '@ world ymin',FMIN
      WRITE(3,'(A,F8.3)') '@ world ymax',FMAX
      WRITE(3,'(A)') '@ yaxis  ticklabel start type spec'
      WRITE(3,'(A,I4)') '@ yaxis  ticklabel start',5*INT(FMIN/5)
      DO I=0,5+ND
        IF (I.LE.9) WRITE(3,'(A3,I1,A)') '@ s',I,' color 1'
        IF (I.GT.9) WRITE(3,'(A3,I2,A)') '@ s',I,' color 1'
      ENDDO
      WRITE(3,'(A)') '@ xaxis  ticklabel type spec'
      WRITE(3,'(A)') '@ xaxis  tick type spec'
      WRITE(3,'(A,I2)') '@ xaxis  tick spec ',ND+1
      X=0
      DO I=1,ND+1
        WRITE(3,'(A,I2,A2,F8.3)') '@ xaxis  tick ',I-1,', ',X
        WRITE(3,'(A,I2,3A)') '@ xaxis  ticklabel ',I-1,', "',
     &                      TEXT(I)(1:ITEXT(TEXT(I),10)),'"'
        X=X+STEP(I)
      ENDDO
      WRITE(3,'(A)') '@ yaxis  tick major 5'
      WRITE(3,'(A)') '@ yaxis  tick minor 1'
      WRITE(3,'(A)') '@ yaxis  label "Energy (eV)"'
C-- styles for SBAND.xmgr
      IF (NPROJ.EQ.1) THEN
       II=1
       III=2
       DO I=ND,ND+5
         IF (I.LE.9) THEN 
           WRITE(3,'(A3,I1,A)') '@ s',I,' linestyle 0'
           WRITE(3,'(A3,I1,A)') '@ s',I,' symbol size 0.7'
           WRITE(3,'(A3,I1,A,I1)') '@ s',I,' symbol fill ',II
           WRITE(3,'(A3,I1,A,I1)') '@ s',I,' symbol ',III
         ELSE
           WRITE(3,'(A3,I2,A)') '@ s',I,' linestyle 0'
           WRITE(3,'(A3,I2,A)') '@ s',I,' symbol size 0.7'
           WRITE(3,'(A3,I2,A,I1)') '@ s',I,' symbol fill ',II
           WRITE(3,'(A3,I2,A,I1)') '@ s',I,' symbol ',III
         ENDIF
         IF (II.EQ.1) THEN
           II=2
         ELSE
           II=1
         ENDIF
         IF (II.EQ.1) THEN
           IF (III.EQ.3) III=5
           IF (III.EQ.2) III=3
         ENDIF
       ENDDO
      ENDIF
C-- increase the linewidth - trust me, (usually) it looks better
      WRITE(3,'(A)') '@ ps linewidth begin 2'
      CLOSE(3)

C-- now for the BBAND.xmgr - if we have a projection
      IF (NPROJ.GT.1) THEN
       OPEN(3,FILE='BBAND.xmgr',STATUS='UNKNOWN')
       DO NB=1,NT
         WRITE(3,'(A)') '@ type xy'
         X=0
         DO I=1,ND
          DO M=1,NN
           XP(M)=X+(M-1)*STEP(I)/(NN-1)
           WRITE(3,'(2F10.4)') XP(M),FEM(M+(I-1)*NN,2*NB-1)-EFERMI-SIGMA
          ENDDO
          X=X+STEP(I)
         ENDDO
         DO I=ND,1,-1
          X=X-STEP(I)
          DO M=NN,1,-1
           XP(M)=X+(M-1)*STEP(I)/(NN-1)
           WRITE(3,'(2F10.4)') XP(M),FEM(M+(I-1)*NN,2*NB)-EFERMI+SIGMA
          ENDDO
         ENDDO
         WRITE(3,'(2F10.4)') 0,FEM(1,2*NB-1)-EFERMI-SIGMA
         WRITE(3,'(A)') '&'
       ENDDO
C-- grid lines once again, because they should be in front
       X=0
       DO I=1,ND-1
         X=X+STEP(I)
         WRITE(3,'(A)') '@ type xy'
         WRITE(3,'(2F12.4)') X,-100.
         WRITE(3,'(2F12.4)') X,100.
         WRITE(3,'(A)') '&'
       ENDDO
       WRITE(3,'(A)') '@ type xy'
       WRITE(3,'(2F12.4)') 0.,0.
       WRITE(3,'(2F12.4)') X+STEP(ND),0.
       WRITE(3,'(A)') '&'
       DO I=ND,5+2*ND+NT
         IF (I.LE.9) THEN
           WRITE(3,'(A3,I1,A)') '@ s',I,' color 1'
         ELSE IF (I.LE.99) THEN
           WRITE(3,'(A3,I2,A)') '@ s',I,' color 1'
         ELSE 
           WRITE(3,'(A3,I3,A)') '@ s',I,' color 1'
         ENDIF
       ENDDO
       DO I=ND,ND-1+NT
         IF (I.LE.9) THEN
           WRITE(3,'(A3,I1,A)') '@ s',I,' fill 1'
           WRITE(3,'(A3,I1,A)') '@ s',I,' fill with color'
           WRITE(3,'(A3,I1,A)') '@ s',I,' fill color 7'
           WRITE(3,'(A3,I1,A)') '@ s',I,' symbol 0'
           WRITE(3,'(A3,I1,A)') '@ s',I,' linestyle 0'
         ELSE IF (I.LE.99) THEN
           WRITE(3,'(A3,I2,A)') '@ s',I,' fill 1'
           WRITE(3,'(A3,I2,A)') '@ s',I,' fill with color'
           WRITE(3,'(A3,I2,A)') '@ s',I,' fill color 7'
           WRITE(3,'(A3,I2,A)') '@ s',I,' symbol 0'
           WRITE(3,'(A3,I2,A)') '@ s',I,' linestyle 0'
         ELSE
           WRITE(3,'(A3,I3,A)') '@ s',I,' fill 1'
           WRITE(3,'(A3,I3,A)') '@ s',I,' fill with color'
           WRITE(3,'(A3,I3,A)') '@ s',I,' fill color 7'
           WRITE(3,'(A3,I3,A)') '@ s',I,' symbol 0'
           WRITE(3,'(A3,I3,A)') '@ s',I,' linestyle 0'
         ENDIF
       ENDDO

C-- styles for SBAND.xmgr
       II=1
       III=2
       DO I=NT+2*ND,2*ND+NT+5
         IF (I.LE.9) THEN
           WRITE(3,'(A3,I1,A)') '@ s',I,' linestyle 0'
           WRITE(3,'(A3,I1,A)') '@ s',I,' symbol size 0.7'
           WRITE(3,'(A3,I1,A,I1)') '@ s',I,' symbol fill ',II
           WRITE(3,'(A3,I1,A,I1)') '@ s',I,' symbol ',III
         ELSE IF (I.LE.99) THEN
           WRITE(3,'(A3,I2,A)') '@ s',I,' linestyle 0'
           WRITE(3,'(A3,I2,A)') '@ s',I,' symbol size 0.7'
           WRITE(3,'(A3,I2,A,I1)') '@ s',I,' symbol fill ',II
           WRITE(3,'(A3,I2,A,I1)') '@ s',I,' symbol ',III
         ELSE
           WRITE(3,'(A3,I3,A)') '@ s',I,' linestyle 0'
           WRITE(3,'(A3,I3,A)') '@ s',I,' symbol size 0.7'
           WRITE(3,'(A3,I3,A,I1)') '@ s',I,' symbol fill ',II
           WRITE(3,'(A3,I3,A,I1)') '@ s',I,' symbol ',III
         ENDIF
         IF (II.EQ.1) THEN
           II=2
         ELSE
           II=1
         ENDIF
         IF (II.EQ.1) THEN
           IF (III.EQ.3) III=5
           IF (III.EQ.2) III=3
         ENDIF
       ENDDO
       CLOSE(3)
      ENDIF

      GOTO 1000

  100 WRITE(*,*)'ERROR: InitFile BANDCTR oder EIGENVAL existiert nicht'
      STOP     
 1000 CONTINUE
      END     

      INTEGER FUNCTION  ITEXT(TEXT,N)
      CHARACTER*(*) TEXT
      DO 10 I=N,1,-1
        IF (TEXT(I:I).NE.' ') GOTO 20
   10 CONTINUE
   20 CONTINUE
      ITEXT=I
      RETURN
      END

