cbess00_cpp.txt

BESSEL PROGRAMS IN C/C

      CSI = P2I;
      for (I=1; I<=K; I++) {
        A1 = FKS - FK;
        AK = (FKS+FK)/(A1+FHS);
        RAK = 2.0/(FK+CONER);
        CBR = (FK+ZR)*RAK;
        CBI = ZI*RAK;
        PTR = P2R;
        PTI = P2I;
        P2R = (PTR*CBR-PTI*CBI-P1R)*AK;
        P2I = (PTI*CBR+PTR*CBI-P1I)*AK;
        P1R = PTR;
        P1I = PTI;
        CSR = CSR + P2R;
        CSI = CSI + P2I;
        FKS = A1 - FK + CONER;
        FK = FK - CONER;
      }
/*----------------------------------------------------------------------
!     COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER
!     SCALING
!---------------------------------------------------------------------*/
      TM = ZABS(CSR,CSI);
      PTR = 1.0/TM;
      S1R = P2R*PTR;
      S1I = P2I*PTR;
      CSR = CSR*PTR;
      CSI = -CSI*PTR;
      ZMLT(COEFR, COEFI, S1R, S1I, &STR, &STI);
      ZMLT(STR, STI, CSR, CSI, &S1R, &S1I);
      if (INU > 0 || N > 1) goto e200;
      ZDR = ZR;
      ZDI = ZI;
      if (IFLAG == 1) goto e270;
      goto e240;
/*----------------------------------------------------------------------
!     COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING
!---------------------------------------------------------------------*/
e200: TM = ZABS(P2R,P2I);
      PTR = 1.0/TM;
      P1R = P1R*PTR;
      P1I = P1I*PTR;
      P2R = P2R*PTR;
      P2I = -P2I*PTR;
      ZMLT(P1R, P1I, P2R, P2I, &PTR, &PTI);
      STR = DNU + 0.5 - PTR;
      STI = -PTI;
      ZDIV(STR, STI, ZR, ZI, &STR, &STI);
      STR = STR + 1.0;
      ZMLT(STR, STI, S1R, S1I, &S2R, &S2I);
/*----------------------------------------------------------------------
!     FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH
!     SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
!---------------------------------------------------------------------*/
e210: STR = DNU + 1.0;
      CKR = STR*RZR;
      CKI = STR*RZI;
      if (N == 1)  INU -= 1;
      if (INU > 0) goto e220;
      if (N > 1) goto e215;
      S1R = S2R;
      S1I = S2I;
e215: ZDR = ZR;
      ZDI = ZI;
      if (IFLAG == 1) goto e270;
      goto e240;
e220: INUB = 1;
      if (IFLAG == 1) goto e261;
e225: P1R = CSRR[KFLAG];
      ASCLE = BRY[KFLAG];
      for (I=INUB; I<=INU; I++) {
        STR = S2R;
        STI = S2I;
        S2R = CKR*STR - CKI*STI + S1R;
        S2I = CKR*STI + CKI*STR + S1I;
        S1R = STR;
        S1I = STI;
        CKR = CKR + RZR;
        CKI = CKI + RZI;
        if (KFLAG >= 3) goto e230;
        P2R = S2R*P1R;
        P2I = S2I*P1R;
        STR = ABS(P2R);
        STI = ABS(P2I);
        P2M = DMAX(STR,STI);
        if (P2M <= ASCLE) goto e230;
        KFLAG++;
        ASCLE = BRY[KFLAG];
        S1R = S1R*P1R;
        S1I = S1I*P1R;
        S2R = P2R;
        S2I = P2I;
        STR = CSSR[KFLAG];
        S1R = S1R*STR;
        S1I = S1I*STR;
        S2R = S2R*STR;
        S2I = S2I*STR;
        P1R = CSRR[KFLAG];
e230:;}
      if (N != 1) goto e240;
      S1R = S2R;
      S1I = S2I;
e240: STR = CSRR[KFLAG];
      YR[1] = S1R*STR;
      YI[1] = S1I*STR;
      if (N == 1) {
		vmfree(vmblock);  
	    return;
      }
      YR[2] = S2R*STR;
      YI[2] = S2I*STR;
      if (N == 2) {
		vmfree(vmblock);  
	    return;
      }
      KK = 2;
