IMPLICIT REAL*8 (A-H,O-Z)
      external func
      COMMON/DICUT/ICUT
      DIMENSION NRECPT(3,30),RECTAX(4,30)
C
      REAL*8 XFU(3,2),WF(3)
      REAL*8 XX(3),XC(3),DR(3)
      real*8 func
      common/enrspl/EM(366145),PN(3),p(3),n(3),js
      INTEGER M(3),NN(3),NK(3)
      DIMENSION CO(3,2,100),RKG(4,20)
      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
      WRITE(6,601) A, B, C
  601 FORMAT(' A=',F9.5,' B=',F9.5,' C=',F9.5)
      READ(1,*) CA,CB,CC
      CALL TSLATC(A,B,C,CA,CB,CC0)
      NB=6
      DO 10 J=1,NB
      READ(31) IXMAX,NB,NNEE
      READ(31) (EM(K),K=1,NNEE)
      WRITE(6,*) IXMAX,NB,NNEE
   10 CONTINUE
      PN(1)=1.0
      PN(2)=1.0
      PN(3)=1.0
      n(1)=IXMAX-1
      n(2)=n(1)
      n(3)=n(1)
      JS=3
      DO 51 I=1,3
      P(I)=Pn(I)/FLOAT(N(I))
   51 CONTINUE
      XFU(1,1)=0.0
      XFU(2,1)=0.0
      XFU(3,1)=-1.0D0
      XFU(1,2)=2.0D0
      XFU(2,2)=2.0D0
      XFU(3,2)=1.0D0
      SIZE=120.0
      WF(1)=SIZE
      WF(2)=WF(1)
      WF(3)=WF(1)
      IEOH=1
      AF=30.0
      EF=0.41714
      BT=0.0
      GM=-20.0
      ED=150.0
      EL=1500.0
      CALL AYPSTR(96)
      CALL AYORIG(30.0,30.0)
      CALL TPERSP(XFU,WF,AF,BT,GM,ED,EL)
      CALL TPTRAC(1)   
      CALL TPSCHG(IEOH)
      CALL TPHILP(EF)
      CALL TPHILD(0.0,1,6,5,0)
      CALL LINEWD(1.6)
      CALL TPCIHL(3)
      CALL LINEWD(0.1)
      CALL TPCIHL(0)
      CALL TPHILD(0.0,0,0,0,0)
      DO 11 I=1,3
      NN(I)=129
      NK(I)=1
   11 CONTINUE
      ICUT=0
      CALL TPFXYZ(FUNC,EF,NN,NK,1)
      CALL LINEWD(1.3)
         ICUT=1
         DR(1)=1.0D0
         XC(1)=0.5D0
         DR(2)=1.0D0
         XC(2)=0.5D0
         DR(3)=1.0D0
         XC(3)=0.5D0
      CALL TPSECT(FUNC,EF,DR,XC,0,1,257,257,4,4)
      CALL TPHILD(EF,1,4,4,0)
      CALL TPCIHL(2)
      CALL LINEWD(1.4)
        ICUT=2
         DR(1)=1.0D0
         XC(1)=0.5D0
         DR(2)=1.0D0
         XC(2)=0.5D0
         DR(3)=-1.0D0
         XC(3)=-0.5D0
      CALL TPSECT(FUNC,EF,DR,XC,0,0,257,257,4,4)
      CALL LINEWD(1.5)
        ICUT=0
        XC(1)=1.0
        XC(2)=1.0
        XC(3)=0.0
        DR(1)=1.0
        DR(2)=1.0
        DR(3)=0.0D0
      CALL TPSECT(FUNC,EF,DR,XC,0,0,129,129,4,4)
      CALL AYPEND
      STOP
      END
      FUNCTION FUNCD(XA)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION XA(3)
         WW=FUNC(XA)
         FUNCD=WW
      RETURN
      END
      FUNCTION LFUNC(XA)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/DICUT/ICUT
      DIMENSION XA(3)
      LFUNC=-1
      IF(ICUT.EQ.0) THEN
         KX=1000000.0*XA(1)
         KY=1000000.0*XA(2)
         KZ=1000000.0*XA(3)
         CALL TSKFBZ(KX,KY,KZ,1000000,IND)
         IF(IND.NE.0.AND.XA(3).GT.0.0) LFUNC=1
      ELSE IF(ICUT.EQ.1) THEN
         KX=1000005.0*XA(1)
         KY=1000005.0*XA(2)
         KZ=1000005.0*XA(3)
         CALL TSKFBZ(KX,KY,KZ,1000000,IND)
         IF(IND.NE.0.AND.XA(3).GT.0.0) LFUNC=1
         IF(ABS(XA(2)-XA(1)).GT.1.0D0) LFUNC=1
         IF(IND.EQ.0.AND.XA(3).LT.0.0) LFUNC=1
      ELSE IF(ICUT.EQ.2) THEN
         IF(ABS(XA(2)-XA(1)).GT.1.0D0) LFUNC=1
         IF(XA(3).GT.0.0) LFUNC=1
      END IF   
      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   
      REAL*8 function FUNC(XX)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 KB 
C      DIMENSION EM(NN),PM(3),N(3),XX(3),M(3),PN(3)
      common/enrspl/em(366145),pm(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.PM(I)) GO TO 1                                               
      KB(I)=KB(I)-PM(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                                                 
      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)                                                      
      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     *** 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                                                      
      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