Attribute VB_Name = "Module5" '************************************************************************** '* Procedures and Functions used By programs TZBESJ, TZBESK, TZBESY * '* (Evalute Bessel Functions with complex argument, 1st to 3rd kind) * '* ---------------------------------------------------------------------- * '* Reference: AMOS, DONALD E., SANDIA NATIONAL LABORATORIES, 1983. * '* * '* Visual Basic Release By J-P Moreau, Paris (07/15/2005). * '* (www.jpmoreau.fr) * '************************************************************************** 'MODULE CBESS1 DefDbl A-H, O-Z DefInt I-N 'LIST OF SUBROUTINES: ' ZBINU(ZR, ZI, FNU, KODE, N, CYR(), CYI(), NZ, RL, FNUL, TOL, ELIM, ALIM) ' ZWRSK(ZRR, ZRI, FNU, KODE, N, NZ, CYR(), CYI(), TOL, ELIM, ALIM) Sub ZBINU(ZR, ZI, FNU, KODE, N, CYR(), CYI(), NZ, RL, FNUL, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZBINU '***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZAIRY,ZBIRY ' ZBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE '***ROUTINES CALLED ZABS,ZASYI,ZBUNI,ZMLRI,ZSERI,ZUOIK,ZWRSK '***END PROLOGUE ZBINU } 'Labels: 10,20,30,40,50,60,70,80, 100,110,120,130, 150 ' AZ, DFNU, ZEROI, ZEROR: double ' I, INW, NLAST, NN, NUI, NW: Integer 'Working zones for subroutine ZWRSK Dim CWR(10), CWI(10) ZEROR = 0#: ZEROI = 0# NZ = 0 AZ = ZABS(ZR, ZI) NN = N DFNU = FNU + 1# * (N - 1) If AZ <= 2# Then GoTo 10 If AZ * AZ * 0.25 > DFNU + 1# Then GoTo 20 '----------------------------------------------------------------------- ' POWER SERIES '----------------------------------------------------------------------- 10 ZSERI ZR, ZI, FNU, KODE, NN, CYR(), CYI(), NW, TOL, ELIM, ALIM INW = Abs(NW) NZ = NZ + INW NN = NN - INW If NN = 0 Then GoTo 150 If NW >= 0 Then GoTo 120 DFNU = FNU + 1# * (NN - 1) 20: If AZ < RL Then GoTo 40 If DFNU <= 1# Then GoTo 30 If AZ + AZ < DFNU * DFNU Then GoTo 50 '----------------------------------------------------------------------- ' ASYMPTOTIC EXPANSION FOR LARGE Z '----------------------------------------------------------------------- 30 ZASYI ZR, ZI, FNU, KODE, NN, CYR(), CYI(), NW, RL, TOL, ELIM, ALIM If NW < 0 Then GoTo 130 GoTo 120 40 If DFNU <= 1# Then GoTo 70 '----------------------------------------------------------------------- ' OVERFLOW AND UNDERFLOW TEST ON I SEQUENCE FOR MILLER ALGORITHM '----------------------------------------------------------------------- 50 ZUOIK ZR, ZI, FNU, KODE, 1, NN, CYR(), CYI(), NW, TOL, ELIM, ALIM If NW < 0 Then GoTo 130 NZ = NZ + NW NN = NN - NW If NN = 0 Then GoTo 150 DFNU = FNU + 1# * (NN - 1) If DFNU > FNUL Then GoTo 110 If AZ > FNUL Then GoTo 110 60 If AZ > RL Then GoTo 80 '----------------------------------------------------------------------- ' MILLER ALGORITHM NORMALIZED BY THE SERIES '----------------------------------------------------------------------- 70 ZMLRI ZR, ZI, FNU, KODE, NN, CYR(), CYI(), NW, TOL If NW < 0 Then GoTo 130 GoTo 120 '----------------------------------------------------------------------- ' MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN '----------------------------------------------------------------------- '----------------------------------------------------------------------- ' OVERFLOW TEST ON K FUNCTIONS USED IN WRONSKIAN '----------------------------------------------------------------------- 80 ZUOIK ZR, ZI, FNU, KODE, 2, 2, CYR(), CYI(), NW, TOL, ELIM, ALIM If NW >= 0 Then GoTo 100 NZ = NN For I = 1 To NN CYR(I) = ZEROR CYI(I) = ZEROI Next I GoTo 150 100 If NW > 0 Then GoTo 130 ZWRSK ZR, ZI, FNU, KODE, NN, CYR(), CYI(), NW, CWR(), CWI(), TOL, ELIM, ALIM If NW < 0 Then GoTo 130 GoTo 120 '----------------------------------------------------------------------- ' INCREMENT FNU+NN-1 UP TO FNUL, COMPUTE AND RECUR BACKWARD '----------------------------------------------------------------------- 110 NUI = Int(FNUL - DFNU) + 1 NUI = IMAX(NUI, 0) ZBUNI ZR, ZI, FNU, KODE, NN, CYR(), CYI(), NW, NUI, NLAST, FNUL, TOL, ELIM, ALIM If NW < 0 Then GoTo 130 NZ = NZ + NW If NLAST = 0 Then GoTo 120 NN = NLAST GoTo 60 120 GoTo 150 130 NZ = -1 If NW = -2 Then NZ = -2 150 End Sub 'ZBINU Sub ZWRSK(ZRR, ZRI, FNU, KODE, N, CYR(), CYI(), NZ, CWR(), CWI(), TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZWRSK '***REFER TO ZBESI,ZBESK ' ' ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z) >= 0 BY ' NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN ' '***ROUTINES CALLED D1MACH,ZBKNU,ZRATI,ZABS '***END PROLOGUE ZWRSK ' COMPLEX CINU,CSCL,CT,CY,C1,C2,RCT,ST,Y,ZR 'Labels: 10,20,30,50, 100 ' ACT, ACW, ASCLE, CINUI, CINUR, CSCLR, CTI, CTR, C1I, C1R, ' C2I, C2R, PTI, PTR, RACT, STI, STR: Double ' I, NW: Integer '----------------------------------------------------------------------- ' I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS ' Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE ' WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU. '----------------------------------------------------------------------- Dim STR As Double NZ = 0 ZBKNU ZRR, ZRI, FNU, KODE, 2, CYR(), CYI(), NW, TOL, ELIM, ALIM If NW <> 0 Then GoTo 50 ZRATI ZRR, ZRI, FNU, N, CYR(), CYI(), TOL '----------------------------------------------------------------------- ' RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z), ' R(FNU+J-1,Z)=Y(J), J=1,...,N '----------------------------------------------------------------------- CINUR = 1# CINUI = 0# If KODE = 1 Then GoTo 10 CINUR = Cos(ZRI) CINUI = Sin(ZRI) '----------------------------------------------------------------------- ' ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH ' THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE ' SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT ' THE RESULT IS ON SCALE. '----------------------------------------------------------------------- 10 ACW = ZABS(CWR(2), CWI(2)) ASCLE = 1000 * D1MACH(1) / TOL CSCLR = 1# If ACW > ASCLE Then GoTo 20 CSCLR = 1# / TOL GoTo 30 20 ASCLE = 1# / ASCLE If ACW < ASCLE Then GoTo 30 CSCLR = TOL 30 C1R = CWR(1) * CSCLR C1I = CWI(1) * CSCLR C2R = CWR(2) * CSCLR C2I = CWI(2) * CSCLR STR = CYR(1) STI = CYI(1) '----------------------------------------------------------------------- ' CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0D0/CABS(CT) PREVENTS ' UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT) '----------------------------------------------------------------------- PTR = STR * C1R - STI * C1I PTI = STR * C1I + STI * C1R PTR = PTR + C2R PTI = PTI + C2I CTR = ZRR * PTR - ZRI * PTI CTI = ZRR * PTI + ZRI * PTR ACT = ZABS(CTR, CTI) RACT = 1# / ACT CTR = CTR * RACT CTI = -CTI * RACT PTR = CINUR * RACT PTI = CINUI * RACT CINUR = PTR * CTR - PTI * CTI CINUI = PTR * CTI + PTI * CTR CYR(1) = CINUR * CSCLR CYI(1) = CINUI * CSCLR If N = 1 Then GoTo 100 For I = 2 To N PTR = STR * CINUR - STI * CINUI CINUI = STR * CINUI + STI * CINUR CINUR = PTR STR = CYR(I) STI = CYI(I) CYR(I) = CINUR * CSCLR CYI(I) = CINUI * CSCLR Next I GoTo 100 50 NZ = -1 If NW = -2 Then NZ = -2 100 End Sub 'ZWRSK 'end of file Cbess1.bas