tzbesk_cpp.txt

BESSEL PROGRAMS IN C/C

!         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
!         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
!         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
!         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
!         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
!         OR -PI/2+P.
!
!***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
!                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
!                 COMMERCE, 1955.
!
!               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
!                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
!
!               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
!                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983.
!
!               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
!                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
!                 1018, MAY, 1985
!
!               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
!                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
!                 MATH. SOFTWARE, 1986
!
!***ROUTINES CALLED  ZACON,ZBKNU,ZBUNK,ZUOIK,ZABS,I1MACH,D1MACH
!***END PROLOGUE  ZBESK
!
!     COMPLEX CY, Z */
//Labels: e50,e60,e70,e80,e90,e100,e180,e200,e260

      REAL AA, ALIM, ALN, ARG, AZ, DIG, ELIM, FN, FNUL, RL, R1M5, TOL, UFL, BB;
      int K, K1, K2, MR, NN, NUF, NW;

//***FIRST EXECUTABLE STATEMENT ZBESK
      *IERR = 0;
      *NZ=0;
      if (ZI == 0.0 && ZR == 0.0)  *IERR=1;
      if (FNU < 0.0)  *IERR=1;
      if (KODE < 1 || KODE > 2)  *IERR=1;
      if (N < 1) *IERR=1;
      if (*IERR != 0) return;    //bad parameter(s)
      NN = N;
/*----------------------------------------------------------------------
!     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
!     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
!     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
!     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
!     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
!     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
!     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
!     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
!     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
!---------------------------------------------------------------------*/
      TOL = DMAX(D1MACH(4),1e-18);
      K1 = I1MACH(15);
      K2 = I1MACH(16);
      R1M5 = D1MACH(5);
      K = IMIN(ABS(K1),ABS(K2));
      ELIM = 2.303*(K*R1M5-3.0);
      K1 = I1MACH(14) - 1;
      AA = R1M5*K1;
      DIG = DMIN(AA,18.0);
      AA = AA*2.303;
      ALIM = ELIM + DMAX(-AA,-41.45);
      FNUL = 10.0 + 6.0*(DIG-3.0);
      RL = 1.2*DIG + 3.0;
/*----------------------------------------------------------------------
!     TEST FOR PROPER RANGE
!---------------------------------------------------------------------*/
      AZ = ZABS(ZR,ZI);
      FN = FNU + 1.0*(NN-1);
      AA = 0.5/TOL;
      BB=0.5*I1MACH(9);
      AA = DMIN(AA,BB);
      if (AZ > AA) goto e260;
      if (FN > AA) goto e260;
      AA = SQRT(AA);
      if (AZ > AA)  *IERR=3;
      if (FN > AA)  *IERR=3;
/*----------------------------------------------------------------------
!     OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
!---------------------------------------------------------------------*/
      UFL = EXP(-ELIM);
      if (AZ < UFL) goto e180;
      if (FNU > FNUL) goto e80;
      if (FN <= 1.0) goto e60;
      if (FN > 2.0)  goto e50;
      if (AZ > TOL)  goto e60;
      ARG = 0.5*AZ;
      ALN = -FN*log(ARG);
      if (ALN > ELIM) goto e180;
      goto e60;
e50:  ZUOIK(ZR, ZI, FNU, KODE, 2, NN, CYR, CYI, &NUF, TOL, ELIM, ALIM);
      if (NUF < 0) goto e180;
      *NZ = *NZ + NUF;
      NN = NN - NUF;
/*----------------------------------------------------------------------
!     HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON return FROM CUOIK
!     IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
!---------------------------------------------------------------------*/
      if (NN == 0) goto e100;
e60:  if (ZR < 0.0) goto e70;
/*----------------------------------------------------------------------
!     RIGHT HALF PLANE COMPUTATION, REAL(Z).GE.0.
!---------------------------------------------------------------------*/
      ZBKNU(ZR, ZI, FNU, KODE, NN, CYR, CYI, &NW, TOL, ELIM, ALIM);
      if (NW < 0) goto e200;
      *NZ=NW;
      return;
/*----------------------------------------------------------------------
!     LEFT HALF PLANE COMPUTATION
!     PI/2.LT.ARG(Z).LE.PI AND -PI.LT.ARG(Z).LT.-PI/2.
!---------------------------------------------------------------------*/
e70:  if (*NZ != 0) goto e180;
      MR = 1;
      if (ZI < 0.0)  MR = -1;
      ZACON(ZR, ZI, FNU, KODE, MR, NN, CYR, CYI, &NW, RL, FNUL, TOL, ELIM, ALIM);
      if (NW < 0) goto e200;
      *NZ=NW;
      return;
/*----------------------------------------------------------------------
!     UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
!---------------------------------------------------------------------*/
e80:  MR = 0;
      if (ZR >= 0.0) goto e90;
      MR = 1;
      if (ZI < 0.0)  MR = -1;
e90:  ZBUNK(ZR, ZI, FNU, KODE, MR, NN, CYR, CYI, &NW, TOL, ELIM, ALIM);
      if (NW < 0) goto e200;
      *NZ = *NZ + NW;
      return;
e100: if (ZR < 0.0) goto e180;
      return;
e180: *NZ = 0;
      *IERR=2;
      return;
e200: if (NW == -1) goto e180;
      *NZ=0;
      *IERR=5;
      return;
e260: *NZ=0;
      *IERR=4;
} //ZBESK()


void main()  {

  n=5;

//memory allocation for cyr, cyi
  vmblock = vminit();  
  cyr = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0); //index 0 not used
  cyi = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);

  if (! vmcomplete(vmblock)) {
    LogError ("No Memory", 0, __FILE__, __LINE__);
    return;
  } 	  

  zr=1.0; zi=2.0;

  ZBESK(zr,zi,0.0,1,n,cyr,cyi,&nz,&ierr);

  printf("\n");
  for (i=1; i<=n; i++) {
    printf(" zr(%d) = %10.6f\n", i-1, cyr[i]);
    printf(" zi(%d) = %10.6f\n", i-1, cyi[i]);
  }
  printf(" NZ = %d\n", nz);
  printf(" Error code: %d\n\n", ierr);
  vmfree(vmblock);
  
}

// end of file tzbesk.cpp