      IMPLICIT REAL*8 (A-H,O-Z)
      EXTERNAL FUNC
      COMMON/SPG2/IL,NG,IG(48),JV(2,3,48)
      COMMON/BZPL/INDD,IEOH
      DIMENSION NRECPT(3,30),RECTAX(4,30)
C
      REAL*8 XFU(3,2),WF(3),AF,BT,GM,ED,EL,EF,X1,X2,Y1,Y2,Z1,Z2
      REAL*8 EM,PN,P
      REAL*8 XX(3),XC(3),DR(3)
      COMMON/ENRSPL/EM(1900000),PN(3),P(3),N(3),JS
      INTEGER NN(3),NK(3)
      DIMENSION CO(3,2,100)
      DIMENSION JB(2,3)
      READ(1,*)
      READ(1,*) IL,NGEN,INV
      CALL TSPACE(IL)
C     CALL TSPHDS
C     CALL TSOPDS
      DO 1 I=1,NGEN
      READ(1,*) JA,((JB(J,K),J=1,2),K=1,3)
      CALL TSGENR(JA,JB)
    1 CONTINUE
      CALL TSPGRP(INV)  
      CALL TSPGDS       
      READ(1,*) A,B,C
      READ(1,*) CA,CB,CC
      CALL TSLATC(A,B,C,CA,CB,CC0)
      CALL TSBZEG(NRECPT,RECTAX,NRP,CO,NLIN)
      WRITE(6,901)
  901 FORMAT(/' RECIPROCAL LATTICE VECTOR FOR THE B.Z. FORMATION'
     &       /' UNIT IS A*, X-AXIS//A* Y-AXIS IS IN A* B* PLANE')
      WRITE(6,902)
  902 FORMAT('  NO ','   A* B* C*',
     &     '   RECTANGULAR COORDINATE',7X,'LENGTH')
      DO 41 J=1,NRP
      WRITE(6,900) J,(NRECPT(I,J),I=1,3),(RECTAX(I,J),I=1,4)
  900 FORMAT(I4,2X,3I3,4F10.5)
   41 CONTINUE
      WRITE(6,*) ' BRILLOUIN ZONE EDGES' 
      DO 31 J=1,NLIN
      WRITE(6,903) J,((CO(K,I,J),K=1,3),I=1,2)
  903 FORMAT(1H ,I4,' FROM(',3F9.5,') TO(',3F9.5,')')
   31 CONTINUE
      IF(IL.EQ.0) B=SQRT(3.0D0)*A
      READ(5,*) NB,ISPIN
      READ(21) NB1,NB2,NSPIN
	  WRITE(6,*) NB1,NB2,NSPIN
      READ(21) PN,N,JS
	  WRITE(6,*) PN
	  WRITE(6,*) N,JS
      DO 10 JSPIN=1,NSPIN
      DO 10 J=NB1,NB2
      READ(21) IXMAX,NBF,ISP,NNEE
	  IF(NNEE.GT.1900000) THEN
	      WRITE(6,*) ' NNEE=',NNEE
		  STOP ' NNEE IS TOO LARAGE'
      END IF 
      READ(21) (EM(K),K=1,NNEE)
      WRITE(6,*) IXMAX,ISP,NBF,NNEE
	  IF(NBF.EQ.NB.AND.ISP.EQ.ISPIN) GO TO 111
   10 CONTINUE
  111 CONTINUE
       IF(EM(1).EQ.0.0) STOP
      DO 51 I=1,3
      P(I)=PN(I)/DFLOAT(N(I))
   51 CONTINUE
      DO 52 I=0,8
	    XX(1)=I/8.0D0
		XX(2)=0.0
		XX(3)=0.0
		WAEK=FUNC(XX)
		WRITE(6,6000) XX,WAEK
 6000   FORMAT(4F10.5)
   52 CONTINUE
      IF(WAEK.LT.0.001D0) STOP
      XFU(1,1)=-PN(1)
      XFU(2,1)=-PN(2)
      XFU(3,1)=-PN(3)
      XFU(1,2)=PN(1)
      XFU(2,2)=PN(2)
      XFU(3,2)=PN(3)
C
C    CO AND RKG ARE GINEN IN UNIT OF A*
C
C      SIZE=80.0
      READ(5,*) SIZE
	  WASIZE=XFU(1,2)-XFU(1,1)
      WF(1)=SIZE
      WF(2)=SIZE*(A/B)*(XFU(2,2)-XFU(2,1))/WASIZE
      WF(3)=SIZE*(A/C)*(XFU(3,2)-XFU(3,1))/WASIZE
