text:
! **********************************************************************************
! *          ROUTINES FOR COMPUTING THE SPHERICAL VECTOR WAVE FUNCTIONS            *
! *    -----------------------------------------------------------------------     *
! *    Partial list of subroutines:                                                *
! *      MN,                  MN_infinit,      MN_DS,          MN_poles_COMP,      *
! *      MN_poles_COMP1,      MN_DS_COMP,      MN_DS_COMP1,    MN_poles_LAY,       *
! *      MN_DS_LAY,           MN_complete,     MN_infinit_complete,                *
! *      MN_infinit_reflection_complete,       MN_left_right,  MN_anisotrop,       *
! *      MN_m_anisotrop.                                                           *
! **********************************************************************************
subroutine MN (index, z, theta, m, Nrank, Nmax, MV, NV)
!-----------------------------------------------------------------------------------
! The routine computes the localized vector spherical wave functions (excepting    !
! the factor exp(j*m*phi)) for the azimuthal mode m.                               !
!                                                                                !
! Input parameters:                                                                !
! - index (integer) - if index = 1, the regular vector wave functions are          !
!   computed. If index = 3, the radiating vector wave functions are computed.      !
! - z (complex,) - z = k * r, where k is the wave number and r is the module of    !
!   the position vector.                                                           !
! - theta (real) - zenith angle.                                                   !
! - m (integer) - azimuthal mode.                                                  !
! - Nrank (integer) - maximum expansion order.                                     !
! - Nmax (integer) - specifies the dimension of the vectors.                       !
!   For m = 0, Nmax = Nrank, while for m > 0, Nmax = Nrank - |m| + 1.              !
!                                                                                !
! Output parameters:                                                               !
! - MV, NV (complex arrays) - vector spherical wave functions.                     !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: index, m, Nrank, Nmax
  real(O)    :: theta
  complex(O) :: z, MV(3,Nmax), NV(3,Nmax)
!
  integer    :: k, n
  real(O)    :: nm, mr, nr, fp, ft, fl
  complex(O) :: zc
  real(O),allocatable    :: Pnm(:), dPnm(:), pinm(:), taunm(:)
  complex(O),allocatable :: jh(:), jhd(:)
!
  allocate (jh(0:Nrank), jhd(0:Nrank))
  allocate (Pnm(0:Nrank), dPnm(0:Nrank), pinm(0:Nrank), taunm(0:Nrank))
  zc = z
  if (abs(zc) < MachEps) zc = cmplx(MachEps,MachEps,O)
  if (index == 1) then
    call besel_j (zc, Nrank, jh, jhd)
  else if (index == 3) then
    call besel_h (zc, Nrank, jh, jhd)
  endif
  call leg_normalized (theta, abs(m), Nrank, Pnm, dPnm, pinm, taunm)
  mr = real(m,O)
  do k = 1, Nmax
    if (m == 0) then
      n = k
    else
      n = abs(m) + k - 1
    end if
    nr = real(n * (n + 1),O)
    nm = 1._O / sqrt(2._O * nr)
    fp = mr * pinm(n) * nm
    ft = taunm(n) * nm
    fl = nr * Pnm(n) * nm
    MV(1,k) =   zero
    MV(2,k) =   jh(n)  * im * fp
    MV(3,k) = - jh(n)  * ft
    NV(1,k) =   jh(n)  * fl / zc
    NV(2,k) =   jhd(n) * ft / zc
    NV(3,k) =   jhd(n) * im * fp / zc
  end do
  deallocate (jh, jhd)
  deallocate (Pnm, dPnm, pinm, taunm)
end subroutine MN
!***********************************************************************************
subroutine MN_infinit (theta, phi, m, Nrank, Nmax, Minf, Ninf)
!-----------------------------------------------------------------------------------
! The routine computes the localized vector spherical wave functions at infinity   !
! (excepting the factor exp(j*k*R)/(k*R)) for the azimuthal mode m.                !
!                                                                                !
! Input parameters:                                                                !
! - theta, phi (real variables) - polar angles.                                    !
! - m (integer) - azimuthal mode.                                                  !
! - Nrank (integer) - maximum expansion order.                                     !
! - Nmax (integer) - specifies the dimension of the vectors.                       !
!   For m = 0, Nmax = Nrank, while for m > 0, Nmax = Nrank - |m| + 1.              !
!                                                                                !
! Output parameters:                                                               !
! - Minf, Ninf (complex arrays) - vector spherical wave functions at infinity.     !
!-----------------------------------------------------------------------------------
  use parameters
  implicit none
  integer    :: m, Nrank, Nmax
  real(O)    :: theta, phi
  complex(O) :: Minf(3,Nmax), Ninf(3,Nmax)
!
  integer    :: k, n
  real(O)    :: arg, nm, mr
  complex(O) :: fact, factc, factt, factp
  real(O),allocatable :: Pnm(:), dPnm(:), pinm(:), taunm(:)
!
  allocate (Pnm(0:Nrank), dPnm(0:Nrank), pinm(0:Nrank), taunm(0:Nrank))
  call leg_normalized (theta, abs(m), Nrank, Pnm, dPnm, pinm, taunm)
  mr   = real(m,O)
  arg  = mr * phi
  fact = exp(im * arg)
  do k = 1, Nmax
    if (m == 0) then
      n = k
    else
      n = abs(m) + k - 1
    end if
    nm    = real(2 * n * (n + 1),O)
    nm    = 1._O / sqrt(nm)
    factc = fact  * (-im)**(n + 1) * nm
    factp = factc * mr * pinm(n)
    factt = factc * taunm(n)
    Minf(1,k) =   zero
    Minf(2,k) =   im * factp
    Minf(3,k) = - factt
    Ninf(1,k) =   zero
    Ninf(2,k) =   im * factt
    Ninf(3,k) = - factp
  end do
  deallocate (Pnm, dPnm, pinm, taunm)
end subroutine MN_infinit
!***********************************************************************************
subroutine MN_DS (index, k, r, theta, zRe, zIm, m, Nrank, MV, NV)
!-----------------------------------------------------------------------------------
! The routine computes the distributed vector spherical wave functions (excepting  !
! the factor exp(j*m*phi)) for the azimuthal mode m.                               !
!                                                                                !
! Input parameters:                                                                !
! - index (integer) - if index = 1, the regular vector wave functions are          !
!   computed. If index = 3, the radiating vector wave functions are computed.      !
! - k (complex) - wave number.                                                     !
! - r (real) - module of the position vector.                                      !
! - theta (real) - zenith angle.                                                   !
! - zRe, zIm (real arrays) - coordinates of distributed sources in the complex     !
!   plane.                                                                         !
! - m (integer) - azimuthal mode.                                                  !
! - Nrank (integer) - maximum expansion order or the number of distributed         !
!   sources.                                                                       !
!                                                                                !
! Output parameters:                                                               !
! - MV, NV (complex arrays) - distributed vector spherical wave functions.         !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: index, m, Nrank
  real(O)    :: r, theta, zRe(Nrank), zIm(Nrank)
  complex(O) :: k, MV(3,Nrank), NV(3,Nrank)
!
  integer    :: p, n, Nbes
  real(O)    :: ro, z, nm, mr, nr, cth, sth, dz
  complex(O) :: RR, sint, cost, argJ, sinc, cosc, Pmm, pimm, taumm, fp, ft, fl,     &
                factp, factt, factl, roc, dzc
  complex(O),allocatable :: jh(:), jhd(:)
!
  if (m == 0) then
    n = 1
  else
    n = abs(m)
  end if
  Nbes = n + 1
  mr   = real(m,O)
  nr   = real(n * (n + 1),O)
  nm   = 1._O / sqrt(2._O * nr)
  allocate (jh(0:Nbes), jhd(0:Nbes))
  sth = sin(theta)
  cth = cos(theta)
  ro  = r * sth
  z   = r * cth
  roc = cmplx(ro,0.0,O)
  do p = 1, Nrank
    dz  = z - zRe(p)
    dzc = cmplx(dz,0.0,O)
    RR  = sqrt(roc * roc + (dzc - im * zIm(p))**2)
    if (abs(RR) < MachEps) RR = cmplx(MachEps,MachEps,O)
    sint = ro / RR
    cost = (dzc - im * zIm(p)) / RR
    argJ = k * RR
    if (index == 1) then
      call besel_j (argJ, Nbes, jh, jhd)
    else if (index == 3) then
      call besel_h (argJ, Nbes, jh, jhd)
    end if
    call P_mm (sint, cost, abs(m), Pmm)
    call pi_mm (sint, abs(m), pimm)
    call tau_mm (sint, cost, abs(m), taumm)
    sinc = sth * cost - cth * sint
    cosc = cth * cost + sth * sint
    fp   = im * mr * pimm * nm
    ft   = taumm * nm
    fl   = nr * Pmm * nm
    factp   =   jh(n)  * fp
    MV(1,p) =   factp  * sinc
    MV(2,p) =   factp  * cosc
    MV(3,p) = - jh(n)  * ft
    factl   =   jh(n)  * fl
    factt   =   jhd(n) * ft
    NV(1,p) =   (   factl * cosc + factt * sinc) / argJ
    NV(2,p) =   ( - factl * sinc + factt * cosc) / argJ
    NV(3,p) =   jhd(n) * fp / argJ
  end do
  deallocate (jh, jhd)
end subroutine MN_DS
!***********************************************************************************
subroutine MN_poles_COMP (index, kk, r, theta, m, Npart, zpart, Nrankp, Nmax, MV, NV)
!-----------------------------------------------------------------------------------
! The routine computes the localized vector spherical wave functions (excepting    !
! the factor exp(j*m*phi)) for a composite particle and the azimuthal mode m.      !
!                                                                                !
! Input parameters:                                                                !
! - index (integer) - if index = 1, the regular vector wave functions are          !
!   computed. If index = 3, the radiating vector wave functions are computed.      !
! - kk (complex) - wave number.                                                    !
! - r (real) - module of the position vector.                                      !
! - theta (real) - zenith  angle.                                                  !
! - m (integer) - azimuthal mode.                                                  !
! - Npart (integer) - number of homogeneous regions.                               !
! - zpart (real array) -  axial coordinates of the region origins.                 !
! - Nrankp (integer array) - maximum expansion orders of the regions.              !
! - Nmax (integer) - specifies the dimension of the vectors.                       !
!                                                                                !
! Output parameters:                                                               !
! - MV, NV (complex arrays) - vector spherical wave functions.                     !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: index, m, Npart, Nmax, Nrankp(Npart)
  real(O)    :: r, theta, zpart(Npart)
  complex(O) :: kk, MV(3,Nmax), NV(3,Nmax)
!
  integer    :: k, n, p, Nstart, Nmaxp, Nbes
  real(O)    :: r1, theta1, sinc, cosc, nm, mr, nr, fp, ft, fl, cth
  complex(O) :: zz, nvr, nvtheta, nvphi, mvtheta
  real(O),allocatable    :: Pnm(:), dPnm(:), pinm(:), taunm(:)
  complex(O),allocatable :: jh(:), jhd(:)
!
  cth  = cos(theta)
  Nbes = 1
  do p = 1, Npart
    if (Nrankp(p) > Nbes) then
      Nbes = Nrankp(p)
    end if
  end do
  allocate (jh(0:Nbes), jhd(0:Nbes))
  allocate (Pnm(0:Nbes), dPnm(0:Nbes), pinm(0:Nbes), taunm(0:Nbes))
  mr = real(m,O)
  Nstart = 0
  do p = 1, Npart
    if (Nrankp(p) >= abs(m)) then
      if (m == 0) then
        Nmaxp = Nrankp(p)
      else
        Nmaxp = Nrankp(p) - abs(m) + 1
      end if
      r1 = sqrt(r**2 + zpart(p)**2 - 2._O * r * zpart(p) * cth)
      if (r1 < MachEps) r1 = MachEps
      theta1 = acos((r * cth - zpart(p)) / r1)
      zz    = kk * r1
      if (index == 1) then
        call besel_j (zz, Nrankp(p), jh, jhd)
      else if (index == 3) then
        call besel_h (zz, Nrankp(p), jh, jhd)
      end if
      call leg_normalized (theta1, abs(m), Nrankp(p), Pnm, dPnm, pinm, taunm)
      sinc = sin(theta - theta1)
      cosc = cos(theta - theta1)
      do k = 1, Nmaxp
        if (m == 0) then
          n = k
        else
          n = abs(m) + k - 1
        end if
        nr = real(n * (n + 1),O)
        nm = 1._O / sqrt(2._O * nr)
        fp = mr * pinm(n) * nm
        ft = taunm(n) * nm
        fl = nr * Pnm(n) * nm
        mvtheta = jh(n) * im * fp
        MV(1,Nstart+k) =    mvtheta * sinc
        MV(2,Nstart+k) =    mvtheta * cosc
        MV(3,Nstart+k) = - jh(n)   * ft
        nvr     = fl * jh(n) / zz
        nvtheta = jhd(n) * ft / zz
        nvphi   = jhd(n) * im * fp / zz
        NV(1,Nstart+k) =   nvr * cosc + nvtheta * sinc
        NV(2,Nstart+k) = - nvr * sinc + nvtheta * cosc
        NV(3,Nstart+k) =   nvphi
      end do
      Nstart = Nstart + Nmaxp
    endif
  end do
  deallocate (jh, jhd)
  deallocate (Pnm, dPnm, pinm, taunm)
