/*************************************************************************
*    Procedures and Functions used By programs TZBESJ, TZBESK, TZBESY    *
*    (Evalute Bessel Functions with complex argument, 1st to 3rd kind)   *
* ---------------------------------------------------------------------- *
* Reference:  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES, 1983.       *
*                                                                        *
*                         C++ Release By J-P Moreau, Paris (07/10/2005). *
*************************************************************************/
#include <basis.h>
#include <vmblock.h>

#include "complex.h"

  //Headers of functions used below
  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);

  void ZUNK1(REAL ZR, REAL ZI, REAL FNU, int KODE, int MR, int N, REAL *YR, 
             REAL *YI, int *NZ, REAL TOL, REAL ELIM, REAL ALIM);

  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);

  void ZUNIK(REAL, REAL, REAL, int, int, REAL, int, REAL *, REAL *, REAL *, REAL *, 
	         REAL *, REAL *, REAL *, REAL *, REAL *, REAL *);

  void ZUCHK(REAL, REAL, int *, REAL, REAL);

  void ZS1S2(REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, int *, REAL, REAL, int *);

  void ZUNHJ(REAL, REAL, REAL, int,  REAL, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, 
	         REAL *, REAL *, REAL *, REAL *, REAL *, REAL *);

  void ZAIRY(REAL, REAL, int, int, REAL *, REAL *, int *, int *);

  void ZBINU(REAL, REAL, REAL, int, int, REAL *, REAL *, int *, REAL, REAL, REAL, 
	         REAL, REAL); 

  void ZBKNU(REAL, REAL, REAL, int, int, REAL *, REAL *, int *, REAL, 
	         REAL, REAL);


void ZBUNK(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  ZBUNK
!***REFER TO  ZBESK,ZBESH
!
!     ZBUNK COMPUTES THE K BESSEL FUNCTION FOR FNU.GT.FNUL.
!     ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR K(FNU,Z)
!     IN ZUNK1 AND THE EXPANSION FOR H(2,FNU,Z) IN ZUNK2
!
!***ROUTINES CALLED  ZUNK1,ZUNK2
!***END PROLOGUE  ZBUNK
!     COMPLEX Y,Z */
//Label: e10

      REAL AX, AY;

      *NZ = 0;
      AX = ABS(ZR)*1.7321;
      AY = ABS(ZI);
      if (AY > AX) goto e10;
/*----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR K(FNU,Z) FOR LARGE FNU APPLIED IN
!     -PI/3 <= ARG(Z) <= PI/3
!---------------------------------------------------------------------*/
      ZUNK1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, ALIM);
      return;
/*----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR H(2,FNU,Z*EXP(M*HPI)) FOR LARGE FNU
!     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
!     AND HPI=PI/2
!---------------------------------------------------------------------*/
e10:  ZUNK2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, ALIM);
} //ZBUNK()


