cbess2_cpp.txt

BESSEL PROGRAMS IN C/C

      S1I = CYI[1];
      S2R = CYR[2];
      S2I = CYI[2];
      C1R = CSRR[IFLAG];
      ASCLE = BRY[IFLAG];
      K = ND - 2;
      FN = 1.0*K;
      for (I=3; I<=ND; I++) {
        C2R = S2R;
        C2I = S2I;
        S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I);
        S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R);
        S1R = C2R;
        S1I = C2I;
        C2R = S2R*C1R;
        C2I = S2I*C1R;
        YR[K] = C2R;
        YI[K] = C2I;
        K = K - 1;
        FN = FN - 1.0;
        if (IFLAG >= 3) goto e90;
        STR = ABS(C2R);
        STI = ABS(C2I);
        C2M = DMAX(STR,STI);
        if (C2M <= ASCLE) goto e90;
        IFLAG++;
        ASCLE = BRY[IFLAG];
        S1R = S1R*C1R;
        S1I = S1I*C1R;
        S2R = C2R;
        S2I = C2I;
        S1R = S1R*CSSR[IFLAG];
        S1I = S1I*CSSR[IFLAG];
        S2R = S2R*CSSR[IFLAG];
        S2I = S2I*CSSR[IFLAG];
        C1R = CSRR[IFLAG];
e90: ;} // i loop
e100: vmfree(vmblock);
	  return;
/*----------------------------------------------------------------------
!     SET UNDERFLOW AND UPDATE PARAMETERS
!---------------------------------------------------------------------*/
e110: if (RS1 > 0.0) goto e120;
      YR[ND] = ZEROR;
      YI[ND] = ZEROI;
      *NZ = *NZ + 1;
      ND--;
      if (ND == 0) goto e100;
      ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, &NUF, TOL, ELIM, ALIM);
      if (NUF < 0) goto e120;
      ND -= NUF;
      *NZ = *NZ + NUF;
      if (ND == 0) goto e100;
      FN = FNU + 1.0*(ND-1);
      if (FN >= FNUL) goto e30;
      *NLAST = ND;
	  vmfree(vmblock);
      return;
e120: *NZ = -1;
	  vmfree(vmblock);
      return;
e130: if (RS1 > 0.0) goto e120;
      *NZ = N;
      for (I=1; I<=N; I++) {
        YR[I] = ZEROR;
        YI[I] = ZEROI;
      }
	  vmfree(vmblock);
} //ZUNI1()


void ZUNI2(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI,
		   int *NZ, int *NLAST, REAL FNUL, REAL TOL, REAL ELIM, REAL ALIM) {
/***BEGIN PROLOGUE  ZUNI2
!***REFER TO  ZBESI,ZBESK
!
!     ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF
!     UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I
!     OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.
!
!     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
!     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
!     NLAST <> 0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
!     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1 < FNUL.
!     Y(I)=CZERO FOR I=NLAST+1,N.
!
!***ROUTINES CALLED  ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,ZABS
!***END PROLOGUE  ZUNI2
!     COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS,
!     CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN */
//Labels: e10,e20,e30,e40,e50,e60,e70,e80,e100,e110,e120,e130,e140,e150

      REAL AARG, AIC, AII, AIR, ANG, APHI, ARGI, ARGR, ASCLE, ASUMI, ASUMR,
      BSUMI, BSUMR, CIDI, CONEI, CONER, CRSC, CSCL, C1R, C2I, C2M, C2R,
      DAII, DAIR, FN, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI, RZR, STI,
      STR, S1I, S1R, S2I, S2R, ZBI, ZBR, ZEROI, ZEROR, ZETA1I, ZETA1R,
      ZETA2I, ZETA2R, ZNI, ZNR;
      int I, IFLAG, IN0, INU, J, K, NAI, ND, NDAI, NN, NUF, NW, IDUM;
      REAL *BRY, *CIPR, *CIPI, *CSSR, *CSRR, *CYR, *CYI;
	  void *vmblock = NULL;

//    Initialize BRY, ..., CYR, CYI
      vmblock = vminit(); 
	  BRY  = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0);
	  CIPR = (REAL *) vmalloc(vmblock, VEKTOR,  5, 0);
      CIPI = (REAL *) vmalloc(vmblock, VEKTOR,  5, 0);
	  CSSR = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0);
      CSRR = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0);
      CYR  = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      CYI  = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      if (! vmcomplete(vmblock)) {
        LogError ("No Memory", 0, __FILE__, __LINE__);
        return;
	  } 	 

      ZEROR=0.0; ZEROI=0.0; CONER=1.0; CONEI=0.0;
      CIPR[1]= 1.0; CIPI[1]=0.0; CIPR[2]=0.0; CIPI[2]= 1.0;
      CIPR[3]=-1.0; CIPI[3]=0.0; CIPR[4]=0.0; CIPI[4]=-1.0;

      HPI=1.57079632679489662; AIC=1.265512123484645396;

      *NZ = 0;
      ND = N;
      *NLAST = 0;