end subroutine MN_poles_COMP
!***********************************************************************************
subroutine MN_poles_COMP1 (index, ipart, kk, r, theta, m, Npart, zpart, Nrankp,     &
           Nmax, MV, NV)
!-----------------------------------------------------------------------------------
! The routine computes the localized vector spherical wave functions (excepting    !
! the factor exp(j*m*phi)) for the homogeneous region ipart (which is a part of    !
! the composite particle) and the azimuthal mode m.                                !
!                                                                                !
! Input parameters:                                                                !
! - index (integer) - if index = 1, the regular vector wave functions are          !
!   computed. If index = 3, the radiating vector wave functions are computed.      !
! - ipart (integer) - index of the homogeneous region.                             !
! - kk (complex) - wave number.                                                    !
! - r (real) - module of the position vector.                                      !
! - theta (real) - zenith angle.                                                   !
! - m (integer) - azimuthal mode.                                                  !
! - Npart (integer) - number of homogeneous regions.                               !
! - zpart (real array) - axial coordinates of the region origins.                  !
! - Nrankp (integer array) - maximum expansion orders of the regions.              !
! - Nmax (integer) - specifies the dimension of the vectors.                       !
!                                                                                !
! Output parameters:                                                               !
! - MV, NV (complex arrays) - vector spherical wave functions.                     !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: index, m, Npart, Nmax, Nrankp(Npart), ipart
  real(O)    :: r, theta, zpart(Npart)
  complex(O) :: kk, MV(3,Nmax), NV(3,Nmax)
!
  integer    :: k, n, p, Nstart, Nmaxp, Nbes
  real(O)    :: r1, theta1, sinc, cosc, nm, mr, nr, fp, ft, fl, cth
  complex(O) :: zz, nvr, nvtheta, nvphi, mvtheta
  real(O),allocatable    :: Pnm(:), dPnm(:), pinm(:), taunm(:)
  complex(O),allocatable :: jh(:), jhd(:)
!
  cth  = cos(theta)
  Nbes = 1
  do p = 1, Npart
    if (Nrankp(p) > Nbes) then
      Nbes = Nrankp(p)
    end if
  end do
  allocate (jh(0:Nbes), jhd(0:Nbes))
  allocate (Pnm(0:Nbes), dPnm(0:Nbes), pinm(0:Nbes), taunm(0:Nbes))
  mr = real(m,O)
  do k = 1, 3
    do p = 1, Nmax
      MV(k,p) = zero
      NV(k,p) = zero
    end do
  end do
  Nstart = 0
  do p = 1, Npart
    if (Nrankp(p) >= abs(m)) then
      if (m == 0) then
        Nmaxp = Nrankp(p)
      else
        Nmaxp = Nrankp(p) - abs(m) + 1
      end if
      if (p == ipart) then
        r1 = sqrt(r**2 + zpart(p)**2 - 2._O * r * zpart(p) * cth)
        if (r1 < MachEps) r1 = MachEps
        theta1 = acos((r * cth - zpart(p)) / r1)
        zz     = kk * r1
        if (index == 1) then
          call besel_j (zz, Nrankp(p), jh, jhd)
        else if (index == 3) then
          call besel_h (zz, Nrankp(p), jh, jhd)
        end if
        call leg_normalized (theta1, abs(m), Nrankp(p), Pnm, dPnm, pinm, taunm)
        sinc = sin(theta - theta1)
        cosc = cos(theta - theta1)
        do k = 1, Nmaxp
          if (m == 0) then
            n = k
          else
            n = abs(m) + k - 1
          end if
          nr = real(n * (n + 1),O)
          nm = 1._O / sqrt(2._O * nr)
          fp = mr * pinm(n) * nm
          ft = taunm(n) * nm
          fl = nr * Pnm(n) * nm
          mvtheta = jh(n) * im * fp
          MV(1,Nstart+k) =   mvtheta * sinc
          MV(2,Nstart+k) =   mvtheta * cosc
          MV(3,Nstart+k) = - jh(n)  * ft
          nvr     = fl * jh(n) / zz
          nvtheta = jhd(n) * ft / zz
          nvphi   = jhd(n) * im * fp / zz
          NV(1,Nstart+k) =   nvr * cosc + nvtheta * sinc
          NV(2,Nstart+k) = - nvr * sinc + nvtheta * cosc
          NV(3,Nstart+k) =   nvphi
        end do
      end if
      Nstart = Nstart + Nmaxp
    endif
  end do
  deallocate (jh, jhd)
  deallocate (Pnm, dPnm, pinm, taunm)
end subroutine MN_poles_COMP1
!***********************************************************************************
subroutine MN_DS_COMP (index, k, r, theta, Npart, Nrankpmax, Nrankp, zRe, zIm, m,   &
           Nrank, MV, NV)
!-----------------------------------------------------------------------------------
! The routine computes the distributed vector spherical wave functions (excepting  !
! the factor exp(j*m*phi)) for a composite particle and the azimuthal mode m.      !
!                                                                                !
! Input parameters:                                                                !
! - index (integer) - if index = 1, the regular vector wave functions are          !
!   computed. If index = 3, the radiating vector wave functions are computed.      !
! - k (complex) - wave number.                                                     !
! - r (real) - module of the position vector.                                      !
! - theta (real) - zenith angle.                                                   !
! - m (integer) - azimuthal mode.                                                  !
! - Npart (integer) - number of homogeneous regions.                               !
! - zRe, zIm (real array) - coordinates of the region origins in the complex       !
!   plane.                                                                         !
! - Nrankpmax (integer) - specifies the dimension of zRe and zIm.                  !
! - Nrankp (integer array) - maximum expansion orders of the homogeneous regions.  !
! - Nrank (integer) - specifies the dimension of the vectors.                      !
!                                                                                !
! Output parameters:                                                               !
! - MV, NV (complex arrays) - vector spherical wave functions.                     !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: index, Npart, Nrankpmax, Nrankp(Npart), m, Nrank
  real(O)    :: r, theta, zRe(Npart,Nrankpmax), zIm(Npart,Nrankpmax)
  complex(O) :: k, MV(3,Nrank), NV(3,Nrank)
!
  integer    :: p, n, Nbes, Nstart, Nrankpl, ipart
  real(O)    :: ro, z, nm, mr, nr, sth, cth, dz
  complex(O) :: RR, sint, cost, argJ, sinc, cosc, Pmm, pimm, taumm, fp, ft, fl,     &
                factp, factt, factl, roc, dzc
  complex(O),allocatable :: jh(:), jhd(:)
!
  if (m == 0) then
    n = 1
  else
    n = abs(m)
  end if
  Nbes = n + 1
  mr   = real(m,O)
  nr   = real(n * (n + 1),O)
  nm   = 1._O / sqrt(2._O * nr)
  allocate (jh(0:Nbes), jhd(0:Nbes))
  sth = sin(theta)
  cth = cos(theta)
  ro  = r * sth
  z   = r * cth
  roc = cmplx(ro,0.0,O)
  Nstart = 0
  do ipart = 1, Npart
    Nrankpl = Nrankp(ipart)
    do p = 1, Nrankpl
      dz  = z - zRe(ipart,p)
      dzc = cmplx(dz,0.0,O)
      RR  = sqrt(roc * roc + (dzc - im * zIm(ipart,p))**2)
      if (abs(RR) < MachEps) RR = cmplx(MachEps,MachEps,O)
      sint = ro / RR
      cost = (dzc - im * zIm(ipart,p)) / RR
      argJ = k * RR
      if (index == 1) then
        call besel_j (argJ, Nbes, jh, jhd)
      else if (index == 3) then
        call besel_h (argJ, Nbes, jh, jhd)
      end if
      call P_mm (sint, cost, abs(m), Pmm)
      call pi_mm (sint, abs(m), pimm)
      call tau_mm (sint, cost, abs(m), taumm)
      sinc  = sth * cost - cth * sint
      cosc  = cth * cost + sth * sint
      fp    = im * mr * pimm * nm
      ft    = taumm * nm
      fl    = nr * Pmm * nm
      factp = jh(n) * fp
      MV(1,Nstart+p) =   factp * sinc
      MV(2,Nstart+p) =   factp * cosc
      MV(3,Nstart+p) = - jh(n) * ft
      factt = jhd(n) * ft
      factl = jh(n)  * fl
      NV(1,Nstart+p) = (  factl * cosc + factt * sinc) / argJ
      NV(2,Nstart+p) = (- factl * sinc + factt * cosc) / argJ
      NV(3,Nstart+p) = jhd(n) * fp / argJ
    end do
    Nstart = Nstart + Nrankpl
  end do
  deallocate (jh, jhd)
end subroutine MN_DS_COMP
!***********************************************************************************
subroutine MN_DS_COMP1 (index, ipart, k, r, theta, Npart, Nrankpmax, Nrankp, zRe,   &
           zIm, m, Nrank, MV, NV)
!-----------------------------------------------------------------------------------
! The routine computes the localized vector spherical wave functions (excepting    !
! the factor exp(j*m*phi)) for the homogeneous region ipart (which is a part of    !
! a composite paticle) and the azimuthal mode m.                                   !
!                                                                                !
! Input parameters:                                                                !
! - index (integer) - if index = 1, the regular vector wave functions are          !
!   computed. If index = 3, the radiating vector wave functions are computed.      !
! - ipart (integer) - index of the homogeneous region.                             !
! - k (complex) - wave number.                                                     !
! - r (real) - module of the position vector.                                      !
! - theta (real) - zenith angle.                                                   !
! - m (integer) - azimuthal mode.                                                  !
! - Npart (integer) - number of homogeneous regions.                               !
! - zRe, zIm (real array) - coordinates of the region origins in the complex       !
!   plane.                                                                         !
! - Nrankpmax (integer) - specifies the dimension of zRe and zIm.                  !
! - Nrankp (integer array) - maximum expansion orders of the homogeneous regions.  !
! - Nrank (integer) - specifies the dimension of the vectors.                      !
!                                                                                !
! Output parameters:                                                               !
! - MV, NV (complex arrays) - vector spherical wave functions.                     !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: index, ipart, Npart, Nrankpmax, Nrankp(Npart), m, Nrank
  real(O)    :: r, theta, zRe(Npart,Nrankpmax), zIm(Npart,Nrankpmax)
  complex(O) :: k, MV(3,Nrank), NV(3,Nrank)
!
  integer    :: l, p, n, Nbes, Nstart, Nrankpl, iipart
  real(O)    :: ro, z, nm, mr, nr, sth, cth, dz
  complex(O) :: RR, sint, cost, argJ, sinc, cosc, Pmm, pimm, taumm, fp, ft, fl,     &
                factp, factl, factt, roc, dzc
  complex(O),allocatable :: jh(:), jhd(:)
!
  if (m == 0) then
    n = 1
  else
    n = abs(m)
  end if
  Nbes = n + 1
  mr   = real(m,O)
  nr   = real(n * (n + 1),O)
  nm   = 1._O / sqrt(2._O * nr)
  allocate (jh(0:Nbes), jhd(0:Nbes))
  do l = 1, 3
    do p = 1, Nrank
      MV(l,p) = zero
      NV(l,p) = zero
    end do
  end do
  sth = sin(theta)
  cth = cos(theta)
  ro  = r * sth
  z   = r * cth
  roc = cmplx(ro,0.0,O)
  Nstart = 0
  do iipart = 1, Npart
    Nrankpl = Nrankp(iipart)
    if (iipart == ipart) then
      do p = 1, Nrankpl
        dz  = z - zRe(ipart,p)
        dzc = cmplx(dz,0.0,O)
        RR  = sqrt(roc * roc + (dzc - im * zIm(ipart,p))**2)
        if (abs(RR) < MachEps) RR = cmplx(MachEps,MachEps,O)
        sint = ro / RR
        cost = (dzc - im * zIm(ipart,p)) / RR
        argJ = k * RR
        if (index == 1) then
          call besel_j (argJ, Nbes, jh, jhd)
        else if (index == 3) then
          call besel_h (argJ, Nbes, jh, jhd)
        end if
        call P_mm (sint, cost, abs(m), Pmm)
        call pi_mm (sint, abs(m), pimm)
        call tau_mm (sint, cost, abs(m), taumm)
        sinc  = sth * cost - cth * sint
        cosc  = cth * cost + sth * sint
        fp    = im * mr * pimm * nm
        ft    = taumm * nm
        fl    = nr * Pmm * nm
        factp = jh(n) * fp
        MV(1,Nstart+p) =   factp * sinc
        MV(2,Nstart+p) =   factp * cosc
        MV(3,Nstart+p) = - jh(n) * ft
        factt = jhd(n) * ft
        factl = jh(n)  * fl
        NV(1,Nstart+p) = (  factl * cosc + factt * sinc) / argJ
        NV(2,Nstart+p) = (- factl * sinc + factt * cosc) / argJ
        NV(3,Nstart+p) = jhd(n) * fp / argJ
      end do
    endif
    Nstart = Nstart + Nrankpl
  end do
  deallocate (jh, jhd)
