!Copyright (c) 2023, Kristopher G. Klein and Daniel Verscharen
!All rights reserved.
!
!This source code is licensed under the BSD-style license found in the
!LICENSE file in the root directory of this source tree.
!
!===============================================================================
!I                                                                             I
!I                              A  L  P  S                                     I
!I                     Arbitrary Linear Plasma Solver                          I
!I                                                                             I
!I                              Version 1.0                                    I
!I                                                                             I
!I  Kristopher Klein   (kgklein@arizona.edu)                                   I
!I  Daniel Verscharen  (d.verscharen@ucl.ac.uk)                                I
!I                                                                             I
!===============================================================================

module alps_nhds
!! Module including the NHDS implementation for bi-Maxwellian/cold-plasma reference cases.
!! The original NHDS code can be found under github.com/danielver02/NHDS.
  implicit none

  private :: calc_ypsilon, besselI, BESSI, BESSI0, BESSI1, WOFZ, dispfunct, calc_chi_cold
  public :: calc_chi

contains

  ! This file is part of NHDS
  ! Copyright (C) 2024 Daniel Verscharen (d.verscharen@ucl.ac.uk)
  !All rights reserved.
  !
  !Redistribution and use in source and binary forms, with or without
  !modification, are permitted provided that the following conditions are met:
  !
  !1. Redistributions of source code must retain the above copyright notice, this
  !   list of conditions and the following disclaimer.
  !2. Redistributions in binary form must reproduce the above copyright notice,
  !   this list of conditions and the following disclaimer in the documentation
  !   and/or other materials provided with the distribution.
  !
  !THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  !ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  !WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  !DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
  !ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  !(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  !LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
  !ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  !(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  !SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  !
  !The views and conclusions contained in the software and documentation are those
  !of the authors and should not be interpreted as representing official policies,
  !either expressed or implied, of the NHDS project.

  !2025-05-28: We are multiplying all of the elements of chi by kperp^2 d_p^2, to make the solution more stable.

  subroutine calc_chi(chi,chi_low,j,kz,kperp,x)
  !! Subroutine that calculates the susceptibility of species j based on NHDS.
  use alps_var, only : bMnmaxs, bMBessel_zeros, bMbetas, bMalphas,bMpdrifts
  use alps_var, only : ms, qs, ns
  use alps_var, only : kperp_norm
  implicit none

  double complex, intent(out) :: chi(3,3)
  !! Susceptibility tensor of species j.

  double complex, intent(out) :: chi_low(3,3,-1:1)
  !! Susceptibility tensor of species j.

  integer, intent(in) :: j
  !! Index for species.

  double precision, intent(in) :: kz
  !! Normalised parallel wavenumber.

  double precision, intent(in) :: kperp
  !! Normalised perpendicular wavenumber.

  double complex, intent(in) :: x
  !! Normalised complex frequency.

  double complex :: Y(3,3)
  !! Y-tensor according to Stix.

  double complex :: Y0(3,3)
  !! n=0 contribution to Y-tensor according to Stix.

  double complex :: Y1(3,3)
  !! n=+ 1 contribution to Y-tensor according to Stix.

  double complex :: Yn1(3,3)
  !! n=- 1 contribution to Y-tensor according to Stix.

  double complex :: Ynew(3,3)
  !! Iteration of Y-tensor according to Stix.

  double precision :: z
  !! Argument of dispersion function.

  double precision :: Omega
  !! Normalised gyro-frequency.

  double precision :: ell
  !! Normalised inertial length.

  double precision :: vtherm
  !! Normalised thermal speed.

  double precision :: Bessel_zero
  !! Limit of Bessel-function calculation.

  double precision :: vdrift
  !! Normalised drift speed of bi-Maxwellian/cold population.

  integer :: n
  !! Index of sum over Bessel functions.

  integer :: i
  !! Index running over tensor.

  integer :: k
  !! Index running over tensor.

  integer :: nmaxrun
  !! Running to maximum index of n for Bessel functions.

  integer :: nmax
  !! Maximum index of n for Bessel functions.

  logical :: Bessel_run
  !! Check whether maximum n has been achieved.


  ! Check if you can use the cold-plasma dispersion relation:
  if (bMbetas(j).EQ.0.d0) then

     call calc_chi_cold(chi,j,kz,kperp,x)
     !No LD, TTD, or CD contribution from a cold species
     chi_low=cmplx(0.d0,0.d0,kind(1.d0))
  else

  Omega=qs(j)/ms(j)
  ell=sqrt(ms(j)/(ns(j)*qs(j)*qs(j)))
  vtherm=sqrt(bMbetas(j)/(ns(j)*ms(j)))
  vdrift=bMpdrifts(j)/ms(j)
  nmax=bMnmaxs(j)
  Bessel_zero=bMBessel_zeros(j)

  z=0.5d0*(kperp*vtherm/Omega)*(kperp*vtherm/Omega)*bMalphas(j)

  do i=1,3
   do k=1,3
      Y(i,k)=0.d0
      Y0(i,k)=0.d0
      Y1(i,k)=0.d0
      Yn1(i,k)=0.d0
   enddo
  enddo


  !determine maximum n for Besselfunction:
  nmaxrun=nmax
  n=0
  Bessel_run=.TRUE.
  do while (Bessel_run)
    if ((n.GE.nmax).OR.(besselI(n,z).LT.Bessel_zero)) then
       nmaxrun=n
       Bessel_run=.FALSE.
    endif
    n=n+1
  enddo


  do n=-nmaxrun,nmaxrun
     call calc_ypsilon(Ynew,j,n,kz,kperp,x)
     do i=1,3
      do k=1,3
        ! Remember that the Bessel functions give I_n(z)*exp(-z), so the factor exp(-z) is absorbed.
         Y(i,k)=Y(i,k)+Ynew(i,k)
         if (n== 0) Y0(i,k)=Ynew(i,k)
         if (n== 1) Y1(i,k)=Y1(i,k)+Ynew(i,k)
         if (n==-1) Yn1(i,k)=Yn1(i,k)+Ynew(i,k)
      enddo
     enddo
  enddo

  chi(1,1)=Y(1,1)/(ell*ell)
  chi(1,2)=Y(1,2)/(ell*ell)
  chi(1,3)=Y(1,3)/(ell*ell)
  chi(2,1)=Y(2,1)/(ell*ell)
  chi(2,2)=Y(2,2)/(ell*ell)
  chi(2,3)=Y(2,3)/(ell*ell)
  chi(3,1)=Y(3,1)/(ell*ell)
  chi(3,2)=Y(3,2)/(ell*ell)
  if (kperp_norm) then
     chi(3,3)=2.d0*x*vdrift/(ell*ell*kz*vtherm*vtherm*bMalphas(j))+Y(3,3)/(ell*ell)
  else
     chi(3,3)=kperp*kperp*2.d0*x*vdrift/(ell*ell*kz*vtherm*vtherm*bMalphas(j))+Y(3,3)/(ell*ell)
  endif

  n=0
  chi_low(1,1,n)=Y0(1,1)/(ell*ell)
  chi_low(1,2,n)=Y0(1,2)/(ell*ell)
  chi_low(1,3,n)=Y0(1,3)/(ell*ell)
  chi_low(2,1,n)=Y0(2,1)/(ell*ell)
  chi_low(2,2,n)=Y0(2,2)/(ell*ell)
  chi_low(2,3,n)=Y0(2,3)/(ell*ell)
  chi_low(3,1,n)=Y0(3,1)/(ell*ell)
  chi_low(3,2,n)=Y0(3,2)/(ell*ell)
  if (kperp_norm) then
     chi_low(3,3,n)=Y0(3,3)/(ell*ell)+2.d0*x*vdrift/(ell*ell*kz*vtherm*vtherm*bMalphas(j))
  else
     chi_low(3,3,n)=Y0(3,3)/(ell*ell)+kperp*kperp*2.d0*x*vdrift/(ell*ell*kz*vtherm*vtherm*bMalphas(j))
  endif

  n=1
  chi_low(1,1,n)=Y1(1,1)/(ell*ell)
  chi_low(1,2,n)=Y1(1,2)/(ell*ell)
  chi_low(1,3,n)=Y1(1,3)/(ell*ell)
  chi_low(2,1,n)=Y1(2,1)/(ell*ell)
  chi_low(2,2,n)=Y1(2,2)/(ell*ell)
  chi_low(2,3,n)=Y1(2,3)/(ell*ell)
  chi_low(3,1,n)=Y1(3,1)/(ell*ell)
  chi_low(3,2,n)=Y1(3,2)/(ell*ell)
  chi_low(3,3,n)=Y1(3,3)/(ell*ell)

  n=-1
  chi_low(1,1,n)=Yn1(1,1)/(ell*ell)
  chi_low(1,2,n)=Yn1(1,2)/(ell*ell)
  chi_low(1,3,n)=Yn1(1,3)/(ell*ell)
  chi_low(2,1,n)=Yn1(2,1)/(ell*ell)
  chi_low(2,2,n)=Yn1(2,2)/(ell*ell)
  chi_low(2,3,n)=Yn1(2,3)/(ell*ell)
  chi_low(3,1,n)=Yn1(3,1)/(ell*ell)
  chi_low(3,2,n)=Yn1(3,2)/(ell*ell)
  chi_low(3,3,n)=Yn1(3,3)/(ell*ell)

  endif

  end subroutine







  subroutine calc_ypsilon(Y,j,n,kz,kperp,x)
  !! Calculates the Y-tensor according to Stix for a bi-Maxwelling, using the NHDS calculation.
  use alps_var, only : bMbetas, bMalphas,bMpdrifts
  use alps_var, only : ms, qs, ns
  use alps_var, only : kperp_norm
  implicit none

  double complex, parameter ::  uniti=(0.d0,1.d0)
  !! Imaginary unit.

  double complex, intent(out) :: Y(3,3)
  !! Y-tensor as defined by Stix.

  integer, intent(in) :: j
  !! Index for species.

  integer, intent(in) :: n
  !! Index of sum over Bessel functions.

  double precision, intent(in) :: kz
  !! Normalised parallel wavenumber.

  double precision, intent(in) :: kperp
  !! Normalised perpendicular wavenumber.

  double complex, intent(in) :: x
  !! Normalised complex frequency.

  double complex :: zeta
  !! Argument of dispersion function.

  double complex :: An
  !! An function according to Stix.

  double complex :: Bn
  !! Bn function according to Stix.

  double complex :: resfac
  !! Resonance factor.

  double precision :: BInz
  !! Modified Bessel function, evaluated.

  double precision :: z
  !! Argument of Bessel function.

  double precision :: zp
  !! Argument of Bessel function, multiplied by kperp^2

  double precision :: Omega
  !! Normalised gyro-frequency.

  double precision :: ell
  !! Normalised inertial length.

  double precision :: vtherm
  !! Normalised thermal speed.

  double precision :: vdrift
  !! Normalised drift speed of bi-Maxwellian/cold population.

  double precision :: dBInzdz
  !! Derivative of Bessel function.

  logical :: kpos
  !! Check whether kpar is positive.

  kpos=.TRUE.
  if (kz.LT.0.d0) kpos=.FALSE.

  Omega=qs(j)/ms(j)
  ell=sqrt(ms(j)/(ns(j)*qs(j)*qs(j)))
  vtherm=sqrt(bMbetas(j)/(ns(j)*ms(j)))
  vdrift=bMpdrifts(j)/ms(j)


  zeta=(x-kz*vdrift-1.d0*n*Omega)/(kz*vtherm)
  resfac=x-kz*vdrift-1.d0*n*Omega
  z=0.5d0*(kperp*vtherm/Omega)*(kperp*vtherm/Omega)*bMalphas(j)
  !zp=z*kperp^2 d_ref^2
  zp=0.5d0*(vtherm/Omega)*(vtherm/Omega)*bMalphas(j)

  An=(bMalphas(j)-1.d0)
  An=An+(1.d0/(kz*vtherm)) *( bMalphas(j)*resfac + 1.d0*n*Omega)*dispfunct(zeta,kpos)

  Bn=(bMalphas(j)*(x-1.d0*n*Omega)-(kz*vdrift-1.d0*n*Omega))/kz
  Bn=Bn+((x-1.d0*n*Omega)*(bMalphas(j)*resfac+1.d0*n*Omega)/(kz*kz*vtherm) )*dispfunct(zeta,kpos)

  if (n.GE.0) then
     BInz=1.d0*besselI(n,z)
     dBInzdz = 5.d-1*(besselI(n+1,z)+besselI(n-1,z))     
     !Changing definition of derivative to avoid division by z
     !dBInzdz=besselI(n+1,z)+1.d0*n*BInz/z
  else
     BInz=1.d0*besselI(-n,z)
     dBInzdz = 5.d-1*(besselI(n+1,z)+besselI(n-1,z))     
     !Changing definition of derivative to avoid division by z
     !dBInzdz=besselI(-n-1,z)+1.d0*n*BInz/z
  endif


  if (kperp_norm) then
     ! The tensor in Stix's (10-57)
     Y(1,1)=1.d0*(n*n)*BInz*An/z
     Y(1,2)=-uniti*n*(BInz-dBInzdz)*An
     Y(1,3)=kperp*n*BInz*Bn/(Omega*z)
     Y(2,1)=uniti*n*(BInz-dBInzdz)*An
     Y(2,2)=(1.d0*(n*n)*BInz/z+2.d0*z*BInz-2.d0*z*dBInzdz)*An
     Y(2,3)=uniti*kperp*(BInz-dBInzdz)*Bn/Omega
     Y(3,1)=kperp*BInz*n*Bn/(Omega*z)
     Y(3,2)=-uniti*kperp*(BInz-dBInzdz)*Bn/Omega
     Y(3,3)=2.d0*(x-1.d0*n*Omega)*BInz*Bn/(kz*vtherm*vtherm*bMalphas(j))
  else
     ! The tensor in Stix's (10-57), multiplied by kperp^2 d_ref^2
     Y(1,1)=1.d0*(n*n)*BInz*An/zp
     Y(1,2)=-uniti*n*(BInz-dBInzdz)*An*kperp*kperp
     Y(1,3)=kperp*n*BInz*Bn/(Omega*zp)
     Y(2,1)=uniti*n*(BInz-dBInzdz)*An*kperp*kperp
     Y(2,2)=(1.d0*(n*n)*BInz/zp+kperp*kperp*2.d0*z*BInz-kperp*kperp*2.d0*z*dBInzdz)*An
     Y(2,3)=uniti*kperp*kperp*kperp*(BInz-dBInzdz)*Bn/Omega
     Y(3,1)=kperp*BInz*n*Bn/(Omega*zp)
     Y(3,2)=-uniti*kperp*kperp*kperp*(BInz-dBInzdz)*Bn/Omega
     Y(3,3)=kperp*kperp*2.d0*(x-1.d0*n*Omega)*BInz*Bn/(kz*vtherm*vtherm*bMalphas(j))
  endif

  end subroutine



  subroutine calc_chi_cold(chi,j,kz,kperp,x)
  !! Subroutine that calculates the susceptibility of species j based on the cold-plasma dispersion relation based on the paper Verscharen & Chandran, ApJ 764, 88, 2013.
  use alps_var, only : bMbetas,bMpdrifts
  use alps_var, only : ms, qs, ns
  use alps_var, only : kperp_norm
  implicit none

  double complex, intent(out) :: chi(3,3)
  !! Susceptibility tensor of species j.

  double complex, parameter ::  uniti=(0.d0,1.d0)
  !! Imaginary unit.

  integer, intent(in) :: j
  !! Index for species.

  double precision, intent(in) :: kz
  !! Normalised parallel wavenumber.

  double precision, intent(in) :: kperp
  !! Normalised perpendicular wavenumber.

  double complex, intent(in) :: x
  !! Normalised complex frequency.

  double precision :: Omega
  !! Normalised gyro-frequency.

  double precision :: ell
  !! Normalised inertial length.

  double precision :: vtherm
  !! Normalised thermal speed.

  double precision :: vdrift
  !! Normalised drift speed of bi-Maxwellian/cold population.

  double complex :: dispR
  !! Susceptibility element R as defined by Verscharen & Chandran

  double complex :: dispL
  !! Susceptibility element L as defined by Verscharen & Chandran

  double complex :: dispP
  !! Susceptibility element P as defined by Verscharen & Chandran

  double complex :: dispM
  !! Susceptibility element M as defined by Verscharen & Chandran

  double complex :: dispJ
  !! Susceptibility element J as defined by Verscharen & Chandran



  Omega=qs(j)/ms(j)
  ell=sqrt(ms(j)/(ns(j)*qs(j)*qs(j)))
  vtherm=sqrt(bMbetas(j)/(ns(j)*ms(j)))
  vdrift=bMpdrifts(j)/ms(j)

  dispR=-(1.d0/(ell*ell))*(x-kz*vdrift)/(x-kz*vdrift+Omega)
  dispL=-(1.d0/(ell*ell))*(x-kz*vdrift)/(x-kz*vdrift-Omega)

  dispP=(x*x/((x-kz*vdrift)**2)) + ((kperp*vdrift)**2/((x-kz*vdrift)**2-Omega**2))
  dispP=-(1.d0/(ell*ell))*dispP

  dispJ=-(1.d0/(ell*ell))*kperp*vdrift*(x-kz*vdrift)/((x-kz*vdrift)**2-Omega**2)

  dispM=uniti*(1.d0/(ell*ell))*kperp*vdrift*Omega/((x-kz*vdrift)**2-Omega**2)


  if (kperp_norm) then
     chi(1,1)=(dispR+dispL)/2.d0
     chi(1,2)=-uniti*(dispR-dispL)/2.d0
     chi(1,3)=dispJ
     chi(2,1)=uniti*(dispR-dispL)/2.d0
     chi(2,2)=(dispR+dispL)/2.d0
     chi(2,3)=dispM
     chi(3,1)=dispJ
     chi(3,2)=-dispM
     chi(3,3)=dispP
  else
     write(*,*) 'ERROR: check kperp norm for cold plasma distribution!!'
     stop
  endif
     
  end subroutine





double precision function besselI(n,x)
  !! Calculates the modified Bessel function of argument x and order n.
  implicit none

  double precision, intent(in) :: x
  !! Argument of the modified Bessel function.

  integer, intent(in) :: n
  !! Order of the modified Bessel function.

  if (n.LT.0) then
  	besselI=BESSI(-n,x)
  else
  	besselI=BESSI(n,x)
  endif

  end function


!
! Calculate the dispersion function:
!
double complex function dispfunct(zeta,kpos)
  !! Calculates the dispersion function based on the complex error function.
  implicit none

  double complex, parameter ::  uniti=(0.d0,1.d0)
  !! Imaginary init.

  double complex, intent(in) :: zeta
  !! Argument of dispersion function.

  logical, intent(in) :: kpos
  !! Check whether kpar is positive.

  logical :: flag
  !! Check whether [[WOFZ(subroutine)]] executes successfully.

  double precision :: U
  !! Real part of the output for [[WOFZ(subroutine)]].

  double precision :: V
  !! Imaginary part of the output for [[WOFZ(subroutine)]].

  double precision :: XI
  !! Real part of the argument for [[WOFZ(subroutine)]].

  double precision :: YI
  !! Imaginary part of the argument for [[WOFZ(subroutine)]].

  double precision :: M_PI = 4.d0*atan(1.d0)
  !! Pi.

  XI=1.d0*real(zeta)
  YI=1.d0*aimag(zeta)

  if (kpos) then
     call WOFZ(XI,YI,U,V,flag)
     dispfunct = uniti*sqrt(M_PI)*(U+uniti*V)
  else
     call WOFZ(-XI,-YI,U,V,flag)
     dispfunct = -uniti*sqrt(M_PI)*(U+uniti*V)
  endif
  end function







  SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
    !!  Given a complex number Z = (XI,YI), this subroutine computes
    !!  the value of the Faddeeva function W(Z) = exp(-Z**2)*erfc(-I*Z),
    !!  where erfc is the complex complementary error function and I
    !! is the imaginary unit.

    !Based on
        !G.P.M. Poppe, C.M.J. Wijers, "More Efficient Computation of the Complex Error-Function",
        !ACM Trans. Math. Software 16, 47 (1990).

        !      ALGORITHM 680, COLLECTED ALGORITHMS FROM ACM.
        !      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
        !      VOL. 16, NO. 1, PP. 47


  !  THE ACCURACY OF THE ALGORITHM FOR Z IN THE 1ST AND 2ND QUADRANT
  !  IS 14 SIGNIFICANT DIGITS; IN THE 3RD AND 4TH IT IS 13 SIGNIFICANT
  !  DIGITS OUTSIDE A CIRCULAR REGION WITH RADIUS 0.126 AROUND A ZERO
  !  OF THE FUNCTION.
  !  ALL REAL VARIABLES IN THE PROGRAM ARE DOUBLE PRECISION.

  !  THE CODE CONTAINS A FEW COMPILER-DEPENDENT PARAMETERS :
  !     RMAXREAL = THE MAXIMUM VALUE OF RMAXREAL EQUALS THE ROOT OF
  !                RMAX = THE LARGEST NUMBER WHICH CAN STILL BE
  !                IMPLEMENTED ON THE COMPUTER IN DOUBLE PRECISION
  !                FLOATING-POINT ARITHMETIC
  !     RMAXEXP  = LN(RMAX) - LN(2)
  !     RMAXGONI = THE LARGEST POSSIBLE ARGUMENT OF A DOUBLE PRECISION
  !                GONIOMETRIC FUNCTION (DCOS, DSIN, ...)
  !  THE REASON WHY THESE PARAMETERS ARE NEEDED AS THEY ARE DEFINED WILL
  !  BE EXPLAINED IN THE CODE BY MEANS OF COMMENTS
  !
  !
  !  PARAMETER LIST
  !     XI     = REAL      PART OF Z
  !     YI     = IMAGINARY PART OF Z
  !     U      = REAL      PART OF W(Z)
  !     V      = IMAGINARY PART OF W(Z)
  !     FLAG   = AN ERROR FLAG INDICATING WHETHER OVERFLOW WILL
  !              OCCUR OR NOT; TYPE LOGICAL;
  !              THE VALUES OF THIS VARIABLE HAVE THE FOLLOWING
  !              MEANING :
  !              FLAG=.FALSE. : NO ERROR CONDITION
  !              FLAG=.TRUE.  : OVERFLOW WILL OCCUR, THE ROUTINE
  !                             BECOMES INACTIVE
  !  XI, YI      ARE THE INPUT-PARAMETERS
  !  U, V, FLAG  ARE THE OUTPUT-PARAMETERS
  !
  !  FURTHERMORE THE PARAMETER FACTOR EQUALS 2/SQRT(PI)
  !
  !  THE ROUTINE IS NOT UNDERFLOW-PROTECTED BUT ANY VARIABLE CAN BE
  !  PUT TO 0 UPON UNDERFLOW;
  !

        IMPLICIT DOUBLE PRECISION (A-H, O-Z)
        IMPLICIT INTEGER (I-N)

        LOGICAL A, B, FLAG
  PARAMETER (FACTOR=1.12837916709551257388D0,RMAXREAL = 0.5D+154,RMAXEXP  = 708.503061461606D0,&
  &RMAXGONI = 3.53711887601422D+15)

        FLAG = .FALSE.

        XABS = DABS(XI)
        YABS = DABS(YI)
        X    = XABS/6.3
        Y    = YABS/4.4

  !     THE FOLLOWING IF-STATEMENT PROTECTS
  !     QRHO = (X**2 + Y**2) AGAINST OVERFLOW
  !
        IF ((XABS.GT.RMAXREAL).OR.(YABS.GT.RMAXREAL)) GOTO 100

        QRHO = X**2 + Y**2

        XABSQ = XABS**2
        XQUAD = XABSQ - YABS**2
        YQUAD = 2*XABS*YABS

       A     = QRHO.LT.0.085264D0

        IF (A) THEN

  !  IF (QRHO.LT.0.085264D0) THEN THE FADDEEVA-FUNCTION IS EVALUATED
  !  USING A POWER-SERIES (ABRAMOWITZ/STEGUN, EQUATION (7.1.5), P.297)
  !  N IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED
  !  ACCURACY

          QRHO  = (1-0.85*Y)*DSQRT(QRHO)
          N     = IDNINT(6 + 72*QRHO)
          J     = 2*N+1
          XSUM  = 1.0/J
          YSUM  = 0.0D0
          DO 10 I=N, 1, -1
            J    = J - 2
            XAUX = (XSUM*XQUAD - YSUM*YQUAD)/I
            YSUM = (XSUM*YQUAD + YSUM*XQUAD)/I
            XSUM = XAUX + 1.0/J
   10     CONTINUE
          U1   = -FACTOR*(XSUM*YABS + YSUM*XABS) + 1.0
          V1   =  FACTOR*(XSUM*XABS - YSUM*YABS)
          DAUX =  DEXP(-XQUAD)
          U2   =  DAUX*DCOS(YQUAD)
          V2   = -DAUX*DSIN(YQUAD)

          U    = U1*U2 - V1*V2
          V    = U1*V2 + V1*U2

        ELSE

  !  IF (QRHO.GT.1.O) THEN W(Z) IS EVALUATED USING THE LAPLACE
  !  CONTINUED FRACTION
  !  NU IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED
  !  ACCURACY
  !
  !  IF ((QRHO.GT.0.085264D0).AND.(QRHO.LT.1.0)) THEN W(Z) IS EVALUATED
  !  BY A TRUNCATED TAYLOR EXPANSION, WHERE THE LAPLACE CONTINUED FRACTION
  !  IS USED TO CALCULATE THE DERIVATIVES OF W(Z)
  !  KAPN IS THE MINIMUM NUMBER OF TERMS IN THE TAYLOR EXPANSION NEEDED
  !  TO OBTAIN THE REQUIRED ACCURACY
  !  NU IS THE MINIMUM NUMBER OF TERMS OF THE CONTINUED FRACTION NEEDED
  !  TO CALCULATE THE DERIVATIVES WITH THE REQUIRED ACCURACY

          IF (QRHO.GT.1.0) THEN
            H    = 0.0D0
            KAPN = 0
            QRHO = DSQRT(QRHO)
            NU   = IDINT(3 + (1442/(26*QRHO+77)))
          ELSE
            QRHO = (1-Y)*DSQRT(1-QRHO)
            H    = 1.88*QRHO
            H2   = 2*H
            KAPN = IDNINT(7  + 34*QRHO)
            NU   = IDNINT(16 + 26*QRHO)
          ENDIF

          B = (H.GT.0.0)

          IF (B) QLAMBDA = H2**KAPN

          RX = 0.0
          RY = 0.0
          SX = 0.0
          SY = 0.0

          DO 11 N=NU, 0, -1
            NP1 = N + 1
            TX  = YABS + H + NP1*RX
            TY  = XABS - NP1*RY
            C   = 0.5/(TX**2 + TY**2)
            RX  = C*TX
            RY  = C*TY
            IF ((B).AND.(N.LE.KAPN)) THEN
              TX = QLAMBDA + SX
              SX = RX*TX - RY*SY
              SY = RY*TX + RX*SY
              QLAMBDA = QLAMBDA/H2
            ENDIF
   11     CONTINUE

          IF (H.EQ.0.0) THEN
            U = FACTOR*RX
            V = FACTOR*RY
          ELSE
            U = FACTOR*SX
            V = FACTOR*SY
          END IF

          IF (YABS.EQ.0.0) U = DEXP(-XABS**2)

        END IF

  !  EVALUATION OF W(Z) IN THE OTHER QUADRANTS

        IF (YI.LT.0.0) THEN
          IF (A) THEN
            U2    = 2*U2
            V2    = 2*V2
          ELSE
            XQUAD =  -XQUAD

  !         THE FOLLOWING IF-STATEMENT PROTECTS 2*EXP(-Z**2)
  !        AGAINST OVERFLOW

            IF ((YQUAD.GT.RMAXGONI).OR. (XQUAD.GT.RMAXEXP)) GOTO 100

            W1 =  2*DEXP(XQUAD)
            U2  =  W1*DCOS(YQUAD)
            V2  = -W1*DSIN(YQUAD)
          END IF

          U = U2 - U
          V = V2 - V
          IF (XI.GT.0.0) V = -V
        ELSE
          IF (XI.LT.0.0) V = -V
        END IF

        RETURN

    100 FLAG = .TRUE.
        RETURN

        END





  ! ----------------------------------------------------------------------
        double precision FUNCTION BESSI(N,X)
        !! Function to calculate the first kind modified Bessel function of integer order N
        !! for any real X.
        ! -------------------------------------------------------------------- *
        !   Reference: From Numath Library By Tuan Dang Trong in Fortran 77.   *
        !                                                                      *
        !                               F90 Release 1.1 By J-P Moreau, Paris.  *
        !                                                                      *
        !   Version 1.1: corected value of P4 in BESSIO (P4=1.2067492 and not  *
        !                1.2067429) Aug. 2011.                                 *

        IMPLICIT DOUBLE PRECISION (A-H, O-Z)
        IMPLICIT INTEGER (I-N)
  !     This subroutine calculates the first kind modified Bessel function
  !     of integer order N, for any REAL X. We use here the classical
  !     recursion formula, when X > N. For X < N, the Miller's algorithm
  !     is used to avoid overflows.
  !     REFERENCE:
  !     C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
  !     MATHEMATICAL TABLES, VOL.5, 1962.
  !
  ! This function calculates I_n(x)*exp(-x) instead of I_n(x)


        DOUBLE PRECISION :: X,TOX,BIM,BI,BIP
        INTEGER,PARAMETER :: IACC = 40
        INTEGER,PARAMETER :: IBIGNO = maxexponent(x)/2
        IF (N.EQ.0) THEN
        BESSI = BESSI0(X)
        RETURN
        ENDIF
        IF (N.EQ.1) THEN
        BESSI = BESSI1(X)
        RETURN
        ENDIF
        IF(X.EQ.0.D0) THEN
        BESSI=0.D0
        RETURN
        ENDIF
        TOX = 2.D0/X
        BIP = 0.D0
        BI  = 1.D0
        BESSI = 0.D0
        M = 2*((N+INT(SQRT(FLOAT(IACC*N)))))
        DO 12 J = M,1,-1
        BIM = BIP+DFLOAT(J)*TOX*BI
        BIP = BI
        BI  = BIM
        IF (EXPONENT(BI).GT.IBIGNO) THEN
        BI  = scale(BI,-IBIGNO)
        BIP = scale(BIP,-IBIGNO)
        BESSI = scale(BESSI,-IBIGNO)
        ENDIF
        IF (J.EQ.N) BESSI = BIP
     12 CONTINUE
        BESSI = BESSI0(X)*(BESSI/BI)
        RETURN
        END




        double precision FUNCTION BESSI0(X)
        !! Auxiliary Bessel functions for N=0, N=1
        ! This function calculates I_0(x)*exp(-x) instead of I_0(x)
        REAL *8 X,Y,P1,P2,P3,P4,P5,P6,P7,  &
        Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX
        DATA P1,P2,P3,P4,P5,P6,P7/1.D0,3.5156229D0,3.0899424D0,1.2067492D0,  &
        0.2659732D0,0.360768D-1,0.45813D-2/
        DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1, &
        0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1,  &
        0.2635537D-1,-0.1647633D-1,0.392377D-2/
        AX=ABS(X)
        IF(AX.LT.3.75D0) THEN
        Y=(X/3.75D0)**2
        BESSI0=(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))*EXP(-AX)
        ELSE
        Y=3.75D0/AX
        BX=1.D0/SQRT(AX)
        AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))
        BESSI0=AX*BX
        ENDIF
        RETURN
        END



        double precision FUNCTION BESSI1(X)
        !! Modified Bessel function of order 1.
        ! This function calculates I_1(x)*exp(-x) instead of I_1(x)
        REAL *8 X,Y,P1,P2,P3,P4,P5,P6,P7,  &
        Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX
        DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0,  &
        0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
        DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1, &
        -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1, &
        -0.2895312D-1,0.1787654D-1,-0.420059D-2/
        AX=ABS(X)
        IF(AX.LT.3.75D0) THEN
          Y=(X/3.75D0)**2
          BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))*EXP(-AX)
        ELSE
          Y=3.75D0/AX
          BX=1.D0/SQRT(AX)
          AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))
          BESSI1=AX*BX
        ENDIF
        RETURN
      END





end module alps_nhds