e250: KK++;
      if (KK > N) {
		vmfree(vmblock);  
	    return;
      }
      P1R = CSRR[KFLAG];
      ASCLE = BRY[KFLAG];
      for (I=KK; I<=N; I++) {
        P2R = S2R;
        P2I = S2I;
        S2R = CKR*P2R - CKI*P2I + S1R;
        S2I = CKI*P2R + CKR*P2I + S1I;
        S1R = P2R;
        S1I = P2I;
        CKR = CKR + RZR;
        CKI = CKI + RZI;
        P2R = S2R*P1R;
        P2I = S2I*P1R;
        YR[I] = P2R;
        YI[I] = P2I;
        if (KFLAG >= 3) goto e260;
        STR = ABS(P2R);
        STI = ABS(P2I);
        P2M = DMAX(STR,STI);
        if (P2M <= ASCLE) goto e260;
        KFLAG = KFLAG + 1;
        ASCLE = BRY[KFLAG];
        S1R = S1R*P1R;
        S1I = S1I*P1R;
        S2R = P2R;
        S2I = P2I;
        STR = CSSR[KFLAG];
        S1R = S1R*STR;
        S1I = S1I*STR;
        S2R = S2R*STR;
        S2I = S2I*STR;
        P1R = CSRR[KFLAG];
e260:;} //I loop
	  vmfree(vmblock);
      return;
/*----------------------------------------------------------------------
!     IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
!---------------------------------------------------------------------*/
e261: HELIM = 0.5*ELIM;
      ELM = EXP(-ELIM);
      CELMR = ELM;
      ASCLE = BRY[1];
      ZDR = ZR;
      ZDI = ZI;
      IC = -1;
      J = 2;
      for (I=1; I<=INU; I++) {
        STR = S2R;
        STI = S2I;
        S2R = STR*CKR-STI*CKI+S1R;
        S2I = STI*CKR+STR*CKI+S1I;
        S1R = STR;
        S1I = STI;
        CKR = CKR+RZR;
        CKI = CKI+RZI;
        AS = ZABS(S2R,S2I);
        ALAS = log(AS);
        P2R = -ZDR+ALAS;
        if (P2R < -ELIM) goto e263;
        ZLOG(S2R,S2I,&STR,&STI,&IDUM);
        P2R = -ZDR+STR;
        P2I = -ZDI+STI;
        P2M = EXP(P2R)/TOL;
        P1R = P2M*COS(P2I);
        P1I = P2M*SIN(P2I);
        ZUCHK(P1R,P1I,&NW,ASCLE,TOL);
        if (NW != 0) goto e263;
        J = 3 - J;
        CYR[J] = P1R;
        CYI[J] = P1I;
        if (IC == I-1) goto e264;
        IC = I;
        goto e262;
e263:   if (ALAS < HELIM) goto e262;
        ZDR = ZDR-ELIM;
        S1R = S1R*CELMR;
        S1I = S1I*CELMR;
        S2R = S2R*CELMR;
        S2I = S2I*CELMR;
e262:;} // I loop
      if (N != 1) goto e270;
      S1R = S2R;
      S1I = S2I;
      goto e270;
e264: KFLAG = 1;
      INUB = I+1;
      S2R = CYR[J];
      S2I = CYI[J];
      J = 3 - J;
      S1R = CYR[J];
      S1I = CYI[J];
      if (INUB <= INU) goto e225;
      if (N != 1) goto e240;
      S1R = S2R;
      S1I = S2I;
      goto e240;
e270: YR[1] = S1R;
      YI[1] = S1I;
      if (N == 1) goto e280;
      YR[2] = S2R;
      YI[2] = S2I;
e280: ASCLE = BRY[1];
      ZKSCL(ZDR,ZDI,FNU,N,YR,YI,NZ,&RZR,&RZI,ASCLE,TOL,ELIM);
      INU = N - (*NZ);
      if (INU <= 0) {
		vmfree(vmblock);  
	    return;
      }
      KK = *NZ + 1;
      S1R = YR[KK];
      S1I = YI[KK];
      YR[KK] = S1R*CSRR[1];
      YI[KK] = S1I*CSRR[1];
      if (INU == 1) {
		vmfree(vmblock);  
	    return;
      }
      KK = *NZ + 2;
      S2R = YR[KK];
      S2I = YI[KK];
      YR[KK] = S2R*CSRR[1];
      YI[KK] = S2I*CSRR[1];
      if (INU == 2) {
		vmfree(vmblock);  
	    return;
      }
      T2 = FNU + 1.0*(KK-1);
      CKR = T2*RZR;
      CKI = T2*RZI;
      KFLAG = 1;
      goto e250;
/*----------------------------------------------------------------------
!     SCALE BY DEXP(Z), IFLAG = 1 CASES
!---------------------------------------------------------------------*/
e290: KODED = 2;
      IFLAG = 1;
      KFLAG = 2;
      goto e120;
/*----------------------------------------------------------------------
!     FNU=HALF ODD INTEGER CASE, DNU=-0.5
!---------------------------------------------------------------------*/
e300: S1R = COEFR;
      S1I = COEFI;
      S2R = COEFR;
      S2I = COEFI;
      goto e210;
e310: *NZ=-2;
	  vmfree(vmblock);  
} //ZBKNU()


