*NZ = 0;
      FMR = 1.0*MR;
      SGN = -SIGN(PI,FMR);
/*----------------------------------------------------------------------
!     CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
!---------------------------------------------------------------------*/
      CSGNI = SGN;
      if (YY <= 0.0)  CSGNI = -CSGNI;
      IFN = INU + N - 1;
      ANG = FNF*SGN;
      CSPNR = COS(ANG);
      CSPNI = SIN(ANG);
      if ((IFN % 2) == 0) goto e190;
      CSPNR = -CSPNR;
      CSPNI = -CSPNI;
/*----------------------------------------------------------------------
!     CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS
!     COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
!     QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
!     CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
!---------------------------------------------------------------------*/
e190: CSR = SAR*CSGNI;
      CSI = CAR*CSGNI;
      IN0 = (IFN % 4) + 1;
      C2R = CIPR[IN0];
      C2I = CIPI[IN0];
      STR = CSR*C2R + CSI*C2I;
      CSI = -CSR*C2I + CSI*C2R;
      CSR = STR;
      ASC = BRY[1];
      IUF = 0;
      KK = N;
      KDFLG = 1;
      IB = IB - 1;
      IC = IB - 1;
      for (K=1; K<=N; K++) {
        FN = FNU + 1.0*(KK-1);
/*----------------------------------------------------------------------
!     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
!     FUNCTION ABOVE
!---------------------------------------------------------------------*/
        if (N > 2) goto e175;
e172:   PHIDR = PHIR[J];
        PHIDI = PHII[J];
        ARGDR = ARGR[J];
        ARGDI = ARGI[J];
        ZET1DR = ZETA1R[J];
        ZET1DI = ZETA1I[J];
        ZET2DR = ZETA2R[J];
        ZET2DI = ZETA2I[J];
        ASUMDR = ASUMR[J];
        ASUMDI = ASUMI[J];
        BSUMDR = BSUMR[J];
        BSUMDI = BSUMI[J];
        J = 3 - J;
        goto e210;
e175:   if (KK == N && IB < N) goto e210;
        if (KK == IB || KK == IC) goto e172;
        ZUNHJ(ZNR, ZNI, FN, 0, TOL, &PHIDR, &PHIDI, &ARGDR,
              &ARGDI, &ZET1DR, &ZET1DI, &ZET2DR, &ZET2DI, &ASUMDR,
              &ASUMDI, &BSUMDR, &BSUMDI);
e210:   if (KODE == 1) goto e220;
        STR = ZBR + ZET2DR;
        STI = ZBI + ZET2DI;
        RAST = FN/ZABS(STR,STI);
        STR = STR*RAST*RAST;
        STI = -STI*RAST*RAST;
        S1R = -ZET1DR + STR;
        S1I = -ZET1DI + STI;
        goto e230;
e220:   S1R = -ZET1DR + ZET2DR;
        S1I = -ZET1DI + ZET2DI;
/*----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!---------------------------------------------------------------------*/
e230:   RS1 = S1R;
        if (ABS(RS1) > ELIM) goto e280;
        if (KDFLG == 1)  IFLAG = 2;
        if (ABS(RS1) < ALIM) goto e240;
/*----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!---------------------------------------------------------------------*/
        APHI = ZABS(PHIDR,PHIDI);
        AARG = ZABS(ARGDR,ARGDI);
        RS1 = RS1 + log(APHI) - 0.25*log(AARG) - AIC;
        if (ABS(RS1) > ELIM) goto e280;
        if (KDFLG == 1)  IFLAG = 1;
        if (RS1 < 0.0) goto e240;
        if (KDFLG == 1)  IFLAG = 3;
e240:   ZAIRY(ARGDR, ARGDI, 0, 2, &AIR, &AII, &NAI, &IDUM);
        ZAIRY(ARGDR, ARGDI, 1, 2, &DAIR, &DAII, &NDAI, &IDUM);
        STR = DAIR*BSUMDR - DAII*BSUMDI;
        STI = DAIR*BSUMDI + DAII*BSUMDR;
        STR = STR + (AIR*ASUMDR-AII*ASUMDI);
        STI = STI + (AIR*ASUMDI+AII*ASUMDR);
        PTR = STR*PHIDR - STI*PHIDI;
        PTI = STR*PHIDI + STI*PHIDR;
        S2R = PTR*CSR - PTI*CSI;
        S2I = PTR*CSI + PTI*CSR;
        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 e250;
        ZUCHK(S2R, S2I, &NW, BRY[1], TOL);
        if (NW == 0) goto e250;
        S2R = ZEROR;
        S2I = ZEROI;
e250:   if (YY <= 0.0)  S2I = -S2I;
        CYR[KDFLG] = S2R;
        CYI[KDFLG] = S2I;
        C2R = S2R;
        C2I = S2I;
        S2R = S2R*CSRR[IFLAG];
        S2I = S2I*CSRR[IFLAG];
/*----------------------------------------------------------------------
!     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
!---------------------------------------------------------------------*/
        S1R = YR[KK];
        S1I = YI[KK];
        if (KODE == 1) goto e270;
        ZS1S2(&ZRR, &ZRI, &S1R, &S1I, &S2R, &S2I, &NW, ASC, ALIM, &IUF);
        *NZ = *NZ + NW;
e270:   YR[KK] = S1R*CSPNR - S1I*CSPNI + S2R;
        YI[KK] = S1R*CSPNI + S1I*CSPNR + S2I;
        KK = KK - 1;
        CSPNR = -CSPNR;
        CSPNI = -CSPNI;
        STR = CSI;
        CSI = -CSR;
        CSR = STR;
        if (C2R != 0.0 || C2I != 0.0) goto e255;
        KDFLG = 1;
        goto e290;
e255:   if (KDFLG == 2) goto e295;
        KDFLG = 2;
        goto e290;
e280:   if (RS1 > 0.0) goto e320;
        S2R = ZEROR;
        S2I = ZEROI;
        goto e250;
e290:;} // K loop 
      K = N;
