C     THIS  PROGRAM  EVALUATES ALL THE  EVEN  CORRELATION  FUNCTIONS
C     UP TO X50 BY MEANS OF PFFAFIANS U1, U6, U7, U11AND NONLINEAR
C     IDENTITIES  FOR  X8,  X11, X20.
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      EXTERNAL  FCT,  FCT1,  FCT2,  FCT3,  FCT4,  FCT5,  FCT6
      EXTERNAL  FCT7, FCT8, FCT9
C      open(unit=17,file="tempout.dat",status="unknown")
      OPEN(UNIT=17, FILE="TEMPOUT.DAT", STATUS="UNKNOWN")
C     L  MUST  BE  SPECIFIED.
      PI=3.14159  26535  89793 D0
      TWOPIR=0.5D0/PI
      XL=-PI
      XU=PI
      MQ=2
      AXU=(XU-XL)/DFLOAT(MQ)
      AKC=0.65847 89484  624084D0
      AKC=0.5D0*DLOG(2.D0+DSQRT(3.D0))
      ALC=0.27465 30721 670275D0
      ALC=0.25D0*DLOG(2.D0*DCOSH(2.D0*AKC)-1.D0)
      WRITE(17,1000) ALC,AKC
1000  FORMAT(15X, 'ALC  = ',F18.16,2X,'AKC  = ',F18.16)
      AL=ALC
      DO 20  I=1, 10
C     THIS PART OF PROGRAM DEFINES THE TERMS NEEDED TO EVALUATE
C     THE CORRELATION FUNCTIONS U1, U6, U7 USING THE PFAFFIAN
C     TECHNIQUE.
      XL1=XL
      XU1=AXU+XL
      V=DTANH(AL)
      SA=8.D0*V*V*V+(1.D0+V*V)**3
      SB=-2.D0*V*(1.D0-V*V)**2
      VV=V*V
      VVV=VV*V
C      T11=(1.D0-VV)*(4.d0*VV+V*(1.D0+VV)**2)
      XR=0.5D0*(DEXP(4.D0*AL)+1.D0)
      AK=0.5D0*DLOG(XR+DSQRT(XR**2-1.D0))
      T1=DTANH(AK)
      T3=DTANH(3.D0*AK)
      AN=0.25D0*(T3+T1)
      BN=0.25D0*(T3-3.D0*T1)
      ASQ=AN*AN
      BSQ=BN*BN
      AB=AN*BN
      A=2.D0*V*(1.D0+V**2)
      C=(1.D0-V)**2
      CC=C*C
      B=C*VV
      F=B-C
      D=B+C
      WDI=1.D0-V
      WSU=1.D0+V
      VDI=V*WDI
      VSU=V*WSU
      TOR=VVV-VV+V+V
      TAR=VV+VV-V+1.D0
      FVC=4.D0*V*C
      EVC=1.D0/(8.D0*VV*CC)
      DVC=2.D0*V*C
      APF=A*A+F*F
      DMF=D*D-F*F
      TAD=2.D0*A*D
      CONS=WDI*TWOPIR
C     THIS PART OF PROGRAM EVALUATES THE INTEGRALS
C     APPEARING IN EACH OF THE ELEMENTS OF THE PFAFFIAN.
      G=1.D0-V-V+6.D0*VV-2.D0*VVV+VV*VV
      YT=0.D0
      Y1T=0.D0
      Y2T=0.D0
      RI1=0.D0
      RI2=0.D0
      RI3=0.D0
      RI4=0.D0
      Y7T=0.D0
      Y8T=0.D0
      Y9T=0.D0
      DO  10  NA=1, MQ
      CALL PFAFF(XL1,XU1,FCT,Y)
      CALL PFAFF(XL1,XU1,FCT1,Y1)
      CALL PFAFF(XL1,XU1,FCT2,Y2)
      CALL PFAFF(XL1,XU1,FCT3,Y3)
      CALL PFAFF(XL1,XU1,FCT4,Y4)
      CALL PFAFF(XL1,XU1,FCT5,Y5)
      CALL PFAFF(XL1,XU1,FCT6,Y6)
      YT=YT+Y
      Y1T=Y1T+Y1
      Y2T=Y2T+Y2
      RI1=Y3+RI1
      RI2=Y4+RI2
      RI3=-Y5/FVC+RI3
      RI4=-Y6/FVC+RI4
      XL1=XU1
      XU1=XU1+AXU
 10   CONTINUE
      FRS=-VDI*RI1-TOR*RI2+TAR*RI3-VDI*RI4
      FRT=V*VDI*RI1+VSU*(RI2+RI3)-WDI*RI4
      FLS=-WDI*RI1+VSU*(RI2+RI3)+V*VDI*RI4
      FLT=-VDI*RI1+TAR*RI2-TOR*RI3-VDI*RI4
      U1=YT*TWOPIR
      APL=Y1T*TWOPIR
      AMI=Y2T*TWOPIR
      U7=U1*U1-APL*AMI
      A12=-U1
      A13=-CONS*FRS
      A14=CONS*FRT
      A23=CONS*FLS
      A24=CONS*FLT
      A34=U1
      U6=-(A12*A34-A13*A24+A14*A23)