end subroutine MN_DS_COMP1
!***********************************************************************************
subroutine MN_poles_LAY (index, kk, r, theta, m, ipart, Npart, Nrankp, Nmaxpmax,    &
           Nmaxp, zpart, MV, NV)
!-----------------------------------------------------------------------------------
! The routine computes the localized vector spherical wave functions (excepting    !
! the factor exp(j*m*phi)) for the layer ipart and the azimuthal mode m.           !
!                                                                                !
! Input Parameters:                                                                !
! - index (integer) - if index = 1, the regular vector wave functions are          !
!   computed. If index = 3, the radiating vector wave functions are computed.      !
! - kk (complex) - wave number.                                                    !
! - r (real) - module of the position vector.                                      !
! - theta (real) - zenith angle.                                                   !
! - m (integer) - azimuthal mode.                                                  !
! - ipart (integer) - index of the layer.                                          !
! - Npart (integer) - number of layers.                                            !
! - zpart (real array) - axial coordinates of the layer origins.                   !
! - Nrankp (integer array) - maximum expansion orders of the layers.               !
! - Nmaxp (integer array) - real dimensions of the vector functions corresponding  !
!   to a specific layer.                                                           !
! - Nmaxpmax (integer) - specifies the physical dimension of the vectors.          !
!                                                                                !
! Output parameters:                                                               !
! - MV, NV (complex arrays) - vector spherical wave functions.                     !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: index, m, ipart, Npart, Nrankp(Npart), Nmaxpmax, Nmaxp(Npart)
  real(O)    :: r, theta, zpart(Npart)
  complex(O) :: kk, MV(3,Nmaxpmax), NV(3,Nmaxpmax)
!
  integer    :: l, k, n, p, Nbes
  real(O)    :: r1, theta1, sinc, cosc, nm, mr, nr, fp, ft, fl, cth
  complex(O) :: zz, nvr, nvtheta, nvphi, mvtheta
  real(O),allocatable    :: Pnm(:), dPnm(:), pinm(:), taunm(:)
  complex(O),allocatable :: jh(:), jhd(:)
!
  cth = cos(theta)
  do l = 1, 3
    do p = 1, Nmaxpmax
      MV(l,p) = zero
      NV(l,p) = zero
    end do
  end do
  Nbes = Nrankp(ipart)
  allocate (jh(0:Nbes), jhd(0:Nbes))
  allocate (Pnm(0:Nbes), dPnm(0:Nbes), pinm(0:Nbes), taunm(0:Nbes))
  mr = real(m,O)
  r1 = sqrt(r**2 + zpart(ipart)**2 - 2._O * r * zpart(ipart) * cth)
  if (r1 < MachEps) r1 = MachEps
  theta1 = acos((r * cth - zpart(ipart)) / r1)
  zz    = kk * r1
  if (index == 1) then
    call besel_j (zz, Nrankp(ipart), jh, jhd)
  else if (index == 3) then
    call besel_h (zz, Nrankp(ipart), jh, jhd)
  end if
  call leg_normalized (theta1, abs(m), Nrankp(ipart), Pnm, dPnm, pinm, taunm)
  sinc = sin(theta - theta1)
  cosc = cos(theta - theta1)
  do k = 1, Nmaxp(ipart)
    if (m == 0) then
      n = k
    else
      n = abs(m) + k - 1
    end if
    nr = real(n * (n + 1),O)
    nm = 1._O / sqrt(2._O * nr)
    fp = mr * pinm(n) * nm
    ft = taunm(n) * nm
    fl = nr * Pnm(n) * nm
    mvtheta = jh(n) * im * fp
    MV(1,k) =   mvtheta * sinc
    MV(2,k) =   mvtheta * cosc
    MV(3,k) = - jh(n)  * ft
    nvr     =   fl * jh(n) / zz
    nvtheta =   jhd(n) * ft / zz
    nvphi   =   jhd(n) * im * fp / zz
    NV(1,k) =   nvr * cosc + nvtheta * sinc
    NV(2,k) = - nvr * sinc + nvtheta * cosc
    NV(3,k) =   nvphi
  end do
  deallocate (jh, jhd)
  deallocate (Pnm, dPnm, pinm, taunm)
end subroutine MN_poles_LAY
!***********************************************************************************
subroutine MN_DS_LAY (index, k, r, theta, ipart, Npart, Nrankpmax, Nrankp,          &
           zRe, zIm, m, MV, NV)
!-----------------------------------------------------------------------------------
! The routine computes the distributed vector spherical wave functions (excepting  !
! the factor exp(j*m*phi)) for the layer ipart and the azimuthal mode m.           !
!                                                                                !
! Input parameters:                                                                !
! - index (integer) - if index = 1, the regular vector wave functions are          !
!   computed. If index = 3, the radiating vector wave functions are computed.      !
! - k (complex) - wave number.                                                     !
! - r (real) - module of the position vector.                                      !
! - theta (real) - zenith angle.                                                   !
! - m (integer) - azimuthal mode.                                                  !
! - ipart (integer) - index of the layer.                                          !
! - Npart (integer) - number of layers.                                            !
! - zRe, zIm (real array) - coordinates of the region origins in the complex       !
!   plane.                                                                         !
! - Nrankp (integer array) - maximum expansion orders of the layers.               !
! - Nrankpmax (integer) - specifies the physical dimension of the vectors.         !
!                                                                                !
! Output parameters:                                                               !
! - MV, NV (complex arrays) - vector spherical wave functions.                     !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: index, ipart, Npart, Nrankpmax, Nrankp(Npart), m
  real(O)    :: r, theta
  real(O)    :: zRe(Npart,Nrankpmax), zIm(Npart,Nrankpmax)
  complex(O) :: k, MV(3,Nrankpmax), NV(3,Nrankpmax)
!
  integer    :: p, n, Nbes, Nrankpl, l
  real(O)    :: ro, z, nm, mr, nr, sth, cth, dz
  complex(O) :: Pmm, pimm, taumm, RR, sint, cost, argJ, sinc, cosc, fp, ft, fl,     &
                factp, factl, factt, roc, dzc
  complex(O),allocatable :: jh(:), jhd(:)
!
  if (m == 0) then
    n = 1
  else
    n = abs(m)
  end if
  Nbes = n + 1
  mr   = real(m,O)
  nr   = real(n * (n + 1),O)
  nm   = 1._O / sqrt(2._O * nr)
  allocate (jh(0:Nbes), jhd(0:Nbes))
  do l = 1, 3
    do p = 1, Nrankpmax
      MV(l,p) = zero
      NV(l,p) = zero
    end do
  end do
  sth = sin(theta)
  cth = cos(theta)
  ro  = r * sth
  z   = r * cth
  roc = cmplx(ro,0.0,O)
  Nrankpl = Nrankp(ipart)
  do p = 1, Nrankpl
    dz  = z - zRe(ipart,p)
    dzc = cmplx(dz,0.0,O)
    RR  = sqrt(roc * roc + (dzc - im * zIm(ipart,p))**2)
    if (abs(RR) < MachEps) RR = cmplx(MachEps,MachEps,O)
    sint = ro / RR
    cost = (dzc - im * zIm(ipart,p)) / RR
    argJ = k * RR
    if (index == 1) then
      call besel_j (argJ, Nbes, jh, jhd)
    else if (index == 3) then
      call besel_h (argJ, Nbes, jh, jhd)
    end if
    call P_mm (sint, cost, abs(m), Pmm)
    call pi_mm (sint, abs(m), pimm)
    call tau_mm (sint, cost, abs(m), taumm)
    sinc  = sth * cost - cth * sint
    cosc  = cth * cost + sth * sint
    fp    = im * mr * pimm * nm
    ft    = taumm * nm
    fl    = nr * Pmm * nm
    factp = jh(n) * fp
    MV(1,p) =   factp * sinc
    MV(2,p) =   factp * cosc
    MV(3,p) = - jh(n) * ft
    factt   = jhd(n) * ft
    factl   = jh(n) * fl
    NV(1,p) = (  factl * cosc + factt * sinc) / argJ
    NV(2,p) = (- factl * sinc + factt * cosc) / argJ
    NV(3,p) = jhd(n) * fp / argJ
  end do
  deallocate (jh, jhd)
end subroutine MN_DS_LAY
!***********************************************************************************
subroutine MN_complete (index, z, theta, phi, Mrank, Nrank, Nmax, plus, sym, MV, NV)
!-----------------------------------------------------------------------------------
! The routine computes the localized vector spherical wave functions for the       !
! azimuthal modes m = 0,1,...,Mrank.                                               !
!                                                                                !
! Input parameters:                                                                !
! - index (integer) - if index = 1, the regular vector wave functions are          !
!   computed. If index = 3, the radiating vector wave functions are computed.      !
! - z (complex) - z = k*r, where k is the wave number and r is the module of the   !
!   position vector.                                                               !
! - theta, phi (real variables) - polar angles.                                    !
! - Mrank (integer) - number of azimuthal modes.                                   !
! - Nrank (integer) - maximum expansion order.                                     !
! - Nmax (integer) - specifies the dimension of the vectors.                       !
! - plus (logical) -  if plus = true the vector functions are computed for the     !
!   azimuthal modes m = 0,+1,-1,...,+Mrank,-Mrank. If plus = false the vector      !
!   functions are computed for the azimuthal modes m = 0,-1,+1,...,-Mrank,+Mrank.  !
! - sym (logical) - if sym = true the azimuthal dependence exp(j*m*phi) is         !
!   omitted. If sym = false, this dependence is taken into account.                !
!                                                                                !
! Output parameters:                                                               !
! - MV, NV (complex arrays) - vector spherical wave functions.                     !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: index, Mrank, Nrank, Nmax
  real(O)    :: theta, phi
  complex(O) :: z, MV(3,Nmax), NV(3,Nmax)
  logical    :: plus, sym
!
  integer    :: k, m, l, ml, N0, n
  real(O)    :: nr, nm, mlr, ft, fl, arg
  complex(O) :: zc, fact, factp, factt, factl
  real(O),allocatable    :: Pnm(:), dPnm(:), pinm(:), taunm(:)
  complex(O),allocatable :: jh(:), jhd(:)
!
  allocate (jh(0:Nrank), jhd(0:Nrank))
  allocate (Pnm(0:Nrank), dPnm(0:Nrank), pinm(0:Nrank), taunm(0:Nrank))
  zc = z
  if (abs(zc) < MachEps) zc = cmplx(MachEps,MachEps,O)
  if (index == 1) then
    call besel_j (zc, Nrank, jh, jhd)
  else if (index == 3) then
    call besel_h (zc, Nrank, jh, jhd)
  end if
  do m = 0, Mrank
    call leg_normalized (theta, m, Nrank, Pnm, dPnm, pinm, taunm)
    if (m == 0) then
      do k = 1, Nrank
        n  = k
        nr = real(n * (n + 1),O)
        nm = 1._O / sqrt(2._O * nr)
        ft = taunm(n) * nm
        fl = nr * Pnm(n) * nm
        MV(1,k) =   zero
        MV(2,k) =   zero
        MV(3,k) = - jh(n)  * ft
        NV(1,k) =   jh(n)  * fl / zc
        NV(2,k) =   jhd(n) * ft / zc
        NV(3,k) =   zero
      end do
    else
      N0 = Nrank + (m - 1) * (2 * Nrank - m + 2)
      ml = m
      if (.not. plus) ml = - m
      do l = 1, 2
        mlr = real(ml,O)
        if (sym) then
          fact = one
        else
          arg  = mlr * phi
          fact = exp(im * arg)
        end if
        do k = 1, Nrank - m + 1
          n     = m + k - 1
          nr    = real(n * (n + 1),O)
          nm    = 1._O / sqrt(2._O * nr)
          factp = im * mlr * pinm(n) * nm * fact
          factt = taunm(n) * nm * fact
          factl = nr * Pnm(n) * nm * fact
          MV(1,N0+k) =   zero
          MV(2,N0+k) =   jh(n)  * factp
          MV(3,N0+k) = - jh(n)  * factt
          NV(1,N0+k) =   jh(n)  * factl / zc
          NV(2,N0+k) =   jhd(n) * factt / zc
          NV(3,N0+k) =   jhd(n) * factp / zc
        end do
        N0 =   N0 + Nrank - m + 1
        ml = - ml
      end do
    end if
  end do
  deallocate (jh, jhd)
  deallocate (Pnm, dPnm, pinm, taunm)