e295: IL = N - K;
      if (IL == 0) { 
		vmfree(vmblock);  
	    return;
      }
/*----------------------------------------------------------------------
!     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
!     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
!     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
!---------------------------------------------------------------------*/
      S1R = CYR[1];
      S1I = CYI[1];
      S2R = CYR[2];
      S2I = CYI[2];
      CSR = CSRR[IFLAG];
      ASCLE = BRY[IFLAG];
      FN = 1.0*(INU+IL);
      for (I=1; I<=IL; I++) {
        C2R = S2R;
        C2I = S2I;
        S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I);
        S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R);
        S1R = C2R;
        S1I = C2I;
        FN = FN - 1.0;
        C2R = S2R*CSR;
        C2I = S2I*CSR;
        CKR = C2R;
        CKI = C2I;
        C1R = YR[KK];
        C1I = YI[KK];
        if (KODE == 1) goto e300;
        ZS1S2(&ZRR, &ZRI, &C1R, &C1I, &C2R, &C2I, &NW, ASC, ALIM, &IUF);
        *NZ = *NZ + NW;
e300:   YR[KK] = C1R*CSPNR - C1I*CSPNI + C2R;
        YI[KK] = C1R*CSPNI + C1I*CSPNR + C2I;
        KK--;
        CSPNR = -CSPNR;
        CSPNI = -CSPNI;
        if (IFLAG >= 3) goto e310;
        C2R = ABS(CKR);
        C2I = ABS(CKI);
        C2M = DMAX(C2R,C2I);
        if (C2M <= ASCLE) goto e310;
        IFLAG++;
        ASCLE = BRY[IFLAG];
        S1R = S1R*CSR;
        S1I = S1I*CSR;
        S2R = CKR;
        S2I = CKI;
        S1R = S1R*CSSR[IFLAG];
        S1I = S1I*CSSR[IFLAG];
        S2R = S2R*CSSR[IFLAG];
        S2I = S2I*CSSR[IFLAG];
        CSR = CSRR[IFLAG];
e310:;} // I loop
	  vmfree(vmblock);
      return;
e320: *NZ = -1;
} //ZUNK2()


