tzbesk_cpp.txt

BESSEL PROGRAMS IN C/C

/****************************************************************
* EVALUATE A K-BESSEL FUNCTION OF COMPLEX ARGUMENT (THIRD KIND) *
* ------------------------------------------------------------- *
* SAMPLE RUN:                                                   *
* (Evaluate K0 to K4 for argument Z=(1.0,2.0) ).                *
*                                                               *
* zr(0) =  -0.242345                                            *
* zi(0) =  -0.176267                                            *
* zr(1) =  -0.300362                                            *
* zi(1) =  -0.151186                                            *
* zr(2) =  -0.483439                                            *
* zi(2) =   0.003548                                            *
* zr(3) =  -0.681436                                            *
* zi(3) =   0.625155                                            *
* zr(4) =   0.199208                                            *
* zi(4) =   2.389181                                            *
* NZ = 0                                                        *
* Error code: 0                                                 *
*                                                               *
* ------------------------------------------------------------- *
* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77.   *
*                                                               *
*                        C++ Release 1.0 By J-P Moreau, Paris   *
*****************************************************************
Note: To link with: CBess0,CBess00,CBess1,CBess2,CBess3,Complex,
                    Basis_r, Vmblock.
--------------------------------------------------------------- */ 
#include <basis.h>
#include <vmblock.h>

#include "complex.h"


    REAL zr,zi;
    REAL *cyr, *cyi; //pointers to vectors of size n+1
    int  i,ierr,n,nz;
    void *vmblock = NULL;

	void ZUOIK(REAL, REAL, REAL, int, int, int, REAL *, REAL *, 
		       int *, REAL, REAL, REAL);

	void ZBKNU(REAL, REAL, REAL, int, int, REAL *, REAL *, int *, REAL, 
	           REAL, REAL);

	void ZACON(REAL, REAL, REAL, int, int, int, REAL *, REAL *, int *, 
		       REAL, REAL, REAL, REAL, REAL);

	void ZBUNK(REAL, REAL, REAL, int, int, int, REAL *, REAL *, int *, 
		       REAL, REAL, REAL);



void ZBESK(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *CYR, 
           REAL *CYI, int *NZ, int *IERR)  {
/***BEGIN PROLOGUE  ZBESK
!***DATE WRITTEN   830501   (YYMMDD)  (Original Fortran Version).
!***REVISION DATE  830501   (YYMMDD)
!***CATEGORY NO.  B5K
!***KEYWORDS  K-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION,
!             MODIFIED BESSEL FUNCTION OF THE SECOND KIND,
!             BESSEL FUNCTION OF THE THIRD KIND
!***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
!***PURPOSE  TO COMPUTE K-BESSEL FUNCTIONS OF COMPLEX ARGUMENT
!***DESCRIPTION
!
!                      ***A DOUBLE PRECISION ROUTINE***
!
!         ON KODE=1, CBESK COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
!         BESSEL FUNCTIONS CY(J)=K(FNU+J-1,Z) FOR REAL, NONNEGATIVE
!         ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z.NE.CMPLX(0.0,0.0)
!         IN THE CUT PLANE -PI < ARG(Z) <= PI. ON KODE=2, CBESK
!         returnS THE SCALED K FUNCTIONS,
!
!         CY(J)=EXP(Z)*K(FNU+J-1,Z) , J=1,...,N,
!
!         WHICH REMOVE THE EXPONENTIAL BEHAVIOR IN BOTH THE LEFT AND
!         RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND
!         NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL
!         FUNCTIONS (REF. 1).
!
!         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
!           ZR,ZI  - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0),
!                    -PI.LT.ARG(Z).LE.PI
!           FNU    - ORDER OF INITIAL K FUNCTION, FNU.GE.0.0D0
!           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
!           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
!                    KODE= 1  returnS
!                             CY(I)=K(FNU+I-1,Z), I=1,...,N
!                        = 2  returnS
!                             CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
!
!         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
!           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
!                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
!                    CY(I)=K(FNU+I-1,Z), I=1,...,N OR
!                    CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
!                    DEPENDING ON KODE
!           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW.
!                    NZ= 0   , NORMAL return
!                    NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE
!                              TO UNDERFLOW, CY(I)=CMPLX(0.0D0,0.0D0),
!                              I=1,...,N WHEN X >= 0. WHEN X < 0,
!                              NZ STATES ONLY THE NUMBER OF UNDERFLOWS
!                              IN THE SEQUENCE.
!
!           IERR   - ERROR FLAG
!                    IERR=0, NORMAL return - COMPUTATION COMPLETED
!                    IERR=1, INPUT ERROR   - NO COMPUTATION
!                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU IS
!                            TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
!                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
!                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
!                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
!                            ACCURACY
!                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
!                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
!                            CANCE BY ARGUMENT REDUCTION
!                    IERR=5, ERROR              - NO COMPUTATION,
!                            ALGORITHM TERMINATION CONDITION NOT MET
!
!***LONG DESCRIPTION
!
!         EQUATIONS OF THE REFERENCE ARE IMPLEMENTED FOR SMALL ORDERS
!         DNU AND DNU+1.0 IN THE RIGHT HALF PLANE X.GE.0.0. FORWARD
!         RECURRENCE GENERATES HIGHER ORDERS. K IS CONTINUED TO THE LEFT
!         HALF PLANE BY THE RELATION
!
!         K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z)
!         MP=MR*PI*I, MR=+1 OR -1, RE(Z) > 0, I^2=-1
!
!         WHERE I(FNU,Z) IS THE I BESSEL FUNCTION.
!
!         FOR LARGE ORDERS, FNU > FNUL, THE K FUNCTION IS COMPUTED
!         BY MEANS OF ITS UNIFORM ASYMPTOTIC EXPANSIONS.
!
!         FOR NEGATIVE ORDERS, THE FORMULA
!
!                       K(-FNU,Z) = K(FNU,Z)
!
!         CAN BE USED.
!
!         CBESK ASSUMES THAT A SIGNIFICANT DIGIT SINH(X) FUNCTION IS
!         AVAILABLE.
!
!         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
!         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
!         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
!         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
!         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
!         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
!         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
!         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
!         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
!         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
!         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
!         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
!         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
!         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
!         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
!         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
!         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
!         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
!         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
!
!         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
!         BESSEL FUNCTION CAN BE EXPRESSED BY P*10^S WHERE P=MAX(UNIT
!         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10^S REPRE-
!         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
!         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
!         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
!         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
!         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
!         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
!         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10^K LARGER
!         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
!         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
!         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER