/****************************************************************
* EVALUATE A Y-BESSEL FUNCTION OF COMPLEX ARGUMENT (SECOND KIND)*
* ------------------------------------------------------------- *
* SAMPLE RUN:                                                   *
* (Evaluate Y0 to Y4 for argument Z=(1.0,2.0) ).                *
*                                                               *
* zr(0) =   1.367419                                            *
* zi(0) =   1.521507                                            *
* zr(1) =  -1.089470                                            *
* zi(1) =   1.314951                                            *
* zr(2) =  -0.751245                                            *
* zi(2) =  -0.123950                                            *
* zr(3) =   0.290153                                            *
* zi(3) =  -0.212119                                            *
* zr(4) =   0.590344                                            *
* zi(4) =  -0.826960                                            *
* NZ = 0                                                        *
* Error code: 0                                                 *
*                                                               *
* ------------------------------------------------------------- *
* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77.   *
*                                                               *
*                        C++ Release 1.0 By J-P Moreau, Paris   *      
*****************************************************************
Note: to link with CBess0,CBess00,CBess1,CBess2,CBess3,Complex.
-------------------------------------------------------------- */
#include <basis.h>
#include <vmblock.h>

#include "complex.h"
//Uses Complex, CBess0, CBess00, CBess1, CBess2, CBess3,
       Basis_r, Vmblock.

    REAL zr, zi;
    REAL *cyr, *cyi, *cwr, *cwi;
    int i, ierr, n, nz;

	void *vmblock = NULL;

    //Headers of functions used below
    void ZBESH(REAL ZR, REAL ZI, REAL FNU, int KODE, int M, int N,
               REAL *CYR, REAL *CYI, int *NZ, int *IERR);

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

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

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

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


