cbess0_cpp.txt

BESSEL PROGRAMS IN C/C

        if (PP < TOL)  IAS = 1;
e60:    if (IBS == 1) goto e90;
        SUMBR = ZEROR;
        SUMBI = ZEROI;
        for (K=1; K<=KMAX; K++) {
          M = L2 + K;
          SUMBR = SUMBR + PR[K]*BETA[M];
          SUMBI = SUMBI + XPI[K]*BETA[M];
          if (AP[K] < ATOL) goto e80;
        }
e80:    *BSUMR = *BSUMR + SUMBR*PP;
        *BSUMI = *BSUMI + SUMBI*PP;
        if (PP < BTOL)  IBS = 1;
e90:    if (IAS == 1 && IBS == 1) goto e110;
        L1 += 30;
        L2 += 30;
      } //IS loop
e110: *ASUMR = *ASUMR + CONER;
      PP = RFNU*RFN13;
      *BSUMR = *BSUMR*PP;
      *BSUMI = *BSUMI*PP;
e120: vmfree(vmblock);
	  return;
/*----------------------------------------------------------------------
!     CABS(W2) > 0.25
!---------------------------------------------------------------------*/
e130: ZSQRT(W2R, W2I, &WR, &WI);
      if (WR < 0.0)  WR = 0.0;
      if (WI < 0.0)  WI = 0.0;
      STR = CONER + WR;
      STI = WI;
      ZDIV(STR, STI, ZBR, ZBI, &ZAR, &ZAI);
      ZLOG(ZAR, ZAI, &ZCR, &ZCI, &IDUM);
      if (ZCI < 0.0)  ZCI = 0.0;
      if (ZCI > HPI) ZCI = HPI;
      if (ZCR < 0.0) ZCR = 0.0;
      ZTHR = (ZCR-WR)*1.5;
      ZTHI = (ZCI-WI)*1.5;
      *ZETA1R = ZCR*FNU;
      *ZETA1I = ZCI*FNU;
      *ZETA2R = WR*FNU;
      *ZETA2I = WI*FNU;
      AZTH = ZABS(ZTHR,ZTHI);
      ANG = THPI;
      if (ZTHR >= 0.0 && ZTHI < 0.0) goto e140;
      ANG = HPI;
      if (ZTHR == 0.0) goto e140;
      ANG = atan(ZTHI/ZTHR);
      if (ZTHR < 0.0) ANG += GXPI;
