!---------------------------------------------------------------------*/
        APHI = ZABS(PHIDR,PHIDI);
        RS1 = RS1 + log(APHI);
        if (ABS(RS1) > ELIM) goto e260;
        if (KDFLG == 1) IFLAG = 1;
        if (RS1 < 0.0) goto e220;
        if (KDFLG == 1) IFLAG = 3;
e220:   STR = PHIDR*SUMDR - PHIDI*SUMDI;
        STI = PHIDR*SUMDI + PHIDI*SUMDR;
        S2R = -CSGNI*STI;
        S2I = CSGNI*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 e230;
        ZUCHK(S2R, S2I, &NW, BRY[1], TOL);
        if (NW == 0) goto e230;
        S2R = ZEROR;
        S2I = ZEROI;
e230:   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 e250;
        ZS1S2(&ZRR, &ZRI, &S1R, &S1I, &S2R, &S2I, &NW, ASC, ALIM, &IUF);
        *NZ = *NZ + NW;
e250:   YR[KK] = S1R*CSPNR - S1I*CSPNI + S2R;
        YI[KK] = CSPNR*S1I + CSPNI*S1R + S2I;
        KK = KK - 1;
        CSPNR = -CSPNR;
        CSPNI = -CSPNI;
        if (C2R != 0.0  || C2I != 0.0)  goto e255;
        KDFLG = 1;
        goto e270;
e255:   if (KDFLG == 2) goto e275;
        KDFLG = 2;
        goto e270;
e260:   if (RS1 > 0.0) goto e300;
        S2R = ZEROR;
        S2I = ZEROI;
        goto e230;
e270:;} // K loop
      K = N;
e275: 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 e280;
        ZS1S2(&ZRR, &ZRI, &C1R, &C1I, &C2R, &C2I, &NW, ASC, ALIM, &IUF);
        *NZ = *NZ + NW;
e280:   YR[KK] = C1R*CSPNR - C1I*CSPNI + C2R;
        YI[KK] = C1R*CSPNI + C1I*CSPNR + C2I;
        KK--;
        CSPNR = -CSPNR;
        CSPNI = -CSPNI;
        if (IFLAG >= 3) goto e290;
        C2R = ABS(CKR);
        C2I = ABS(CKI);
        C2M = DMAX(C2R,C2I);
        if (C2M <= ASCLE) goto e290;
        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];
e290:;} // I loop
	  vmfree(vmblock);
      return;
e300: *NZ = -1;
      vmfree(vmblock);
} //ZUNK1()


void ZUNK2(REAL ZR, REAL ZI, REAL FNU, int KODE, int MR, int N, REAL *YR, 
           REAL *YI, int *NZ, REAL TOL, REAL ELIM, REAL ALIM)  {
/***BEGIN PROLOGUE  ZUNK2
!***REFER TO  ZBESK
!
!     ZUNK2 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
!     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
!     UNIFORM ASYMPTOTIC EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN)
!     WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR
!     -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT
!     HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC-
!     ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
!     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
!
!***ROUTINES CALLED  ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,ZABS
!***END PROLOGUE  ZUNK2
!     COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC,
!     CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ,
!     S1,S2,Y,Z,ZB,ZETA1,ZETA1D,ZETA2,ZETA2D,ZN,ZR */
//Labels: e10,e20,e30,e40,e50,e60,e70,e80,e85,e90,e100,e105,e120,e130,e172,
//	      e175,e180,e190,e210,e220,e230,e240,e250,e255,e270,e280,e290,e295,
//	      e300,e310,e320

      REAL AARG, AIC, AII, AIR, ANG, APHI, ARGDI, ARGDR, ASC, ASCLE, ASUMDI,
      ASUMDR, BSUMDI, BSUMDR,CAR, CKI, CKR, CONEI, CONER, CRSC, CR1I, CR1R,
      CR2I, CR2R, CSCL, CSGNI, CSI, CSPNI, CSPNR, CSR, C1I, C1R, C2I, C2M,
      C2R, DAII, DAIR, FMR, FN, FNF, HPI, PHIDI, PHIDR, PTI, PTR, RAST,
      RAZR, RS1, RZI, RZR, SAR, SGN, STI, STR, S1I, S1R, S2I, S2R, YY,
      ZBI, ZBR, ZEROI, ZEROR, ZET1DI, ZET1DR, ZET2DI, ZET2DR, ZNI, ZNR,
      ZRI, ZRR;
      int I, IB, IFLAG, IFN, IL, IN0, INU, IUF, K, KDFLG, KFLAG, KK, NAI,
      NDAI, NW, IDUM, J, IPARD, IC;
      REAL *BRY, *ASUMR, *ASUMI, *BSUMR, *BSUMI, *PHIR, *PHII, *ARGR, *ARGI, *ZETA1R,
      *ZETA1I, *ZETA2R, *ZETA2I, *CYR, *CYI, *CIPR, *CIPI, *CSSR, *CSRR;

	  void *vmblock = NULL;

//*** First executable statement ZUNK2

      //initialize pointers to vectors
      vmblock = vminit();  
      BRY   = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0); //index 0 not used
      ASUMR = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      ASUMI = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      BSUMR = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      BSUMI = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      PHIR  = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      PHII  = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      ZETA1R = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      ZETA1I = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      ZETA2R = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      ZETA2I = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      CYR  = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      CYI  = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      CIPR = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      CIPI = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      ARGR = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      ARGI = (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;
	  } 	

      printf(" ZUNK2:allocations Ok.\n"); getchar();

      ZEROR=0.0; ZEROI=0.0; CONER=1.0; CONEI=0.0;
      CR1R= 1.0;  CR1I= 1.73205080756887729;
      CR2R=-0.5; CR2I=-8.66025403784438647E-01;
      HPI=1.57079632679489662; 
      AIC=1.26551212348464539;
 
      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;

      KDFLG = 1;
      *NZ = 0;
