'****************************************************************** '* Complex Bessel Function of the 1st Kind of integer order * '* -------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Complex Bessel Function of the 1st Kind of integer order * '* * '* Input complex argument (real imaginary): 1,2 * '* Input integer order: 1 * '* * '* Function value: 1.2918478193366672 1.010488365160812 * '* * '* * '* Basic Release 1.2 By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '* -------------------------------------------------------------- * '* Release 1.1: Corrected bug in Subroutine 500 (ATN replaced by * '* Sub 400 - 11/10/2005. * '* Release 1.2: ZPower replaced by IZpower (integer exponant). * '* Limitations increased in CDIV, MAXK=20 * '* Last Modified 09/21/2009. * '****************************************************************** DEFDBL A-H, O-Z DEFINT I-N MAXK = 20 HALF = .5 ONE = 1# FPF = 5.5 'integer nu order of complex Bessel DIM z(2), z1(2) 'Complex numbers 'auxiliary complex numbers DIM zz(2), zz1(2), zz2(2), zz3(2) DIM a(2), b(2), c(2), sum(2), tmp(2), tmp1(2) CLS PRINT PRINT " Complex Bessel Function of the 1st Kind of integer order" PRINT INPUT " Input complex argument (real imaginary): ", z(1), z(2) INPUT " Input integer order: ", nu PRINT GOSUB 2000 'call CBESSJ(z,nu,z1) PRINT " Function value: "; z1(1); " "; z1(2) PRINT END 'of main program 600 'complex division a=b/c D = c(1) * c(1) + c(2) * c(2) IF D > 1D+20 THEN a(1) = 0#: a(2) = 0#: RETURN END IF IF D < 9.999999999999999D-21 THEN RETURN a(1) = (b(1) * c(1) + b(2) * c(2)) / D a(2) = (b(2) * c(1) - b(1) * c(2)) / D RETURN 700 'complex multiplication a=b*c a(1) = b(1) * c(1) - b(2) * c(2) a(2) = b(1) * c(2) + b(2) * c(1) RETURN 910 'zz = zz1^iex IF iex = 0 THEN zz(1) = 1# zz(2) = 0# ELSEIF iex = 1 THEN zz(1) = zz1(1) zz(2) = zz1(2) ELSE tmp1(1) = zz1(1): tmp1(2) = zz1(2) FOR i = 2 TO iex 'CMUL(tmp1,zz1,tmp) tmp(1) = tmp1(1) * zz1(1) - tmp1(2) * zz1(2) tmp(2) = tmp1(1) * zz1(2) + tmp1(2) * zz1(1) tmp1(1) = tmp(1) tmp1(2) = tmp(2) NEXT i zz(1) = tmp(1) zz(2) = tmp(2) END IF RETURN 950 'Fact(k) f = 1# FOR i = 2 TO k f = f * i NEXT i Fact = f RETURN '******************************************* '* FUNCTION GAMMA(X) * '* --------------------------------------- * '* Returns the value of Gamma(x) in double * '* precision as EXP(LN(GAMMA(X))) for X>0. * '******************************************* 1000 'Gamma(xx) DIM cof(6) cof(1) = 76.18009173# cof(2) = -86.50532033# cof(3) = 24.01409822# cof(4) = -1.231739516# cof(5) = .00120858003# cof(6) = -.00000536382# stp = 2.50662827465# x = xx - ONE tmp = x + FPF tmp = (x + HALF) * LOG(tmp) - tmp ser = 1# FOR j = 1 TO 6 x = x + ONE ser = ser + cof(j) / x NEXT j Gamma = EXP(tmp + LOG(stp * ser)) RETURN 'main complex Bessel of first kind subroutine 2000 'CBESSJ(z, nu, z1) '-------------------------------------------------- ' inf. (-z^2/4)^k ' Jnu(z) = (z/2)^nu x Sum ------------------ ' k=0 k! x Gamma(nu+k+1) ' (nu must be >= 0). Here k=15. '-------------------------------------------------- 'calculate (z/2)^nu in tmp b(1) = z(1): b(2) = z(2): c(1) = 2#: c(2) = 0#: GOSUB 600 zz1(1) = a(1): zz1(2) = a(2) 'zz1=z/2 iex = nu: GOSUB 910 temp(1) = zz(1): temp(2) = zz(2) 'temp=(z/2)^nu sum(1) = 0#: sum(2) = 0# 'calculate Sum FOR k = 0 TO MAXK 'calculate (-z^2/4)^k b(1) = z(1): b(2) = z(2): c(1) = z(1): c(2) = z(2): GOSUB 700 b(1) = -a(1): b(2) = -a(2) 'b=-z^2 c(1) = 4#: c(2) = 0#: GOSUB 600 zz1(1) = a(1): zz1(2) = a(2) 'zz1=(-z^2/4) iex = k: GOSUB 910 b(1) = zz(1): b(2) = zz(2) 'b=(-z^2/4)^k 'divide by k! GOSUB 950 'Fact = k! c(1) = Fact: c(2) = 0#: GOSUB 600 b(1) = a(1): b(2) = a(2) 'divide by Gamma(nu+k+1) xx = 1# * (nu + k + 1): GOSUB 1000 c(1) = Gamma: c(2) = 0#: GOSUB 600 'actualize sum sum(1) = sum(1) + a(1) sum(2) = sum(2) + a(2) NEXT k 'multiply (z/2)^nu by sum 2010 b(1) = temp(1): b(2) = temp(2) c(1) = sum(1): c(2) = sum(2): GOSUB 700 z1(1) = a(1): z1(2) = a(2) RETURN ' end of file cbessj.bas