void ZKSCL(REAL ZRR, REAL ZRI, REAL FNU, int N, REAL *YR, REAL *YI, int *NZ,
           REAL *RZR, REAL *RZI, REAL ASCLE, REAL TOL, REAL ELIM)  {
/***BEGIN PROLOGUE  ZKSCL
!***REFER TO  ZBESK
!
!     SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
!     ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
!     RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
!
!***ROUTINES CALLED  ZUCHK,ZABS,ZLOG
!***END PROLOGUE  ZKSCL
!     COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM */
//Labels: e10, e20, e25, e30, e40, e45

      REAL ACS, AS, CKI, CKR, CSI, CSR, FN, STR, S1I, S1R, S2I, S2R,
      ZEROI, ZEROR, ZDR, ZDI, CELMR, ELM, HELIM, ALAS;
      int I, IC, IDUM, KK, NN, NW;
      REAL *CYR, *CYI;
      void *vmblock = NULL;

//    Initialize 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;
	  }
	  
      ZEROR=0.0; ZEROI=0.0;
      *NZ = 0;
      IC = 0;
      NN = IMIN(2,N);
      for (I=1; I<=NN; I++) {
        S1R = YR[I];
        S1I = YI[I];
        CYR[I] = S1R;
        CYI[I] = S1I;
        AS = ZABS(S1R,S1I);
        ACS = -ZRR + log(AS);
        *NZ = *NZ + 1;
        YR[I] = ZEROR;
        YI[I] = ZEROI;
        if (ACS < -ELIM) goto e10;
        ZLOG(S1R, S1I, &CSR, &CSI, &IDUM);
        CSR = CSR - ZRR;
        CSI = CSI - ZRI;
        STR = EXP(CSR)/TOL;
        CSR = STR*COS(CSI);
        CSI = STR*SIN(CSI);
        ZUCHK(CSR, CSI, &NW, ASCLE, TOL);
        if (NW != 0) goto e10;
        YR[I] = CSR;
        YI[I] = CSI;
        IC = I;
        *NZ = *NZ - 1;
e10: ;}
      if (N == 1) {
		vmfree(vmblock);  
	    return;
      }
      if (IC > 1) goto e20;
      YR[1] = ZEROR;
      YI[1] = ZEROI;
      *NZ = 2;
e20:  if (N == 2) {
		vmfree(vmblock);  
	    return;
      }
      if (*NZ == 0) {
		vmfree(vmblock);  
	    return;
      }
      FN = FNU + 1.0;
      CKR = FN*(*RZR);
      CKI = FN*(*RZI);
      S1R = CYR[1];
      S1I = CYI[1];
      S2R = CYR[2];
      S2I = CYI[2];
      HELIM = 0.5*ELIM;
      ELM = EXP(-ELIM);
      CELMR = ELM;
      ZDR = ZRR;
      ZDI = ZRI;

//    FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
//    S2 GETS LARGER THAN EXP(ELIM/2).

      for (I=3; I<=N; I++) {
        KK = I;
        CSR = S2R;
        CSI = S2I;
        S2R = CKR*CSR - CKI*CSI + S1R;
        S2I = CKI*CSR + CKR*CSI + S1I;
        S1R = CSR;
        S1I = CSI;
        CKR = CKR + (*RZR);
        CKI = CKI + (*RZI);
        AS = ZABS(S2R,S2I);
        ALAS = log(AS);
        ACS = -ZDR + ALAS;
        *NZ = *NZ + 1;
        YR[I] = ZEROR;
        YI[I] = ZEROI;
        if (ACS < -ELIM) goto e25;
        ZLOG(S2R, S2I, &CSR, &CSI, &IDUM);
        CSR = CSR - ZDR;
        CSI = CSI - ZDI;
        STR = EXP(CSR)/TOL;
        CSR = STR*COS(CSI);
        CSI = STR*SIN(CSI);
        ZUCHK(CSR, CSI, &NW, ASCLE, TOL);
        if (NW != 0) goto e25;
        YR[I] = CSR;
        YI[I] = CSI;
        *NZ = *NZ - 1;
        if (IC == KK-1) goto e40;
        IC = KK;
        goto e30;
e25:    if (ALAS < HELIM) goto e30;
        ZDR = ZDR - ELIM;
        S1R = S1R*CELMR;
        S1I = S1I*CELMR;
        S2R = S2R*CELMR;
        S2I = S2I*CELMR;
e30:  ;}
      *NZ = N;
      if (IC == N)  *NZ = N-1;
      goto e45;
e40:  *NZ = KK - 2;
e45:  for (I=1; I<=*NZ; I++) {
        YR[I] = ZEROR;
        YI[I] = ZEROI;
      }
	  vmfree(vmblock);  
} //ZKSCL()

// end of file Cbess00.cpp