MINICERN/mathlib/gen/divon/quad.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:03:25  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       REAL FUNCTION QUAD(D,DEG,LOWER,UPPER,FUN)
  11       INTEGER D,DEG
  12       REAL LOWER(D),UPPER(D),FUN
  13       EXTERNAL FUN
  14       INTEGER D0,DP1,TWOD,PT,NPTS,TIMES,KMAX
  15       REAL X(10,202),Y(10),Z(10),R
  16       DOUBLE PRECISION SUM0,SUM1,SUM2,B0,B1,B2
  17       REAL PI,TWOPI
  18       DATA PI,TWOPI /3.14159E+0,6.28319E+0/
  19       DATA D0 /0/
  20       IF(D.EQ.D0) GOTO 130
  21       D0=D
  22       SQ23=SQRT(2.0E+0/3.0E+0)
  23       SQ13I=1.0E+0/SQRT(3.0E+0)
  24       DP1=D+1
  25       TWOD=2*D
  26       KMAX=D/2
  27       DO 20 K=1,KMAX
  28       DO 10 J=1,DP1
  29       I=J-1
  30       X(2*K-1,J)=SQ23*COS(TWOPI*I*K/DP1)
  31       X(2*K,J)=SQ23*SIN(TWOPI*I*K/DP1)
  32  10   CONTINUE
  33  20   CONTINUE
  34       IF(MOD(D,2).NE.1) GOTO 40
  35       DO 30 I=1,DP1
  36       X(D,I)=-SQ13I*(-1)**I
  37  30   CONTINUE
  38  40   DO 60 K=1,KMAX
  39       DO 50 I=1,TWOD
  40       X(2*K-1,I+DP1)=SQ23*COS((2*K-1)*I*PI/D)
  41       X(2*K,I+DP1)=SQ23*SIN((2*K-1)*I*PI/D)
  42  50   CONTINUE
  43  60   CONTINUE
  44       IF(MOD(D,2).NE.1) GOTO 80
  45       DO 70 I=1,TWOD
  46       X(D,I+DP1)=SQ13I*(-1)**I
  47  70   CONTINUE
  48  80   I=1
  49       GOTO 100
  50  90   I=I+1
  51  100  IF((I).GT.(DP1+TWOD)) GOTO 120
  52       DO 110 J=1,D
  53       X(J,I)=0.5E+0*(X(J,I)+1.0E+0)
  54  110  CONTINUE
  55       GOTO 90
  56  120  R=0.5E+0*SQRT(3.0E+0/5.0E+0)
  57       B0=(25.0E+0*D**2-115.0E+0*D+162.0E+0)/162.0E+0
  58       B1=(70.0E+0-25.0E+0*D)/162.0E+0
  59       B2=25.0E+0/324.0E+0
  60  130  IF(DEG.NE.2) GOTO 140
  61       PT=0
  62       NPTS=DP1
  63       GOTO 180
  64  140  IF(DEG.NE.3) GOTO 150
  65       PT=DP1
  66       NPTS=TWOD
  67       GOTO 180
  68  150  IF(DEG.NE.5) GOTO 160
  69       GOTO 210
  70  160  WRITE(6,170) DEG
  71  170  FORMAT(' ILLEGAL DEGREE =',G13.5)
  72       STOP
  73  180  QUAD=0.0E+0
  74       DO 200 I=1,NPTS
  75       DO 190 J=1,D
  76       Y(J)=(UPPER(J)-LOWER(J))*X(J,I+PT)+LOWER(J)
  77  190  CONTINUE
  78       QUAD=QUAD+FUN(D,Y)
  79  200  CONTINUE
  80       QUAD=QUAD/NPTS
  81       RETURN
  82  210  SUM0=0.0E+0
  83       SUM1=SUM0
  84       SUM2=SUM1
  85       DO 220 J=1,D
  86       Y(J)=0.5E+0*(UPPER(J)+LOWER(J))
  87  220  CONTINUE
  88       SUM0=FUN(D,Y)
  89       DO 230 J=1,D
  90       Z(J)=(UPPER(J)-LOWER(J))*R
  91  230  CONTINUE
  92       DO 240 I=1,D
  93       Y(I)=Y(I)+Z(I)
  94       SUM1=SUM1+FUN(D,Y)
  95       Y(I)=Y(I)-2.0E+0*Z(I)
  96       SUM1=SUM1+FUN(D,Y)
  97       Y(I)=Y(I)+Z(I)
  98  240  CONTINUE
  99       DO 270 I=2,D
 100       Y(I)=Y(I)+Z(I)
 101       JMAX=I-1
 102       DO 260 TIMES=1,2
 103       DO 250 J=1,JMAX
 104       Y(J)=Y(J)+Z(J)
 105       SUM2=SUM2+FUN(D,Y)
 106       Y(J)=Y(J)-2.0E+0*Z(J)
 107       SUM2=SUM2+FUN(D,Y)
 108       Y(J)=Y(J)+Z(J)
 109  250  CONTINUE
 110       Y(I)=Y(I)-2.0E+0*Z(I)
 111  260  CONTINUE
 112       Y(I)=Y(I)+3.0E+0*Z(I)
 113  270  CONTINUE
 114       QUAD=B0*SUM0+B1*SUM1+B2*SUM2
 115       RETURN
 116       END