Mercurial > octave-nkf

annotate libcruft/amos/cbesk.f @ 12312:b10ea6efdc58 release-3-4-x ss-3-3-91

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
version is now 3.3.91
author John W. Eaton <jwe@octave.org>
date Mon, 31 Jan 2011 08:36:58 -0500
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 CBESK(Z, FNU, KODE, N, CY, NZ, IERR)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
2 C***BEGIN PROLOGUE CBESK
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
3 C***DATE WRITTEN 830501 (YYMMDD)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
4 C***REVISION DATE 890801 (YYMMDD)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
5 C***CATEGORY NO. B5K
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
6 C***KEYWORDS K-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
7 C MODIFIED BESSEL FUNCTION OF THE SECOND KIND,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
8 C BESSEL FUNCTION OF THE THIRD KIND
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
9 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
10 C***PURPOSE TO COMPUTE K-BESSEL FUNCTIONS OF COMPLEX ARGUMENT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
11 C***DESCRIPTION
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
12 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
13 C ON KODE=1, CBESK COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
14 C BESSEL FUNCTIONS CY(J)=K(FNU+J-1,Z) FOR REAL, NONNEGATIVE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
15 C ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z.NE.CMPLX(0.0,0.0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
16 C IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESK
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
17 C RETURNS THE SCALED K FUNCTIONS,
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 C CY(J)=EXP(Z)*K(FNU+J-1,Z) , J=1,...,N,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
20 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
21 C WHICH REMOVE THE EXPONENTIAL BEHAVIOR IN BOTH THE LEFT AND
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
22 C RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
23 C NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
24 C FUNCTIONS (REF. 1).
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
25 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
26 C INPUT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
27 C Z - Z=CMPLX(X,Y),Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
28 C FNU - ORDER OF INITIAL K FUNCTION, FNU.GE.0.0E0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
29 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
30 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
31 C KODE= 1 RETURNS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
32 C CY(I)=K(FNU+I-1,Z), I=1,...,N
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
33 C = 2 RETURNS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
34 C CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
35 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
36 C OUTPUT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
37 C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
38 C VALUES FOR THE SEQUENCE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
39 C CY(I)=K(FNU+I-1,Z), I=1,...,N OR
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
40 C CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
41 C DEPENDING ON KODE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
42 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
43 C NZ= 0 , NORMAL RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
44 C NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
45 C DUE TO UNDERFLOW, CY(I)=CMPLX(0.0,0.0),
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
46 C I=1,...,N WHEN X.GE.0.0. WHEN X.LT.0.0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
47 C NZ STATES ONLY THE NUMBER OF UNDERFLOWS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
48 C IN THE SEQUENCE.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
49 C IERR - ERROR FLAG
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
50 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
51 C IERR=1, INPUT ERROR - NO COMPUTATION
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
52 C IERR=2, OVERFLOW - NO COMPUTATION, FNU+N-1 IS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
53 C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
54 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
55 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
56 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
57 C ACCURACY
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
58 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
59 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
60 C CANCE BY ARGUMENT REDUCTION
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
61 C IERR=5, ERROR - NO COMPUTATION,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
62 C ALGORITHM TERMINATION CONDITION NOT MET
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
63 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
64 C***LONG DESCRIPTION
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
65 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
66 C EQUATIONS OF THE REFERENCE ARE IMPLEMENTED FOR SMALL ORDERS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
67 C DNU AND DNU+1.0 IN THE RIGHT HALF PLANE X.GE.0.0. FORWARD
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
68 C RECURRENCE GENERATES HIGHER ORDERS. K IS CONTINUED TO THE LEFT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
69 C HALF PLANE BY THE RELATION
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
70 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
71 C K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
72 C MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
73 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
74 C WHERE I(FNU,Z) IS THE I BESSEL FUNCTION.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
75 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
76 C FOR LARGE ORDERS, FNU.GT.FNUL, THE K FUNCTION IS COMPUTED
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
77 C BY MEANS OF ITS UNIFORM ASYMPTOTIC EXPANSIONS.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
78 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
79 C FOR NEGATIVE ORDERS, THE FORMULA
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
80 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
81 C K(-FNU,Z) = K(FNU,Z)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
82 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
83 C CAN BE USED.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
84 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
85 C CBESK ASSUMES THAT A SIGNIFICANT DIGIT SINH(X) FUNCTION IS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
86 C AVAILABLE.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
87 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
88 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
89 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
90 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
91 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
92 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
93 C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
94 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
95 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
96 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
97 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
98 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
99 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
100 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
101 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
102 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
103 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
104 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
105 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
106 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
107 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
108 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
109 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
110 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
111 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
112 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
113 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
114 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
115 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
116 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
117 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
118 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
119 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
120 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
121 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
122 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
123 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
124 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
125 C OR -PI/2+P.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
126 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
127 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
128 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
129 C COMMERCE, 1955.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
130 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
131 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
132 C BY D. E. AMOS, SAND83-0083, MAY, 1983.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
133 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
134 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
135 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
136 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
137 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
138 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
139 C 1018, MAY, 1985
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
140 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
141 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
142 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
143 C MATH. SOFTWARE, 1986
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
144 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
145 C***ROUTINES CALLED CACON,CBKNU,CBUNK,CUOIK,I1MACH,R1MACH
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
146 C***END PROLOGUE CBESK
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
147 C
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
148 COMPLEX CY, Z
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
149 REAL AA, ALIM, ALN, ARG, AZ, DIG, ELIM, FN, FNU, FNUL, RL, R1M5,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
150 * TOL, UFL, XX, YY, R1MACH, BB
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
151 INTEGER IERR, K, KODE, K1, K2, MR, N, NN, NUF, NW, NZ, I1MACH
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
152 DIMENSION CY(N)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
153 C***FIRST EXECUTABLE STATEMENT CBESK
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
154 IERR = 0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
155 NZ=0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
156 XX = REAL(Z)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
157 YY = AIMAG(Z)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
158 IF (YY.EQ.0.0E0 .AND. XX.EQ.0.0E0) IERR=1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
159 IF (FNU.LT.0.0E0) IERR=1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
160 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
161 IF (N.LT.1) IERR=1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
162 IF (IERR.NE.0) RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
163 NN = N
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
164 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
165 C SET PARAMETERS RELATED TO MACHINE CONSTANTS.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
166 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
167 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
168 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
169 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
170 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
171 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
172 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
173 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
174 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
175 TOL = AMAX1(R1MACH(4),1.0E-18)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
176 K1 = I1MACH(12)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
177 K2 = I1MACH(13)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
178 R1M5 = R1MACH(5)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
179 K = MIN0(IABS(K1),IABS(K2))
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
180 ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
181 K1 = I1MACH(11) - 1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
182 AA = R1M5*FLOAT(K1)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
183 DIG = AMIN1(AA,18.0E0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
184 AA = AA*2.303E0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
185 ALIM = ELIM + AMAX1(-AA,-41.45E0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
186 FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
187 RL = 1.2E0*DIG + 3.0E0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
188 AZ = CABS(Z)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
189 FN = FNU + FLOAT(NN-1)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
190 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
191 C TEST FOR RANGE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
192 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
193 AA = 0.5E0/TOL
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
194 BB=FLOAT(I1MACH(9))*0.5E0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
195 AA=AMIN1(AA,BB)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
196 IF(AZ.GT.AA) GO TO 210
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
197 IF(FN.GT.AA) GO TO 210
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
198 AA=SQRT(AA)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
199 IF(AZ.GT.AA) IERR=3
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
200 IF(FN.GT.AA) IERR=3
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
201 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
202 C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
203 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
204 C UFL = EXP(-ELIM)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
205 UFL = R1MACH(1)*1.0E+3
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
206 IF (AZ.LT.UFL) GO TO 180
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
207 IF (FNU.GT.FNUL) GO TO 80
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
208 IF (FN.LE.1.0E0) GO TO 60
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
209 IF (FN.GT.2.0E0) GO TO 50
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
210 IF (AZ.GT.TOL) GO TO 60
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
211 ARG = 0.5E0*AZ
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
212 ALN = -FN*ALOG(ARG)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
213 IF (ALN.GT.ELIM) GO TO 180
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
214 GO TO 60
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
215 50 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
216 CALL CUOIK(Z, FNU, KODE, 2, NN, CY, NUF, TOL, ELIM, ALIM)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
217 IF (NUF.LT.0) GO TO 180
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
218 NZ = NZ + NUF
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
219 NN = NN - NUF
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
220 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
221 C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
222 C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
223 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
224 IF (NN.EQ.0) GO TO 100
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
225 60 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
226 IF (XX.LT.0.0E0) GO TO 70
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
227 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
228 C RIGHT HALF PLANE COMPUTATION, REAL(Z).GE.0.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
229 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
230 CALL CBKNU(Z, FNU, KODE, NN, CY, NW, TOL, ELIM, ALIM)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
231 IF (NW.LT.0) GO TO 200
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
232 NZ=NW
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
233 RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
234 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
235 C LEFT HALF PLANE COMPUTATION
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
236 C PI/2.LT.ARG(Z).LE.PI AND -PI.LT.ARG(Z).LT.-PI/2.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
237 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
238 70 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
239 IF (NZ.NE.0) GO TO 180
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
240 MR = 1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
241 IF (YY.LT.0.0E0) MR = -1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
242 CALL CACON(Z, FNU, KODE, MR, NN, CY, NW, RL, FNUL, TOL, ELIM,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
243 * ALIM)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
244 IF (NW.LT.0) GO TO 200
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
245 NZ=NW
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
246 RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
247 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
248 C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
249 C-----------------------------------------------------------------------
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
250 80 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
251 MR = 0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
252 IF (XX.GE.0.0E0) GO TO 90
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
253 MR = 1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
254 IF (YY.LT.0.0E0) MR = -1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
255 90 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
256 CALL CBUNK(Z, FNU, KODE, MR, NN, CY, NW, TOL, ELIM, ALIM)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
257 IF (NW.LT.0) GO TO 200
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
258 NZ = NZ + NW
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
259 RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
260 100 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
261 IF (XX.LT.0.0E0) GO TO 180
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
262 RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
263 180 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
264 NZ = 0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
265 IERR=2
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
266 RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
267 200 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
268 IF(NW.EQ.(-1)) GO TO 180
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
269 NZ=0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
270 IERR=5
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
271 RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
272 210 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
273 NZ=0
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
274 IERR=4
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
275 RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
276 END