C
      READ(5,*) IEOH
      READ(5,*) EF
	  READ(5,*) NN,NK
	  WRITE(6,*) NN,NK
	  READ(5,*) AF
	  WRITE(6,*) ' ALPHA=',AF
C      AF=40.0
      BT=0.0
      GM=-20.0
      ED=150.0
      EL=1500.0
      CALL AYPSTR(97)
      CALL AYORIG(40.0,40.0)
	  WRITE(6,*) ' GM=',GM
      CALL TPERSP(XFU,WF,AF,BT,GM,ED,EL)
      CALL TPTRAC(1)   
      CALL TPSCHG(IEOH)
      CALL TPHILP(EF)
      CALL TPCIHL(3)
      INDD=0
      CALL TPHILD(0.0,1,6,5,0)
      CALL LINEWD(1.6)
      DO 42 I=1,NLIN
        IF(IL.GE.1) THEN
          X1=CO(1,1,I)
          Y1=CO(2,1,I)*(B/A)
          Z1=CO(3,1,I)*(C/A)
          X2=CO(1,2,I)
          Y2=CO(2,2,I)*(B/A)
          Z2=CO(3,2,I)*(C/A)
        ELSE IF(IL.EQ.0) THEN
          XK1=CO(1,1,I)
		  XK2=CO(1,2,I)
		  YK1=CO(2,1,I)
		  YK2=CO(2,2,I)
		  ZK1=CO(3,1,I)
		  ZK2=CO(3,2,I)
		  CALL GETKVC(XK1,YK1,ZK1,AK1,BK1,CK1)
		  CALL GETKVC(XK2,YK2,ZK2,AK2,BK2,CK2)
		  X1=AK1+BK1
		  X2=AK2+BK2
		  Y1=BK1-AK1
		  Y2=BK2-AK2
		  WRITE(6,610) AK1,BK1,CK1,AK2,BK2,CK2
		  Z1=CK1
		  Z2=CK2
		END IF
	  WRITE(6,610) X1,Y1,Z1,X2,Y2,Z2
  610 FORMAT(3F5.2,' TO',3F5.2)
      CALL TPLINE(X1,Y1,Z1,X2,Y2,Z2,30)
   42 CONTINUE
      CALL TPHILD(0.0,0,0,0,0)
      CALL LINEWD(1.5)
      CALL TPCIHL(0)
      CALL TPLINE(-1.0D0,0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,30)
      CALL TPLINE(0.0D0,-1.0D0,0.0D0,0.0D0,1.0D0,0.0D0,30)
      CALL TPLINE(0.0D0,0.0D0,-1.0D0,0.0D0,0.0D0,1.0D0,30)
      CALL LINEWD(0.1)
      CALL TPFXYZ(FUNC,EF,NN,NK,0)
      CALL LINEWD(0.11)
      DO 43 IRP=1,NRP
      IF(IL.GE.1) THEN
        DO 44 I=1,3
          DR(I)=NRECPT(I,IRP)
          XC(I)=0.5D0*NRECPT(I,IRP)
   44   CONTINUE
      ELSE IF(IL.EQ.0) THEN
         DR(1)=(NRECPT(1,IRP)+NRECPT(2,IRP))
         DR(2)=(NRECPT(2,IRP)-NRECPT(1,IRP))
         DR(3)=NRECPT(3,IRP)
         XC(1)=0.5D0*DR(1)
         XC(2)=0.5D0*DR(2)
         XC(3)=0.5D0*DR(3)
      END IF
         CALL TPSECT(FUNC,EF,DR,XC,0,1,241,241,2,2)
   43 CONTINUE
      CALL LINEWD(0.52)
      DO 46 I=1,3
        IF(I.EQ.1) THEN
           DR(1)=1.0D0
           DR(2)=0.0D0
           DR(3)=0.0D0
        ELSE IF(I.EQ.2) THEN
            DR(2)=1.0D0
            DR(1)=0.0D0
            DR(3)=0.0D0
        ELSE IF(I.EQ.3) THEN
            DR(1)=0.0D0
		    DR(2)=0.0D0
			DR(3)=1.0D0
		END IF
	     D0 45 I=1.3
		 XC(I)=0.0D0
   45    CONTINUE
        CALL TPSECT(FUNC,EF,DR,XC,0,1,241,241,2,2)
   46 CONTINUE
      CALL AYPEND
      STOP
      END