C    THEFOLLOWING PROGRAM CALCULATES U11 BY MEANS OF PFFAFIAN
      XL1=XL
      XU1=AXU+XL
      DO 11 KA=1, MQ
      CALL PFAFF(XL1,XU1,FCT7,Y7)
      CALL PFAFF(XL1,XU1,FCT8,Y8)
      CALL PFAFF(XL1,XU1,FCT9,Y9)
      Y7T=Y7T+Y7
      Y8T=Y8T+Y8
      Y9T=Y9T+Y9
      XL1=XU1
      XU1=XU1+AXU
 11   CONTINUE
      W2=Y7T*TWOPIR
      W3=Y8T*TWOPIR
      W4=Y9T*TWOPIR
      Z=1.D0-V*V
      ST1=-U1*(-U1*U1+Z*Z*W2*W3+Z*Z*W4*W4)
      ST2=-Z*Z*W4*(W4*U1-Z*W4*W3-Z*W2*W4)
      ST3=Z*Z*W2*(-W3*U1+Z*W4*W4-Z*W2*W2)
      ST4=-Z*Z*W3*(Z*W3*W3-Z*W4*W4+U1*W2)
      ST5=-Z*Z*W4*(-Z*W3*W4-Z*W4*W2+W4*U1)
      U11=ST1+ST2+ST3+ST4+ST5
      X82=U11
C     U6  FROM  STEPHENSON, J. Math. Phys.  5, 1009(1964);
C     ibid  7, 1123(1966);\, ibid 11, 413(1970)
      U6=-(A12*A34-A13*A24+A14*A23)
C     IDENTITY (17)
      X1=AN+(2.D0*AN+BN)*U1
C     GKB  NONLINEAR  IDENTITIES, Physica 93A,  138(1978)
C     X8=X1**2+X1*X3-X2**2,  		(NL-1)
C     X11=X2*(X1+X1-X3),           		(NL-2)
C     X20=X1**2+X2**2-X3**2 		(NL-3)
C     X11=AN*(X2+X2+X33)+BN*X5 	LINEAR IDENTITY (68)
C     X5=U6,  X6=U7,  X33=U2,  X34=U3,  X35=U4,  X81=U5
C     X2=U1,  X37=U8,  X38=U9, X41=U10,  X82=U11
      X2=U1
C     IDENTITY (115)
      X7=3.D0*AN*U1+BN
      X5=U6
      X6=U7
C    LINEAR IDENTITIES (1), (68) AND NONLINEAR IDENTITY (NL-2)
C    ARE SOLVED FOR X3, X11 AND U2
      U2=(2.D0*U1*(X1-AN*U1-AN)-(BN+AN*U1)*U6)/(AN+BN*U1)
      X33=U2
      X3=AN*(2.D0*X2+X5)+BN*X33
      X11=AN*(2.D0*X2+U2)+BN*X5
C     NONLINEAR IDENTITY (NL-1)
      X8=X1**2+X1*X3-X2**2
C     NONLINEAR IDENTITY (NL-3)
      X20=X1**2+X2**2-X3**2
C     IDENTITIES (5), (76)
      X4=((AN*X2-BN*X8)-(ASQ-AB)*(X1+X3))/(ASQ-BSQ)
      X17=(X2-AN*(X1+X3+X4))/BN
C     IDENTITY (89)
      X9=AN*(X1+X4+X17)+BN*X3
C     IDENTITY (68)
C     IDENTITY (72)
      X12=AN*(X2+X5+X33)+BN*X2
C     IDENTITY (160)
      X10=AN*(X3+X4+X17)+BN*X1