/*----------------------------------------------------------------------
!     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
!     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
!     EXP(ALIM)=EXP(ELIM)*TOL
!---------------------------------------------------------------------*/
      CSCL = 1.0/TOL;
      CRSC = TOL;
      CSSR[1] = CSCL;
      CSSR[2] = CONER;
      CSSR[3] = CRSC;
      CSRR[1] = CRSC;
      CSRR[2] = CONER;
      CSRR[3] = CSCL;
      BRY[1] = 1000*D1MACH(1)/TOL;
/*----------------------------------------------------------------------
!     ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
!---------------------------------------------------------------------*/
      ZNR = ZI;
      ZNI = -ZR;
      ZBR = ZR;
      ZBI = ZI;
      CIDI = -CONER;
      INU = (int) floor(FNU);
      ANG = HPI*(FNU-1.0*INU);
      C2R = COS(ANG);
      C2I = SIN(ANG);
      IN0 = INU + N - 1;
      IN0 = (IN0 % 4) + 1;
      STR = C2R*CIPR[IN0] - C2I*CIPI[IN0];
      C2I = C2R*CIPI[IN0] + C2I*CIPR[IN0];
      C2R = STR;
      if (ZI > 0.0) goto e10;
      ZNR = -ZNR;
      ZBI = -ZBI;
      CIDI = -CIDI;
      C2I = -C2I;
/*----------------------------------------------------------------------
!     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
!---------------------------------------------------------------------*/
e10:  FN = DMAX(FNU,1.0);
      ZUNHJ(ZNR, ZNI, FN, 1, TOL, &PHIR, &PHII, &ARGR, &ARGI, &ZETA1R,
            &ZETA1I, &ZETA2R, &ZETA2I, &ASUMR, &ASUMI, &BSUMR, &BSUMI);
      if (KODE == 1) goto e20;
      STR = ZBR + ZETA2R;
      STI = ZBI + ZETA2I;
      RAST = FN/ZABS(STR,STI);
      STR = STR*RAST*RAST;
      STI = -STI*RAST*RAST;
      S1R = -ZETA1R + STR;
      S1I = -ZETA1I + STI;
      goto e30;
e20:  S1R = -ZETA1R + ZETA2R;
      S1I = -ZETA1I + ZETA2I;
e30:  RS1 = S1R;
      if (ABS(RS1) > ELIM) goto e150;
