IMPLICIT REAL*8 (A-H,O-Z)
external func
COMMON/BZPL/KG(3,20),NG,IEOH,INDC,IBR,ICUT
common/enrspl/EM(366145),PN(3),p(3),n(3),js
DIMENSION NRECPT(3,30),RECTAX(4,30),CO(3,2,100)
DIMENSION JB(2,3)
REAL*8 XFU(3,2),WF(3)
REAL*8 XX(3),XC(3),DR(3)
INTEGER M(3),NN(3),NK(3)
DATA XFU/-1.0,-1.0,-1.0,1.0,1.0,1.0/
DATA WF/120.0,120.0,120.0/
DATA AF,BT,GM,ED,EL/10.0,0.0,-20.0,150.0,1500.0/
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)
CALL TSBZEG(NRECPT,RECTAX,NRP,CO,NLIN)
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
C
CALL AYPSTR(51)
CALL AYORIG(30.0,30.0)
CALL TPERSP(XFU,WF,AF,BT,GM,ED,EL)
CALL TPTRAC(1)
CALL TPSCHG(IEOH)
EF=0.41714
CALL TPHILP(EF)
CALL TPHILD(0.0,1,6,5,0)
CALL LINEWD(1.6)
CALL TPCIHL(3)
DO 42 I=1,NLIN
X1=CO(1,1,I)
Y1=CO(2,1,I)
Z1=CO(3,1,I)
X2=CO(1,2,I)
Y2=CO(2,2,I)
Z2=CO(3,2,I)
CALL TPLINE(X1,Y1,Z1,X2,Y2,Z2,30)
42 CONTINUE
C
CALL LINEWD(2.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 TPLINE(0.0,0.0,0.0,-0.75D0,-0.75D0,0.0,30)
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
CALL TPFXYZ(FUNC,EF,NN,NK,1)
CALL LINEWD(1.2)
XC(1)=0.0
XC(2)=0.0
XC(3)=0.0
DR(1)=0.0
DR(2)=0.0
DR(3)=1.0D0
CALL TPSECT(FUNC,EF,DR,XC,0,1,129,129,4,4)
DR(1)=1.0D0
DR(2)=-1.0D0
DR(3)=0.0
CALL TPSECT(FUNC,EF,DR,XC,0,1,129,129,4,4)
DR(1)=1.0D0
DR(2)=0.0
DR(3)=0.0
CALL TPSECT(FUNC,EF,DR,XC,0,1,129,129,4,4)
CALL LINEWD(1.1)
DO 12 IX=1,2
DO 12 IY=1,2
DO 12 IZ=1,2
DR(1)=1.0D0*(3-2*IX)
XC(1)=0.5D0*DR(1)
DR(2)=1.0D0*(3-2*IY)
XC(2)=0.5D0*DR(2)
DR(3)=1.0D0*(3-2*IZ)
XC(3)=0.5D0*DR(3)
CALL TPSECT(FUNC,EF,DR,XC,0,1,257,257,4,4)
12 CONTINUE
CALL AYPEND
STOP
END
FUNCTION FUNCD(XA)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XA(3)
FUNCD=1.0
KX=1000000.0*XA(1)-1000000
KY=1000000.0*XA(2)-1000000
KZ=1000000.0*XA(3)-1000000
CALL TSKFBZ(KX,KY,KZ,1000000,IND)
IF(IND.NE.0) FUNCD=-1.0
RETURN
END
FUNCTION LFUNC(XA)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XA(3)
LFUNC=-1
KX=1000000.0*XA(1)
KY=1000000.0*XA(2)
KZ=1000000.0*XA(3)
CALL TSKFBZ(KX,KY,KZ,1000000,IND)
IF(IND.EQ.0) LFUNC=1
IF(XA(1).LT.0.0.AND.XA(2).LT.XA(1)-1.0D-4
& .AND.XA(3).GT.0.0) 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
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