MINICERN/mathlib/gen/c/c344cod2.inc

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:09  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 *
  10 * c344cod2.inc
  11 *
  12       D=TC(8)*ND+TC(9)
  13       Q=0
  14       R0=AZ
  15       IF(NMAX .GT. 0) THEN
  16        V=0.5*D/NMAX
  17        IF(V .LE. 10) THEN
  18         T=TC(1)
  19         DO 6 J = 2,6
  20     6   T=V*T+TC(J)
  21        ELSE
  22         U=LOG(V)-TC(7)
  23         T=V/(U*(1+(TC(7)-LOG(U))/(1+U)))
  24        END IF
  25        Q=NMAX*T
  26       END IF
  27       T=1
  28       IF(YA .LT. D) THEN
  29        V=TC(10)*(D-YA)/R0
  30        IF(V .LE. 10) THEN
  31         T=TC(1)
  32         DO 7 J = 2,6
  33     7   T=V*T+TC(J)
  34        ELSE
  35         U=LOG(V)-TC(7)
  36         T=V/(U*(1+(TC(7)-LOG(U))/(1+U)))
  37        END IF
  38       END IF
  39       NU=1+MAX(Q,TC(11)*R0*T)
  40
  41       A2=A+A
  42       MU=-1
  43     2 MU=MU+1
  44       AL=1
  45       DO 3 N = 1,NU
  46       XN=N
  47     3 AL=AL*(XN+A2)/(XN+1)
  48       C=(-I)**NU
  49       R=0
  50       S=0
  51       DO 4 N = NU,1,-1
  52       XN=N
  53       A2N=2*(A+XN)
  54       R=ZA/(A2N-ZA*R)
  55       IF(N .LE. NMAX) RR(N-1)=R
  56       AL=AL*(XN+1)/(XN+A2)
  57       S=R*(A2N*AL*C+S)
  58     4 C=I*C
  59       CB(0)=F/(1+S)
  60       DO 10 N = 0,NMAX-1
  61    10 CB(N+1)=RR(N)*CB(N)
  62       IF(Y .LT. 0) THEN
  63        DO 11 N = 0,NMAX
  64    11  CB(N)=GCONJG(CB(N))
  65       END IF
  66       DO 12 N = 0,NMAX
  67       IF(ABS(CB(N)-CBA(N)) .GT. EPS*ABS(CB(N))) THEN
  68        DO 13 M = 0,NMAX
  69    13  CBA(M)=CB(M)
  70        NU=NU+5
  71        IF(MU .LE. 50) GO TO 2
  72        WRITE(ERRTXT,104) Z,A
  73        CALL MTLPRT(NAME,'C344.4',ERRTXT)
  74        RETURN
  75       END IF
  76    12 CONTINUE
  77       RETURN
  78   101 FORMAT('ILLEGAL ARGUMENT Z = ',1P,2D15.8)
  79   102 FORMAT('ILLEGAL ORDER A = ',1P,D15.8)
  80   103 FORMAT('ILLEGAL NMAX = ',I5)
  81   104 FORMAT('NO CONVERGENCE, Z = ',1P,2D13.6,' A = ',D13.6,
  82      1       ' TRY SMALLER ND')
  83       END