cbess0_cpp.txt

BESSEL PROGRAMS IN C/C

/*************************************************************************
*            Functions used By programs TZBESJ, TZBESK, TZBESY           *
*    (Evalute Bessel Functions with complex argument, 1st to 3rd kind)   *
* ---------------------------------------------------------------------- *
* Reference:  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES, 1983.       *
*                                                                        *
*                         C++ Release By J-P Moreau, Paris (07/01/2005). *
*************************************************************************/
// Functions using only module COMPLEX.cpp
#include <math.h>
#include <basis.h>
#include <vmblock.h>

#include "complex.h"


REAL Sign(REAL a, REAL b) {
  if (b <0.0) return (-ABS(a));
  else return ABS(a);
}


void ZUCHK(REAL YR, REAL YI, int *NZ, REAL ASCLE, REAL TOL) {
/***BEGIN PROLOGUE  ZUCHK
!***REFER TO ZSERI,ZUOIK,ZUNK1,ZUNK2,ZUNI1,ZUNI2,ZKSCL
!
!   Y ENTERS AS A SCALED QUANTITY WHOSE MAGNITUDE IS GREATER THAN
!   EXP(-ALIM)=ASCLE=1000*D1MACH(1)/TOL. THE TEST IS MADE TO SEE
!   IF THE MAGNITUDE OF THE REAL OR IMAGINARY PART WOULD UNDERFLOW
!   WHEN Y IS SCALED (BY TOL) TO ITS PROPER VALUE. Y IS ACCEPTED
!   IF THE UNDERFLOW IS AT LEAST ONE PRECISION BELOW THE MAGNITUDE
!   OF THE LARGEST COMPONENT; OTHERWISE THE PHASE ANGLE DOES NOT HAVE
!   ABSOLUTE ACCURACY AND AN UNDERFLOW IS ASSUMED.
!
!***ROUTINES CALLED  (NONE)
!***END PROLOGUE  ZUCHK
!
!     COMPLEX Y */
REAL SS, ST, WR, WI;
      *NZ = 0;
      WR = ABS(YR);
      WI = ABS(YI);
      ST = DMIN(WR,WI);
      if (ST > ASCLE) return;
      SS = DMAX(WR,WI);
      ST = ST/TOL;
      if (SS < ST) *NZ = 1;
} //ZUCHK()