C FUNCTION   BZIORO ====*====3====*====4====*====5====*====6====*====7
C
C   HIDDEN LINE FUNCTION FOR B.Z.PLOT
C
C                 1988.10.18 :  A. YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
C      FUNCTION BZIORO(XA)
      FUNCTION FUNCD(XA)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/SPG2/IL,NG,IG(48),JV(2,3,48)
      COMMON/BZPL/INDD,IEOH
      DIMENSION XA(3)
      FUNCD=1.0*IEOH
      IF(INDD.EQ.0) THEN
        IF(IL.GE.1) THEN
         KX=1000000.0*XA(1)
         KY=1000000.0*XA(2)
         KZ=1000000.0*XA(3)
        ELSE IF(IL.EQ.0) THEN
         KX=1000000.0D0*0.5D0*(XA(1)-XA(2))
         KY=1000000.0D0*0.5D0*(XA(1)+XA(2))
         KZ=1000000.0D0*XA(3)
        END IF
         CALL TSKFBZ(KX,KY,KZ,1000000,IND)
         IF(IND.NE.0) FUNCD=-1.0*IEOH
      ELSE IF(INDD.EQ.1) THEN
         WW=FUNC(XA)
         FUNCD=WW
      END IF
      RETURN
      END
      FUNCTION LFUNC(XA)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/SPG2/IL,NG,IG(48),JV(2,3,48)
      DIMENSION XA(3)
      LFUNC=-1
        IF(IL.GE.1) THEN
         KX=1000000.0*XA(1)
         KY=1000000.0*XA(2)
         KZ=1000000.0*XA(3)
        ELSE IF(IL.EQ.0) THEN
         KX=1000000.0*0.5D0*(XA(1)-XA(2))
         KY=1000000.0*0.5D0*(XA(1)+XA(2))
         KZ=1000000.0*XA(3)
        END IF
         CALL TSKFBZ(KX,KY,KZ,1000000,IND)
         IF(IND.EQ.0) LFUNC=1
C		WRITE(6,*) ' IL=',IL,XA,KX,KY,KZ,IND	
      IF(XA(1).LT.0.0D0.AND.XA(2).LT.0.0D0.AND.XA(3).GT.0.0D0)
     &          LFUNC=1
      RETURN
      END
C      FUNCTION TLINTP(XX,JS)                                        
C     *** THREE DIMENSIONAL INTERPOLATION BY A LINEAR FORMULA                  
C     *** PERIODICITY AND F(X,Y,Z)=F(ABS(X),ABS(Y),ABS(Z))                     
C     *** ARE ASSUMED.                                                         
C     *** IS=1  ORTHORHOMBIC SYMMETRY                                          
C     *** IS=2  TETRAGONAL   SYMMETRY                                          
C     *** IS=3  CUBIC        SYMMETRY                                          
      FUNCTION FUNC(XX)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 KB 
      COMMON/ENRSPL/EM(1900000),PN(3),P(3),N(3),IS
      DIMENSION IFAC(3)
      DIMENSION A(4),IPER(3,6),FF(2,2,2),F(8),KB(3)                            
      EQUIVALENCE(F(1),FF(1,1,1))                                              
      DATA IPER/1,2,3, 1,3,2, 3,1,2, 2,1,3, 2,3,1, 3,2,1/                      
      DATA IFAC/1,2,4/
      data niold,ipold/0,0/
      save niold,ipold
      REAL*8 XX(3)
C      DO 53 I=1,3
C     PM(I)=PN(I)                                                              
C     N(I)=M(I)                                                                
C  53 CONTINUE                                                                 
C                                                                              
C     *** KB IN THE FIRST ZONE *** 
C      write(6,*) xx,ipold,niold                                            
      DO 1 I=1,3                                                               
      KB(I)=XX(I)                                                              
    2 KB(I)=ABS(KB(I))                                                         
      IF(KB(I).LE.PN(I)) GO TO 1                                               
      KB(I)=KB(I)-PN(I)*2.0                                                    
      GO TO 2                                                                  
    1 CONTINUE                                                                 
      IF(IS.EQ.1) GO TO 4                                                      
C     *** KB(1).GE.KB(2).GE.KB(3) ***                                          
      IF(KB(1).GE.KB(2)) GO TO 3                                               
      W=KB(1)                                                                  
      KB(1)=KB(2)                                                              
      KB(2)=W                                                                  
    3 CONTINUE                                                                 
      IF(IS.EQ.2) GO TO 4                                                      
      IF(KB(2).GE.KB(3)) GO TO 4                                               
      W=KB(2)                                                                  
      KB(2)=KB(3)                                                              
      KB(3)=W                                                                  
      IF(KB(1).GE.KB(2)) GO TO 4                                               
      W=KB(1)                                                                  
      KB(1)=KB(2)                                                              
      KB(2)=W