e140: PP = pow(AZTH,EX2);
      ANG = ANG*EX2;
      ZETAR = PP*cos(ANG);
      ZETAI = PP*sin(ANG);
      if (ZETAI < 0.0)  ZETAI = 0.0;
      *ARGR = ZETAR*FN23;
      *ARGI = ZETAI*FN23;
      ZDIV(ZTHR, ZTHI, ZETAR, ZETAI, &RTZTR, &RTZTI);
      ZDIV(RTZTR, RTZTI, WR, WI, &ZAR, &ZAI);
      TZAR = ZAR + ZAR;
      TZAI = ZAI + ZAI;
      ZSQRT(TZAR, TZAI, &STR, &STI);
      *PHIR = STR*RFN13;
      *PHII = STI*RFN13;
      if (IPMTR == 1) goto e120;
      RAW = 1.0/SQRT(AW2);
      STR = WR*RAW;
      STI = -WI*RAW;
      TFNR = STR*RFNU*RAW;
      TFNI = STI*RFNU*RAW;
      RAZTH = 1.0/AZTH;
      STR = ZTHR*RAZTH;
      STI = -ZTHI*RAZTH;
      RZTHR = STR*RAZTH*RFNU;
      RZTHI = STI*RAZTH*RFNU;
      ZCR = RZTHR*AR[2];
      ZCI = RZTHI*AR[2];
      RAW2 = 1.0/AW2;
      STR = W2R*RAW2;
      STI = -W2I*RAW2;
      T2R = STR*RAW2;
      T2I = STI*RAW2;
      STR = T2R*C[2] + C[3];
      STI = T2I*C[2];
      UPR[2] = STR*TFNR - STI*TFNI;
      UXPI[2] = STR*TFNI + STI*TFNR;
      *BSUMR = UPR[2] + ZCR;
      *BSUMI = UXPI[2] + ZCI;
      *ASUMR = ZEROR;
      *ASUMI = ZEROI;
      if (RFNU < TOL) goto e220;
      PRZTHR = RZTHR;
      PRZTHI = RZTHI;
      PTFNR = TFNR;
      PTFNI = TFNI;
      UPR[1] = CONER;
      UXPI[1] = CONEI;
      PP = 1.0;
      BTOL = TOL*(ABS(*BSUMR)+ABS(*BSUMI));
      KS = 0;
      KP1 = 2;
      L = 3;
      IAS = 0;
      IBS = 0;
      LR=2;
      while (LR<=12) {
        LRP1 = LR + 1;
/*----------------------------------------------------------------------
!     COMPUTE TWO ADDITIONAL CR, DR, AND UP FOR TWO MORE TERMS IN
!     NEXT SUMA AND SUMB
!---------------------------------------------------------------------*/
        for (K=LR; K<=LRP1; K++) {
          KS++;
          KP1++;
          L++;
          ZAR = C[L];
          ZAI = ZEROI;
          for (J=2; J<=KP1; J++) {
            L++;
            STR = ZAR*T2R - T2I*ZAI + C[L];
            ZAI = ZAR*T2I + ZAI*T2R;
            ZAR = STR;
          }
          STR = PTFNR*TFNR - PTFNI*TFNI;
          PTFNI = PTFNR*TFNI + PTFNI*TFNR;
          PTFNR = STR;
          UPR[KP1] = PTFNR*ZAR - PTFNI*ZAI;
          UXPI[KP1] = PTFNI*ZAR + PTFNR*ZAI;
          CRR[KS] = PRZTHR*BR[KS+1];
          CRI[KS] = PRZTHI*BR[KS+1];
          STR = PRZTHR*RZTHR - PRZTHI*RZTHI;
          PRZTHI = PRZTHR*RZTHI + PRZTHI*RZTHR;
          PRZTHR = STR;
          DRR[KS] = PRZTHR*AR[KS+2];
          DRI[KS] = PRZTHI*AR[KS+2];
		} //K loop
        PP = PP*RFNU2;
        if (IAS == 1) goto e180;
        SUMAR = UPR[LRP1];
        SUMAI = UXPI[LRP1];
        JU = LRP1;
        for (JR=1; JR<=LR; JR++) {
          JU = JU - 1;
          SUMAR = SUMAR + CRR[JR]*UPR[JU] - CRI[JR]*UXPI[JU];
          SUMAI = SUMAI + CRR[JR]*UXPI[JU] + CRI[JR]*UPR[JU];
        }
        *ASUMR = *ASUMR + SUMAR;
        *ASUMI = *ASUMI + SUMAI;
        TEST = ABS(SUMAR) + ABS(SUMAI);
        if (PP < TOL && TEST < TOL)  IAS = 1;
e180:   if (IBS == 1) goto e200;
        SUMBR = UPR[LR+2] + UPR[LRP1]*ZCR - UXPI[LRP1]*ZCI;
        SUMBI = UXPI[LR+2] + UPR[LRP1]*ZCI + UXPI[LRP1]*ZCR;
        JU = LRP1;
        for (JR=1; JR<=LR; JR++) {
          JU = JU - 1;
          SUMBR = SUMBR + DRR[JR]*UPR[JU] - DRI[JR]*UXPI[JU];
          SUMBI = SUMBI + DRR[JR]*UXPI[JU] + DRI[JR]*UPR[JU];
        }
        *BSUMR = *BSUMR + SUMBR;
        *BSUMI = *BSUMI + SUMBI;
        TEST = ABS(SUMBR) + ABS(SUMBI);
        if (PP < BTOL && TEST < BTOL)  IBS = 1;
e200:   if (IAS == 1 && IBS == 1) goto e220;
        LR += 2;
      } //while LR
e220: *ASUMR = *ASUMR + CONER;
      STR = -(*BSUMR)*RFN13;
      STI = -(*BSUMI)*RFN13;
      ZDIV(STR, STI, RTZTR, RTZTI, BSUMR, BSUMI);
      goto e120;
} //ZUNHJ()


