← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/special/amos/zbuni.f

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/amos/zbuni.f
 
1
      SUBROUTINE ZBUNI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NUI, NLAST,
 
2
     * FNUL, TOL, ELIM, ALIM)
 
3
C***BEGIN PROLOGUE  ZBUNI
 
4
C***REFER TO  ZBESI,ZBESK
 
5
C
 
6
C     ZBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z).GT.
 
7
C     FNUL AND FNU+N-1.LT.FNUL. THE ORDER IS INCREASED FROM
 
8
C     FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING
 
9
C     ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z)
 
10
C     ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2
 
11
C
 
12
C***ROUTINES CALLED  ZUNI1,ZUNI2,AZABS,D1MACH
 
13
C***END PROLOGUE  ZBUNI
 
14
C     COMPLEX CSCL,CSCR,CY,RZ,ST,S1,S2,Y,Z
 
15
      DOUBLE PRECISION ALIM, AX, AY, CSCLR, CSCRR, CYI, CYR, DFNU,
 
16
     * ELIM, FNU, FNUI, FNUL, GNU, RAZ, RZI, RZR, STI, STR, S1I, S1R,
 
17
     * S2I, S2R, TOL, YI, YR, ZI, ZR, AZABS, ASCLE, BRY, C1R, C1I, C1M,
 
18
     * D1MACH
 
19
      INTEGER I, IFLAG, IFORM, K, KODE, N, NL, NLAST, NUI, NW, NZ
 
20
      DIMENSION YR(N), YI(N), CYR(2), CYI(2), BRY(3)
 
21
      NZ = 0
 
22
      AX = DABS(ZR)*1.7321D0
 
23
      AY = DABS(ZI)
 
24
      IFORM = 1
 
25
      IF (AY.GT.AX) IFORM = 2
 
26
      IF (NUI.EQ.0) GO TO 60
 
27
      FNUI = DBLE(FLOAT(NUI))
 
28
      DFNU = FNU + DBLE(FLOAT(N-1))
 
29
      GNU = DFNU + FNUI
 
30
      IF (IFORM.EQ.2) GO TO 10
 
31
C-----------------------------------------------------------------------
 
32
C     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
 
33
C     -PI/3.LE.ARG(Z).LE.PI/3
 
34
C-----------------------------------------------------------------------
 
35
      CALL ZUNI1(ZR, ZI, GNU, KODE, 2, CYR, CYI, NW, NLAST, FNUL, TOL,
 
36
     * ELIM, ALIM)
 
37
      GO TO 20
 
38
   10 CONTINUE
 
39
C-----------------------------------------------------------------------
 
40
C     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
 
41
C     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
 
42
C     AND HPI=PI/2
 
43
C-----------------------------------------------------------------------
 
44
      CALL ZUNI2(ZR, ZI, GNU, KODE, 2, CYR, CYI, NW, NLAST, FNUL, TOL,
 
45
     * ELIM, ALIM)
 
46
   20 CONTINUE
 
47
      IF (NW.LT.0) GO TO 50
 
48
      IF (NW.NE.0) GO TO 90
 
49
      STR = AZABS(CYR(1),CYI(1))
 
50
C----------------------------------------------------------------------
 
51
C     SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
 
52
C----------------------------------------------------------------------
 
53
      BRY(1)=1.0D+3*D1MACH(1)/TOL
 
54
      BRY(2) = 1.0D0/BRY(1)
 
55
      BRY(3) = BRY(2)
 
56
      IFLAG = 2
 
57
      ASCLE = BRY(2)
 
58
      CSCLR = 1.0D0
 
59
      IF (STR.GT.BRY(1)) GO TO 21
 
60
      IFLAG = 1
 
61
      ASCLE = BRY(1)
 
62
      CSCLR = 1.0D0/TOL
 
63
      GO TO 25
 
64
   21 CONTINUE
 
65
      IF (STR.LT.BRY(2)) GO TO 25
 
66
      IFLAG = 3
 
67
      ASCLE=BRY(3)
 
68
      CSCLR = TOL
 
69
   25 CONTINUE
 
70
      CSCRR = 1.0D0/CSCLR
 
71
      S1R = CYR(2)*CSCLR
 
72
      S1I = CYI(2)*CSCLR
 
73
      S2R = CYR(1)*CSCLR
 
74
      S2I = CYI(1)*CSCLR
 
75
      RAZ = 1.0D0/AZABS(ZR,ZI)
 
76
      STR = ZR*RAZ
 
77
      STI = -ZI*RAZ
 