void ZUNK1(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  ZUNK1
!***REFER TO  ZBESK
!
!     ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
!     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
!     UNIFORM ASYMPTOTIC EXPANSION.
!     MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
!     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
!
!***ROUTINES CALLED  ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,ZABS
!***END PROLOGUE  ZUNK1
!     COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO,
!     C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR */
//Labels: e10,e20,e30,e40,e50,e60,e70,e75,e80,e90,e95,e100,e120,e160,e170,
//        e172,e175,e180,e200,e210,e220,e230,e250,e255,e260,e270,e275,e280,
//        e290,e300

      REAL ANG, APHI, ASC, ASCLE, CKI, CKR, CONEI, CONER, CRSC, CSCL, CSGNI,
      CSPNI, CSPNR, CSR, C1I, C1R, C2I, C2M, C2R, FMR, FN, FNF, PHIDI,
      PHIDR, RAST, RAZR, RS1, RZI, RZR, SGN, STI, STR, SUMDI, SUMDR,
      S1I, S1R, S2I, S2R, ZEROI, ZEROR, ZET1DI, ZET1DR, ZET2DI, ZET2DR,
      ZRI, ZRR;
      int I, IB, IFLAG, IFN, II, IL, INU, IUF, K, KDFLG, KFLAG, KK, NW, INITD,
      IC, IPARD, J, M;
      REAL *BRY, *SUMR, *SUMI, *ZETA1R, *ZETA1I, *ZETA2R, *ZETA2I, *CYR, *CYI;
      int  *INIT;          //size 0..3
      REAL **CWRKR, **CWRKI; //size 0..16, 0..3     
	  REAL *CSSR, *CSRR, *PHIR, *PHII;
      REAL *TMPR, *TMPI;   //size 0..16
      void *vmblock = NULL;

//***  First executable statement ZUNK1

      //initialize pointers to vectors
      vmblock = vminit();  
      BRY  = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0); //index 0 not used
      SUMR = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      SUMI = (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);
      INIT = (int *) vmalloc(vmblock, VEKTOR,  4, 0);
      CWRKR = (REAL **) vmalloc(vmblock, MATRIX,  17, 4);
      CWRKI = (REAL **) vmalloc(vmblock, MATRIX,  17, 4);
      TMPR = (REAL *) vmalloc(vmblock, VEKTOR,  17, 0);
      TMPI = (REAL *) vmalloc(vmblock, VEKTOR,  17, 0);
      CSSR = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0);
      CSRR = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0);
      PHIR = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      PHII = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);

      if (! vmcomplete(vmblock)) {
        LogError ("No Memory", 0, __FILE__, __LINE__);
        return;
	  }
	  
	  printf(" ZUNK1:allocations Ok.\n"); getchar();

      ZEROR=0.0; ZEROI=0.0; CONER=1.0; CONEI=0.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*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:  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);
        INIT[J] = 0;
        for (II=1; II<=16; II++) {
          TMPR[II]=CWRKR[II][J];
          TMPI[II]=CWRKI[II][J];
        }
        ZUNIK(ZRR, ZRI, FN, 2, 0, TOL, INIT[J], &PHIR[J], &PHII[J],
              &ZETA1R[J], &ZETA1I[J], &ZETA2R[J], &ZETA2I[J], &SUMR[J], &SUMI[J],
              TMPR, TMPI);
        for (II=1; II<=16; II++) {
          CWRKR[II][J]=TMPR[II];
          CWRKI[II][J]=TMPI[II];
        }
        if (KODE == 1) goto e20;
        STR = ZRR + ZETA2R[J];
        STI = ZRI + ZETA2I[J] ;
        RAST = FN/ZABS(STR,STI);
        STR = STR*RAST*RAST;
        STI = -STI*RAST*RAST;
        S1R = ZETA1R[J] - STR;
        S1I = ZETA1I[J] - STI;
        goto e30;
e20:    S1R = ZETA1R[J] - ZETA2R[J];
        S1I = ZETA1I[J] - ZETA2I[J];
e30:    RS1 = S1R;
/*----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!---------------------------------------------------------------------*/
        if (ABS(RS1) > ELIM) goto e60;
        if (KDFLG == 1)  KFLAG = 2;
        if (ABS(RS1) < ALIM)  goto e40;
/*----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!---------------------------------------------------------------------*/
        APHI = ZABS(PHIR[J],PHII[J]);
        RS1 = RS1 + log(APHI);
        if (ABS(RS1) > ELIM) goto e60;
        if (KDFLG == 1) KFLAG = 1;
        if (RS1 < 0.0)  goto e40;
        if (KDFLG == 1) KFLAG = 3;
/*----------------------------------------------------------------------
!     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
!     EXPONENT EXTREMES
!---------------------------------------------------------------------*/
e40:    S2R = PHIR[J]*SUMR[J] - PHII[J]*SUMI[J];
        S2I = PHIR[J]*SUMI[J] + PHII[J]*SUMR[J];
        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 e50;
        ZUCHK(S2R, S2I, &NW, BRY[1], TOL);
        if (NW != 0) goto e60;
e50:    CYR[KDFLG] = S2R;
        CYI[KDFLG] = S2I;
        YR[I] = S2R*CSRR[KFLAG];
        YI[I] = S2I*CSRR[KFLAG];
        if (KDFLG == 2) goto e75;
        KDFLG = 2;
        goto e70;
e60:    if (RS1 > 0.0) goto e300;
/*----------------------------------------------------------------------
!     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!---------------------------------------------------------------------*/
        if (ZR < 0.0) goto e300;
        KDFLG = 1;
        YR[I]=ZEROR;
        YI[I]=ZEROI;
        NZ=NZ+1;
        if (I == 1) goto e70;
        if (YR[I-1] == ZEROR && YI[I-1] == ZEROI) goto e70;
        YR[I-1]=ZEROR;
        YI[I-1]=ZEROI;
        *NZ=*NZ+1;
e70: ;} // I loop
      I = N;