REAL DGAMLN(REAL Z, int *IERR) {
/***BEGIN PROLOGUE  DGAMLN
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  830501   (YYMMDD)
!***CATEGORY NO.  B5F
!***KEYWORDS  GAMMA FUNCTION,LOGARITHM OF GAMMA FUNCTION
!***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
!***PURPOSE  TO COMPUTE THE LOGARITHM OF THE GAMMA FUNCTION
!***DESCRIPTION
!
!               **** A REAL PRECISION ROUTINE ****
!         DGAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR
!         Z > 0.  THE ASYMPTOTIC EXPANSION IS USED TO GENERATE VALUES
!         GREATER THAN ZMIN WHICH ARE ADJUSTED BY THE RECURSION
!         G(Z+1)=Z*G(Z) FOR Z <= ZMIN.  THE FUNCTION WAS MADE AS
!         PORTABLE AS POSSIBLE BY COMPUTIMG ZMIN FROM THE NUMBER OF BASE
!         10 DIGITS IN A WORD, RLN=MAX(-ALOG10(R1MACH(4)),0.5E-18)
!         LIMITED TO 18 DIGITS OF (RELATIVE) ACCURACY.
!
!         SINCE INTEGER ARGUMENTS ARE COMMON, A TABLE LOOK UP ON 100
!         VALUES IS USED FOR SPEED OF EXECUTION.
!
!     DESCRIPTION OF ARGUMENTS
!
!         INPUT      Z IS D0UBLE PRECISION
!           Z      - ARGUMENT, Z > 0
!
!         OUTPUT      DGAMLN IS REAL PRECISION
!           DGAMLN  - NATURAL LOG OF THE GAMMA FUNCTION AT Z <> 0
!           IERR    - ERROR FLAG
!                     IERR=0, NORMAL RETURN, COMPUTATION COMPLETED
!                     IERR=1, Z <= 0.0,      NO COMPUTATION
!
!
!***REFERENCES  COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
!                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
!***ROUTINES CALLED  I1MACH,D1MACH (of module COMPLEX).
!***END PROLOGUE  DGAMLN */
//Labels: e10,e20,e40,e50,e70

      REAL CON, FLN, FZ, RLN, S, TLG, TRM, TST, T1, WDTOL, ZDMY,
      ZINC, ZM, ZMIN, ZP, ZSQ;
      int I, I1M, K, MZ, NZ;
      REAL *CF;
      REAL *GLN;

	  void *vmblock = NULL;

//    LNGAMMA(N), N=1,100
      vmblock = vminit();  
      GLN = (REAL *) vmalloc(vmblock, VEKTOR,  101, 0); //index 0 not used
      CF  = (REAL *) vmalloc(vmblock, VEKTOR,  23, 0); 
      if (! vmcomplete(vmblock)) {
        LogError ("No Memory", 0, __FILE__, __LINE__);
        return 0;
	  } 
      GLN[1]= 0.00000000000000000;     GLN[2]= 0.00000000000000000;
      GLN[3]= 6.93147180559945309e-01; GLN[4]= 1.79175946922805500;
      GLN[5]= 3.17805383034794562;     GLN[6]= 4.78749174278204599;
      GLN[7]= 6.57925121201010100;     GLN[8]= 8.52516136106541430;
      GLN[9]= 1.06046029027452502e+01; GLN[10]=1.28018274800814696e+01;
      GLN[11]=1.51044125730755153e+01; GLN[12]=1.75023078458738858e+01;
      GLN[13]=1.99872144956618861e+01; GLN[14]=2.25521638531234229e+01;
      GLN[15]=2.51912211827386815e+01; GLN[16]=2.78992713838408916e+01;
      GLN[17]=3.06718601060806728e+01; GLN[18]=3.35050734501368889e+01;
      GLN[19]=3.63954452080330536e+01; GLN[20]=3.93398841871994940e+01;
      GLN[21]=4.23356164607534850e+01; GLN[22]=4.53801388984769080e+01;

      GLN[23]=4.84711813518352239e+01; GLN[24]=5.16066755677643736e+01;
      GLN[25]=5.47847293981123192e+01; GLN[26]=5.80036052229805199e+01;
      GLN[27]=6.12617017610020020e+01; GLN[28]=6.45575386270063311e+01;
      GLN[29]=6.78897431371815350e+01; GLN[30]=7.12570389671680090e+01;
      GLN[31]=7.46582363488301644e+01; GLN[32]=7.80922235533153106e+01;
      GLN[33]=8.15579594561150372e+01; GLN[34]=8.50544670175815174e+01;
      GLN[35]=8.85808275421976788e+01; GLN[36]=9.21361756036870925e+01;
      GLN[37]=9.57196945421432025e+01; GLN[38]=9.93306124547874269e+01;
      GLN[39]=1.02968198614513813e+02; GLN[40]=1.06631760260643459e+02;
      GLN[41]=1.10320639714757395e+02; GLN[42]=1.14034211781461703e+02;
      GLN[43]=1.17771881399745072e+02; GLN[44]=1.21533081515438634e+02;

      GLN[45]=1.25317271149356895e+02; GLN[46]=1.29123933639127215e+02;
      GLN[47]=1.32952575035616310e+02; GLN[48]=1.36802722637326368e+02;
      GLN[49]=1.40673923648234259e+02; GLN[50]=1.44565743946344886e+02;
      GLN[51]=1.48477766951773032e+02; GLN[52]=1.52409592584497358e+02;
      GLN[53]=1.56360836303078785e+02; GLN[54]=1.60331128216630907e+02;
      GLN[55]=1.64320112263195181e+02; GLN[56]=1.68327445448427652e+02;
      GLN[57]=1.72352797139162802e+02; GLN[58]=1.76395848406997352e+02;
      GLN[59]=1.80456291417543771e+02; GLN[60]=1.84533828861449491e+02;
      GLN[61]=1.88628173423671591e+02; GLN[62]=1.92739047287844902e+02;
      GLN[63]=1.96866181672889994e+02; GLN[64]=2.01009316399281527e+02;
      GLN[65]=2.05168199482641199e+02; GLN[66]=2.09342586752536836e+02;

      GLN[67]=2.13532241494563261e+02; GLN[68]=2.17736934113954227e+02;
      GLN[69]=2.21956441819130334e+02; GLN[70]=2.26190548323727593e+02;
      GLN[71]=2.30439043565776952e+02; GLN[72]=2.34701723442818268e+02;
      GLN[73]=2.38978389561834323e+02; GLN[74]=2.43268849002982714e+02;
      GLN[75]=2.47572914096186884e+02; GLN[76]=2.51890402209723194e+02;
      GLN[77]=2.56221135550009525e+02; GLN[78]=2.60564940971863209e+02;
      GLN[79]=2.64921649798552801e+02; GLN[80]=2.69291097651019823e+02;
      GLN[81]=2.73673124285693704e+02; GLN[82]=2.78067573440366143e+02;
      GLN[83]=2.82474292687630396e+02; GLN[84]=2.86893133295426994e+02;
      GLN[85]=2.91323950094270308e+02; GLN[86]=2.95766601350760624e+02;
      GLN[87]=3.00220948647014132e+02; GLN[88]=3.04686856765668715e+02;

      GLN[89]=3.09164193580146922e+02; GLN[90]=3.13652829949879062e+02;
      GLN[91]=3.18152639620209327e+02; GLN[92]=3.22663499126726177e+02;
      GLN[93]=3.27185287703775217e+02; GLN[94]=3.31717887196928473e+02;
      GLN[95]=3.36261181979198477e+02; GLN[96]=3.40815058870799018e+02;
      GLN[97]=3.45379407062266854e+02; GLN[98]=3.49954118040770237e+02;
      GLN[99]=3.54539085519440809e+02; GLN[100]=3.59134205369575399e+02;

//    COEFFICIENTS OF ASYMPTOTIC EXPANSION
      CF[1]= 8.33333333333333333e-02; CF[2]= -2.77777777777777778e-03;
      CF[3]= 7.93650793650793651e-04; CF[4]= -5.95238095238095238e-04;
      CF[5]= 8.41750841750841751e-04; CF[6]= -1.91752691752691753e-03;
      CF[7]= 6.41025641025641026e-03; CF[8]= -2.95506535947712418e-02;
      CF[9]= 1.79644372368830573e-01; CF[10]=-1.39243221690590112e+00;
      CF[11]=1.34028640441683920e+01; CF[12]=-1.56848284626002017e+02;
      CF[13]=2.19310333333333333e+03; CF[14]=-3.61087712537249894e+04;
      CF[15]=6.91472268851313067e+05; CF[16]=-1.52382215394074162e+07;
      CF[17]=3.82900751391414141e+08; CF[18]=-1.08822660357843911e+10;
      CF[19]=3.47320283765002252e+11; CF[20]=-1.23696021422692745e+13;
      CF[21]=4.88788064793079335e+14; CF[22]=-2.13203339609193739e+16;

//    LN(2*PI)
      CON=1.83787706640934548;

      *IERR=0;
      if (Z <= 0.0) goto e70;
      if (Z > 101.0) goto e10;
      NZ = (int) floor(Z);
      FZ = Z - 1.0*NZ;
      if (FZ > 0.0) goto e10;
      if (NZ > 100) goto e10;
	  vmfree(vmblock);
      return (GLN[NZ]);

e10:  WDTOL = D1MACH(4);
      WDTOL = DMAX(WDTOL,0.5e-18);
      I1M = I1MACH(14);
      RLN = D1MACH(5)*I1M;
      FLN = DMIN(RLN,20.0);
      FLN = DMAX(FLN,3.0);
      FLN = FLN - 3.0;
      ZM = 1.8000 + 0.3875*FLN;
      MZ = (int) (floor(ZM) + 1);
      ZMIN = 1.0*MZ;
      ZDMY = Z;
      ZINC = 0.0;
      if (Z >= ZMIN) goto e20;
      ZINC = ZMIN - 1.0*NZ;
      ZDMY = Z + ZINC;
e20:  ZP = 1.0/ZDMY;
      T1 = CF[1]*ZP;
      S = T1;
      if (ZP < WDTOL) goto e40;
      ZSQ = ZP*ZP;
      TST = T1*WDTOL;
      for (K=2; K<=22; K++) {
        ZP = ZP*ZSQ;
        TRM = CF[K]*ZP;
        if (ABS(TRM) < TST) goto e40;
        S += TRM;
      }
e40:  if (ZINC != 0.0) goto e50;
      TLG = log(Z);
	  vmfree(vmblock);
      return (Z*(TLG-1.0) + 0.5*(CON-TLG) + S);
e50:  ZP = 1.0;
      NZ = (int) floor(ZINC);
      for (I=1; I<=NZ; I++)  ZP *= (Z+1.0*(I-1));
      TLG = log(ZDMY);
	  vmfree(vmblock);
      return (ZDMY*(TLG-1.0) - log(ZP) + 0.5*(CON-TLG) + S);

e70:  *IERR=1;
      vmfree(vmblock);
      return 0;
} //DGAMLN()