C     IDENTITIES (23), (49) ARE SOLVED FOR X13 AND X28
      X13=((AN*X8-BN*X20)-(ASQ-AB)*(X1+X11))/(ASQ-BSQ)
      X28=(X8-AN*(X1+X11+X13))/BN
C IDENTY (69) YIELDS X34
      X34=(X4-AN*(X2+X5+X6))/BN
C IDENTITY (70)
      X14=AN*(X2+X6+X34)+BN*X5
C IDENTITY (100)
      X15=AN*(X5+X6+X34)+BN*X2
C IDENTITY (60)
      X39=(X17-AN*(X2+X10)-BN*X10)/AN
C IDENTITY (56)
      X21=AN*(X1+X13+X28)+BN*X11
C IDENTITY (24)
      X31=(X9-AN*(X1+X12+X14))/BN
C IDENTITY (18)
      X37=(X17-AN*(X2+X33+X34))/BN
C IDENTITY (54)
      X24=AN*(X1+X14+X31)+BN*X12
C IDENTITY (146)
      X22=(X11-AN*(X2+X20))/(AN+BN)
C IDENTITY (144)
      X25=(X14-AN*(X6+X24))/(AN+BN)
C IDENTITY (145)
      X26=(X13-AN*(X5+X21)-BN*X22)/AN
C IDENTITY (148)
      X23=(X12-AN*(X5+X21)-BN*X26)/AN
C IDENTITY (154)
      X16=AN*(X6+2.D0*X25)+BN*X24
C IDENTITY (142)
      X18=AN*(X2+X23+X26)+BN*X21
C IDENTITY (141)
      X19=AN*(X2+X22+X22)+BN*X20
C IDENTITY (27)
      X36=(X10-AN*(X4+X18+X19))/BN
C IDENTITY (40)
      X27=AN*(X6+X34+X37)+BN*X34
C IDENTITY (46)
      X40=(X25-AN*(X3+X16)-BN*X18)/AN
C IDENTITY (20)
      X41=(X19-AN*(X2+X34+X34))/BN
C IDENTITY (39)
      X32=AN*(X2+X34+X41)+BN*X34
C IDENTITY (129)
      X42=AN*(X10+X10+X39)+BN*X2
C IDENTITY (121)
      X46=AN*(X11+X13+X28)+BN*X1
C IDENTITY (125)
      X47=AN*(X18+X19+X36)+BN*X4
C IDENTITY (123)
      X49=AN*(X16+X18+X40)+BN*X3
C IDENTITY (122)
      X48=AN*(X12+X14+X31)+BN*X1
C IDENTITY (119)
      X50=AN*(X34+X34+X37)+BN*X6
C IDENTITY (177)
      X43=AN*(X20+X22+X22)+BN*X2
C IDENTITY (178)
      X44=AN*(X21+X23+X26)+BN*X5
C IDENTITY (180)
      X45=AN*(X24+X25+X25)+BN*X6
C IDENTITY (6)
      X51=(X10-AN*(X7+X12+X15))/BN
C IDENTITY (92)
      X58=(X43-BN*X26-AN*(X9+X22))/AN
C IDENTITY (93)
      X59=(X44-BN*X25-AN*(X8+X23))/AN
C IDENTITY (11)
      X74=(X39-AN*(X17+X36+X40))/BN
C IDENTITY (104)
C      X63=AN*(X12+X15+X51)+BN*X7
C IDENTITY (118)
      X53=AN*(X34+X34+X41)+BN*X2
C IDENTITY (135)
C      X60=(X33-AN*(X17+X36+X40))/BN
C IDENTITY (138)
C      X61=(X34-AN*(X15+X31+X40))/BN
C IDENTITY (22)
      X81=0.5D0*(X74-AN*X37-BN*X82)/AN
C IDENTITY (3)
      X38=(X36-AN*(X33+X37)-BN*X81)/AN
C IDENTITY (19)
      X35=(X18-AN*(X5+X34)-BN*X38)/AN
C IDENTITY (41)
      X30=AN*(X5+X35+X38)+BN*X34