C     *** FIND THE BLOCK WHERE THE POINT IS ***                                
    4 IX=KB(1)/P(1)                                                            
      IF(IX.EQ.N(1)) IX=N(1)-1                                                 
      IY=KB(2)/P(2)                                                            
      IF(IY.EQ.N(2)) IY=N(2)-1                                                 
      IZ=KB(3)/P(3)                                                            
      IF(IZ.EQ.N(3)) IZ=N(3)-1
C      WRITE(6,*) ' in FUNC',XX
C	  WRITE(6,*) KB
C      WRITE(6,*) IX,IY,IZ	  
      IF(IS.EQ.1) NI=(IX*(N(2)+1)+IY)*(N(3)+1)+IZ+1                            
      IF(IS.EQ.2) NI=((IX*(IX+1))/2+IY)*(N(3)+1)+IZ+1                          
      IF(IS.EQ.3) NI=(IX*(IX+1)*(IX+2))/6                                      
     1  +(IY*(IY+1))/2+IZ+1                                                    
C     *** IF THE BLOCK IS THE SAME, GO TO 40 ***                               
      IF(NI.EQ.NIOLD) GO TO 40                                                 
C     *** PICK UP THE ENERGIES OF 8 POINTS ***                                 
      NK=NI                                                                    
      IF(IS.EQ.1) MX=(N(2)+1)*(N(3)+1)                                         
      IF(IS.EQ.2) MX=(IX+1)*(N(3)+1)                                           
      IF(IS.EQ.3) MX=((IX+1)*(IX+2))/2                                         
      DO 41 KX=1,2                                                             
      DO 42 KY=1,2                                                             
      DO 43 KZ=1,2                                                             
      FF(KX,KY,KZ)=EM(NK)
C      WRITE(6,*) ' NK=',NK,EM(NK)	  
      NK=NK+1                                                                  
   43 CONTINUE                                                                 
      IF(IS.NE.3) NK=NK-2+(N(3)+1)                                             
      IF(IS.EQ.3) NK=NK-2+IY+KY                                                
   42 CONTINUE                                                                 
      IF(KX.LT.2) NK=NI+MX                                                     
   41 CONTINUE                                                                 
   40 CONTINUE
C      WRITE(6,*) F   
C     *** FIND THE TETRAHEDRON ***                                             
      X=KB(1)-P(1)*FLOAT(IX)                                                   
      Y=KB(2)-P(2)*FLOAT(IY)                                                   
      Z=KB(3)-P(3)*FLOAT(IZ)                                                   
      IF(X.GE.Y) GO TO 11                                                      
      IF(Y.GE.Z) GO TO 12
C      WRITE(6,*) X,Y,Z	  
      IP=6                                                                     
      GO TO 13                                                                 
   12 IF(X.GE.Z) GO TO 24                                                      
      IP=5                                                                     
      GO TO 13                                                                 
   24 IP=4                                                                     
      GO TO 13                                                                 
   11 IF(Y.GE.Z) GO TO 21                                                      
      IF(X.GE.Z) GO TO 22                                                      
      IP=3                                                                     
      GO TO 13                                                                 
   22 IP=2                                                                     
      GO TO 13                                                                 
   21 IP=1                                                                     
C     *** IF THE TETRAHEDRON IS THE SAME, GO TO 10 ***                         
   13 IF(NIOLD.NE.NI) GO TO 14                                                 
      IF(IPOLD.EQ.IP) GO TO 10                                                 
   14 IPOLD=IP                                                                 
      NIOLD=NI                                                                 
C     *** DETERMINE THE COEFFICIENTS OF LINEAR FORM ***                        
      K1=IPER(1,IP)                                                            
      K2=IPER(2,IP)                                                            
      K3=IPER(3,IP)                                                            
      I1=IFAC(K1)                                                              
      I2=IFAC(K2)                                                              
      I3=IFAC(K3)                                                              
      A(4)=F(1)                                                                
      JA=I1+1                                                                  
      A(K1)=(F(JA)-F(1))/P(K1)                                                 
      JB=I1+I2+1                                                               
      A(K2)=(F(JB)-F(JA))/P(K2)                                                
      A(K3)=(F(8)-F(JB))/P(K3)                                                 
   10 FUNC=X*A(1)+Y*A(2)+Z*A(3)+A(4)                                         
      RETURN                                                                   
      END                                                                      
      SUBROUTINE CLOCKM(ITIME)
C
C       FOR UNIX
C
      DIMENSION TERY(2)
C
C     etime for AIX FORTRAN
      CALL ETIME_(TERY)
C     etime for SUN FORTARAN
C     CALL ETIME(TERY)
      ITIME=TERY(1)*1000
      RETURN
      END
      SUBROUTINE PTIME(T)
      REAL*8 T
      CALL CLOCKM(IT)
      T=IT/3600000.0
      RETURN
      END
