github.com/gopherd/gonum@v0.0.4/mathext/internal/amos/amoslib/zbesy.f (about) 1 SUBROUTINE ZBESY(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, CWRKR, CWRKI, 2 * IERR) 3 C***BEGIN PROLOGUE ZBESY 4 C***DATE WRITTEN 830501 (YYMMDD) 5 C***REVISION DATE 890801 (YYMMDD) 6 C***CATEGORY NO. B5K 7 C***KEYWORDS Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT, 8 C BESSEL FUNCTION OF SECOND KIND 9 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 10 C***PURPOSE TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT 11 C***DESCRIPTION 12 C 13 C ***A DOUBLE PRECISION ROUTINE*** 14 C 15 C ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 16 C BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE 17 C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE 18 C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESY RETURNS THE SCALED 19 C FUNCTIONS 20 C 21 C CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z) 22 C 23 C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND 24 C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION 25 C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS 26 C (REF. 1). 27 C 28 C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION 29 C ZR,ZI - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0), 30 C -PI.LT.ARG(Z).LE.PI 31 C FNU - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0D0 32 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 33 C KODE= 1 RETURNS 34 C CY(I)=Y(FNU+I-1,Z), I=1,...,N 35 C = 2 RETURNS 36 C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N 37 C WHERE Y=AIMAG(Z) 38 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 39 C CWRKR, - DOUBLE PRECISION WORK VECTORS OF DIMENSION AT 40 C CWRKI AT LEAST N 41 C 42 C OUTPUT CYR,CYI ARE DOUBLE PRECISION 43 C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS 44 C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE 45 C CY(I)=Y(FNU+I-1,Z) OR 46 C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N 47 C DEPENDING ON KODE. 48 C NZ - NZ=0 , A NORMAL RETURN 49 C NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO 50 C UNDERFLOW (GENERALLY ON KODE=2) 51 C IERR - ERROR FLAG 52 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 53 C IERR=1, INPUT ERROR - NO COMPUTATION 54 C IERR=2, OVERFLOW - NO COMPUTATION, FNU IS 55 C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH 56 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE 57 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT 58 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE 59 C ACCURACY 60 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- 61 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- 62 C CANCE BY ARGUMENT REDUCTION 63 C IERR=5, ERROR - NO COMPUTATION, 64 C ALGORITHM TERMINATION CONDITION NOT MET 65 C 66 C***LONG DESCRIPTION 67 C 68 C THE COMPUTATION IS CARRIED OUT BY THE FORMULA 69 C 70 C Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I 71 C 72 C WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z) 73 C AND H(2,FNU,Z) ARE CALCULATED IN CBESH. 74 C 75 C FOR NEGATIVE ORDERS,THE FORMULA 76 C 77 C Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU) 78 C 79 C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD 80 C INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE 81 C POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)* 82 C SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS 83 C NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A 84 C LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM 85 C CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, 86 C WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF 87 C ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z). 88 C 89 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 90 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS 91 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. 92 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN 93 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG 94 C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS 95 C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. 96 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS 97 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS 98 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE 99 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS 100 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 101 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION 102 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION 103 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN 104 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT 105 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS 106 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. 107 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. 108 C 109 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 110 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 111 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 112 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 113 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 114 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 115 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 116 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 117 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 118 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 119 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 120 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 121 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 122 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 123 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 124 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 125 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 126 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 127 C OR -PI/2+P. 128 C 129 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 130 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 131 C COMMERCE, 1955. 132 C 133 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 134 C BY D. E. AMOS, SAND83-0083, MAY, 1983. 135 C 136 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 137 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 138 C 139 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 140 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 141 C 1018, MAY, 1985 142 C 143 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 144 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 145 C MATH. SOFTWARE, 1986 146 C 147 C***ROUTINES CALLED ZBESH,I1MACH,D1MACH 148 C***END PROLOGUE ZBESY 149 C 150 C COMPLEX CWRK,CY,C1,C2,EX,HCI,Z,ZU,ZV 151 DOUBLE PRECISION CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2R, 152 * ELIM, EXI, EXR, EY, FNU, HCII, STI, STR, TAY, ZI, ZR, DEXP, 153 * D1MACH, ASCLE, RTOL, ATOL, AA, BB, TOL 154 INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH 155 DIMENSION CYR(N), CYI(N), CWRKR(N), CWRKI(N) 156 C***FIRST EXECUTABLE STATEMENT ZBESY 157 IERR = 0 158 NZ=0 159 IF (ZR.EQ.0.0D0 .AND. ZI.EQ.0.0D0) IERR=1 160 IF (FNU.LT.0.0D0) IERR=1 161 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 162 IF (N.LT.1) IERR=1 163 IF (IERR.NE.0) RETURN 164 HCII = 0.5D0 165 CALL ZBESH(ZR, ZI, FNU, KODE, 1, N, CYR, CYI, NZ1, IERR) 166 IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170 167 CALL ZBESH(ZR, ZI, FNU, KODE, 2, N, CWRKR, CWRKI, NZ2, IERR) 168 IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170 169 NZ = MIN0(NZ1,NZ2) 170 IF (KODE.EQ.2) GO TO 60 171 DO 50 I=1,N 172 STR = CWRKR(I) - CYR(I) 173 STI = CWRKI(I) - CYI(I) 174 CYR(I) = -STI*HCII 175 CYI(I) = STR*HCII 176 50 CONTINUE 177 RETURN 178 60 CONTINUE 179 TOL = DMAX1(D1MACH(4),1.0D-18) 180 K1 = I1MACH(15) 181 K2 = I1MACH(16) 182 K = MIN0(IABS(K1),IABS(K2)) 183 R1M5 = D1MACH(5) 184 C----------------------------------------------------------------------- 185 C ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT 186 C----------------------------------------------------------------------- 187 ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) 188 EXR = DCOS(ZR) 189 EXI = DSIN(ZR) 190 EY = 0.0D0 191 TAY = DABS(ZI+ZI) 192 IF (TAY.LT.ELIM) EY = DEXP(-TAY) 193 IF (ZI.LT.0.0D0) GO TO 90 194 C1R = EXR*EY 195 C1I = EXI*EY 196 C2R = EXR 197 C2I = -EXI 198 70 CONTINUE 199 NZ = 0 200 RTOL = 1.0D0/TOL 201 ASCLE = D1MACH(1)*RTOL*1.0D+3 202 DO 80 I=1,N 203 C STR = C1R*CYR(I) - C1I*CYI(I) 204 C STI = C1R*CYI(I) + C1I*CYR(I) 205 C STR = -STR + C2R*CWRKR(I) - C2I*CWRKI(I) 206 C STI = -STI + C2R*CWRKI(I) + C2I*CWRKR(I) 207 C CYR(I) = -STI*HCII 208 C CYI(I) = STR*HCII 209 AA = CWRKR(I) 210 BB = CWRKI(I) 211 ATOL = 1.0D0 212 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 75 213 AA = AA*RTOL 214 BB = BB*RTOL 215 ATOL = TOL 216 75 CONTINUE 217 STR = (AA*C2R - BB*C2I)*ATOL 218 STI = (AA*C2I + BB*C2R)*ATOL 219 AA = CYR(I) 220 BB = CYI(I) 221 ATOL = 1.0D0 222 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 85 223 AA = AA*RTOL 224 BB = BB*RTOL 225 ATOL = TOL 226 85 CONTINUE 227 STR = STR - (AA*C1R - BB*C1I)*ATOL 228 STI = STI - (AA*C1I + BB*C1R)*ATOL 229 CYR(I) = -STI*HCII 230 CYI(I) = STR*HCII 231 IF (STR.EQ.0.0D0 .AND. STI.EQ.0.0D0 .AND. EY.EQ.0.0D0) NZ = NZ 232 * + 1 233 80 CONTINUE 234 RETURN 235 90 CONTINUE 236 C1R = EXR 237 C1I = EXI 238 C2R = EXR*EY 239 C2I = -EXI*EY 240 GO TO 70 241 170 CONTINUE 242 NZ = 0 243 RETURN 244 END