/*************************************************************************
*            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/05/2005). *
*************************************************************************/
#include <math.h>
#include <basis.h>
#include <vmblock.h>

#include "complex.h"

  //headers of functions used below
  void ZUNI1(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI,
	         int *NZ, int *NLAST, REAL FNUL, REAL TOL, REAL ELIM, REAL ALIM);

  void ZUNI2(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI,
	         int *NZ, int *NLAST, REAL FNUL, 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 ZUOIK(REAL, REAL, REAL, int, int, int, REAL *, REAL *, int *, REAL, REAL, REAL);

  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 ZBUNI(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI, 
           int *NZ, int NUI, int *NLAST, REAL FNUL, REAL TOL, REAL ELIM, REAL ALIM) {
/***BEGIN PROLOGUE  ZBUNI
!***REFER TO  ZBESI,ZBESK
!
!     ZBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z) >
!     FNUL AND FNU+N-1 < FNUL. THE ORDER IS INCREASED FROM
!     FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING
!     ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z)
!     ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2.
!
!***ROUTINES CALLED  ZUNI1,ZUNI2,ZABS,D1MACH
!***END PROLOGUE  ZBUNI
!     COMPLEX CSCL,CSCR,CY,RZ,ST,S1,S2,Y,Z */
//Labels: e10,e20,e21,e25,e30,e40,e50,e60,e70,e80,e90

      REAL AX, AY, CSCLR, CSCRR, DFNU, FNUI, GNU, RAZ, RZI, RZR, STI, STR,
      S1I, S1R, S2I, S2R, ASCLE, C1R, C1I, C1M;
      int I, IFLAG, IFORM, K, NL, NW;
      REAL BRY[4];

      *NZ = 0;
      AX = ABS(ZR)*1.7321;
      AY = ABS(ZI);
      IFORM = 1;
      if (AY > AX) IFORM = 2;
      if (NUI == 0) goto e60;
      FNUI = 1.0*NUI;
      DFNU = FNU + 1.0*(N-1);
      GNU = DFNU + FNUI;
      if (IFORM == 2) goto e10;
/*----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
!     -PI/3 <= ARG(Z) <= PI/3
!---------------------------------------------------------------------*/
      ZUNI1(ZR, ZI, GNU, KODE, 2, YR, YI, &NW, NLAST, FNUL, TOL, ELIM, ALIM);
      goto e20;
/*----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
!     APPLIED IN PI/3 < ABS(ARG(Z)) <= PI/2 WHERE M=+I OR -I
!     AND HPI=PI/2
!---------------------------------------------------------------------*/
e10:  ZUNI2(ZR, ZI, GNU, KODE, 2, YR, YI, &NW, NLAST, FNUL, TOL, ELIM, ALIM);
e20:  if (NW < 0) goto e50;
      if (NW != 0) goto e90;
      STR = ZABS(YR[1],YI[1]);
/*---------------------------------------------------------------------
!     SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
!--------------------------------------------------------------------*/
      BRY[1]=1000*D1MACH(1)/TOL;
      BRY[2] = 1.0/BRY[1];
      BRY[3] = BRY[2];
      IFLAG = 2;
      ASCLE = BRY[2];
      CSCLR = 1.0;
      if (STR > BRY[1]) goto e21;
      IFLAG = 1;
      ASCLE = BRY[1];
      CSCLR = 1.0/TOL;
      goto e25;
e21:  if (STR < BRY[2]) goto e25;
      IFLAG = 3;
      ASCLE=BRY[3];
      CSCLR = TOL;
e25:  CSCRR = 1.0/CSCLR;
      S1R = YR[2]*CSCLR;
      S1I = YI[2]*CSCLR;
      S2R = YR[1]*CSCLR;
      S2I = YI[1]*CSCLR;
      RAZ = 1.0/ZABS(ZR,ZI);
      STR = ZR*RAZ;
      STI = -ZI*RAZ;
      RZR = (STR+STR)*RAZ;
      RZI = (STI+STI)*RAZ;
      for (I=1; I<=NUI; I++) {
        STR = S2R;
        STI = S2I;
        S2R = (DFNU+FNUI)*(RZR*STR-RZI*STI) + S1R;
        S2I = (DFNU+FNUI)*(RZR*STI+RZI*STR) + S1I;
        S1R = STR;
        S1I = STI;
        FNUI = FNUI - 1.0;
        if (IFLAG >= 3) goto e30;
        STR = S2R*CSCRR;
        STI = S2I*CSCRR;
        C1R = ABS(STR);
        C1I = ABS(STI);
        C1M = DMAX(C1R,C1I);
        if (C1M <= ASCLE) goto e30;
        IFLAG++;
        ASCLE = BRY[IFLAG];
        S1R = S1R*CSCRR;
        S1I = S1I*CSCRR;
        S2R = STR;
        S2I = STI;
        CSCLR = CSCLR*TOL;
        CSCRR = 1.0/CSCLR;
        S1R = S1R*CSCLR;
        S1I = S1I*CSCLR;
        S2R = S2R*CSCLR;
        S2I = S2I*CSCLR;
e30: ;}
      YR[N] = S2R*CSCRR;
      YI[N] = S2I*CSCRR;
      if (N == 1) return;
      NL = N - 1;
      FNUI = 1.0*NL;
      K = NL;
      for (I=1; I<=NL; I++) {
        STR = S2R;
        STI = S2I;
        S2R = (FNU+FNUI)*(RZR*STR-RZI*STI) + S1R;
        S2I = (FNU+FNUI)*(RZR*STI+RZI*STR) + S1I;
        S1R = STR;
        S1I = STI;
        STR = S2R*CSCRR;
        STI = S2I*CSCRR;
        YR[K] = STR;
        YI[K] = STI;
        FNUI = FNUI - 1.0;
        K--;
        if (IFLAG >= 3) goto e40;
        C1R = ABS(STR);
        C1I = ABS(STI);
        C1M = DMAX(C1R,C1I);
        if (C1M <= ASCLE) goto e40;
        IFLAG++;
        ASCLE = BRY[IFLAG];
        S1R = S1R*CSCRR;
        S1I = S1I*CSCRR;
        S2R = STR;
        S2I = STI;
        CSCLR = CSCLR*TOL;
        CSCRR = 1.0/CSCLR;
        S1R = S1R*CSCLR;
        S1I = S1I*CSCLR;
        S2R = S2R*CSCLR;
        S2I = S2I*CSCLR;
e40: ;} // I loop
      return;
e50:  *NZ = -1;
      if (NW == -2) *NZ=-2;
      return;
e60:  if (IFORM == 2) goto e70;
/*----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
!     -PI/3 <= ARG(Z) <= PI/3
!---------------------------------------------------------------------*/
      ZUNI1(ZR, ZI, FNU, KODE, N, YR, YI, &NW, NLAST, FNUL, TOL, ELIM, ALIM);
      goto e80;
/*----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
!     APPLIED IN PI/3 < ABS(ARG(Z)) <= PI/2 WHERE M=+I OR -I
!     AND HPI=PI/2
!---------------------------------------------------------------------*/
e70:  ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, &NW, NLAST, FNUL, TOL, ELIM, ALIM);
e80:  if (NW < 0) goto e50;
      *NZ = NW;
      return;
e90:  *NLAST = N;
} //ZBUNI()