void ZACON(REAL ZR, REAL ZI, REAL FNU, int KODE, int MR, int N, REAL *YR, REAL *YI, int *NZ,
           REAL RL, REAL FNUL, REAL TOL, REAL ELIM, REAL ALIM)  {
/***BEGIN PROLOGUE  ZACON
!***REFER TO  ZBESK,ZBESH
!
!     ZACON APPLIES THE ANALYTIC CONTINUATION FORMULA
!
!         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
!                 MP=PI*MR*CMPLX(0.0,1.0)
!
!     TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
!     HALF Z PLANE
!
!***ROUTINES CALLED  ZBINU,ZBKNU,ZS1S2,D1MACH,ZABS,ZMLT
!***END PROLOGUE  ZACON
!     COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST,
!     S1,S2,Y,Z,ZN */
//Labels: e10,e20,e30,e40,e50,e60,e70,e80,e90

      REAL ARG, ASCLE, AS2, AZN, BSCLE, CKI, CKR, CONEI, CONER, CPN, CSCL, CSCR,
      CSGNI, CSGNR, CSPNI, CSPNR, CSR, C1I, C1M, C1R, C2I, C2R, FMR, FN,
      PTI, PTR, RAZN, RZI, RZR, SC1I, SC1R, SC2I, SC2R, SGN, SPN, STI, STR,
      S1I, S1R, S2I, S2R, YY, ZEROI, ZEROR, ZNI, ZNR;
      int I, INU, IUF, KFLAG, NN, NW;
      REAL *CYR, *CYI, *CSSR, *CSRR, *BRY;
      void *vmblock = NULL;

//*** First executable statement ZACON

      //initialize pointers to vectors
      vmblock = vminit();  
      BRY   = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0); //index 0 not used
      CYR  = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      CYI  = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      CSSR = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0);
      CSRR = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0);

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

      ZEROR=0.0; ZEROI=0.0; CONER=1.0; CONEI=0.0;
      NZ = 0;
      ZNR = -ZR;
      ZNI = -ZI;
      NN = N;
      ZBINU(ZNR, ZNI, FNU, KODE, NN, YR, YI, &NW, RL, FNUL, TOL, ELIM, ALIM);
      if (NW < 0) goto e90;
/*----------------------------------------------------------------------
!     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
!---------------------------------------------------------------------*/
      NN = IMIN(2,N);
      ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, &NW, TOL, ELIM, ALIM);
      if (NW != 0) goto e90;
      S1R = CYR[1];
      S1I = CYI[1];
      FMR = 1.0*MR;
      SGN = -SIGN(PI,FMR);
      CSGNR = ZEROR;
      CSGNI = SGN;
      if (KODE == 1) goto e10;
      YY = -ZNI;
      CPN = COS(YY);
      SPN = SIN(YY);
      ZMLT(CSGNR, CSGNI, CPN, SPN, &CSGNR, &CSGNI);
/*----------------------------------------------------------------------
!     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
!     WHEN FNU IS LARGE
!---------------------------------------------------------------------*/
e10:  INU = (int) floor(FNU);
      ARG = (FNU-1.0*INU)*SGN;
      CPN = COS(ARG);
      SPN = SIN(ARG);
      CSPNR = CPN;
      CSPNI = SPN;
      if ((INU % 2) == 0) goto e20;
      CSPNR = -CSPNR;
      CSPNI = -CSPNI;
e20:  IUF = 0;
      C1R = S1R;
      C1I = S1I;
      C2R = YR[1];
      C2I = YI[1];
      ASCLE = 1000.0*D1MACH(1)/TOL;
      if (KODE == 1) goto e30;
      ZS1S2(&ZNR, &ZNI, &C1R, &C1I, &C2R, &C2I, &NW, ASCLE, ALIM, &IUF);
      *NZ = *NZ + NW;
      SC1R = C1R;
      SC1I = C1I;
e30:  ZMLT(CSPNR, CSPNI, C1R, C1I, &STR, &STI);
      ZMLT(CSGNR, CSGNI, C2R, C2I, &PTR, &PTI);
      YR[1] = STR + PTR;
      YI[1] = STI + PTI;
      if (N == 1) {
		vmfree(vmblock);  
	    return;
      }
      CSPNR = -CSPNR;
      CSPNI = -CSPNI;
      S2R = CYR[2];
      S2I = CYI[2];
      C1R = S2R;
      C1I = S2I;
      C2R = YR[2];
      C2I = YI[2];
      if (KODE == 1) goto e40;
      ZS1S2(&ZNR, &ZNI, &C1R, &C1I, &C2R, &C2I, &NW, ASCLE, ALIM, &IUF);
      *NZ = *NZ + NW;
      SC2R = C1R;
      SC2I = C1I;