end subroutine MN_complete
!***********************************************************************************
subroutine MN_infinit_complete (theta, phi, Mrank, Nrank, Nmax, plus, Minf, Ninf)
!-----------------------------------------------------------------------------------
! The routine computes the localized vector spherical wave functions at infinity   !
! (excepting the factor exp(j*k*R)/(k*R)) for the azimuthal modes                  !
!  m = 0,1,...,Mrank.                                                              !
!                                                                                !
! Input parameters:                                                                !
! - theta, phi (real variables) - polar angles.                                    !
! - Mrank (integer) - number of azimuthal modes.                                   !
! - Nrank (integer) - maximum expansion order.                                     !
! - Nmax (integer) - specifies the dimension of the vectors.                       !
! - plus (logical) - if plus = true the vector functiosn are computed for the      !
!   azimuthal modes m = 0,+1,-1,...,+Mrank,-Mrank. If plus = false, the vector     !
!   functions are computed for the azimuthal modes m = 0,-1,+1,...,-Mrank,+Mrank.  !
!                                                                                !
! Output parameters:                                                               !
! - Minf, Ninf (complex arrays) - vector spherical wave functions.                 !
!-----------------------------------------------------------------------------------
  use parameters
  implicit none
  integer    :: Mrank, Nrank, Nmax
  real(O)    :: theta, phi
  complex(O) :: Minf(3,Nmax), Ninf(3,Nmax)
  logical    :: plus
!
  integer    :: k, m, l, ml, N0, n
  real(O)    :: nm, mlr, arg
  complex(O) :: fact, factc, factt, factp
  real(O),allocatable :: Pnm(:), dPnm(:), pinm(:), taunm(:)
!
  allocate (Pnm(0:Nrank), dPnm(0:Nrank), pinm(0:Nrank), taunm(0:Nrank))
  do m = 0, Mrank
    call leg_normalized (theta, m, Nrank, Pnm, dPnm, pinm, taunm)
    if (m == 0) then
      do k = 1, Nrank
        n     = k
        nm    = real(2 * n * (n + 1),O)
        nm    = 1._O / sqrt(nm)
        factc = (-im)**(n + 1) * nm
        factt = factc * taunm(n)
        Minf(1,k) =   zero
        Minf(2,k) =   zero
        Minf(3,k) = - factt
        Ninf(1,k) =   zero
        Ninf(2,k) =   im * factt
        Ninf(3,k) =   zero
      end do
    else
      N0 = Nrank + (m - 1) * (2 * Nrank - m + 2)
      ml = m
      if (.not. plus) ml = - m
      do l = 1, 2
        mlr  = real(ml,O)
        arg  = mlr * phi
        fact = exp(im * arg)
        do k = 1, Nrank - m + 1
          n     = m + k - 1
          nm    = real(2 * n * (n + 1),O)
          nm    = 1._O / sqrt(nm)
          factc = fact * (-im)**(n + 1) * nm
          factp = factc * mlr * pinm(n)
          factt = factc * taunm(n)
          Minf(1,N0+k) =   zero
          Minf(2,N0+k) =   im * factp
          Minf(3,N0+k) = - factt
          Ninf(1,N0+k) =   zero
          Ninf(2,N0+k) =   im * factt
          Ninf(3,N0+k) = - factp
        end do
        N0 =   N0 + Nrank - m + 1
        ml = - ml
      end do
    end if
  end do
  deallocate (Pnm, dPnm, pinm, taunm)
end subroutine  MN_infinit_complete
!***********************************************************************************
subroutine MN_infinit_reflection_complete (wavenumber, z0, ind_ref, theta, phi,     &
           Mrank, Nrank, Nmax, Minf, Ninf)
!-----------------------------------------------------------------------------------
! The routine computes the reflected localized vector spherical wave functions     !
! at infinity (excepting the factor exp(j*k*R)/(k*R)) and the azimuthal modes      !
! m = 0,+1,-1,...,+Mrank,-Mrank.                                                   !
!                                                                                !
! Input parameters:                                                                !
! - wavenumber (real) - wave number in the ambient medium.                         !
! - z0 (real) - axial coordinate of the substrate.                                 !
! - ind_ref (complex) - relative refractive index of the substrate.                !
! - theta, phi (real variables) - polar angles.                                    !
! - Mrank (integer) - number of azimuthal modes.                                   !
! - Nrank (integer) - maximum expansion order.                                     !
! - Nmax (integer) - specifies the dimension of the vectors.                       !
!                                                                                !
! Output parameters:                                                               !
! - Minf, Ninf (complex arrays) - vector spherical wave functions.                 !
!-----------------------------------------------------------------------------------
  use parameters
  implicit none
  integer    :: Mrank, Nrank, Nmax, k, m, N0, ml, l, n
  real(O)    :: theta, theta0, phi, wavenumber, z0, arga, mlr, nm
  complex(O) :: ind_ref, Rpar, Rperp, Tpar, Tperp, cost, sint, cost0, fact,         &
                facta, factr, factc, factt, factp, arg, Minf(3,Nmax), Ninf(3,Nmax)
  real(O),allocatable :: Pnm(:), dPnm(:), pinm(:), taunm(:)
!
  allocate (Pnm(0:Nrank), dPnm(0:Nrank), pinm(0:Nrank), taunm(0:Nrank))
  theta0 = Pi - theta
  cost0 = cmplx(cos(theta0),0.0,O)
  call Fresnel_aer_sub (cost0, ind_ref, Rpar, Rperp, Tpar, Tperp, cost, sint)
  arg   = 2._O * cost0 * wavenumber * z0
  factr = exp(im * arg)
  do m = 0, Mrank
    call leg_normalized (theta0, m, Nrank, Pnm, dPnm, pinm, taunm)
    if (m == 0) then
      do k = 1, Nrank
        n     = k
        nm    = real(2 * n * (n + 1),O)
        nm    = 1._O / sqrt(nm)
        factc = factr * (-im)**(n + 1) * nm
        factt = factc * taunm(n)
        Minf(1,k) =   zero
        Minf(2,k) =   zero
        Minf(3,k) = - factt * Rperp
        Ninf(1,k) =   zero
        Ninf(2,k) =   im * factt * Rpar
        Ninf(3,k) =   zero
      end do
    else
      N0 = Nrank + (m - 1) * (2 * Nrank - m + 2)
      ml = m
      do l = 1, 2
        mlr   = real(ml,O)
        arga  = mlr * phi
        facta = exp(im * arga)
        fact  = facta * factr
        do k = 1, Nrank - m + 1
          n     = m + k - 1
          nm    = real(2 * n * (n + 1),O)
          nm    = 1._O / sqrt(nm)
          factc = fact * (-im)**(n + 1) * nm
          factp = factc * mlr * pinm(n)
          factt = factc * taunm(n)
          Minf(1,N0+k) =   zero
          Minf(2,N0+k) =   im * factp * Rpar
          Minf(3,N0+k) = - factt * Rperp
          Ninf(1,N0+k) =   zero
          Ninf(2,N0+k) =   im * factt * Rpar
          Ninf(3,N0+k) = - factp * Rperp
        end do
        N0 =   N0 + Nrank - m + 1
        ml = - ml
      end do
    end if
  end do
  deallocate (Pnm, dPnm, pinm, taunm)
end subroutine MN_infinit_reflection_complete
!***********************************************************************************
subroutine MN_left_right (Nmax, mvl, nvl, mvr, nvr, mv, nv)
  use parameters
  implicit none
  integer    :: Nmax, i, j
  complex(O) :: mvl(3,Nmax), nvl(3,Nmax), mvr(3,Nmax), nvr(3,Nmax), mv(3,Nmax),     &
                nv(3,Nmax)
!
  do i = 1, 3
    do j = 1, Nmax
      mv(i,j) = mvl(i,j) + nvl(i,j)
      nv(i,j) = mvr(i,j) - nvr(i,j)
    end do
  end do
end subroutine MN_left_right
!***********************************************************************************
subroutine MN_anisotrop (k1, ind_ref, ind_refZ, rBIG, thetaBIG, phiBIG, alphaP,    &
           betaP, gammaP, Mrank, Nrank, Nmax, Nbeta, Xe, Ye, Xh, Yh)
!-----------------------------------------------------------------------------------
! The routine computes the system of vector wave functions (vector quasispherical  !
! wave functions) for uniaxial anisotropic particles.                              !
!                                                                                !
! Input parameters:                                                                !
! - k1 (complex) - wave number.                                                    !
! - ind_ref, ind_refZ (complex variables) - relative refractive indices of the     !
!   uniaxial anisotropic particle.                                                 !
! - rBIG, thetaBIG, phiBIG (real avriables) - spherical coordinates of the         !
!   generic point in the particle coordinate system.                               !
! - alphaP, betaP, gammaP (real variables) - Euler angles specifying the           !
!   orientation of the principal coordinate system with respect to the particle    !
!   coordinate system.                                                             !
! - Mrank (integer) - number of azimuthal modes.                                   !
! - Nrank (integer) - maximum expansion order.                                     !
! - Nmax (integer) - specifies the dimension of the vectors.                       !
! - Nbeta (integer) - number of integration points.                                !
!                                                                                !
! Output parameters:                                                               !
! - Xe, Ye, Xh, Yh (complex arrays) - vector wave functions for the uniaxial       !
!   anisotropic particle.                                                          !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: Mrank, Nrank, Nmax, Nbeta
  real(O)    :: rBIG, thetaBIG, phiBIG, alphaP, betaP, gammaP
  complex(O) :: k1, ind_ref, ind_refZ, Xe(3,Nmax), Ye(3,Nmax), Xh(3,Nmax), Yh(3,Nmax)
!
  integer    :: i, k, p, m, n, N0, ml, l
  real(O)    :: xBIG, yBIG, zBIG, x, y, z, r, theta, phi, abeta, bbeta, beta, fact, &
                Pi2, ct, st, cp, sp, cth, sth, cph, sph, cb, sb, nm, mr, mrp, mrm,  &
                arg, argp, argm, s, mlr, stcp, stsp, sthcph, sthsph, cthcph, cthsph
  complex(O) :: f1, f2, x1, x2, y1, y2, fcp, fsp, f1s, ex1, ex2, f12p, f12m, o1, o2,&
                s1, c1, s2, c2, ft, fp, fJ1, fJ2, fJ1p, fJ1m, fJ2p, fJ2m, Imx1,     &
                Imx2, ImCx1, ImCx2, ImSx1, ImSx2, Xe1, Xe2, Xe3, Ye1, Ye2, Ye3,     &
                Xh1, Xh2, Xh3, Yh1, Yh2, Yh3, ind_refZc, ind_ref2, ind_refZ2,       &
                k1rsth, k1rcth, Pi2im, Piimp, Piimm
  real(O),allocatable    :: Pnm(:), dPnm(:), pinm(:), taunm(:), xbeta(:), wbeta(:)
  complex(O),allocatable :: J1(:), J2(:)
!
  Pi2  = 2._O * Pi
  ct   = cos(thetaBIG)
  st   = sin(thetaBIG)
  cp   = cos(phiBIG)
  sp   = sin(phiBIG)
  stcp = st * cp
  stsp = st * sp
  xBIG = rBIG * stcp
  yBIG = rBIG * stsp
  zBIG = rBIG * ct
  call T_cartesian_global_local (xBIG, yBIG, zBIG, alphaP, betaP, gammaP, x, y, z)
  call T_cartesian_spherical (x, y, z, r, theta, phi)
  cth = cos(theta)
  sth = sin(theta)
  cph = cos(phi)
  sph = sin(phi)
  sthcph = sth * cph
  sthsph = sth * sph
  cthcph = cth * cph
  cthsph = cth * sph
  fcp = Pi2 * im * cph
  fsp = Pi2 * im * sph
  k1rsth = k1 * r * sth
  k1rcth = k1 * r * cth
  do i = 1, 3
    do k = 1, Nmax
      Xe(i,k) = zero
      Ye(i,k) = zero
      Xh(i,k) = zero
      Yh(i,k) = zero
    end do
  end do
  allocate (xbeta(Nbeta), wbeta(Nbeta))
  abeta = 0._O
  bbeta = Pi
  call Gauss_Legendre (abeta, bbeta, Nbeta, wbeta, xbeta)
  allocate (J1(0:Mrank+1), J2(0:Mrank+1))
  allocate (Pnm(0:Nrank), dPnm(0:Nrank), pinm(0:Nrank), taunm(0:Nrank))
  ind_refZc = ind_refZ
  if (abs(ind_refZc) < MachEps) ind_refZc = cmplx(MachEps,MachEps,O)
  ind_ref2  = ind_ref * ind_ref
  ind_refZ2 = ind_refZc * ind_refZc
  do p = 1, Nbeta
    beta = xbeta(p)
    cb   = cos(beta)
    sb   = sin(beta)
    f1   = (ind_refZ2 * cb**2 + ind_ref2 * sb**2) / ind_refZ2
    f2   = (ind_refZ2 - ind_ref2) * sb * cb / ind_refZ2
    f1s  = sqrt(f1)
    x1   = k1rsth * sb
    x2   = x1 / f1s
    y1   = k1rcth * cb
    y2   = y1 / f1s
    ex1  = exp(im * y1)
    ex2  = exp(im * y2)
    f12p =   f1 * cb + f2 * sb
    f12m = - f1 * sb + f2 * cb
    call bes_J (x1, Mrank + 1, J1)
    call bes_J (x2, Mrank + 1, J2)
    fact = - wbeta(p) * sb / 4._O / Pi
    do m = 0, Mrank
      if (m == 0) then
        call Leg_normalized (beta, m, Nrank, Pnm, dPnm, pinm, taunm)
        Imx1  = Pi2 * J1(0)
        Imx2  = Pi2 * J2(0)
        ImCx1 = fcp * J1(1)
        ImCx2 = fcp * J2(1)
        ImSx1 = fsp * J1(1)
        ImSx2 = fsp * J2(1)
        o1 = Imx1  * ex1
        o2 = Imx2  * ex2
        s1 = ImSx1 * ex1
        c1 = ImCx1 * ex1
        s2 = ImSx2 * ex2
        c2 = ImCx2 * ex2
        do k = 1, Nrank
          n   = k
          nm  = real(2 * n * (n + 1),O)
          nm  = 1._O / sqrt(nm)
          ft  = taunm(n) * nm * (-im)**(n + 1) * fact
          f1  = im * ft
          Xe1 = - f1 * s1
          Xe2 =   f1 * c1
          Xe3 =   zero
          Xe(1,k) = Xe(1,k) + Xe1
          Xe(2,k) = Xe(2,k) + Xe2
          Xe(3,k) = Xe(3,k) + Xe3