void ZUNI1(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI,
           int *NZ, int *NLAST, REAL FNUL, REAL TOL, REAL ELIM, REAL ALIM)  {
/***BEGIN PROLOGUE  ZUNI1
!***REFER TO  ZBESI,ZBESK
!
!     ZUNI1 COMPUTES I(FNU,Z)  BY MEANS OF THE UNIFORM ASYMPTOTIC
!     EXPANSION FOR I(FNU,Z) IN -PI/3 <= ARG Z <= PI/3.
!
!     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
!     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
!     NLAST <> 0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
!     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1 < FNUL.
!     Y(I)=CZERO FOR I=NLAST+1,N.
!
!***ROUTINES CALLED  ZUCHK,ZUNIK,ZUOIK,D1MACH,ZABS
!***END PROLOGUE  ZUNI1
!     COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1,
!     S2,Y,Z,ZETA1,ZETA2 */
//Labels: e10,e20,e30,e40,e50,e60,e70,e90,e100,e110,e120,e130

      REAL APHI, ASCLE, CONEI, CONER, CRSC, CSCL, C1R, C2I, C2M, C2R, FN,
      PHII, PHIR, RAST, RS1, RZI, RZR, STI, STR, SUMI, SUMR, S1I, S1R,
      S2I, S2R, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R;
      int I, IFLAG, INIT, K, M, ND, NN, NUF, NW;
      REAL CWRKR[17], CWRKI[17];
      REAL *BRY, *CSSR, *CSRR, *CYR, *CYI;
      void *vmblock = NULL;

//    Initialize BRY, CSSR, CSRR, CYR, CYI
      vmblock = vminit(); 
	  BRY  = (REAL *) vmalloc(vmblock, VEKTOR,  4, 0);
	  CSSR = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      CSRR = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      CYR  = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 0);
      CYI  = (REAL *) vmalloc(vmblock, VEKTOR,  N+1, 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;
      ND = N;
      *NLAST = 0;
/*----------------------------------------------------------------------
!     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
!     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
!     EXP(ALIM)=EXP(ELIM)*TOL
!---------------------------------------------------------------------*/
      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;
/*----------------------------------------------------------------------
!     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
!---------------------------------------------------------------------*/
      FN = DMAX(FNU,1.0);
      INIT = 0;
      ZUNIK(ZR, ZI, FN, 1, 1, TOL, INIT, &PHIR, &PHII, &ZETA1R,
            &ZETA1I, &ZETA2R, &ZETA2I, &SUMR, &SUMI, CWRKR, CWRKI);
      if (KODE == 1) goto e10;
      STR = ZR + ZETA2R;
      STI = ZI + ZETA2I;
      RAST = FN/ZABS(STR,STI);
      STR = STR*RAST*RAST;
      STI = -STI*RAST*RAST;
      S1R = -ZETA1R + STR;
      S1I = -ZETA1I + STI;
      goto e20;
e10:  S1R = -ZETA1R + ZETA2R;
      S1I = -ZETA1I + ZETA2I;
e20:  RS1 = S1R;
      if (ABS(RS1) > ELIM) goto e130;
e30:  NN = IMIN(2,ND);
      for (I=1; I<=NN; I++) {
        FN = FNU + 1.0*(ND-I);
        INIT = 0;
        ZUNIK(ZR, ZI, FN, 1, 0, TOL, INIT, &PHIR, &PHII, &ZETA1R,
              &ZETA1I, &ZETA2R, &ZETA2I, &SUMR, &SUMI, CWRKR, CWRKI);
        if (KODE == 1) goto e40;
        STR = ZR + ZETA2R;
        STI = ZI + ZETA2I;
        RAST = FN/ZABS(STR,STI);
        STR = STR*RAST*RAST;
        STI = -STI*RAST*RAST;
        S1R = -ZETA1R + STR;
        S1I = -ZETA1I + STI + ZI;
        goto e50;
e40:    S1R = -ZETA1R + ZETA2R;
        S1I = -ZETA1I + ZETA2I;
/*----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!---------------------------------------------------------------------*/
e50:    RS1 = S1R;
        if (ABS(RS1) > ELIM) goto e110;
        if (I == 1)  IFLAG = 2;
        if (ABS(RS1) < ALIM) goto e60;
/*----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!---------------------------------------------------------------------*/
        APHI = ZABS(PHIR,PHII);
        RS1 = RS1 + log(APHI);
        if (ABS(RS1) > ELIM) goto e110;
        if (I == 1)  IFLAG = 1;
        if (RS1 < 0.0) goto e60;
        if (I == 1)  IFLAG = 3;
/*----------------------------------------------------------------------
!     SCALE S1 IF CABS(S1) < ASCLE
!---------------------------------------------------------------------*/
e60:    S2R = PHIR*SUMR - PHII*SUMI;
        S2I = PHIR*SUMI + PHII*SUMR;
        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 e70;
        ZUCHK(S2R, S2I, &NW, BRY[1], TOL);
        if (NW != 0) goto e110;
e70:    CYR[I] = S2R;
        CYI[I] = S2I;
        M = ND - I + 1;
        YR[M] = S2R*CSRR[IFLAG];
        YI[M] = S2I*CSRR[IFLAG];
      } // i loop
      if (ND <= 2) goto e100;
      RAST = 1.0/ZABS(ZR,ZI);
      STR = ZR*RAST;
      STI = -ZI*RAST;
      RZR = (STR+STR)*RAST;
      RZI = (STI+STI)*RAST;
      BRY[2] = 1.0/BRY[1];
      BRY[3] = D1MACH(2);
      S1R = CYR[1];