void ZBESY(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *CYR, REAL *CYI,
           int *NZ, REAL *CWRKR, REAL *CWRKI, int *IERR)  {

/***BEGIN PROLOGUE  ZBESY
!***DATE WRITTEN   830501   (YYMMDD)    (Original Fortran Version).
!***REVISION DATE  830501   (YYMMDD)
!***CATEGORY NO.  B5K
!***KEYWORDS  Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT,
!             BESSEL FUNCTION OF SECOND KIND
!***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
!***PURPOSE  TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT
!***DESCRIPTION
!
!                      ***A DOUBLE PRECISION ROUTINE***
!
!         ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
!         BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE
!         ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
!         -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESY RETURNS THE SCALED
!         FUNCTIONS
!
!         CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z)   I = 1,...,N , Y=AIMAG(Z)
!
!         WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
!         LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
!         ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
!         (REF. 1).
!
!         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
!           ZR,ZI  - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0),
!                    -PI.LT.ARG(Z).LE.PI
!           FNU    - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0D0
!           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
!                    KODE= 1  RETURNS
!                             CY(I)=Y(FNU+I-1,Z), I=1,...,N
!                        = 2  RETURNS
!                             CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N
!                             WHERE Y=AIMAG(Z)
!           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
!           CWRKR, - DOUBLE PRECISION WORK VECTORS OF DIMENSION AT
!           CWRKI    AT LEAST N
!
!         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
!           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
!                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
!                    CY(I)=Y(FNU+I-1,Z)  OR
!                    CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y))  I=1,...,N
!                    DEPENDING ON KODE.
!           NZ     - NZ=0 , A NORMAL RETURN
!                    NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO
!                    UNDERFLOW (GENERALLY ON KODE=2)
!           IERR   - ERROR FLAG
!                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
!                    IERR=1, INPUT ERROR   - NO COMPUTATION
!                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU IS
!                            TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
!                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
!                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
!                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
!                            ACCURACY
!                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
!                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
!                            CANCE BY ARGUMENT REDUCTION
!                    IERR=5, ERROR              - NO COMPUTATION,
!                            ALGORITHM TERMINATION CONDITION NOT MET
!
!***LONG DESCRIPTION
!
!         THE COMPUTATION IS CARRIED OUT BY THE FORMULA
!
!             Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I
!
!         WHERE I^2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z)
!         AND H(2,FNU,Z) ARE CALCULATED IN CBESH.
!
!         FOR NEGATIVE ORDERS,THE FORMULA
!
!             Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
!
!         CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD
!         INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE
!         POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)*
!         SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS
!         NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A
!         LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM
!         CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS,
!         WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF
!         ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z).
!
!         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
!         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
!         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
!         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
!         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
!         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
!         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
!         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
!         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
!         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
!         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
!         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
!         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
!         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
!         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
!         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
!         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
!         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
!         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
!
!         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
!         BESSEL FUNCTION CAN BE EXPRESSED BY P*10^S WHERE P=MAX(UNIT
!         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10^S REPRE-
!         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
!         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
!         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
!         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
!         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
!         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
!         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10^K LARGER
!         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
!         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
!         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
!         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
!         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
!         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
!         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
!         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
!         OR -PI/2+P.
!
!***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
!                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
!                 COMMERCE, 1955.
!
!               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
!                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
!
!               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
!                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
!
!               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
!                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
!                 1018, MAY, 1985
!
!               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
!                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
!                 MATH. SOFTWARE, 1986
!
!***ROUTINES CALLED  ZBESH,I1MACH,D1MACH
!***END PROLOGUE  ZBESY
!
!     COMPLEX CWRK,CY,C1,C2,EX,HCI,Z */
//Labels: e60, e70, e90, e170

      REAL C1I, C1R, C2I, C2R, ELIM, EXI, EXR, EY, HCII, STI, STR, TAY;
      int I, K, K1, K2, NZ1, NZ2;

//***FIRST EXECUTABLE STATEMENT ZBESY
      *IERR = 0;
      *NZ=0;
      if (ZR == 0.0 && ZI == 0.0) *IERR=1;
      if (FNU < 0.0) *IERR=1;
      if (KODE < 1 || KODE > 2) *IERR=1;
      if (N < 1) *IERR=1;
      if (*IERR != 0) return;   //bad parameter(s)
      HCII = 0.5;
      ZBESH(ZR, ZI, FNU, KODE, 1, N, CYR, CYI, &NZ1, IERR);
      if (*IERR != 0 && *IERR != 3) goto e170;
      ZBESH(ZR, ZI, FNU, KODE, 2, N, CWRKR, CWRKI, &NZ2, IERR);
      if (*IERR != 0 && *IERR != 3) goto e170;
      *NZ = IMIN(NZ1,NZ2);
      if (KODE == 2) goto e60;
      for (I=1; I<=N; I++) {
        STR = CWRKR[I] - CYR[I];
        STI = CWRKI[I] - CYI[I];
        CYR[I] = -STI*HCII;
        CYI[I] = STR*HCII;
      }
      return;
e60:  K1 = I1MACH(15);
      K2 = I1MACH(16);
      K = IMIN(ABS(K1),ABS(K2));
/*----------------------------------------------------------------------
!     ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
!---------------------------------------------------------------------*/
      ELIM = 2.303*(K*D1MACH(5)-3.0);
      EXR = COS(ZR);
      EXI = SIN(ZR);
      EY = 0.0;
      TAY = ABS(ZI+ZI);
      if (TAY < ELIM)  EY = EXP(-TAY);
      if (ZI < 0.0)  goto e90;
      C1R = EXR*EY;
      C1I = EXI*EY;
      C2R = EXR;
      C2I = -EXI;
e70:  *NZ = 0;
      for (I=1; I<=N; I++) {
        STR = C1R*CYR[I] - C1I*CYI[I];
        STI = C1R*CYI[I] + C1I*CYR[I];
        STR = -STR + C2R*CWRKR[I] - C2I*CWRKI[I];
        STI = -STI + C2R*CWRKI[I] + C2I*CWRKR[I];
        CYR[I] = -STI*HCII;
        CYI[I] = STR*HCII;
        if (STR == 0.0 && STI == 0.0 && EY == 0.0)  *NZ = *NZ + 1;
      }
      return;
e90:  C1R = EXR;
      C1I = EXI;
      C2R = EXR*EY;
      C2I = -EXI*EY;
      goto e70;
e170: *NZ = 0;
} //ZBESY()


void ZBESH(REAL ZR, REAL ZI, REAL FNU, int KODE, int M, int N,
           REAL *CYR, REAL *CYI, int *NZ, int *IERR)  {
/***BEGIN PROLOGUE  ZBESH
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  830501   (YYMMDD)
!***CATEGORY NO.  B5K
!***KEYWORDS  H-BESSEL FUNCTIONS,BESSEL FUNCTIONS OF COMPLEX ARGUMENT,
!             BESSEL FUNCTIONS OF THIRD KIND,HANKEL FUNCTIONS
!***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
!***PURPOSE  TO COMPUTE THE H-BESSEL FUNCTIONS OF A COMPLEX ARGUMENT
!***DESCRIPTION
!
!                      ***A DOUBLE PRECISION ROUTINE***
!         ON KODE=1, ZBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
!         HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1
!         OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX
!         Z.NE.CMPLX(0.0,0.0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI.
!         ON KODE=2, ZBESH RETURNS THE SCALED HANKEL FUNCTIONS
!
!         CY(I)=EXP(-MM*Z*I)*H(M,FNU+J-1,Z)       MM=3-2*M,   I^2=-1.
!
!         WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER AND
!         LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN THE
!         NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1).
!
!         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
!           ZR,ZI  - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0),
!                    -PT.LT.ARG(Z).LE.PI
!           FNU    - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0D0
!           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
!                    KODE= 1  RETURNS
!                             CY(J)=H(M,FNU+J-1,Z),   J=1,...,N
!                        = 2  RETURNS
!                             CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M))
!                                  J=1,...,N  ,  I^2=-1
!           M      - KIND OF HANKEL FUNCTION, M=1 OR 2
!           N      - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1