C IDENTITY (38)
      X29=AN*(X5+X34+X38)+BN*X35
      WRITE(17,104)AL, AK
      WRITE(17,30)U1,U2,U6,U7,U11
      WRITE(17,31)X1,X2,X3,X4,X5
      WRITE(17,32)X6,X7,X8,X9,X10
      WRITE(17,33)X11,X12,X13,X14,X15
      WRITE(17,34)X16,X17,X18,X19,X20
      WRITE(17,35)X21,X22,X23,X24,X25
      WRITE(17,36)X26,X27,X28,X29,X30
      WRITE(17,37)X31,X32,X33,X34,X35
      WRITE(17,38)X36,X37,X38,X39,X40
      WRITE(17,39)X41,X42,X43,X44,X45
      WRITE(17,40)X46,X47,X48,X49,X50
 30   FORMAT(1X,'  U   = ', 5(2X,F13.10))
 31   FORMAT(1X,' X1   = ',5(2X,F13.10))
 32   FORMAT(1X,' X6   = ',5(2X,F13.10))
 33   FORMAT(1X,' X11  = ',5(2X,F13.10))
 34   FORMAT(1X,' X16  = ',5(2X,F13.10))
 35   FORMAT(1X,' X21  = ',5(2X,F13.10))
 36   FORMAT(1X,' X26  = ',5(2X,F13.10))
 37   FORMAT(1X,' X31  = ',5(2X,F13.10))
 38   FORMAT(1X,' X36  = ',5(2X,F13.10))
 39   FORMAT(1X,' X41  = ',5(2X,F13.10))
 40   FORMAT(1X,' X46  = ',5(2X,F13.10))
104   FORMAT(20X,'AL = ',F13.10,2X,'AK = ',F13.10)
      AL=AL+0.1d0
20    CONTINUE
      STOP
      END
C
      REAL*8 FUNCTION FCT(X)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      FNM=A-D*DCOS(X)
      VFO=DSQRT(APF-TAD*DCOS(X)+DMF*(DCOS(X)**2))
      FCT=FNM/VFO
      RETURN
      END
C
      REAL*8 FUNCTION FCT1(X)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      FNM1=(A-D*DCOS(X))*DCOS(X)-F*DSIN(X)**2
      VF1=DSQRT(APF-TAD*DCOS(X)+DMF*(DCOS(X)**2))
      FCT1=FNM1/VF1
      RETURN
      END
C
      REAL*8 FUNCTION FCT2(X)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      FNM2=(A-D*DCOS(X))*DCOS(X)+F*DSIN(X)**2
      VF2=DSQRT(APF-TAD*DCOS(X)+DMF*(DCOS(X)**2))
      FCT2=FNM2/VF2
      RETURN
      END
C
      REAL*8 FUNCTION FCT3(X)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      VF3=DSQRT(APF-TAD*DCOS(X)+DMF*(DCOS(X)**2))
      FCT3=1.D0/VF3
      RETURN
      END
C
      REAL*8 FUNCTION FCT4(X)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      VF4=DSQRT(APF-TAD*DCOS(X)+DMF*(DCOS(X)**2))
      FCT4=DCOS(X)/VF4
      RETURN
      END
C
      REAL*8 FUNCTION FCT5(X)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      VF5=DSQRT(APF-TAD*DCOS(X)+DMF*(DCOS(X)**2))
      FCT5=1.D0-(G-DVC*DCOS(X))/VF5
      RETURN
      END
C
      REAL*8 FUNCTION FCT6(X)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      VFO=DSQRT(APF-TAD*DCOS(X)+DMF*(DCOS(X)**2))
      VF6=1.D0-(G-DVC*DCOS(X))/VFO
      FCT6=(2.D0*DCOS(X)-1.D0)*VF6
      RETURN
      END
C     (2,2)RD
      REAL*8  FUNCTION  FCT7(X)
      IMPLICIT  REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      AA=SA+SB*DCOS(X)
      BB=SB*(1.D0+DCOS(X))
      CC=-SB*DSIN(X)
      FAC=BB*BB+CC*CC
      F1=DSQRT(AA*AA-BB*BB-CC*CC)
      F2=F1-AA
      XP=-(1.D0+V*V+2.D0*V*V*V)
      Y=V*V*(1.D0+V*V)+2.D0*V
      Z=V*(1.D0-V*V)
      XN1=XP*DCOS(2.D0*X)+Y*DCOS(X)
      XN2=BB*(DCOS(X)+DCOS(2.D0*X))+CC*(DSIN(X)-DSIN(2.D0*X))
      XN3=Z*F2*XN2/FAC
      FCT7=(XN1+XN3)/F1
      RETURN
      END