void ZUOIK(REAL ZR, REAL ZI, REAL FNU, int KODE, int IKFLG, int N, REAL *YR,
		   REAL *YI, int *NUF, REAL TOL, REAL ELIM, REAL ALIM)  {
/***BEGIN PROLOGUE  ZUOIK
!***REFER TO  ZBESI,ZBESK,ZBESH
!
!     ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
!     EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
!     (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
!     WHERE ALIM < ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
!     EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
!     THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
!     MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE
!     EXP(-ELIM)=SMALLEST MACHINE NUMBER*1E+03 AND EXP(-ALIM)=
!     EXP(-ELIM)/TOL.
!
!     IKFLG=1 MEANS THE I SEQUENCE IS TESTED
!          =2 MEANS THE K SEQUENCE IS TESTED
!     NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
!         =-1 MEANS AN OVERFLOW WOULD OCCUR
!     IKFLG=1 AND NUF > 0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
!             THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
!     IKFLG=2 AND NUF = N MEANS ALL Y VALUES WERE SET TO ZERO
!     IKFLG=2 AND 0 < NUF < N NOT CONSIDERED. Y MUST BE SET BY
!             ANOTHER ROUTINE.
!
!***ROUTINES CALLED  ZUCHK,ZUNHJ,ZUNIK,D1MACH,ZABS,ZLOG
!***END PROLOGUE  ZUOIK
!     COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN,ZR */

//Labels: e10,e20,e30,e40,e50,e60,e70,e80,e90,e110,e120,e130,e140,e150,e160,e170,
//	      e180,e190,e200,e210

      REAL AARG, AIC, APHI, ARGI, ARGR, ASUMI, ASUMR, ASCLE, AX, AY,
      BSUMI, BSUMR, CZI, CZR, FNN, GNN, GNU, PHII, PHIR, RCZ,
      STR, STI, SUMI, SUMR, ZBI, ZBR, ZEROI, ZEROR, ZETA1I,
      ZETA1R, ZETA2I, ZETA2R, ZNI, ZNR, ZRI, ZRR;
      int I, IDUM, IFORM, INIT, NN, NW;
      REAL *CWRKR, *CWRKI;
      void *vmblock = NULL;

      ZEROR=0.0; ZEROI=0.0;
//    Initialize CWRKR, CWRKI
      vmblock = vminit();  
      CWRKR = (REAL *) vmalloc(vmblock, VEKTOR,  17, 0); //index 0 not used
      CWRKI = (REAL *) vmalloc(vmblock, VEKTOR,  17, 0);
      if (! vmcomplete(vmblock)) {
        LogError ("No Memory", 0, __FILE__, __LINE__);
        return;
	  }
	  
      AIC=1.265512123484645396;
      *NUF = 0;
      NN = N;
      ZRR = ZR;
      ZRI = ZI;
      if (ZR >= 0.0) goto e10;
      ZRR = -ZR;
      ZRI = -ZI;
e10:  ZBR = ZRR;
      ZBI = ZRI;
      AX = ABS(ZR)*1.7321;
      AY = ABS(ZI);
      IFORM = 1;
      if (AY > AX) IFORM = 2;
      GNU = DMAX(FNU,1.0);
      if (IKFLG == 1) goto e20;
      FNN = 1.0*NN;
      GNN = FNU + FNN - 1.0;
      GNU = DMAX(GNN,FNN);
/*----------------------------------------------------------------------
!     ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
!     REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET
!     THE SIGN OF THE IMAGINARY PART CORRECT.
!---------------------------------------------------------------------*/
e20:  if (IFORM == 2) goto e30;
      INIT = 0;
      ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, &PHIR, &PHII, &ZETA1R, &ZETA1I, &ZETA2R,
            &ZETA2I, &SUMR, &SUMI, CWRKR, CWRKI);
      CZR = -ZETA1R + ZETA2R;
      CZI = -ZETA1I + ZETA2I;
      goto e50;
e30:  ZNR = ZRI;
      ZNI = -ZRR;
      if (ZI > 0.0) goto e40;
      ZNR = -ZNR;
e40:  ZUNHJ(ZNR, ZNI, GNU, 1, TOL, &PHIR, &PHII, &ARGR, &ARGI, &ZETA1R,
            &ZETA1I, &ZETA2R, &ZETA2I, &ASUMR, &ASUMI, &BSUMR, &BSUMI);
      CZR = -ZETA1R + ZETA2R;
      CZI = -ZETA1I + ZETA2I;
      AARG = ZABS(ARGR,ARGI);
e50:  if (KODE == 1) goto e60;
      CZR -= ZBR;
      CZI -= ZBI;
e60:  if (IKFLG == 1) goto e70;
      CZR = -CZR;
      CZI = -CZI;
e70:  APHI = ZABS(PHIR,PHII);
      RCZ = CZR;
/*----------------------------------------------------------------------
!     OVERFLOW TEST
!---------------------------------------------------------------------*/
      if (RCZ > ELIM) goto e210;
      if (RCZ < ALIM) goto e80;
      RCZ += log(APHI);
      if (IFORM == 2)  RCZ -= (0.25*log(AARG) - AIC);
      if (RCZ > ELIM) goto e210;
      goto e130;
/*----------------------------------------------------------------------
!     UNDERFLOW TEST
!---------------------------------------------------------------------*/
e80:  if (RCZ < -ELIM) goto e90;
      if (RCZ > -ALIM) goto e130;
      RCZ += log(APHI);
      if (IFORM == 2) RCZ -= (0.25*log(AARG) - AIC);
      if (RCZ > -ELIM) goto e110;
e90:  for (I=1; I<=NN; I++) {
        YR[I] = ZEROR;
        YI[I] = ZEROI;
      }
      *NUF = NN;
	  vmfree(vmblock);
      return;
e110: ASCLE = 1000*D1MACH(1)/TOL;
      ZLOG(PHIR, PHII, &STR, &STI, &IDUM);
      CZR = CZR + STR;
      CZI = CZI + STI;
      if (IFORM == 1) goto e120;