e40:  ZMLT(CSPNR, CSPNI, C1R, C1I, &STR, &STI);
      ZMLT(CSGNR, CSGNI, C2R, C2I, &PTR, &PTI);
      YR[2] = STR + PTR;
      YI[2] = STI + PTI;
      if (N == 2) {
	    vmfree(vmblock);	  
	    return;
      }
      CSPNR = -CSPNR;
      CSPNI = -CSPNI;
      AZN = ZABS(ZNR,ZNI);
      RAZN = 1.0/AZN;
      STR = ZNR*RAZN;
      STI = -ZNI*RAZN;
      RZR = (STR+STR)*RAZN;
      RZI = (STI+STI)*RAZN;
      FN = FNU + 1.0;
      CKR = FN*RZR;
      CKI = FN*RZI;
/*---------------------------------------------------------------------
!     SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
!---------------------------------------------------------------------*/
      CSCL = 1.0/TOL;
      CSCR = TOL;
      CSSR[1] = CSCL;
      CSSR[2] = CONER;
      CSSR[3] = CSCR;
      CSRR[1] = CSCR;
      CSRR[2] = CONER;
      CSRR[3] = CSCL;
      BRY[1] = ASCLE;
      BRY[2] = 1.0/ASCLE;
      BRY[3] = D1MACH(2);
      AS2 = ZABS(S2R,S2I);
      KFLAG = 2;
      if (AS2 > BRY[1]) goto e50;
      KFLAG = 1;
      goto e60;
e50:  if (AS2 < BRY[2]) goto e60;
      KFLAG = 3;
e60:  BSCLE = BRY[KFLAG];
      S1R = S1R*CSSR[KFLAG];
      S1I = S1I*CSSR[KFLAG];
      S2R = S2R*CSSR[KFLAG];
      S2I = S2I*CSSR[KFLAG];
      CSR = CSRR[KFLAG];
      for (I=3; I<=N; I++) {
        STR = S2R;
        STI = S2I;
        S2R = CKR*STR - CKI*STI + S1R;
        S2I = CKR*STI + CKI*STR + S1I;
        S1R = STR;
        S1I = STI;
        C1R = S2R*CSR;
        C1I = S2I*CSR;
        STR = C1R;
        STI = C1I;
        C2R = YR[I];
        C2I = YI[I];
        if (KODE == 1) goto e70;
        if (IUF < 0) goto e70;
        ZS1S2(&ZNR, &ZNI, &C1R, &C1I, &C2R, &C2I, &NW, ASCLE, ALIM, &IUF);
        *NZ = *NZ + NW;
        SC1R = SC2R;
        SC1I = SC2I;
        SC2R = C1R;
        SC2I = C1I;
        if (IUF != 3) goto e70;
        IUF = -4;
        S1R = SC1R*CSSR[KFLAG];
        S1I = SC1I*CSSR[KFLAG];
        S2R = SC2R*CSSR[KFLAG];
        S2I = SC2I*CSSR[KFLAG];
        STR = SC2R;
        STI = SC2I;
e70:    PTR = CSPNR*C1R - CSPNI*C1I;
        PTI = CSPNR*C1I + CSPNI*C1R;
        YR[I] = PTR + CSGNR*C2R - CSGNI*C2I;
        YI[I] = PTI + CSGNR*C2I + CSGNI*C2R;
        CKR = CKR + RZR;
        CKI = CKI + RZI;
        CSPNR = -CSPNR;
        CSPNI = -CSPNI;
        if (KFLAG >= 3) goto e80;
        PTR = ABS(C1R);
        PTI = ABS(C1I);
        C1M = DMAX(PTR,PTI);
        if (C1M <= BSCLE) goto e80;
        KFLAG++;
        BSCLE = BRY[KFLAG];
        S1R = S1R*CSR;
        S1I = S1I*CSR;
        S2R = STR;
        S2I = STI;
        S1R = S1R*CSSR[KFLAG];
        S1I = S1I*CSSR[KFLAG];
        S2R = S2R*CSSR[KFLAG];
        S2I = S2I*CSSR[KFLAG];
        CSR = CSRR[KFLAG];
e80: ;} // I loop
	  vmfree(vmblock);
      return;
e90:  *NZ = -1;
      if (NW = -2) *NZ=-2;
	  vmfree(vmblock);
} //ZACON()

//end of file Cbess3.cpp