!
          f2  = f12p * ft
          Ye1 = f2   * c2
          Ye2 = f2   * s2
          Ye3 = f12m * ft * o2
          Ye(1,k) = Ye(1,k) + Ye1
          Ye(2,k) = Ye(2,k) + Ye2
          Ye(3,k) = Ye(3,k) + Ye3
!
          f1  =   ft * cb
          Xh1 =   f1 * c1
          Xh2 =   f1 * s1
          Xh3 = - ft * sb * o1
          Xh(1,k) = Xh(1,k) + Xh1
          Xh(2,k) = Xh(2,k) + Xh2
          Xh(3,k) = Xh(3,k) + Xh3
!
          f2  =   f1s * im * ft
          Yh1 = - f2 * s2
          Yh2 =   f2 * c2
          Yh3 =   zero
          Yh(1,k) = Yh(1,k) + Yh1
          Yh(2,k) = Yh(2,k) + Yh2
          Yh(3,k) = Yh(3,k) + Yh3
        end do
      else
        N0    = Nrank + (m - 1) * (2 * Nrank - m + 2)
        mr    = real(m,O)
        mrp   = real(m + 1,O)
        mrm   = real(m - 1,O)
        arg   = mr  * phi
        argp  = mrp * phi
        argm  = mrm * phi
        Pi2im = Pi2 * im**m
        Piimp = Pi  * im**(m + 1)
        Piimm = Pi  * im**(m - 1)
        call Leg_normalized (beta, m, Nrank, Pnm, dPnm, pinm, taunm)
        ml = m
        s  = 1._O
        do l = 1, 2
          mlr   = real(ml,O)
          f1    = exp(s * im * arg)
          fJ1   = f1 * J1(m)
          fJ2   = f1 * J2(m)
          f1    = exp(s * im * argp)
          fJ1p  = f1 * J1(m+1)
          fJ2p  = f1 * J2(m+1)
          f1    = exp(s * im * argm)
          fJ1m  = f1 * J1(m-1)
          fJ2m  = f1 * J2(m-1)
          f1    = Pi2im
          Imx1  = f1 * fJ1
          Imx2  = f1 * fJ2
          f1    = Piimp * fJ1p
          f2    = Piimm * fJ1m
          ImCx1 =   f1 + f2
          ImSx1 = - im * (s * f1 - s * f2)
          f1    = Piimp * fJ2p
          f2    = Piimm * fJ2m
          ImCx2 = f1 + f2
          ImSx2 = - im * (s * f1 - s * f2)
          o1 = Imx1  * ex1
          o2 = Imx2  * ex2
          s1 = ImSx1 * ex1
          c1 = ImCx1 * ex1
          s2 = ImSx2 * ex2
          c2 = ImCx2 * ex2
          do k = 1, Nrank - m + 1
            n   = m + k - 1
            nm  = real(2 * n * (n + 1),O)
            nm  = 1._O / sqrt(nm)
            fp  = mlr * pinm(n) * nm * (-im)**(n + 1) * fact
            ft  = taunm(n) * nm * (-im)**(n + 1) * fact
            f1  = im * ft
            f2  = f12p * fp
            Xe1 = - f1 * s1 + f2 * c2
            Xe2 =   f1 * c1 + f2 * s2
            Xe3 =   f12m * fp * o2
            Xe(1,k+N0) = Xe(1,k+N0) + Xe1
            Xe(2,k+N0) = Xe(2,k+N0) + Xe2
            Xe(3,k+N0) = Xe(3,k+N0) + Xe3
!
            f1  =   im * fp
            f2  =   f12p * ft
            Ye1 = - f1 * s1 + f2 * c2
            Ye2 =   f1 * c1 + f2 * s2
            Ye3 =   f12m * ft * o2
            Ye(1,k+N0) = Ye(1,k+N0) + Ye1
            Ye(2,k+N0) = Ye(2,k+N0) + Ye2
            Ye(3,k+N0) = Ye(3,k+N0) + Ye3
!
            f1  =   ft * cb
            f2  =   f1s * im * fp
            Xh1 =   f1 * c1 - f2 * s2
            Xh2 =   f1 * s1 + f2 * c2
            Xh3 = - ft * sb * o1
            Xh(1,k+N0) = Xh(1,k+N0) + Xh1
            Xh(2,k+N0) = Xh(2,k+N0) + Xh2
            Xh(3,k+N0) = Xh(3,k+N0) + Xh3
!
            f1  =   fp * cb
            f2  =   f1s * im * ft
            Yh1 =   f1 * c1 - f2 * s2
            Yh2 =   f1 * s1 + f2 * c2
            Yh3 = - fp * sb * o1
            Yh(1,k+N0) = Yh(1,k+N0) + Yh1
            Yh(2,k+N0) = Yh(2,k+N0) + Yh2
            Yh(3,k+N0) = Yh(3,k+N0) + Yh3
          end do
          N0 =   N0 + Nrank - m + 1
          ml = - m
          s  = - s
        end do
      end if
    end do
  end do
! --- pass to spherical coordinates ---
  do k = 1, Nmax
    Xe1 = Xe(1,k)
    Xe2 = Xe(2,k)
    Xe3 = Xe(3,k)
    Xe(1,k) =   Xe1 * sthcph + Xe2 * sthsph + Xe3 * cth
    Xe(2,k) =   Xe1 * cthcph + Xe2 * cthsph - Xe3 * sth
    Xe(3,k) = - Xe1 * sph    + Xe2 * cph
!
    Ye1 = Ye(1,k)
    Ye2 = Ye(2,k)
    Ye3 = Ye(3,k)
    Ye(1,k) =   Ye1 * sthcph + Ye2 * sthsph + Ye3 * cth
    Ye(2,k) =   Ye1 * cthcph + Ye2 * cthsph - Ye3 * sth
    Ye(3,k) = - Ye1 * sph    + Ye2 * cph
!
    Xh1 = Xh(1,k)
    Xh2 = Xh(2,k)
    Xh3 = Xh(3,k)
    Xh(1,k) =   Xh1 * sthcph + Xh2 * sthsph + Xh3 * cth
    Xh(2,k) =   Xh1 * cthcph + Xh2 * cthsph - Xh3 * sth
    Xh(3,k) = - Xh1 * sph    + Xh2 * cph
!
    Yh1 = Yh(1,k)
    Yh2 = Yh(2,k)
    Yh3 = Yh(3,k)
    Yh(1,k) =   Yh1 * sthcph + Yh2 * sthsph + Yh3 * cth
    Yh(2,k) =   Yh1 * cthcph + Yh2 * cthsph - Yh3 * sth
    Yh(3,k) = - Yh1 * sph    + Yh2 * cph
  end do
  deallocate (J1, J2, xbeta, wbeta)
  deallocate (Pnm, dPnm, pinm, taunm)
end subroutine MN_anisotrop
!***********************************************************************************
subroutine MN_m_anisotrop (k1, ind_ref, ind_refZ, rBIG, thetaBIG, phiBIG, alphaP,   &
           betaP, gammaP, m, Nrank, Nmax, Nbeta, Xe, Ye, Xh, Yh)

!-----------------------------------------------------------------------------------
! The routine computes the system of vector wave functions (vector quasispherical  !
! wave functions) for uniaxial anisotropic particles and the azimuthal mode m.     !
!                                                                                !
! Input parameters:                                                                !
! - k1 (complex) - wave number.                                                    !
! - ind_ref, ind_refZ (complex variables) - relative refractive indices of the     !
!   uniaxial anisotropic particle.                                                 !
! - rBIG, thetaBIG, phiBIG (real avriables) - spherical coordinates of the         !
!   generic point in the particle coordinate system.                               !
! - alphaP, betaP, gammaP (real variables) - Euler angles specifying the           !
!   orientation of the principal coordinate system with respect to the particle    !
!   coordinate system.                                                             !
! - m (integer) - azimuthal mode. m is postive or negative.                        !
! - Nrank (integer) - maximum expansion order.                                     !
! - Nmax (integer) - specifies the dimension of the vectors.                       !
! - Nbeta (integer) - number of integration points.                                !
!                                                                                !
! Output parameters:                                                               !
! - Xe, Ye, Xh, Yh (complex arrays) - vector wave functions for the uniaxial       !
!   anisotropic particle.                                                          !
!-----------------------------------------------------------------------------------
  use parameters
  use derived_parameters
  implicit none
  integer    :: m, Nrank, Nmax, Nbeta
  real(O)    :: rBIG, thetaBIG, phiBIG, alphaP, betaP, gammaP
  complex(O) :: k1, ind_ref, ind_refZ, Xe(3,Nmax), Ye(3,Nmax), Xh(3,Nmax), Yh(3,Nmax)
!
  integer    :: i, k, p, mabs, n
  real(O)    :: xBIG, yBIG, zBIG, x, y, z, r, theta, phi, abeta, bbeta, beta, fact, &
                Pi2, ct, st, cp, sp, cth, sth, cph, sph, cb, sb, nm, mr, mrp, mrm,  &
                arg, argp, argm, s, mlr, stcp, stsp, sthcph, sthsph, cthcph, cthsph
  complex(O) :: f1, f2, x1, x2, y1, y2, fcp, fsp, f1s, ex1, ex2, f12p, f12m, o1, o2,&
                s1, c1, s2, c2, ft, fp, fJ1, fJ2, fJ1p, fJ1m, fJ2p, fJ2m, Imx1,     &
                Imx2, ImCx1, ImCx2, ImSx1, ImSx2, Xe1, Xe2, Xe3, Ye1, Ye2, Ye3,     &
                Xh1, Xh2, Xh3, Yh1, Yh2, Yh3, ind_refZc, ind_ref2, ind_refZ2,       &
                k1rsth, k1rcth, Pi2im, Piimp, Piimm
  real(O),allocatable    :: Pnm(:), dPnm(:), pinm(:), taunm(:), xbeta(:), wbeta(:)
  complex(O),allocatable :: J1(:), J2(:)
!
  Pi2  = 2._O * Pi
  ct   = cos(thetaBIG)
  st   = sin(thetaBIG)
  cp   = cos(phiBIG)
  sp   = sin(phiBIG)
  stcp = st * cp
  stsp = st * sp
  xBIG = rBIG * stcp
  yBIG = rBIG * stsp
  zBIG = rBIG * ct
  call T_cartesian_global_local (xBIG, yBIG, zBIG, alphaP, betaP, gammaP, x, y, z)
  call T_cartesian_spherical (x, y, z, r, theta, phi)
  cth = cos(theta)
  sth = sin(theta)
  cph = cos(phi)
  sph = sin(phi)
  sthcph = sth * cph
  sthsph = sth * sph
  cthcph = cth * cph
  cthsph = cth * sph
  fcp = Pi2 * im * cph
  fsp = Pi2 * im * sph
  k1rsth = k1 * r * sth
  k1rcth = k1 * r * cth
  do i = 1, 3
    do k = 1, Nmax
      Xe(i,k) = zero
      Ye(i,k) = zero
      Xh(i,k) = zero
      Yh(i,k) = zero
    end do
  end do
  allocate (xbeta(Nbeta), wbeta(Nbeta))
  abeta = 0._O
  bbeta = Pi
  call Gauss_Legendre (abeta, bbeta, Nbeta, wbeta, xbeta)
  mabs  = abs(m)
  mr    = real(mabs,O)
  mrp   = real(mabs + 1,O)
  mrm   = real(mabs - 1,O)
  arg   = mr  * phi
  argp  = mrp * phi
  argm  = mrm * phi
  Pi2im = Pi2 * im**mabs
  Piimp = Pi  * im**(mabs + 1)
  Piimm = Pi  * im**(mabs - 1)
  allocate (J1(0:mabs+1), J2(0:mabs+1))
  allocate (Pnm(0:Nrank), dPnm(0:Nrank), pinm(0:Nrank), taunm(0:Nrank))
  ind_refZc = ind_refZ
  if (abs(ind_refZc) < MachEps) ind_refZc = cmplx(MachEps,MachEps,O)
  ind_ref2  = ind_ref * ind_ref
  ind_refZ2 = ind_refZc * ind_refZc
  do p = 1, Nbeta
    beta = xbeta(p)
    cb   = cos(beta)
    sb   = sin(beta)
    f1   = (ind_refZ2 * cb**2 + ind_ref2 * sb**2) / ind_refZ2
    f2   = (ind_refZ2 - ind_ref2) * sb * cb / ind_refZ2
    f1s  = sqrt(f1)
    x1   = k1rsth * sb
    x2   = x1 / f1s
    y1   = k1rcth * cb
    y2   = y1 / f1s
    ex1  = exp(im * y1)
    ex2  = exp(im * y2)
    f12p =   f1 * cb + f2 * sb
    f12m = - f1 * sb + f2 * cb
    call bes_J (x1, mabs + 1, J1)
    call bes_J (x2, mabs + 1, J2)
    fact = - wbeta(p) * sb / 4._O / Pi
    if (m == 0) then
      call Leg_normalized (beta, 0, Nrank, Pnm, dPnm, pinm, taunm)
      Imx1  = Pi2 * J1(0)
      Imx2  = Pi2 * J2(0)
      ImCx1 = fcp * J1(1)
      ImCx2 = fcp * J2(1)
      ImSx1 = fsp * J1(1)
      ImSx2 = fsp * J2(1)
      o1 = Imx1  * ex1
      o2 = Imx2  * ex2
      s1 = ImSx1 * ex1
      c1 = ImCx1 * ex1
      s2 = ImSx2 * ex2
      c2 = ImCx2 * ex2
      do k = 1, Nmax
        n   = k
        nm  = real(2 * n * (n + 1),O)
        nm  = 1._O / sqrt(nm)
        ft  = taunm(n) * nm * (-im)**(n + 1) * fact
        f1  = im * ft
        Xe1 = - f1 * s1
        Xe2 =   f1 * c1
        Xe3 =   zero
        Xe(1,k) = Xe(1,k) + Xe1
        Xe(2,k) = Xe(2,k) + Xe2
        Xe(3,k) = Xe(3,k) + Xe3