void ZUNIK(REAL ZRR, REAL ZRI, REAL FNU, int IKFLG, int IPMTR, REAL TOL,
           int INIT, REAL *PHIR, REAL *PHII, REAL *ZETA1R, REAL *ZETA1I, REAL *ZETA2R, 
		   REAL *ZETA2I, REAL *SUMR, REAL *SUMI, REAL *CWRKR, REAL *CWRKI)  {
/***BEGIN PROLOGUE  ZUNIK
!***REFER TO  ZBESI,ZBESK
!
!        ZUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC
!        EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2
!        RESPECTIVELY BY:
!
!        W(FNU,ZR) = PHI*EXP(ZETA)*SUM
!
!        WHERE       ZETA=-ZETA1 + ZETA2       OR
!                          ZETA1 - ZETA2
!
!        THE FIRST CALL MUST HAVE INIT=0. SUBSEQUENT CALLS WITH THE
!        SAME ZR AND FNU WILL RETURN THE I OR K FUNCTION ON IKFLG=
!        1 OR 2 WITH NO CHANGE IN INIT. CWRK IS A COMPLEX WORK
!        ARRAY. IPMTR=0 COMPUTES ALL PARAMETERS. IPMTR=1 COMPUTES PHI,
!        ZETA1,ZETA2.
!
!***ROUTINES CALLED  ZDIV,ZLOG,ZSQRT
!***END PROLOGUE  ZUNIK
!     COMPLEX CFN,CON,CONE,CRFN,CWRK,CZERO,P,HI,S,SR,SUM,T,T2,ZETA1,
!     ZETA2,ZN,ZR */
//Labels: e30, e40, e60

	  REAL *C;
      void *vmblock = NULL;

      REAL AC, CONEI, CONER, CRFNI, CRFNR, RFN, SI, SR, SRI, SRR, STI, STR,
      TEST, TI, TR, T2I, T2R, ZEROI, ZEROR, ZNI, ZNR;
      int I, IDUM, J, K, L;
      REAL CON[3]; //index 0 not used

//*** First executable statement ZUNIK

      ZEROR=0.0; ZEROI= 0.0; CONER=1.0; CONEI=0.0;
	  CON[0]=0.0;
      CON[1]=3.98942280401432678e-01;
      CON[2]=1.25331413731550025;

      //Initialize table C
      vmblock = vminit();  
      C = (REAL *) vmalloc(vmblock, VEKTOR,  121, 0); //index 0 not used
      if (! vmcomplete(vmblock)) {
        LogError ("No Memory", 0, __FILE__, __LINE__);
        return;
	  }
      C[1]=  1.00000000000000000e+00; C[2]= -2.08333333333333333e-01;
      C[3]=  1.25000000000000000e-01; C[4]=  3.34201388888888889e-01;
      C[5]= -4.01041666666666667e-01; C[6]=  7.03125000000000000e-02;
      C[7]= -1.02581259645061728e+00; C[8]=  1.84646267361111111e+00;
      C[9]= -8.91210937500000000e-01; C[10]= 7.32421875000000000e-02;
      C[11]= 4.66958442342624743e+00; C[12]=-1.12070026162229938e+01;
      C[13]= 8.78912353515625000e+00; C[14]=-2.36408691406250000e+00;
      C[15]= 1.12152099609375000e-01; C[16]=-2.82120725582002449e+01;
      C[17]= 8.46362176746007346e+01; C[18]=-9.18182415432400174e+01;
      C[19]= 4.25349987453884549e+01; C[20]=-7.36879435947963170e+00;
      C[21]= 2.27108001708984375e-01; C[22]= 2.12570130039217123e+02;
      C[23]=-7.65252468141181642e+02; C[24]= 1.05999045252799988e+03;

      C[25]=-6.99579627376132541e+02; C[26]= 2.18190511744211590e+02;
      C[27]=-2.64914304869515555e+01; C[28]= 5.72501420974731445e-01;
      C[29]=-1.91945766231840700e+03; C[30]= 8.06172218173730938e+03;
      C[31]=-1.35865500064341374e+04; C[32]= 1.16553933368645332e+04;
      C[33]=-5.30564697861340311e+03; C[34]= 1.20090291321635246e+03;
      C[35]=-1.08090919788394656e+02; C[36]= 1.72772750258445740e+00;
      C[37]= 2.02042913309661486e+04; C[38]=-9.69805983886375135e+04;