IMPLICIT REAL*8 (A-H,O-Z)
EXTERNAL FUNC
COMMON/USER1/PI,SQ2,SQ2PI,SQ3
COMMON/USER3/RHO12,RHO13,RHO23
COMMON/USER5/SQ3R
C N = NUMBER OF POINTS OF INTEGRATION FOR ITERATED 2-POINT GAUSS-LEGEND
C NGL = NUMBER OF POINTS FOR GAUSS-LEGENDRE
C NNC = NUMBER OF POINTS FOR NEWTON-COTES
C NX , NY, NZ = NUMBER OF POINTS FOR ITERATED NEWTON COTES
C AX,AY,AZ,BX,BY,BZ = LIMITS OF INTEGRATION
open (unit=10,file='quad.out',status='unknown')
PI=3.1415926538979
SQ2=DSQRT(2.0D0)
SQ3=DSQRT(3.0D0)
SQ3R=1.0D0/SQ3
SQ2PI=DSQRT(2.D0*PI)
RHO12=0.12
RHO13=0.23
RHO23=0.15
AX=-1.0
BX=2.0
AY=-0.8
BY=2.1
AZ=-2.7
BZ=1.96
I=10
CALL TRINT(AX,BX,AY,BY,AZ,BZ,I,FUNC,ANSA,IER)
NGL=11
CALL GAULG3(AX,BX,AY,BY,AZ,BZ,NGL,NGL,NGL,FUNC,ANSG,IER)
NNC=9
C NEWCOT
CALL NEWCT3(AX,BX,AY,BY,AZ,BZ,NNC,NNC,NNC,FUNC,ANSN,IER)
C CLARK3
P1=CLARK3(BX,BY,BZ,RHO12,RHO13,RHO23)
P2=CLARK3(AX,BY,BZ,RHO12,RHO13,RHO23)
P3=CLARK3(BX,BY,AZ,RHO12,RHO13,RHO23)
P4=CLARK3(AX,BY,AZ,RHO12,RHO13,RHO23)
P5=CLARK3(AX,AY,BZ,RHO12,RHO13,RHO23)
P6=CLARK3(AX,AY,AZ,RHO12,RHO13,RHO23)
P7=CLARK3(BX,AY,BZ,RHO12,RHO13,RHO23)
P8=CLARK3(BX,AY,AZ,RHO12,RHO13,RHO23)
C You should verify that if you are integrating a rectangular region in
C 3 dimensions with an algorithm that integrates from a fixed point to
C infinity, you have to integrate over 8 infinite regions and apply formula
C ANSC which follows
ANSC=-P1+P2+P3-P4-P5+P6+P7-P8
P1=DUTT3(BX,BY,BZ,RHO12,RHO13,RHO23)
P2=DUTT3(AX,BY,BZ,RHO12,RHO13,RHO23)
P3=DUTT3(BX,BY,AZ,RHO12,RHO13,RHO23)
P4=DUTT3(AX,BY,AZ,RHO12,RHO13,RHO23)
P5=DUTT3(AX,AY,BZ,RHO12,RHO13,RHO23)
P6=DUTT3(AX,AY,AZ,RHO12,RHO13,RHO23)
P7=DUTT3(BX,AY,BZ,RHO12,RHO13,RHO23)
P8=DUTT3(BX,AY,AZ,RHO12,RHO13,RHO23)
ANSD=-P1+P2+P3-P4-P5+P6+P7-P8
WRITE (10,1003)
1003 FORMAT(///' INTEGRALS '//' TRINT CLARK DUTT 1 GAULG NEWCT'//)
WRITE (10,1004) ANSA,ANSC,ANSD,ANSG,ANSN
1004 FORMAT(1X,5F15.10)
2000 CONTINUE
C Here we do another experiment
AX=-2.0
AY=-2.0
AZ=-2.0
write (10,3005)
3005 format(//' DUTT CLARK')
DO 3000 I=1,5
P6=DUTT3(AX,AY,AZ,RHO12,RHO13,RHO23)
P7=CLARK3(AX,AY,AZ,RHO12,RHO13,RHO23)
WRITE (10,3001) P6,P7
3001 FORMAT(1X,2F10.6)
AX=AX+1.
AY=AY+1.
3000 AZ=AZ+1.
close (10)
STOP
END
FUNCTION FUNC(X,Y,Z)
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/USER1/PI,SQ2,SQ2PI,SQ3
COMMON/USER3/RHO12,RHO13,RHO23
DELTA=1+2.*RHO12*RHO13*RHO23-RHO12**2-RHO13**2-RHO23**2
EXX=((1-RHO23**2)*X**2+(1.-RHO13**2)*Y**2+(1.-RHO12**2)*Z**2
1 +2.*((RHO13*RHO23-RHO12)*X*Y+(RHO12*RHO23-RHO13)*X*Z+
2 (RHO12*RHO13-RHO23)*Y*X))/(2.*DELTA)
FUNC=DEXP(-EXX)/(SQ2PI**3*DSQRT(DELTA))
RETURN
END
******************************************************************************
OUTPUT FILE
******************************************************************************
INTEGRALS
TRINT CLARK DUTT GAULG NEWCT
0.6243602183 0.6071457172 0.6171489376 0.6243495086 0.6239873387
DUTT CLARK
0.935255 0.940097
0.621245 0.612904
0.165022 0.163101
0.010244 0.012209
0.000106 0.000217