MINICERN/mathlib/gen/d/cauchy64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:14  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       FUNCTION DCAUCH(F,A,B,S,EPS)
  12 #endif
  13 #if !defined(CERNLIB_DOUBLE)
  14       FUNCTION RCAUCH(F,A,B,S,EPS)
  15 #endif
  16 #include "gen/imp64.inc"
  17       CHARACTER NAME*(*)
  18       CHARACTER*80 ERRTXT
  19 #if defined(CERNLIB_DOUBLE)
  20       PARAMETER (NAME = 'DCAUCH')
  21 #endif
  22 #if !defined(CERNLIB_DOUBLE)
  23       PARAMETER (NAME = 'RCAUCH')
  24 #endif
  25       EXTERNAL F
  26
  27       DIMENSION X(12),W(12)
  28
  29       PARAMETER (Z1 = 1, HF = Z1/2, CST = 5*Z1/1000)
  30
  31       DATA X( 1) /9.6028985649753623D-1/, W( 1) /1.0122853629037626D-1/
  32       DATA X( 2) /7.9666647741362674D-1/, W( 2) /2.2238103445337447D-1/
  33       DATA X( 3) /5.2553240991632899D-1/, W( 3) /3.1370664587788729D-1/
  34       DATA X( 4) /1.8343464249564980D-1/, W( 4) /3.6268378337836198D-1/
  35       DATA X( 5) /9.8940093499164993D-1/, W( 5) /2.7152459411754095D-2/
  36       DATA X( 6) /9.4457502307323258D-1/, W( 6) /6.2253523938647893D-2/
  37       DATA X( 7) /8.6563120238783174D-1/, W( 7) /9.5158511682492785D-2/
  38       DATA X( 8) /7.5540440835500303D-1/, W( 8) /1.2462897125553387D-1/
  39       DATA X( 9) /6.1787624440264375D-1/, W( 9) /1.4959598881657673D-1/
  40       DATA X(10) /4.5801677765722739D-1/, W(10) /1.6915651939500254D-1/
  41       DATA X(11) /2.8160355077925891D-1/, W(11) /1.8260341504492359D-1/
  42       DATA X(12) /9.5012509837637440D-2/, W(12) /1.8945061045506850D-1/
  43
  44       IF(S .EQ. A .OR. S .EQ. B) THEN
  45        H=0
  46        WRITE(ERRTXT,101) S
  47        CALL MTLPRT(NAME,'D104.1',ERRTXT)
  48       ELSEIF(S .LT. MIN(A,B) .OR. S .GT. MAX(A,B)) THEN
  49 #if defined(CERNLIB_DOUBLE)
  50        H=DGAUSS(F,A,B,EPS)
  51 #endif
  52 #if !defined(CERNLIB_DOUBLE)
  53        H= GAUSS(F,A,B,EPS)
  54 #endif
  55       ELSE
  56        IF(2*S .LE. A+B) THEN
  57 #if defined(CERNLIB_DOUBLE)
  58         H=DGAUSS(F,2*S-A,B,EPS)
  59 #endif
  60 #if !defined(CERNLIB_DOUBLE)
  61         H= GAUSS(F,2*S-A,B,EPS)
  62 #endif
  63         B0=S-A
  64        ELSE
  65 #if defined(CERNLIB_DOUBLE)
  66         H=DGAUSS(F,A,2*S-B,EPS)
  67 #endif
  68 #if !defined(CERNLIB_DOUBLE)
  69         H= GAUSS(F,A,2*S-B,EPS)
  70 #endif
  71         B0=B-S
  72        ENDIF
  73        C=CST/B0
  74        BB=0
  75
  76     1  AA=BB
  77        BB=B0
  78     2  C1=HF*(BB+AA)
  79        C2=HF*(BB-AA)
  80        C3=S+C1
  81        C4=S-C1
  82        S8=0
  83        DO 3 I = 1,4
  84        U=C2*X(I)
  85     3  S8=S8+W(I)*((F(C3+U)+F(C4-U))+(F(C3-U)+F(C4+U)))
  86        S8=C2*S8
  87        S16=0
  88        DO 4 I = 5,12
  89        U=C2*X(I)
  90     4  S16=S16+W(I)*((F(C3+U)+F(C4-U))+(F(C3-U)+F(C4+U)))
  91        S16=C2*S16
  92        IF(ABS(S16-S8) .LE. EPS*(1+ABS(S16))) GO TO 5
  93        BB=C1
  94        IF(1+ABS(C*C2) .NE. 1) GO TO 2
  95        H=0
  96        CALL MTLPRT(NAME,'D104.2','TOO HIGH ACCURACY REQUIRED')
  97        GO TO 9
  98     5  H=H+S16
  99        IF(BB .NE. B0) GO TO 1
 100       END IF
 101 #if defined(CERNLIB_DOUBLE)
 102     9 DCAUCH=H
 103 #endif
 104 #if !defined(CERNLIB_DOUBLE)
 105     9 RCAUCH=H
 106 #endif
 107       RETURN
 108   101 FORMAT('SINGULARITY  S = ',D15.8,' AT END-POINT OF INTERVAL')
 109       END