MINICERN/mathlib/gen/g/vviset.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:48  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE VVISET(RKA,BE2,MODE,XL,XU)
  11
  12       COMMON /G116C1/ H(7),T0,T1,T,OMEGA,A(155),B(155),X0
  13       EXTERNAL G116F1,G116F2
  14       CHARACTER NAME*(*)
  15       CHARACTER*80 ERRTXT
  16       PARAMETER (NAME = 'VVISET')
  17       DIMENSION XP(8),XQ(6)
  18
  19 C     H1 = 5 ln 10 - ln(2/PI**2), h2 = 3 ln 10, h0 = Ln 0.0005...
  20
  21       PARAMETER (PI = 3.14159 265, EU = 0.57721 566)
  22       PARAMETER (PI2 = 2*PI, RPI = 1/PI)
  23       PARAMETER (H1 = 5*2.30258 509-1.59631 259, H2 = 3*2.30258 509)
  24       PARAMETER (H0 =-7.6, PIH = PI/2, EPS = 1E-5)
  25
  26       DATA  XP /9.29,  2.47, 0.89, 0.36, 0.15, 0.07, 0.03, 0.02/
  27       DATA  XQ /0.012, 0.03, 0.08, 0.26, 0.87, 3.83/
  28
  29       IF(RKA .LT. 0.01 .OR. RKA .GT. 10) THEN
  30        WRITE(ERRTXT,101) RKA
  31        CALL MTLPRT(NAME,'G116.1',ERRTXT)
  32       ELSEIF(BE2 .LT. 0 .OR. BE2 .GT. 1) THEN
  33        WRITE(ERRTXT,102) BE2
  34        CALL MTLPRT(NAME,'G116.2',ERRTXT)
  35       ELSE
  36        H(5)=1-BE2*(1-EU)-H0/RKA
  37        H(6)=BE2
  38        H(7)=1-BE2
  39        H4=H0/RKA-(1+BE2*EU)
  40        H5=LOG(RKA)
  41        H6=1/RKA
  42        T0=(H4-H(5)*H5-(H(5)+BE2)*(LOG(H(5))+REXPIN(H(5)))+EXP(-H(5)))/
  43      1     H(5)
  44        DO 1 LP = 1,8
  45        IF(RKA .GE. XP(LP)) GO TO 11
  46     1  CONTINUE
  47        LP=9
  48    11  DO 2 LQ = 1,6
  49        IF(RKA .LE. XQ(LQ)) GO TO 22
  50     2  CONTINUE
  51        LQ=7
  52    22  CALL RZERO(-LP-0.5,LQ-7.5,U,XX,EPS,1000,G116F2)
  53        Q=1/U
  54        T1=H4*Q-H5-(1+BE2*Q)*(LOG(ABS(U))+REXPIN(U))+EXP(-U)*Q
  55
  56        T=T1-T0
  57        OMEGA=PI2/T
  58        H(1)=RKA*(2+BE2*EU)+H1
  59        IF(RKA .GE. 0.07) H(1)=H(1)+H2
  60        H(2)=BE2*RKA
  61        H(3)=H6*OMEGA
  62        H(4)=PIH*OMEGA
  63        CALL RZERO(5.,155.,X0,XX,EPS,1000,G116F1)
  64        N=X0+1
  65
  66        D=RPI*EXP(RKA*(1+BE2*(EU-H5)))
  67        A(N)=0
  68        IF(MODE .EQ. 0) A(N)=RPI*OMEGA
  69        Q=-1
  70        Q2=2
  71        DO 3 K = 1,N-1
  72        L=N-K
  73        X=OMEGA*K
  74        X1=H6*X
  75        C1=LOG(X)-COSINT(X1)
  76        C2=SININT(X1)
  77        C3=SIN(X1)
  78        C4=COS(X1)
  79        XF1=RKA*(BE2*C1-C4)-X*C2
  80        XF2=X*C1+RKA*(C3+BE2*C2)+T0*X
  81        IF(MODE .EQ. 0) THEN
  82         D1=Q*D*OMEGA*EXP(XF1)
  83         A(L)=D1*COS(XF2)
  84         B(L)=-D1*SIN(XF2)
  85        ELSE
  86         D1=Q*D*EXP(XF1)/K
  87         A(L)=D1*SIN(XF2)
  88         B(L)=D1*COS(XF2)
  89         A(N)=A(N)+Q2*A(L)
  90        ENDIF
  91        Q=-Q
  92     3  Q2=-Q2
  93       ENDIF
  94       XL=T0
  95       XU=T1
  96       RETURN
  97   101 FORMAT('KAPPA =   ',E10.3,' - OUT OF RANGE')
  98   102 FORMAT('BETA**2 = ',E10.3,' - OUT OF RANGE')
  99       END