!
        f2  = f12p * ft
        Ye1 = f2   * c2
        Ye2 = f2   * s2
        Ye3 = f12m * ft * o2
        Ye(1,k) = Ye(1,k) + Ye1
        Ye(2,k) = Ye(2,k) + Ye2
        Ye(3,k) = Ye(3,k) + Ye3
!
        f1  =   ft * cb
        Xh1 =   f1 * c1
        Xh2 =   f1 * s1
        Xh3 = - ft * sb * o1
        Xh(1,k) = Xh(1,k) + Xh1
        Xh(2,k) = Xh(2,k) + Xh2
        Xh(3,k) = Xh(3,k) + Xh3
!
        f2  =   f1s * im * ft
        Yh1 = - f2 * s2
        Yh2 =   f2 * c2
        Yh3 =   zero
        Yh(1,k) = Yh(1,k) + Yh1
        Yh(2,k) = Yh(2,k) + Yh2
        Yh(3,k) = Yh(3,k) + Yh3
      end do
    else
      call Leg_normalized (beta, mabs, Nrank, Pnm, dPnm, pinm, taunm)
      if (m > 0) then
        s = 1._O
      else
        s = - 1._O
      end if
      mlr   = real(m,O)
      f1    = exp(s * im * arg)
      fJ1   = f1 * J1(mabs)
      fJ2   = f1 * J2(mabs)
      f1    = exp(s * im * argp)
      fJ1p  = f1 * J1(mabs+1)
      fJ2p  = f1 * J2(mabs+1)
      f1    = exp(s * im * argm)
      fJ1m  = f1 * J1(mabs-1)
      fJ2m  = f1 * J2(mabs-1)
      f1    = Pi2im
      Imx1  = f1 * fJ1
      Imx2  = f1 * fJ2
      f1    = Piimp * fJ1p
      f2    = Piimm * fJ1m
      ImCx1 =   f1 + f2
      ImSx1 = - im * (s * f1 - s * f2)
      f1    = Piimp * fJ2p
      f2    = Piimm * fJ2m
      ImCx2 = f1 + f2
      ImSx2 = - im * (s * f1 - s * f2)
      o1 = Imx1  * ex1
      o2 = Imx2  * ex2
      s1 = ImSx1 * ex1
      c1 = ImCx1 * ex1
      s2 = ImSx2 * ex2
      c2 = ImCx2 * ex2
      do k = 1, Nmax
        n   = mabs + k - 1
        nm  = real(2 * n * (n + 1),O)
        nm  = 1._O / sqrt(nm)
        fp  = mlr * pinm(n) * nm * (-im)**(n + 1) * fact
        ft  = taunm(n) * nm * (-im)**(n + 1) * fact
        f1  = im * ft
        f2  = f12p * fp
        Xe1 = - f1 * s1 + f2 * c2
        Xe2 =   f1 * c1 + f2 * s2
        Xe3 =   f12m * fp * o2
        Xe(1,k) = Xe(1,k) + Xe1
        Xe(2,k) = Xe(2,k) + Xe2
        Xe(3,k) = Xe(3,k) + Xe3
!
        f1  =   im * fp
        f2  =   f12p * ft
        Ye1 = - f1 * s1 + f2 * c2
        Ye2 =   f1 * c1 + f2 * s2
        Ye3 =   f12m * ft * o2
        Ye(1,k) = Ye(1,k) + Ye1
        Ye(2,k) = Ye(2,k) + Ye2
        Ye(3,k) = Ye(3,k) + Ye3
!
        f1  =   ft * cb
        f2  =   f1s * im * fp
        Xh1 =   f1 * c1 - f2 * s2
        Xh2 =   f1 * s1 + f2 * c2
        Xh3 = - ft * sb * o1
        Xh(1,k) = Xh(1,k) + Xh1
        Xh(2,k) = Xh(2,k) + Xh2
        Xh(3,k) = Xh(3,k) + Xh3
!
        f1  =   fp * cb
        f2  =   f1s * im * ft
        Yh1 =   f1 * c1 - f2 * s2
        Yh2 =   f1 * s1 + f2 * c2
        Yh3 = - fp * sb * o1
        Yh(1,k) = Yh(1,k) + Yh1
        Yh(2,k) = Yh(2,k) + Yh2
        Yh(3,k) = Yh(3,k) + Yh3
      end do
    end if
  end do
! --- pass to spherical coordinates ---
  do k = 1, Nmax
    Xe1 = Xe(1,k)
    Xe2 = Xe(2,k)
    Xe3 = Xe(3,k)
    Xe(1,k) =   Xe1 * sthcph + Xe2 * sthsph + Xe3 * cth
    Xe(2,k) =   Xe1 * cthcph + Xe2 * cthsph - Xe3 * sth
    Xe(3,k) = - Xe1 * sph    + Xe2 * cph
!
    Ye1 = Ye(1,k)
    Ye2 = Ye(2,k)
    Ye3 = Ye(3,k)
    Ye(1,k) =   Ye1 * sthcph + Ye2 * sthsph + Ye3 * cth
    Ye(2,k) =   Ye1 * cthcph + Ye2 * cthsph - Ye3 * sth
    Ye(3,k) = - Ye1 * sph    + Ye2 * cph
!
    Xh1 = Xh(1,k)
    Xh2 = Xh(2,k)
    Xh3 = Xh(3,k)
    Xh(1,k) =   Xh1 * sthcph + Xh2 * sthsph + Xh3 * cth
    Xh(2,k) =   Xh1 * cthcph + Xh2 * cthsph - Xh3 * sth
    Xh(3,k) = - Xh1 * sph    + Xh2 * cph
!
    Yh1 = Yh(1,k)
    Yh2 = Yh(2,k)
    Yh3 = Yh(3,k)
    Yh(1,k) =   Yh1 * sthcph + Yh2 * sthsph + Yh3 * cth
    Yh(2,k) =   Yh1 * cthcph + Yh2 * cthsph - Yh3 * sth
    Yh(3,k) = - Yh1 * sph    + Yh2 * cph
  end do
  deallocate (J1, J2, xbeta, wbeta)
  deallocate (Pnm, dPnm, pinm, taunm)
end subroutine MN_m_anisotrop

max_stars_repo_name: johnaparker/MiePy

max_issues_repo_name: johnaparker/MiePy

max_forks_repo_name: johnaparker/MiePy


text:
*DECK CQRSL
      SUBROUTINE CQRSL (X, LDX, N, K, QRAUX, Y, QY, QTY, B, RSD, XB,
     +   JOB, INFO)
C***BEGIN PROLOGUE  CQRSL
C***PURPOSE  Apply the output of CQRDC to compute coordinate transfor-
C            mations, projections, and least squares solutions.
C***LIBRARY   SLATEC (LINPACK)
C***CATEGORY  D9, D2C1
C***TYPE      COMPLEX (SQRSL-S, DQRSL-D, CQRSL-C)
C***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX, ORTHOGONAL TRIANGULAR,
C             SOLVE
C***AUTHOR  Stewart, G. W., (U. of Maryland)
C***DESCRIPTION
C
C     CQRSL applies the output of CQRDC to compute coordinate
C     transformations, projections, and least squares solutions.
C     For K .LE. MIN(N,P), let XK be the matrix
C
C            XK = (X(JVPT(1)),X(JVPT(2)), ... ,X(JVPT(K)))
C
C     formed from columns JVPT(1), ... ,JVPT(K) of the original
C     N x P matrix X that was input to CQRDC (if no pivoting was
C     done, XK consists of the first K columns of X in their
C     original order).  CQRDC produces a factored unitary matrix Q
C     and an upper triangular matrix R such that
C
C              XK = Q * (R)
C                       (0)
C
C     This information is contained in coded form in the arrays
C     X and QRAUX.
C
C     On Entry
C
C        X      COMPLEX(LDX,P).
C               X contains the output of CQRDC.
C
C        LDX    INTEGER.
C               LDX is the leading dimension of the array X.
C
C        N      INTEGER.
C               N is the number of rows of the matrix XK.  It must
C               have the same value as N in CQRDC.
C
C        K      INTEGER.
C               K is the number of columns of the matrix XK.  K
C               must not be greater than (N,P), where P is the
C               same as in the calling sequence to CQRDC.
C
C        QRAUX  COMPLEX(P).
C               QRAUX contains the auxiliary output from CQRDC.
C
C        Y      COMPLEX(N)
C               Y contains an N-vector that is to be manipulated
C               by CQRSL.
C
C        JOB    INTEGER.
C               JOB specifies what is to be computed.  JOB has
C               the decimal expansion ABCDE, with the following
C               meaning.
C
C                    If A .NE. 0, compute QY.
C                    If B,C,D, or E .NE. 0, compute QTY.
C                    If C .NE. 0, compute B.
C                    If D .NE. 0, compute RSD .
C                    If E .NE. 0, compute  XB.
C
C               Note that a request to compute B, RSD, or XB
C               automatically triggers the computation of QTY, for
C               which an array must be provided in the calling
C               sequence.
C
C     On Return
C
C        QY     COMPLEX(N).
C               QY contains Q*Y, if its computation has been
C               requested.
C
C        QTY    COMPLEX(N).
C               QTY contains CTRANS(Q)*Y, if its computation has
C               been requested.  Here CTRANS(Q) is the conjugate
C               transpose of the matrix Q.
C
C        B      COMPLEX(K)
C               B contains the solution of the least squares problem
C
C                    minimize NORM2(Y - XK*B),
C
C               if its computation has been requested.  (Note that
C               if pivoting was requested in CQRDC, the J-th
C               component of B will be associated with column JVPT(J)
C               of the original matrix X that was input into CQRDC.)
C
C        RSD    COMPLEX(N).
C               RSD contains the least squares residual Y - XK*B,
C               if its computation has been requested.  RSD is
C               also the orthogonal projection of Y onto the
C               orthogonal complement of the column space of XK.
C
C        XB     COMPLEX(N).
C               XB contains the least squares approximation XK*B,
C               if its computation has been requested.  XB is also
C               the orthogonal projection of Y onto the column space
C               of X.
C
C        INFO   INTEGER.
C               INFO is zero unless the computation of B has
C               been requested and R is exactly singular.  In
C               this case, INFO is the index of the first zero
C               diagonal element of R and B is left unaltered.
C
C     The parameters QY, QTY, B, RSD, and XB are not referenced
C     if their computation is not requested and in this case
C     can be replaced by dummy variables in the calling program.
C     To save storage, the user may in some cases use the same
C     array for different parameters in the calling sequence.  A
C     frequently occurring example is when one wishes to compute
C     any of B, RSD, or XB and does not need Y or QTY.  In this
C     case one may identify Y, QTY, and one of B, RSD, or XB, while
C     providing separate arrays for anything else that is to be
C     computed.  Thus the calling sequence
C
C          CALL CQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO)
C
C     will result in the computation of B and RSD, with RSD
C     overwriting Y.  More generally, each item in the following
C     list contains groups of permissible identifications for
C     a single calling sequence.
C
C          1. (Y,QTY,B) (RSD) (XB) (QY)
C
C          2. (Y,QTY,RSD) (B) (XB) (QY)
C
C          3. (Y,QTY,XB) (B) (RSD) (QY)
C
C          4. (Y,QY) (QTY,B) (RSD) (XB)
C
C          5. (Y,QY) (QTY,RSD) (B) (XB)
C
C          6. (Y,QY) (QTY,XB) (B) (RSD)
C
C     In any group the value returned in the array allocated to
C     the group corresponds to the last member of the group.
C
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
C***ROUTINES CALLED  CAXPY, CCOPY, CDOTC
C***REVISION HISTORY  (YYMMDD)
C   780814  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  CQRSL
      INTEGER LDX,N,K,JOB,INFO
      COMPLEX X(LDX,*),QRAUX(*),Y(*),QY(*),QTY(*),B(*),RSD(*),XB(*)
