'******************************************************************
'*    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