e75:  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 e160;
/*----------------------------------------------------------------------
!     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;
      INITD = 0;
      for (II=1; II<=16; II++) {
        TMPR[II]=CWRKR[II][3];
        TMPI[II]=CWRKI[II][3];
      }
      ZUNIK(ZRR, ZRI, FN, 2, IPARD, TOL, INITD, &PHIDR, &PHIDI, &ZET1DR,
            &ZET1DI, &ZET2DR, &ZET2DI, &SUMDR, &SUMDI, TMPR, TMPI);
      for (II=1; II<=16; II++) {
        CWRKR[II][3]=TMPR[II];
        CWRKI[II][3]=TMPI[II];
      }
      if (KODE == 1) goto e80;
      STR = ZRR + ZET2DR;
      STI = ZRI + ZET2DI;
      RAST = FN/ZABS(STR,STI);
      STR = STR*RAST*RAST;
      STI = -STI*RAST*RAST;
      S1R = ZET1DR - STR;
      S1I = ZET1DI - STI;
      goto e90;
e80:  S1R = ZET1DR - ZET2DR;
      S1I = ZET1DI - ZET2DI;
e90:  RS1 = S1R;
      if (ABS(RS1) > ELIM) goto e95;
      if (ABS(RS1) < ALIM) goto e100;
/*----------------------------------------------------------------------
!     REFINE ESTIMATE AND TEST
!---------------------------------------------------------------------*/
      APHI = ZABS(PHIDR,PHIDI);
      RS1 = RS1+log(APHI);
      if (ABS(RS1) < ELIM) goto e100;
e95:  if (ABS(RS1) > 0.0) goto e300;
/*----------------------------------------------------------------------
!     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!---------------------------------------------------------------------*/
      if (ZR < 0.0) goto e300;
      *NZ = N;
      for (I=1; I<=N; I++) {
        YR[I] = ZEROR;
        YI[I] = ZEROI;
      }
	  vmfree(vmblock);  
      return;
/*----------------------------------------------------------------------
!     FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
!---------------------------------------------------------------------*/
e100: 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 e120;
        STR = ABS(C2R);
        STI = ABS(C2I);
        C2M = DMAX(STR,STI);
        if (C2M <= ASCLE) goto e120;
        KFLAG++;
        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];
e120:;} // I loop
e160: if (MR == 0) {
		vmfree(vmblock);  
	    return;
      }
/*----------------------------------------------------------------------
!     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
!---------------------------------------------------------------------*/
      *NZ = 0;
      FMR = 1.0*MR;
      SGN = -SIGN(PI,FMR);
/*----------------------------------------------------------------------
!     CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
!---------------------------------------------------------------------*/
      CSGNI = SGN;
      INU = (int) floor(FNU);
      FNF = FNU - 1.0*INU;
      IFN = INU + N - 1;
      ANG = FNF*SGN;
      CSPNR = COS(ANG);
      CSPNI = SIN(ANG);
      if ((IFN % 2) == 0) goto e170;
      CSPNR = -CSPNR;
      CSPNI = -CSPNI;
e170: 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
!---------------------------------------------------------------------*/
        M=3;
        if (N > 2) goto e175;
e172:   INITD = INIT[J];
        PHIDR = PHIR[J];
        PHIDI = PHII[J];
        ZET1DR = ZETA1R[J];
        ZET1DI = ZETA1I[J];
        ZET2DR = ZETA2R[J];
        ZET2DI = ZETA2I[J];
        SUMDR = SUMR[J];
        SUMDI = SUMI[J];
        M = J;
        J = 3 - J;
        goto e180;
e175:   if (KK == N && IB < N) goto e180;
        if (KK == IB || KK == IC) goto e172;
        INITD = 0;
e180:   for (II=1; II<=16; II++) {
          TMPR[II]=CWRKR[II][M];
          TMPI[II]=CWRKI[II][M];
        }
        ZUNIK(ZRR, ZRI, FN, 1, 0, TOL, INITD, &PHIDR, &PHIDI,
              &ZET1DR, &ZET1DI, &ZET2DR, &ZET2DI, &SUMDR, &SUMDI,
              TMPR,TMPI);
        for (II=1; II<=16; II++) {
          CWRKR[II][M]=TMPR[II];
          CWRKI[II][M]=TMPI[II];
        }
        if (KODE == 1) goto e200;
        STR = ZRR + ZET2DR;
        STI = ZRI + ZET2DI;
        RAST = FN/ZABS(STR,STI);
        STR = STR*RAST*RAST;
        STI = -STI*RAST*RAST;
        S1R = -ZET1DR + STR;
        S1I = -ZET1DI + STI;
        goto e210;
e200:   S1R = -ZET1DR + ZET2DR;
        S1I = -ZET1DI + ZET2DI;
/*----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!---------------------------------------------------------------------*/
e210:   RS1 = S1R;
        if (ABS(RS1) > ELIM) goto e260;
        if (KDFLG == 1)  IFLAG = 2;
        if (ABS(RS1) < ALIM) goto e220;
/*----------------------------------------------------------------------
!     REFINE  TEST AND SCALE