Mercurial > hg > octave-lyh

annotate libcruft/amos/cwrsk.f @ 11027:a2bdd9180a3c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
__print_parse_opts__.m: Properly crop eps image and fix tests.
author Ben Abbott <bpabbott@mac.com>
date Mon, 27 Sep 2010 19:38:16 -0400
parents 82be108cc558
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
1 SUBROUTINE CWRSK(ZR, FNU, KODE, N, Y, NZ, CW, TOL, ELIM, ALIM)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
2 C***BEGIN PROLOGUE CWRSK
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
3 C***REFER TO CBESI,CBESK
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
4 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
5 C CWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
6 C NORMALIZING THE I FUNCTION RATIOS FROM CRATI BY THE WRONSKIAN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
7 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
8 C***ROUTINES CALLED CBKNU,CRATI,R1MACH
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
9 C***END PROLOGUE CWRSK
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
10 COMPLEX CINU, CSCL, CT, CW, C1, C2, RCT, ST, Y, ZR
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
11 REAL ACT, ACW, ALIM, ASCLE, ELIM, FNU, S1, S2, TOL, YY
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
12 INTEGER I, KODE, N, NW, NZ
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
13 DIMENSION Y(N), CW(2)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
14 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
15 C I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
16 C Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
17 C WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
18 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
19 NZ = 0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
20 CALL CBKNU(ZR, FNU, KODE, 2, CW, NW, TOL, ELIM, ALIM)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
21 IF (NW.NE.0) GO TO 50
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
22 CALL CRATI(ZR, FNU, N, Y, TOL)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
23 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
24 C RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
25 C R(FNU+J-1,Z)=Y(J), J=1,...,N
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
26 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
27 CINU = CMPLX(1.0E0,0.0E0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
28 IF (KODE.EQ.1) GO TO 10
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
29 YY = AIMAG(ZR)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
30 S1 = COS(YY)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
31 S2 = SIN(YY)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
32 CINU = CMPLX(S1,S2)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
33 10 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
34 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
35 C ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
36 C THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
37 C SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
38 C THE RESULT IS ON SCALE.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
39 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
40 ACW = CABS(CW(2))
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
41 ASCLE = 1.0E+3*R1MACH(1)/TOL
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
42 CSCL = CMPLX(1.0E0,0.0E0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
43 IF (ACW.GT.ASCLE) GO TO 20
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
44 CSCL = CMPLX(1.0E0/TOL,0.0E0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
45 GO TO 30
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
46 20 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
47 ASCLE = 1.0E0/ASCLE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
48 IF (ACW.LT.ASCLE) GO TO 30
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
49 CSCL = CMPLX(TOL,0.0E0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
50 30 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
51 C1 = CW(1)*CSCL
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
52 C2 = CW(2)*CSCL
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
53 ST = Y(1)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
54 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
55 C CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0E0/CABS(CT) PREVENTS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
56 C UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
57 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
58 CT = ZR*(C2+ST*C1)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
59 ACT = CABS(CT)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
60 RCT = CMPLX(1.0E0/ACT,0.0E0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
61 CT = CONJG(CT)*RCT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
62 CINU = CINU*RCT*CT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
63 Y(1) = CINU*CSCL
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
64 IF (N.EQ.1) RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
65 DO 40 I=2,N
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
66 CINU = ST*CINU
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
67 ST = Y(I)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
68 Y(I) = CINU*CSCL
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
69 40 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
70 RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
71 50 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
72 NZ = -1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
73 IF(NW.EQ.(-2)) NZ=-2
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
74 RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
75 END