78
      RZR = (STR+STR)*RAZ
 
79
      RZI = (STI+STI)*RAZ
 
80
      DO 30 I=1,NUI
 
81
        STR = S2R
 
82
        STI = S2I
 
83
        S2R = (DFNU+FNUI)*(RZR*STR-RZI*STI) + S1R
 
84
        S2I = (DFNU+FNUI)*(RZR*STI+RZI*STR) + S1I
 
85
        S1R = STR
 
86
        S1I = STI
 
87
        FNUI = FNUI - 1.0D0
 
88
        IF (IFLAG.GE.3) GO TO 30
 
89
        STR = S2R*CSCRR
 
90
        STI = S2I*CSCRR
 
91
        C1R = DABS(STR)
 
92
        C1I = DABS(STI)
 
93
        C1M = DMAX1(C1R,C1I)
 
94
        IF (C1M.LE.ASCLE) GO TO 30
 
95
        IFLAG = IFLAG+1
 
96
        ASCLE = BRY(IFLAG)
 
97
        S1R = S1R*CSCRR
 
98
        S1I = S1I*CSCRR
 
99
        S2R = STR
 
100
        S2I = STI
 
101
        CSCLR = CSCLR*TOL
 
102
        CSCRR = 1.0D0/CSCLR
 
103
        S1R = S1R*CSCLR
 
104
        S1I = S1I*CSCLR
 
105
        S2R = S2R*CSCLR
 
106
        S2I = S2I*CSCLR
 
107
   30 CONTINUE
 
108
      YR(N) = S2R*CSCRR
 
109
      YI(N) = S2I*CSCRR
 
110
      IF (N.EQ.1) RETURN
 
111
      NL = N - 1
 
112
      FNUI = DBLE(FLOAT(NL))
 
113
      K = NL
 
114
      DO 40 I=1,NL
 
115
        STR = S2R
 
116
        STI = S2I
 
117
        S2R = (FNU+FNUI)*(RZR*STR-RZI*STI) + S1R
 
118
        S2I = (FNU+FNUI)*(RZR*STI+RZI*STR) + S1I
 
119
        S1R = STR
 
120
        S1I = STI
 
121
        STR = S2R*CSCRR
 
122
        STI = S2I*CSCRR
 
123
        YR(K) = STR
 
124
        YI(K) = STI
 
125
        FNUI = FNUI - 1.0D0
 
126
        K = K - 1
 
127
        IF (IFLAG.GE.3) GO TO 40
 
128
        C1R = DABS(STR)
 
129
        C1I = DABS(STI)
 
130
        C1M = DMAX1(C1R,C1I)
 
131
        IF (C1M.LE.ASCLE) GO TO 40
 
132
        IFLAG = IFLAG+1
 
133
        ASCLE = BRY(IFLAG)
 
134
        S1R = S1R*CSCRR
 
135
        S1I = S1I*CSCRR
 
136
        S2R = STR
 
137
        S2I = STI
 
138
        CSCLR = CSCLR*TOL
 
139
        CSCRR = 1.0D0/CSCLR
 
140
        S1R = S1R*CSCLR
 
141
        S1I = S1I*CSCLR
 
142
        S2R = S2R*CSCLR
 
143
        S2I = S2I*CSCLR
 
144
   40 CONTINUE
 
145
      RETURN
 
146
   50 CONTINUE
 
147
      NZ = -1
 
148
      IF(NW.EQ.(-2)) NZ=-2
 
149
      RETURN
 
150
   60 CONTINUE
 
151
      IF (IFORM.EQ.2) GO TO 70
 
152
C-----------------------------------------------------------------------
 
153
C     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
 
154
C     -PI/3.LE.ARG(Z).LE.PI/3
 
155
C-----------------------------------------------------------------------
 
156
      CALL ZUNI1(ZR, ZI, FNU, KODE, N, YR, YI, NW, NLAST, FNUL, TOL,
 
157
     * ELIM, ALIM)
 
158
      GO TO 80
 
159
   70 CONTINUE
 
160
C-----------------------------------------------------------------------
 
161
C     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
 
162
C     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
 
163
C     AND HPI=PI/2
 
164
C-----------------------------------------------------------------------
 
165
      CALL ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, NW, NLAST, FNUL, TOL,
 
166
     * ELIM, ALIM)
 
167
   80 CONTINUE
 
168
      IF (NW.LT.0) GO TO 50
 
169
      NZ = NW
 
170
      RETURN
 
171
   90 CONTINUE
 
172
      NLAST = N
 
173
      RETURN
 
174
      END