MINICERN/mathlib/gen/c/nzerfz64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:53  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       FUNCTION NZERFZ(F,ZP,N)
  11 C
  12 #include "gen/impc64.inc"
  13 C
  14       CHARACTER*(*) NAME
  15       PARAMETER(NAME='NZERFZ')
  16 #include "gen/def64.inc"
  17      + W,X,R,EPS,CST,CONST,R1,HF
  18       DIMENSION ZP(*),W(12),X(12)
  19
  20       PARAMETER (PI = 3.14159 26535 89793 24D0)
  21 #if defined(CERNLIB_CRAY)
  22       PARAMETER (I = (0E0,1E0))
  23 #endif
  24 #if !defined(CERNLIB_CRAY)
  25       PARAMETER (I = (0D0,1D0))
  26 #endif
  27       PARAMETER ( CPI2 = 2*PI*I)
  28       PARAMETER (R1 = 1, HF = R1/2, DZ = (1+I)*1D-8)
  29       PARAMETER (EPS = 1D-4, CST = 0.005D0, NFMAX = 200000)
  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       H=0
  45       DO 10 J = 1,N
  46       NF=0
  47       A=ZP(J)
  48       IF(J .LT. N) THEN
  49        B=ZP(J+1)
  50       ELSE
  51        B=ZP(1)
  52       ENDIF
  53       IF(B .EQ. A) GO TO 10
  54       CONST=CST/ABS(B-A)
  55       BB=A
  56     1 AA=BB
  57       BB=B
  58     2 C1=HF*(BB+AA)
  59       C2=HF*(BB-AA)
  60       S8=0
  61       DO 3 K = 1,4
  62       U=C2*X(K)
  63       FPP=F(C1+U+DZ)
  64       FPM=F(C1+U-DZ)
  65       FMP=F(C1-U+DZ)
  66       FMM=F(C1-U-DZ)
  67     3 S8=S8+W(K)*(((FPP-FPM)/DZ)/(FPP+FPM)+((FMP-FMM)/DZ)/(FMP+FMM))
  68       S16=0
  69       DO 4 K = 5,12
  70       U=C2*X(K)
  71       FPP=F(C1+U+DZ)
  72       FPM=F(C1+U-DZ)
  73       FMP=F(C1-U+DZ)
  74       FMM=F(C1-U-DZ)
  75     4 S16=S16+W(K)*(((FPP-FPM)/DZ)/(FPP+FPM)+((FMP-FMM)/DZ)/(FMP+FMM))
  76       S16=C2*S16
  77       NF=NF+48
  78       IF(ABS(S16-C2*S8) .LE. EPS*(1+ABS(S16))) THEN
  79        H=H+S16
  80        IF(BB .NE. B) GO TO 1
  81       ELSE
  82        BB=C1
  83        IF(1+CONST*ABS(C2) .NE. 1 .AND. NF .LE. NFMAX) GO TO 2
  84        R=0
  85        CALL MTLPRT(NAME,'C210.1','PROBLEMS WITH INTEGRATION,'//
  86      +             'POLYGON TOO NEAR TO A ZERO OR SINGULARITY ?')
  87        GO TO 99
  88       ENDIF
  89    10 CONTINUE
  90       R=H/CPI2
  91    99 NZERFZ=NINT(ABS(R))
  92       RETURN
  93       END