C     (0,1)RS
      REAL*8  FUNCTION  FCT8(X)
      IMPLICIT  REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      AA=SA+SB*DCOS(X)
      BB=SB*(1.D0+DCOS(X))
      CC=-SB*DSIN(X)
      F1=DSQRT(AA*AA-BB*BB-CC*CC)
      F2=F1-AA
      FAC=BB*BB+CC*CC
      XP=-(1.D0+V*V+2.D0*V*V*V)
      Y=V*V*(1.D0+V*V)+2.D0*V
      Z=V*(1.D0-V*V)
      AL=F2*(BB*(XP+Z*DCOS(X))+CC*Z*DSIN(X))/FAC
      DL=(BB*BB-CC*CC)*(Y+Z*DCOS(X))-2.D0*Z*BB*CC*DSIN(X)
      XN1=(F2*F2)*DL/(FAC*FAC)
      FCT8=(AL+XN1)/F1
      RETURN
      END
C     (1,2)LD
      REAL*8  FUNCTION  FCT9(X)
      IMPLICIT  REAL*8(A-H,O-Z)
      COMMON A,D,F,G,APF,TAD,DMF,DVC, V, SA, SB
      AA=SA+SB*DCOS(X)
      BB=SB*(1.D0+DCOS(X))
      CC=-SB*DSIN(X)
      FAC=BB*BB+CC*CC
      F1=DSQRT(AA*AA-BB*BB-CC*CC)
      F2=F1-AA
      P=1.D0-V*V
      Q=-V*V*(1.D0-V*V)
      R=-(V+2.D0*V*V+V*V*V)
      XN1=F2*(BB*DCOS(X)+CC*DSIN(X))*P/FAC
      XN2=F2*Q*(BB*DCOS(2.D0*X)-CC*DSIN(2.D0*X))/FAC
      XN3=R*(DCOS(2.D0*X)+DCOS(X))
      FCT9=(XN1+XN2+XN3)/F1
      RETURN
      END
C
      SUBROUTINE PFAFF(XL,XU,FCT,Y)
      IMPLICIT REAL*8(A-H,O-Z)
      A=0.5D0*(XL+XU)
      B=XU-XL
      C=0.49863193092474078D0*B
      Y=0.35093050047350483D-2*(FCT(A+C)+FCT(A-C))
      C=0.49280575577263417D0*B
      Y=Y+0.8137197365452835D-2*(FCT(A+C)+FCT(A-C))
      C=0.48238112779375322D0*B
      Y=Y+0.12696032654631030D-1*(FCT(A+C)+FCT(A-C))
      C=0.46745303796886984D0*B
      Y=Y+0.17136931456510717D-1*(FCT(A+C)+FCT(A-C))
      C=0.44816057788302606D0*B
      Y=Y+0.21417949011113340D-1*(FCT(A+C)+FCT(A-C))
      C=0.42468380686628499D0*B
      Y=Y+0.25499029631188088D-1*(FCT(A+C)+FCT(A-C))
      C=0.39724189798397120D0*B
      Y=Y+0.29342046739267774D-1*(FCT(A+C)+FCT(A-C))
      C=0.36609105937014484D0*B
      Y=Y+0.32911111388180923D-1*(FCT(A+C)+FCT(A-C))
      C=0.33152213346510760D0*B
      Y=Y+0.36172897054424253D-1*(FCT(A+C)+FCT(A-C))
      C=0.29385787862038116D0*B
      Y=Y+0.39096947893535153D-1*(FCT(A+C)+FCT(A-C))
      C=0.25344995446611470D0*B
      Y=Y+0.41655962113473378D-1*(FCT(A+C)+FCT(A-C))
      C=0.21067563806531767D0*B
      Y=Y+0.43826046502201906D-1*(FCT(A+C)+FCT(A-C))
      C=0.16593430114106382D0*B
      Y=Y+0.45586939347881942D-1*(FCT(A+C)+FCT(A-C))
      C=0.11964368112606854D0*B
      Y=Y+0.46922199540402283D-1*(FCT(A+C)+FCT(A-C))
      C=0.7223598079139825D-1*B
      Y=Y+0.47819360039637430D-1*(FCT(A+C)+FCT(A-C))
      C=0.24153832843869158D-1*B
      Y=B*(Y+0.48270044257363900D-1*(FCT(A+C)+FCT(A-C)))
      RETURN
      END
%%%%%%%%%%%%%%  END  OF  PROGRAM  FOR  HONEYCOMB  CORRELATION  FUNCTIONS
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% Appendix 5, Part 3. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