C
      INTEGER I,J,JJ,JU,KP1
      COMPLEX CDOTC,T,TEMP
      LOGICAL CB,CQY,CQTY,CR,CXB
      COMPLEX ZDUM
      REAL CABS1
      CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM))
C***FIRST EXECUTABLE STATEMENT  CQRSL
C
C     SET INFO FLAG.
C
      INFO = 0
C
C     DETERMINE WHAT IS TO BE COMPUTED.
C
      CQY = JOB/10000 .NE. 0
      CQTY = MOD(JOB,10000) .NE. 0
      CB = MOD(JOB,1000)/100 .NE. 0
      CR = MOD(JOB,100)/10 .NE. 0
      CXB = MOD(JOB,10) .NE. 0
      JU = MIN(K,N-1)
C
C     SPECIAL ACTION WHEN N=1.
C
      IF (JU .NE. 0) GO TO 40
         IF (CQY) QY(1) = Y(1)
         IF (CQTY) QTY(1) = Y(1)
         IF (CXB) XB(1) = Y(1)
         IF (.NOT.CB) GO TO 30
            IF (CABS1(X(1,1)) .NE. 0.0E0) GO TO 10
               INFO = 1
            GO TO 20
   10       CONTINUE
               B(1) = Y(1)/X(1,1)
   20       CONTINUE
   30    CONTINUE
         IF (CR) RSD(1) = (0.0E0,0.0E0)
      GO TO 250
   40 CONTINUE
C
C        SET UP TO COMPUTE QY OR QTY.
C
         IF (CQY) CALL CCOPY(N,Y,1,QY,1)
         IF (CQTY) CALL CCOPY(N,Y,1,QTY,1)
         IF (.NOT.CQY) GO TO 70
C
C           COMPUTE QY.
C
            DO 60 JJ = 1, JU
               J = JU - JJ + 1
               IF (CABS1(QRAUX(J)) .EQ. 0.0E0) GO TO 50
                  TEMP = X(J,J)
                  X(J,J) = QRAUX(J)
                  T = -CDOTC(N-J+1,X(J,J),1,QY(J),1)/X(J,J)
                  CALL CAXPY(N-J+1,T,X(J,J),1,QY(J),1)
                  X(J,J) = TEMP
   50          CONTINUE
   60       CONTINUE
   70    CONTINUE
         IF (.NOT.CQTY) GO TO 100
C
C           COMPUTE CTRANS(Q)*Y.
C
            DO 90 J = 1, JU
               IF (CABS1(QRAUX(J)) .EQ. 0.0E0) GO TO 80
                  TEMP = X(J,J)
                  X(J,J) = QRAUX(J)
                  T = -CDOTC(N-J+1,X(J,J),1,QTY(J),1)/X(J,J)
                  CALL CAXPY(N-J+1,T,X(J,J),1,QTY(J),1)
                  X(J,J) = TEMP
   80          CONTINUE
   90       CONTINUE
  100    CONTINUE
C
C        SET UP TO COMPUTE B, RSD, OR XB.
C
         IF (CB) CALL CCOPY(K,QTY,1,B,1)
         KP1 = K + 1
         IF (CXB) CALL CCOPY(K,QTY,1,XB,1)
         IF (CR .AND. K .LT. N) CALL CCOPY(N-K,QTY(KP1),1,RSD(KP1),1)
         IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 120
            DO 110 I = KP1, N
               XB(I) = (0.0E0,0.0E0)
  110       CONTINUE
  120    CONTINUE
         IF (.NOT.CR) GO TO 140
            DO 130 I = 1, K
               RSD(I) = (0.0E0,0.0E0)
  130       CONTINUE
  140    CONTINUE
         IF (.NOT.CB) GO TO 190
C
C           COMPUTE B.
C
            DO 170 JJ = 1, K
               J = K - JJ + 1
               IF (CABS1(X(J,J)) .NE. 0.0E0) GO TO 150
                  INFO = J
                  GO TO 180
  150          CONTINUE
               B(J) = B(J)/X(J,J)
               IF (J .EQ. 1) GO TO 160
                  T = -B(J)
                  CALL CAXPY(J-1,T,X(1,J),1,B,1)
  160          CONTINUE
  170       CONTINUE
  180       CONTINUE
  190    CONTINUE
         IF (.NOT.CR .AND. .NOT.CXB) GO TO 240
C
C           COMPUTE RSD OR XB AS REQUIRED.
C
            DO 230 JJ = 1, JU
               J = JU - JJ + 1
               IF (CABS1(QRAUX(J)) .EQ. 0.0E0) GO TO 220
                  TEMP = X(J,J)
                  X(J,J) = QRAUX(J)
                  IF (.NOT.CR) GO TO 200
                     T = -CDOTC(N-J+1,X(J,J),1,RSD(J),1)/X(J,J)
                     CALL CAXPY(N-J+1,T,X(J,J),1,RSD(J),1)
  200             CONTINUE
                  IF (.NOT.CXB) GO TO 210
                     T = -CDOTC(N-J+1,X(J,J),1,XB(J),1)/X(J,J)
                     CALL CAXPY(N-J+1,T,X(J,J),1,XB(J),1)
  210             CONTINUE
                  X(J,J) = TEMP
  220          CONTINUE
  230       CONTINUE
  240    CONTINUE
  250 CONTINUE
      RETURN
      END

ext: f


program test_read_print
    use mpi_f08
    implicit none
    integer :: n_ranks, my_rank
    integer :: i

    call MPI_Init()
    call MPI_Comm_size(MPI_COMM_WORLD, n_ranks)
    call MPI_Comm_rank(MPI_COMM_WORLD, my_rank)

    read *, i
    print *, "my_rank=", my_rank, ", i=", i

    call MPI_Finalize()
end program test_read_print


text:
      PROGRAM DO
      REAL A(100)
c
      K = 24
      X = 1.1
      Y = 82.2
      DO I = 1,20
         IF ( SIN(X) .GT. 1. ) THEN
         X = X*X
         X = X+X
         ENDIF
         CALL P(X)
         IF ( X .GT. Y ) THEN
            X = X - 0.5
            Y = Y -0.1
         ELSE
            X = X + 0.5
            CALL P(X)
         ENDIF
         A(I) = 1.
         PRINT '(2F10.5)', X, Y
      ENDDO
      END
C
      SUBROUTINE P(Q)
      Q = Q + 1.
c      W = W - 1.
      END


text:
../../src/fp_format_cvt.v
../../src/fp_sum_block.v
../../src/int_sum_block.v
../../src/NV_NVDLA_CDP_DP_bufferin.v
../../src/NV_NVDLA_CDP_DP_cvtin.v
../../src/NV_NVDLA_CDP_DP_cvtout.v
../../src/NV_NVDLA_CDP_DP_INTP_unit.v
../../src/NV_NVDLA_CDP_DP_intp.v
../../src/NV_NVDLA_CDP_DP_LUT_CTRL_unit.v
../../src/NV_NVDLA_CDP_DP_LUT_ctrl.v
../../src/NV_NVDLA_CDP_DP_lut.v
../../src/NV_NVDLA_CDP_DP_MUL_unit.v
../../src/NV_NVDLA_CDP_DP_mul.v
../../src/NV_NVDLA_CDP_DP_nan.v
../../src/NV_NVDLA_CDP_DP_sum.v
../../src/NV_NVDLA_CDP_DP_syncfifo.v
../../src/NV_NVDLA_CDP_dp.v
../../src/NV_NVDLA_CDP_RDMA_cq.v
../../src/NV_NVDLA_CDP_RDMA_eg.v
../../src/NV_NVDLA_CDP_RDMA_ig.v
../../src/NV_NVDLA_CDP_RDMA_REG_dual.v
../../src/NV_NVDLA_CDP_RDMA_REG_single.v
../../src/NV_NVDLA_CDP_RDMA_reg.v
../../src/NV_NVDLA_CDP_rdma.v
../../src/NV_NVDLA_CDP_REG_dual.v
../../src/NV_NVDLA_CDP_REG_single.v
../../src/NV_NVDLA_CDP_reg.v
../../src/NV_NVDLA_CDP_slcg.v
../../src/NV_NVDLA_cdp.v
../../src/NV_NVDLA_CDP_wdma.v

max_stars_repo_name: animeshbchowdhury/robust-pnr-time

max_issues_repo_name: animeshbchowdhury/robust-pnr-time

max_forks_repo_name: animeshbchowdhury/robust-pnr-time