e40:  NN = IMIN(2,ND);
      for (I=1; I<=NN; I++) {
        FN = FNU + 1.0*(ND-I);
        ZUNHJ(ZNR, ZNI, FN, 0, TOL, &PHIR, &PHII, &ARGR, &ARGI, &ZETA1R, 
			  &ZETA1I, &ZETA2R, &ZETA2I, &ASUMR, &ASUMI, &BSUMR, &BSUMI);
        if (KODE == 1) goto e50;
        STR = ZBR + ZETA2R;
        STI = ZBI + ZETA2I;
        RAST = FN/ZABS(STR,STI);
        STR = STR*RAST*RAST;
        STI = -STI*RAST*RAST;
        S1R = -ZETA1R + STR;
        S1I = -ZETA1I + STI + ABS(ZI);
        goto e60;
e50:    S1R = -ZETA1R + ZETA2R;
        S1I = -ZETA1I + ZETA2I;
/*----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!---------------------------------------------------------------------*/
e60:    RS1 = S1R;
        if (ABS(RS1) > ELIM) goto e120;
        if (I == 1)  IFLAG = 2;
        if (ABS(RS1) < ALIM) goto e70;
/*----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!---------------------------------------------------------------------*/
        APHI = ZABS(PHIR,PHII);
        AARG = ZABS(ARGR,ARGI);
        RS1 = RS1 + log(APHI) - 0.25*log(AARG) - AIC;
        if (ABS(RS1) > ELIM) goto e120;
        if (I == 1)  IFLAG = 1;
        if (RS1 < 0.0) goto e70;
        if (I == 1)  IFLAG = 3;
/*----------------------------------------------------------------------
!     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
!     EXPONENT EXTREMES
!---------------------------------------------------------------------*/
e70:    ZAIRY(ARGR, ARGI, 0, 2, &AIR, &AII, &NAI, &IDUM);
        ZAIRY(ARGR, ARGI, 1, 2, &DAIR, &DAII, &NDAI, &IDUM);
        STR = DAIR*BSUMR - DAII*BSUMI;
        STI = DAIR*BSUMI + DAII*BSUMR;
        STR = STR + (AIR*ASUMR-AII*ASUMI);
        STI = STI + (AIR*ASUMI+AII*ASUMR);
        S2R = PHIR*STR - PHII*STI;
        S2I = PHIR*STI + PHII*STR;
        STR = EXP(S1R)*CSSR[IFLAG];
        S1R = STR*COS(S1I);
        S1I = STR*SIN(S1I);
        STR = S2R*S1R - S2I*S1I;
        S2I = S2R*S1I + S2I*S1R;
        S2R = STR;
        if (IFLAG != 1) goto e80;
        ZUCHK(S2R, S2I, &NW, BRY[1], TOL);
        if (NW != 0) goto e120;
e80:    if (ZI <= 0.0)  S2I = -S2I;
        STR = S2R*C2R - S2I*C2I;
        S2I = S2R*C2I + S2I*C2R;
        S2R = STR;
        CYR[I] = S2R;
        CYI[I] = S2I;
        J = ND - I + 1;
        YR[J] = S2R*CSRR[IFLAG];
        YI[J] = S2I*CSRR[IFLAG];
        STR = -C2I*CIDI;
        C2I = C2R*CIDI;
        C2R = STR;
      } // I loop
      if (ND <= 2) goto e110;
      RAZ = 1.0/ZABS(ZR,ZI);
      STR = ZR*RAZ;
      STI = -ZI*RAZ;
      RZR = (STR+STR)*RAZ;
      RZI = (STI+STI)*RAZ;
      BRY[2] = 1.0/BRY[1];
      BRY[3] = D1MACH(2);
      S1R = CYR[1];
      S1I = CYI[1];
      S2R = CYR[2];
      S2I = CYI[2];
      C1R = CSRR[IFLAG];
      ASCLE = BRY[IFLAG];
      K = ND - 2;
      FN = 1.0*K;
      for (I=3; I<=ND; I++) {
        C2R = S2R;
        C2I = S2I;
        S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I);
        S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R);
        S1R = C2R;
        S1I = C2I;
        C2R = S2R*C1R;
        C2I = S2I*C1R;
        YR[K] = C2R;
        YI[K] = C2I;
        K = K - 1;
        FN = FN - 1.0;
        if (IFLAG >= 3) goto e100;
        STR = ABS(C2R);
        STI = ABS(C2I);
        C2M = DMAX(STR,STI);
        if (C2M <= ASCLE) goto e100;
        IFLAG = IFLAG + 1;
        ASCLE = BRY[IFLAG];
        S1R = S1R*C1R;
        S1I = S1I*C1R;
        S2R = C2R;
        S2I = C2I;
        S1R = S1R*CSSR[IFLAG];
        S1I = S1I*CSSR[IFLAG];
        S2R = S2R*CSSR[IFLAG];
        S2I = S2I*CSSR[IFLAG];
        C1R = CSRR[IFLAG];
e100:;}
e110: vmfree(vmblock);
	  return;
e120: if (RS1 > 0.0) goto e140;
/*----------------------------------------------------------------------
!     SET UNDERFLOW AND UPDATE PARAMETERS
!---------------------------------------------------------------------*/
      YR[ND] = ZEROR;
      YI[ND] = ZEROI;
      *NZ = *NZ + 1;
      ND = ND - 1;
      if (ND == 0) goto e110;
      ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, &NUF, TOL, ELIM, ALIM);
      if (NUF < 0) goto e140;
      ND -= NUF;
      *NZ = *NZ + NUF;
      if (ND == 0) goto e110;
      FN = FNU + 1.0*(ND-1);
      if (FN < FNUL) goto e130;
      FN = CIDI;
      J = NUF + 1;
      K = (J % 4) + 1;
      S1R = CIPR[K];
      S1I = CIPI[K];
      if (FN < 0.0)  S1I = -S1I;
      STR = C2R*S1R - C2I*S1I;
      C2I = C2R*S1I + C2I*S1R;
      C2R = STR;
      goto e40;
e130: *NLAST = ND;
      vmfree(vmblock);
      return;
e140: *NZ = -1;
	  vmfree(vmblock);
      return;
e150: if (RS1 < 0.0) goto e140;
      *NZ = N;
      for (I=1; I<=N; I++) {
        YR[I] = ZEROR;
        YI[I] = ZEROI;
      }
	  vmfree(vmblock);
} // ZUNI2()

//end of file cbess2.cpp