/*----------------------------------------------------------------------
!     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
!     THE UNDERFLOW LIMIT
!---------------------------------------------------------------------*/
      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.0*D1MACH(1)/TOL;
      BRY[2] = 1.0/BRY[1];
      BRY[3] = D1MACH(2);
      ZRR = ZR;
      ZRI = ZI;
      if (ZR >= 0.0) goto e10;
      ZRR = -ZR;
      ZRI = -ZI;
e10:  YY = ZRI;
      ZNR = ZRI;
      ZNI = -ZRR;
      ZBR = ZRR;
      ZBI = ZRI;
      INU = (int) floor(FNU);
      FNF = FNU - 1.0*INU;
      ANG = -HPI*FNF;
      CAR = COS(ANG);
      SAR = SIN(ANG);
      C2R = HPI*SAR;
      C2I = -HPI*CAR;
      KK = (INU % 4) + 1;
      STR = C2R*CIPR[KK] - C2I*CIPI[KK];
      STI = C2R*CIPI[KK] + C2I*CIPR[KK];
      CSR = CR1R*STR - CR1I*STI;
      CSI = CR1R*STI + CR1I*STR;
      if (YY > 0.0) goto e20;
      ZNR = -ZNR;
      ZBI = -ZBI;
/*----------------------------------------------------------------------
!     K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
!     QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
!     CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
!---------------------------------------------------------------------*/
e20:  J = 2;
      for (I=1; I<=N; I++) {
/*----------------------------------------------------------------------
!     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
!---------------------------------------------------------------------*/
        J = 3 - J;
        FN = FNU + 1.0*(I-1);
        ZUNHJ(ZNR, ZNI, FN, 0, TOL, &PHIR[J], &PHII[J], &ARGR[J], &ARGI[J],
              &ZETA1R[J], &ZETA1I[J], &ZETA2R[J], &ZETA2I[J], &ASUMR[J],
              &ASUMI[J], &BSUMR[J], &BSUMI[J]);
        if (KODE == 1) goto e30;
        STR = ZBR + ZETA2R[J];
        STI = ZBI + ZETA2I[J];
        RAST = FN/ZABS(STR,STI);
        STR = STR*RAST*RAST;
        STI = -STI*RAST*RAST;
        S1R = ZETA1R[J] - STR;
        S1I = ZETA1I[J] - STI;
        goto e40;
e30:    S1R = ZETA1R[J] - ZETA2R[J];
        S1I = ZETA1I[J] - ZETA2I[J];
/*----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!---------------------------------------------------------------------*/
e40:    RS1 = S1R;
        if (ABS(RS1) > ELIM) goto e70;
        if (KDFLG == 1)  KFLAG = 2;
        if (ABS(RS1) < ALIM) goto e50;
/*----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!---------------------------------------------------------------------*/
        APHI = ZABS(PHIR[J],PHII[J]);
        AARG = ZABS(ARGR[J],ARGI[J]);
        RS1 = RS1 + log(APHI) - 0.25*log(AARG) - AIC;
        if (ABS(RS1) > ELIM) goto e70;
        if (KDFLG == 1)  KFLAG = 1;
        if (RS1 < 0.0) goto e50;
        if (KDFLG == 1)  KFLAG = 3;
/*----------------------------------------------------------------------
!     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
!     EXPONENT EXTREMES
!---------------------------------------------------------------------*/
e50:    C2R = ARGR[J]*CR2R - ARGI[J]*CR2I;
        C2I = ARGR[J]*CR2I + ARGI[J]*CR2R;
        ZAIRY(C2R, C2I, 0, 2, &AIR, &AII, &NAI, &IDUM);
        ZAIRY(C2R, C2I, 1, 2, &DAIR, &DAII, &NDAI, &IDUM);
        STR = DAIR*BSUMR[J] - DAII*BSUMI[J];
        STI = DAIR*BSUMI[J] + DAII*BSUMR[J];
        PTR = STR*CR2R - STI*CR2I;
        PTI = STR*CR2I + STI*CR2R;
        STR = PTR + (AIR*ASUMR[J]-AII*ASUMI[J]);
        STI = PTI + (AIR*ASUMI[J]+AII*ASUMR[J]);
        PTR = STR*PHIR[J] - STI*PHII[J];
        PTI = STR*PHII[J] + STI*PHIR[J];
        S2R = PTR*CSR - PTI*CSI;
        S2I = PTR*CSI + PTI*CSR;
        STR = EXP(S1R)*CSSR[KFLAG];
        S1R = STR*COS(S1I);
        S1I = STR*SIN(S1I);
        STR = S2R*S1R - S2I*S1I;
        S2I = S1R*S2I + S2R*S1I;
        S2R = STR;
        if (KFLAG != 1) goto e60;
        ZUCHK(S2R, S2I, &NW, BRY[1], TOL);
        if (NW != 0) goto e70;
e60:    if (YY <= 0.0)  S2I = -S2I;
        CYR[KDFLG] = S2R;
        CYI[KDFLG] = S2I;
        YR[I] = S2R*CSRR[KFLAG];
        YI[I] = S2I*CSRR[KFLAG];
        STR = CSI;
        CSI = -CSR;
        CSR = STR;
        if (KDFLG == 2) goto e85;
        KDFLG = 2;
        goto e80;
e70:    if (RS1 > 0.0) goto e320;
/*----------------------------------------------------------------------
!     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!---------------------------------------------------------------------*/
        if (ZR < 0.0) goto e320;
        KDFLG = 1;
        YR[I]=ZEROR;
        YI[I]=ZEROI;
        *NZ=*NZ+1;
        STR = CSI;
        CSI =-CSR;
        CSR = STR;
        if (I == 1) goto e80;
        if (YR[I-1] == ZEROR && YI[I-1] == ZEROI) goto e80;
        YR[I-1]=ZEROR;
        YI[I-1]=ZEROI;
        *NZ = *NZ + 1;
e80: ;} // I loop
      I = N;