text: Shearer Tree Service provides pruning, planting, removals, stump grinding and consultation of Tree Services trees.  Davis resident since 1996.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module rlat_imat
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!       NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-----------------------------------------------------------------------
!
! !ROUTINE: rlat_imat - a lookup table of latitude for indexed matrices
!
! !INTERFACE:
!	use rlat_imat[, only : nveclat, veclats]
!
! !DESCRIPTION:
!	Module rlat_imat defines a lookup table of latitude values for
!	indexed matrix.
!
! !EXAMPLES: (to do)
! !BUGS: (to do)
! !SEE ALSO: (to do)
! !REVISION HISTORY:
! 	04Mar96 - J. Guo	- (to do)
!	19Sep96 - J. Guo	- add `save' to all public variables
!_______________________________________________________________________

  use config, only : MXveclat
  implicit none

  private	!  except ...
  public	:: nveclat, veclats
  public	:: MXveclat
  public        :: rlat_imat0

  integer,save :: nveclat		! actual size of veclats(:)
  real,   save :: veclats(MXveclat)	! a lookup table

Contains

  subroutine rlat_imat0
  end subroutine rlat_imat0

end module rlat_imat
!.

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
Davis Plastic Models is a campus student organizations organization that was founded to build plastic models. Although the majority of DPM is centered around Gundam as Gundam kits are the easiest to obtain, we do have a small but passionate community of people who make military kits.

Our goal is to provide a fun and comfortable environment for our builders to come hang out and build their model kits in the presence of their fellow peers. We also provide a service where people who are new to this hobby can seek assistance from our more experienced club members. This allows members to rapidly improve in ways they could have never achieved alone. Members are encouraged to give and receive help on different aspects of plastic models.

Although obtaining model kits and Gunpla may be difficult, Davis Plastic Models has a distributor, http://modelgrade.net/ Model Grade that offers us discounts on various products.

At the beginning of each quarter, we offer heavily discounted beginner Gunpla kits, and we teach people who may be new to the hobby the tips and tricks to being a good builder. Of course, existing members also have access to the group buy, and usually participate in the order as well. However, at the end of the quarter, we do a bulk purchase from Modelgrade. The group buy is discounted as well (usually 1015%) and comes in at the beginning of the following quarter. We use any extra remaining group funds to subsidize purchases as well.

Meetings are in Wellman 115, or whatever room in Wellman Hall that we happen to find open. In the case of three day weekends, well either be in one of the Tercero lounges or somewhere else on campus (wherever we can get in).

During meetings, we just build stuff and have fun. Usually every week well have a short tutorial on a different aspect of building, so newcomers will learn a lot. DPM also has a variety of basic tools (flush cutters, knives, etc) that are free of use to anyone, but paid members have priority if theres a dispute over who gets access. Bring your kits to club and have fun!

Membership is $5 per quarter, and is just to secure the funds to replace club equipment and purchase new stuff for the club. Whatever isnt spent is saved up to a cap, and then used to subsidize for the bulk order or new member kits. Members get priority on club tools, as well as the ability to participate in the bulk orders at the end of the quarter.

Occasionally DPM will also have activities outside of club. For instance, DPM will make occasional trips to Hobbytown in Vacaville to pick up painting supplies.


text:
      subroutine endout
c
c******************************************************************
c                                                                 *
c     Called from subroutine sedout.                              *
c     Summarizes soil loss information for the entire simulation  *
c     period.                                                     *
c                                                                 *
c******************************************************************
c
      include 'pmxpln.inc'
      include 'pmxtls.inc'
      include 'pmxtil.inc'
      include 'pmxhil.inc'
c
c******************************************************************
c                                                                 *
c   Common Blocks                                                 *
c                                                                 *
c******************************************************************
c
      include 'ccliyr.inc'
c
      include 'cstmflg.inc'
c
      include 'cstruc.inc'
c
      include 'csumirr.inc'
c
c*********************************************************************
c                                                                    *
c    sumirr variables udated                                         *
c    irrunt                                                          *
c                                                                    *
c*********************************************************************
c
      include 'csumout.inc'
c
      include 'cunicon.inc'
c     read outopt, units
c
      include 'cupdate.inc'
c
c******************************************************************
c                                                                 *
c   Local Variables                                               *
c     avernf :                                                    *
c     averun :                                                    *
c     unit(2):                                                    *
c     ii     :                                                    *
c     aveirr :                                                    *
c                                                                 *
c******************************************************************
c
      real avernf, averun, aveirr, sumh2o, pctirg,avemlt
      integer ii, ipl, jkl, sumevt
c
      character*5 unit(5)
      data unit /'   mm', 'kg/m2','  in.', 'in/yr','mm/yr'/
c
      ii = 1
c
      mxevnt = 0
      mxmelt = 0
c
      do 10 jkl = 1, nplane
        if (nrunot(jkl).gt.mxevnt) mxevnt = nrunot(jkl)
        if (nmunot(jkl).gt.mxmelt) mxmelt = nmunot(jkl)
   10 continue
c
      if (nirrt.eq.0) then
        write (31,1100) ii, nyear
c       end new
c
        if(outopt.eq.1.and.units.eq.1)then
c
c       write english or metric units for abbreviated average annual output file
c
          write (31,1300) nraint, traint/25.4, unit(3), mxevnt,
     1      trunot(nplane)/25.4, unit(3), mxmelt,
     1      tmunot(nplane)/25.4, unit(3)
        else
          write (31,1350) nraint, traint, mxevnt,
     1      trunot(nplane), mxmelt, tmunot(nplane)
        end if
      else
        write (31,1000) ii, nyear
c
        if(outopt.eq.1.and.units.eq.1)then
c
c       write english or metric units for abbreviated average annual output file
c
          write (31,1200) nraint, traint/25.4, nirrt,
     1      tirrt/25.4, mxevnt, trunot(nplane)/25.4,
     1      mxmelt, tmunot(nplane)/25.4, ncommt
        else
          write (31,1250) nraint, traint, nirrt, tirrt,
     1      mxevnt, trunot(nplane), mxmelt,
     1      tmunot(nplane), ncommt
        end if
      end if
c
      avernf = traint / float(jyear)
      averun = trunot(nplane) / float(jyear)
      avemlt = tmunot(nplane) / float(jyear)
c
      if (nirrt.eq.0) then
c
        if(outopt.eq.1.and.units.eq.1)then
c
c       write english or metric units for abbreviated average annual output file
c
          write (31,1600) jyear,avernf/25.4,unit(4),averun/25.4,unit(4),
     1      avemlt/25.4,unit(4)
        else
        write (31,1650) jyear,avernf,unit(1),averun,unit(1),avemlt,
     1   unit(1)
        end if
      else
        aveirr = tirrt / float(jyear)
c
        if(outopt.eq.1.and.units.eq.1)then
c
c       write english or metric units for abbreviated average annual output file
c
          write (31,1400) jyear,avernf/25.4,unit(4),aveirr/25.4,unit(4),
     1      averun/25.4,unit(4),avemlt/25.4,unit(4)
        else
          write (31,1450) jyear,avernf,unit(1),aveirr,unit(1),averun,
     1      unit(1),avemlt,unit(1)
        end if
c
        do 20 ipl = 1, nplane
          sumevt = nrunot(ipl) + nmunot(ipl)
          sumh2o = trunot(ipl) + tmunot(ipl)
          if (sumevt.gt.0.0) then
            pctirg = ((irrunt(ipl)*1000.0) / sumh2o ) * 100.0
          else
            pctirg = 0.0
          end if
          write (31,1500) ipl, pctirg
   20   continue
        irper=pctirg/100
c
      end if
c
      return
 1000 format (//'I.   RAINFALL, IRRIGATION, AND RUNOFF SUMMARY',/,5x,8(
     1    '-'),2x,10('-'),2x,3('-'),1x,6('-'),1x,7('-'),//,6x,
     1    'total summary: ',' years ',i4,' - ',i4)
 1100 format (//'I.   RAINFALL AND RUNOFF SUMMARY',/,5x,8('-'),1x,3('-'
     1    ),1x,6('-'),1x,7('-'),//,6x,'total summary: ',' years ',i4,
     1    ' - ',i4)
 1200 format (/5x,i5,
     1    ' storms produced                       ',f12.4,1x,
     1    ' in. of precipitation',/,5x,i5,
     1    ' irrigations supplied                  ',f12.4,1x,
     1    ' in. of water',/,5x,i5,
     1    ' rain/irrigation runoff events produced',f12.4,1x,
     1    ' in. of runoff',/,5x,i5,
     1    ' snow melts and/or',/,10x,
     1    '   events during winter produced       ',f12.4,1x,
     1    ' in. of runoff',/,5x,i5,
     1    ' irrigations occurred on days with rainfall',/)
 1250 format (/5x,i5,
     1    ' storms produced                       ',f9.2,
     1    ' mm of precipitation',/,5x,i5,
     1    ' irrigations supplied                  ',f9.2,
     1    ' mm of water',/,5x,i5,
     1    ' rain/irrigation runoff events produced',f9.2,
     1    ' mm of runoff',/,5x,i5,
     1    ' snow melts and/or',/,10x,
     1    '   events during winter produced       ',f9.2,
     1    ' mm of runoff',/,5x,i5,
     1    ' irrigations occurred on days with rainfall',/)
 1300 format(/5x,i5,
     1    ' storms produced                       ',f12.4,1x,a,
     1    ' of precipitation',/,5x,i5,
     1    ' rain storm runoff events produced     ',f12.4,1x,a,
     1    ' of runoff',/,5x,i5,
     1    ' snow melts and/or',/,10x,
     1    '   events during winter produced       ',f12.4,1x,a,
     1    ' of runoff',/)
 1350 format(/5x,i5,
     1    ' storms produced                       ',f9.2,
     1    ' mm of precipitation',/,5x,i5,
     1    ' rain storm runoff events produced     ',f9.2,
     1    ' mm of runoff',/,5x,i5,
     1    ' snow melts and/or',/,10x,
     1    '   events during winter produced       ',f9.2,
     1    ' mm of runoff',/)
 1400 format (6x,'annual averages'/6x,'---------------'//6x,
     1    '  Number of years                              ',
     1    3x,i4,/,6x,
     1    '  Mean annual precipitation                    ',
     1    f10.4,1x,a,/,6x,
     1    '  Mean annual irrigation                       ',
     1    f10.4,1x,a,/,6x,
     1    '  Mean annual runoff from rainfall/irrigation  ',
     1    f10.4,1x,a,/,6x,
     1    '  Mean annual runoff from snow melt',/,6x,
     1    '    and/or events during winter                ',
     1    f10.4,1x,a,//,6x,
     1    '      Overland     Fraction of Mean Annual Runoff'/6x,
     1    '    Flow Element      Attributed to Irrigation,'/6x,
     1    '       Number                     %')
 1450 format (6x,'annual averages'/6x,'---------------'//6x,
     1    '  Number of years                              ',
     1    3x,i4,/,6x,
     1    '  Mean annual precipitation                    ',
     1    f7.2,1x,a,/,6x,
     1    '  Mean annual irrigation                       ',
     1    f7.2,1x,a,/,6x,
     1    '  Mean annual runoff from rainfall/irrigation  ',
     1    f7.2,1x,a,/,6x,
     1    '  Mean annual runoff from snow melt',/,6x,
     1    '    and/or events during winter                ',
     1    f7.2,1x,a,//,6x,
     1    '      Overland     Fraction of Mean Annual Runoff'/6x,
     1    '    Flow Element      Attributed to Irrigation,'/6x,
     1    '       Number                     %')
 1500 format (15x,i2,21x,f5.1)
 1600 format (6x,'annual averages'/6x,'---------------'//6x,
     1    '  Number of years                              ',
     1    3x,i4,/,6x,
     1    '  Mean annual precipitation                    ',
     1    f10.4,1x,a,/,6x,
     1    '  Mean annual runoff from rainfall             ',
     1    f10.4,1x,a,/,6x,
     1    '  Mean annual runoff from snow melt',/,6x,
     1    '    and/or rain storm during winter            ',
     1    f10.4,1x,a/)
 1650 format (6x,'annual averages'/6x,'---------------'//6x,
     1    '  Number of years                              ',
     1    3x,i4,/,6x,
     1    '  Mean annual precipitation                    ',
     1    f7.2,1x,a,/,6x,
     1    '  Mean annual runoff from rainfall             ',
     1    f7.2,1x,a,/,6x,
     1    '  Mean annual runoff from snow melt',/,6x,
     1    '    and/or rain storm during winter            ',
     1    f7.2,1x,a/)
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
      subroutine tcorr(word,latt,pos,master,lfno)
      use tfstk
      use ffs
      use tffitcode
      use ffs_pointer, only:idelc,idvalc,idtypec
      use tfcsi,only:cssetp
      implicit real*8 (a-h,o-z)
      integer*8 latt(nlat),le,le1,ix,i1,i2
      dimension pos(nlat)
      integer*4 master(nlat)
      character*(*) word
      character*20 name
      logical exist,exist1,abb,temat,abbrev,add
      add=.false.
      call peekwd(word,next)
      if(abbrev(word,'A_DD','_'))then
        call cssetp(next)
        add=.true.
      endif
      call tfgetr(cl,1.d0,'DELCOR cor. length',lfno,exist)
      if(.not. exist)then
        go to 9990
      endif
      call tfgetr(dx,1.d0,'DELCOR delx',lfno,exist)
      if(.not. exist)then
        go to 9990
      endif
      call tfgetr(dy,1.d0,'DELCOR dely',lfno,exist)
      if(.not. exist)then
        go to 9990
      endif
      al=pos(nlat)
      if(.not. cell)then
        al=al*2.d0
      endif
      nc=128
      do while(nc .le. 32768 .and. nc*.5d0 .lt. al)
        nc=nc*2
      enddo
      ix=ktaloc(nc*2)
      call tvcorr(rlist(ix),cl,al,nc)
      exist1=.false.
 1    call peekwd(word,next)
      if(word .eq. ' ')then
        go to 9980
      endif
      abb=ifany(word,'*%{}',1) .gt. 0
      exist=.false.
      do 10 i=1,nlat-1
        id=idtypec(i)
        if(id .eq. 1 .or.
     1     (id .gt. 8 .and. id .ne. 20 .and.
     $       id .ne. 31 .and. id .ne. 41))then
          go to 10
        endif
        if(temat(i,name,word))then
          exist=.true.
          call cssetp(next)
          if(id .eq. 41)then
            ie=i
          else
            ie=master(i)
            if(ie .le. 0)then
              go to 10
            endif
          endif
          i0=i
          s=(pos(i0)+pos(ie+1))*.5d0
          ns=mod(int(nc*s/al)+nc/2,nc)+1
          ns1=mod(ns,nc)+1
          w=mod(nc*s/al+nc/2-(ns-1),1.d0)
          i1=ix+(ns-1)*2
          i2=ix+(ns1-1)*2
          delx=(rlist(i1)*(1.d0-w)+rlist(i2)*w)*dx
          dely=(rlist(i1+1)*(1.d0-w)+rlist(i2+1)*w)*dy
          if(id .eq. 20)then
            istep=ie-i0
          else
            istep=1
          endif
          do 110 j=i0,ie,istep
            if(master(j) .ne. 0 .or. id .eq. 41)then
              le=latt(j)+5
              le1=le+1
              dx1=delx
              dy1=dely
              if(id .eq. 2)then
                le=le+4
                le1=le+1
              elseif(id .eq. 20)then
                le=le-2
                le1=le+1
                if(.not. add)then
                  dx1=rlist(idvalc(j)+3)-delx
                  dy1=rlist(idvalc(j)+4)-dely
                else
                  dx1=-delx
                  dy1=-dely
                endif
              elseif(id .eq. 31)then
                le=le-2
                le1=le+1
                if(.not. add)then
                  dx1=rlist(idvalc(j)+13)-delx
                  dy1=rlist(idvalc(j)+14)-dely
                else
                  dx1=-delx
                  dy1=-dely
                endif
              elseif(id .eq. 41)then
                le=latt(j)+15
                le1=latt(j)+17
              endif
              if(add)then
                rlist(le)=rlist(le)+dx1
                rlist(le1)=rlist(le1)+dy1
              else
                rlist(le)=dx1
                rlist(le1)=dy1
              endif
              if(id .eq. 41 .and. i .eq. 1)then
                dxi=rlist(le)
                dyi=rlist(le1)
              endif
            endif
110       continue
          if(.not. abb)then
            go to 11
          endif
        endif
10    continue
11    call tfadjst(latt,pos)
      exist1=exist1 .or. exist
      if(exist)then
        go to 1
      else
        if(abb)then
          call cssetp(next)
          go to 1
        else
          if(.not. exist1)then
            call termes(lfno,'?Undefined element ',word)
          endif
          go to 9980
        endif
      endif
9980  call tfree(ix)
      return
9990  call termes(lfno,
     1'Syntax: DELCOR corlength delx dely element [element1...]',' ')
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD