Attribute VB_Name = "Module2" DefDbl A-H, O-Z DefInt I-N ' --------------------------------------------------------------------- ' Utility subroutines used by any program from Numath library ' with complex numbers z = (zr,zi). ' --------------------------------------------------------------------- ' Reference: From Numath Library By Tuan Dang Trong in Fortran 77 ' [BIBLI 18]. ' ' Visual Basic Release 1.0 By J-P Moreau, Paris. ' (www.jpmoreau.fr) ' --------------------------------------------------------------------- 'MODULE COMPLEX 'LIST OF SUBROUTINES ' Function ZABS(AR1, AI1): Double ' Subroutine ZSQRT(AR, AI, BR, BI) ' Subroutine ZEXP(AR, AI, BR, BI) ' Subroutine ZMLT(AR, AI, BR, BI, CR,CI) ' Subroutine ZDIV(AR, AI, BR, BI, CR,CI) ' Subroutine ZLOG(AR, AI, BR, BI, IERR) ' Function D1MACH(Ir): Double ' Function XI1MACH(I): Double ' Function DMAX(a,b): Double ' Function DMIN(a,b): Double ' Function IMAX(ia,ib):integer ' Function IMIN(ia,ibr):integer ' Function ZABS(ByVal AR1 As Double, ByVal AI1 As Double) '***BEGIN PROLOGUE ZABS '***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY ' ZABS COMPUTES THE ABSOLUTE VALUE OR MAGNITUDE OF A DOUBLE ' PRECISION COMPLEX VARIABLE CMPLX(AR1,AI1) '***ROUTINES CALLED (NONE) '***END PROLOGUE ZABS 'Labels: 510, 520 ' U, V, Q, S: double U = Abs(AR1) V = Abs(AI1) S = U + V '--------------------------------------------------------------------- ' S*1.0 MAKES AN UNNORMALIZED UNDERFLOW ON CDC MACHINES INTO A ' TRUE FLOATING ZERO '--------------------------------------------------------------------} S = S * 1# If S = 0# Then GoTo 520 If U > V Then GoTo 510 Q = U / V ZABS = V * Sqr(1# + Q * Q) GoTo 530 510 Q = V / U ZABS = U * Sqr(1# + Q * Q) GoTo 530 520 ZABS = 0# 530 End Function 'ZABS Sub ZSQRT(ByVal AR, ByVal AI, BR, BI) '***BEGIN PROLOGUE ZSQRT '***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY ' ' DOUBLE PRECISION COMPLEX SQUARE ROOT, B=CSQRT(A) ' '***ROUTINES CALLED ZABS '***END PROLOGUE ZSQRT} 'Labels: 610,620,630,640,650,660,670 ' ZM, DTHETA, DPI, DRT: double DRT = 0.707106781186548 DPI = 3.14159265358979 ZM = ZABS(AR, AI) ZM = Sqr(ZM) If AR = 0# Then GoTo 610 If AI = 0# Then GoTo 620 DTHETA = Atn(AI / AR) If DTHETA <= 0# Then GoTo 640 If AR < 0# Then DTHETA = DTHETA - DPI GoTo 650 610 If AI > 0# Then GoTo 660 If AI < 0# Then GoTo 670 BR = 0# BI = 0# GoTo 680 620 If AR > 0# Then GoTo 630 BR = 0# BI = Sqr(Abs(AR)) GoTo 680 630 BR = Sqr(AR) BI = 0# GoTo 680 640 If AR < 0# Then DTHETA = DTHETA + DPI 650 DTHETA = DTHETA * 0.5 BR = ZM * Cos(DTHETA) BI = ZM * Sin(DTHETA) GoTo 680 660 BR = ZM * DRT BI = ZM * DRT GoTo 680 670 BR = ZM * DRT BI = -ZM * DRT 680 End Sub 'ZSQRT Sub ZEXP(ByVal AR, ByVal AI, BR, BI) '***BEGIN PROLOGUE ZEXP '***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY ' ' DOUBLE PRECISION COMPLEX EXPONENTIAL FUNCTION B=EXP(A) ' '***ROUTINES CALLED (NONE) '***END PROLOGUE ZEXP ' ZM, CA, CB: double ZM = Exp(AR) CA = ZM * Cos(AI) CB = ZM * Sin(AI) BR = CA BI = CB End Sub 'ZEXP Sub ZMLT(AR, AI, BR, BI, CR, CI) '***BEGIN PROLOGUE ZMLT '***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY ' ' DOUBLE PRECISION COMPLEX MULTIPLY, C=A*B. ' '***ROUTINES CALLED (NONE) '***END PROLOGUE ZMLT ' CA, CB: double CA = AR * BR - AI * BI CB = AR * BI + AI * BR CR = CA CI = CB End Sub 'ZMLT Sub ZDIV(AR, AI, BR, BI, CR, CI) '***BEGIN PROLOGUE ZDIV '***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY ' ' DOUBLE PRECISION COMPLEX DIVIDE C=A/B. '***ROUTINES CALLED ZABS '***END PROLOGUE ZDIV ' BM, CA, CB, CC, CD: double BM = 1# / ZABS(BR, BI) CC = BR * BM CD = BI * BM CA = (AR * CC + AI * CD) * BM CB = (AI * CC - AR * CD) * BM CR = CA CI = CB End Sub 'ZDIV Sub ZLOG(AR, AI, BR, BI, IERR) '***BEGIN PROLOGUE ZLOG '***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY ' DOUBLE PRECISION COMPLEX LOGARITHM B=CLOG(A) ' IERR=0,NORMAL RETURN IERR=1, Z=CMPLX(0.0,0.0) '***ROUTINES CALLED ZABS '***END PROLOGUE ZLOG } 'Labels: 1010,1020,1030,1040,1050,1060 'ZM, DTHETA, DPI, DHPI: double DPI = 3.14159265358979 DHPI = 1.5707963267949 IERR = 0 If AR = 0# Then GoTo 1010 If AI = 0# Then GoTo 1020 DTHETA = Atn(AI / AR) If DTHETA <= 0# Then GoTo 1040 If AR < 0# Then DTHETA = DTHETA - DPI GoTo 1050 1010 If AI = 0# Then GoTo 1060 BI = DHPI BR = Log(Abs(AI)) If AI < 0# Then BI = -BI GoTo 1070 1020 If AR > 0# Then GoTo 1030 BR = Log(Abs(AR)) BI = DPI GoTo 1070 1030 BR = Log(AR) BI = 0# GoTo 1070 1040 If AR < 0# Then DTHETA = DTHETA + DPI 1050 ZM = ZABS(AR, AI) BR = Log(ZM) BI = DTHETA GoTo 1070 1060 IERR = 1 1070 End Sub 'ZLOG Function D1MACH(II) '***BEGIN PROLOGUE D1MACH '***DATE WRITTEN 750101 (YYMMDD) '***REVISION DATE 860501 (YYMMDD) '***CATEGORY NO. R1 '***KEYWORDS MACHINE CONSTANTS '***AUTHOR FOX, P. A., (BELL LABS) ' HALL, A. D., (BELL LABS) ' SCHRYER, N. L., (BELL LABS) '***PURPOSE RETURN DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS. '***DESCRIPTION ' D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS ' FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION ' SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED ' AS FOLLOWS, FOR EXAMPLE ' D = D1MACH(II) ' WHERE II=1,...,5. THE (OUTPUT) VALUE OF D ABOVE IS ' DETERMINED BY THE (INPUT) VALUE OF II. THE RESULTS FOR ' VARIOUS VALUES OF II ARE DISCUSSED BELOW. ' DOUBLE-PRECISION MACHINE CONSTANTS ' D1MACH( 1) = B^(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. ' D1MACH( 2) = B^EMAX*(1 - B^(-T)), THE LARGEST MAGNITUDE. ' D1MACH( 3) = B^(-T), THE SMALLEST RELATIVE SPACING. ' D1MACH( 4) = B^(1-T), THE LARGEST RELATIVE SPACING. ' D1MACH( 5) = LOG10(B) '***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A ' PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL ' SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. '***ROUTINES CALLED XERROR '***END PROLOGUE D1MACH Dim DMACH(5) As Double '***FIRST EXECUTABLE STATEMENT D1MACH If II < 1 Or II > 5 Then Form1.Print " D1MACH -- I OUT OF BOUNDS" End If 'For IBM PC or APOLLO DMACH(1) = 2.22559E-308 DMACH(2) = 1.79728E+308 DMACH(3) = 1.11048E-16 DMACH(4) = 2.22096E-16 DMACH(5) = 0.301029995663981 D1MACH = DMACH(II) End Function 'D1MACH} Function XI1MACH(II) '***BEGIN PROLOGUE I1MACH '***DATE WRITTEN 750101 (YYMMDD) '***REVISION DATE 890313 (YYMMDD) '***CATEGORY NO. R1 '***KEYWORDS LIBRARY=SLATEC,TYPE=INTEGER(I1MACH-I),MACHINE CONSTANTS '***AUTHOR FOX, P. A., (BELL LABS) ' HALL, A. D., (BELL LABS) ' SCHRYER, N. L., (BELL LABS) '***PURPOSE Return integer (here double) machine dependent constants. '***DESCRIPTION ' XI1MACH can be used to obtain machine-dependent parameters ' for the local machine environment. It is a function ' subroutine with one (input) argument, and can be called ' as follows, for example ' XK = XI1MACH(II) ' where II=1,...,16. The (output) value of XK above is ' determined by the (input) value of II. The results for ' various values of II are discussed below. ' I/O unit numbers. ' XI1MACH(1) = the standard input unit. ' XI1MACH(2) = the standard output unit. ' XI1MACH(3) = the standard punch unit. ' XI1MACH(4) = the standard error message unit. ' Words. ' XI1MACH(5) = the number of bits per integer storage unit. ' XI1MACH(6) = the number of characters per integer storage unit. ' Integers. ' assume integers are represented in the S-digit, base-A form ' sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) ' where 0 .LE. X(I) .LT. A for I=0,...,S-1. ' XI1MACH(7) = A, the base. ' XI1MACH(8) = S, the number of base-A digits. ' XI1MACH(9) = A^S - 1, the largest magnitude. ' Floating-Point Numbers. ' Assume floating-point numbers are represented in the T-digit, ' base-B form ' sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) ' where 0 .LE. X(I) .LT. B for II=1,...,T, ' 0 .LT. X(1), and EMIN .LE. E .LE. EMAX. ' XI1MACH(10) = B, the base. ' Single-Precision ' XI1MACH(11) = T, the number of base-B digits. ' XI1MACH(12) = EMIN, the smallest exponent E. ' XI1MACH(13) = EMAX, the largest exponent E. ' Double-Precision ' XI1MACH(14) = T, the number of base-B digits. ' XI1MACH(15) = EMIN, the smallest exponent E. ' XI1MACH(16) = EMAX, the largest exponent E. ' The values of XI1MACH(1) - XI1MACH(4) should be checked for ' consistency with the local operating system. '***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A ' PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL ' SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. '***ROUTINES CALLED (NONE) '***END PROLOGUE XI1MACH Dim XIMACH(16) As Double '***FIRST EXECUTABLE STATEMENT XD1MACH If II < 1 Or II > 16 Then Form1.Print " XI1MACH -- I OUT OF BOUNDS" End If 'For IBM PC or APOLLO XIMACH(1) = 5# XIMACH(2) = 6# XIMACH(3) = 0# XIMACH(4) = 0# XIMACH(5) = 32# XIMACH(6) = 4# XIMACH(7) = 2# XIMACH(8) = 31# XIMACH(9) = 2147483647# XIMACH(10) = 2# XIMACH(11) = 24# XIMACH(12) = -125# XIMACH(13) = 127# XIMACH(14) = 53# XIMACH(15) = -1021# XIMACH(16) = 1023# XI1MACH = XIMACH(II) End Function 'XI1MACH Function DMAX(a, b) If a >= b Then DMAX = a Else DMAX = b End If End Function Function DMIN(ByVal a, ByVal b) If a <= b Then DMIN = a Else DMIN = b End If End Function Function IMAX(ia, ib) If ia >= ib Then IMAX = ia Else IMAX = ib End If End Function Function IMIN(ia, ib) If ia <= ib Then IMIN = ia Else IMIN = ib End If End Function 'end of file complex.bas