e85:  RAZR = 1.0/ZABS(ZRR,ZRI);
      STR = ZRR*RAZR;
      STI = -ZRI*RAZR;
      RZR = (STR+STR)*RAZR;
      RZI = (STI+STI)*RAZR;
      CKR = FN*RZR;
      CKI = FN*RZI;
      IB = I + 1;
      if (N < IB) goto e180;
/*----------------------------------------------------------------------
!     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
!     ON UNDERFLOW.
!---------------------------------------------------------------------*/
      FN = FNU + 1.0*(N-1);
      IPARD = 1;
      if (MR != 0)  IPARD = 0;
      ZUNHJ(ZNR, ZNI, FN, IPARD, TOL, &PHIDR, &PHIDI, &ARGDR, &ARGDI,
            &ZET1DR, &ZET1DI, &ZET2DR, &ZET2DI, &ASUMDR, &ASUMDI, &BSUMDR, &BSUMDI);
      if (KODE == 1) goto e90;
      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 e100;
e90:  S1R = ZET1DR - ZET2DR;
      S1I = ZET1DI - ZET2DI;
e100: RS1 = S1R;
      if (ABS(RS1) > ELIM) goto e105;
      if (ABS(RS1) < ALIM) goto e120;
/*----------------------------------------------------------------------
!     REFINE ESTIMATE AND TEST
!---------------------------------------------------------------------*/
      APHI = ZABS(PHIDR,PHIDI);
      RS1 = RS1+log(APHI);
      if (ABS(RS1) < ELIM) goto e120;
e105: if (RS1 > 0.0) goto e320;
/*----------------------------------------------------------------------
!     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!---------------------------------------------------------------------*/
      if (ZR < 0.0) goto e320;
      *NZ = N;
      for (I=1; I<=N; I++) {
        YR[I] = ZEROR;
        YI[I] = ZEROI;
      }
      vmfree(vmblock);
      return;
e120: S1R = CYR[1];
      S1I = CYI[1];
      S2R = CYR[2];
      S2I = CYI[2];
      C1R = CSRR[KFLAG];
      ASCLE = BRY[KFLAG];
      for (I=IB; I<=N; I++) {
        C2R = S2R;
        C2I = S2I;
        S2R = CKR*C2R - CKI*C2I + S1R;
        S2I = CKR*C2I + CKI*C2R + S1I;
        S1R = C2R;
        S1I = C2I;
        CKR = CKR + RZR;
        CKI = CKI + RZI;
        C2R = S2R*C1R;
        C2I = S2I*C1R;
        YR[I] = C2R;
        YI[I] = C2I;
        if (KFLAG >= 3) goto e130;
        STR = ABS(C2R);
        STI = ABS(C2I);
        C2M = DMAX(STR,STI);
        if (C2M <= ASCLE) goto e130;
        KFLAG = KFLAG + 1;
        ASCLE = BRY[KFLAG];
        S1R = S1R*C1R;
        S1I = S1I*C1R;
        S2R = C2R;
        S2I = C2I;
        S1R = S1R*CSSR[KFLAG];
        S1I = S1I*CSSR[KFLAG];
        S2R = S2R*CSSR[KFLAG];
        S2I = S2I*CSSR[KFLAG];
        C1R = CSRR[KFLAG];
e130:;} // I loop
e180: if (MR == 0) {
        vmfree(vmblock);
	    return;
      }
/*----------------------------------------------------------------------
!     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
!---------------------------------------------------------------------*/