! This module contains subroutines related to the dispersion relation ! for EC waves in a plasma. ! ! In particular, ! `colddisp` gives the analytical solution of the cold plasma ! dispersion; ! `warmdisp` solves numerically the warm plasma dispersion, either in ! full or in the weakly relativistic approximation; ! `harmnumber` computes the possible harmonic numbers and can be used ! to limit the number of calls to `warmdisp`. ! ! Note: the module must be initialised by calling `expinit` when using ! the fully relativistic model (argument `model` > 1). module dispersion use const_and_precisions, only : wp_, zero, one, half, two, & im, czero, cunit, pi, sqrt_pi implicit none ! global constants integer, parameter :: npts = 500 real(wp_), parameter :: tmax = 5.0_wp_ real(wp_), parameter :: dtex = 2*tmax/dble(npts) ! global variables real(wp_), dimension(npts+1), save :: ttv, extv private public expinit, colddisp, warmdisp public zetac, harmnumber contains pure function colddisp(X, Y, Npl, sox) result(npr) ! Cold plasma dispersion relation ! ! Given the parallel component of refractive index N∥ ! returns the orthogonal one N⊥. ! ! Reference: IFP-CNR Internal Report FP 05/1 - App. A ! subroutine arguments ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω real(wp_), intent(in) :: X, Y ! cold refractive index: N∥, N⊥ real(wp_), intent(in) :: Npl real(wp_) :: Npr ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X integer, intent(in) :: sox ! local variables real(wp_) :: Y2, Npl2, delta, dX ! The formula is: ! ! N⊥ = √(N² - N∥²) ! ! where ! N² = 1 - X - ½ XY² (1 + N∥² ± Δ)/(1 - X - Y²) ! Δ² = (1 - N∥²)² + 4N∥²(1 - X)/Y² ! Npl2 = Npl**2 Y2 = Y**2 dX = 1 - X delta = sqrt((1 - Npl2)**2 + 4*Npl2*dX/Y2) Npr = sqrt(dx - Npl2 - X*Y2 * (1 + Npl2 + sox*delta)/(dX - Y2)/2) end function colddisp pure subroutine harmnumber(Y, mu, Npl2, weakly, nhmin, nhmax) ! Computes the min/max possible harmonic numbers ! ! The condition for the wave-particle resonance to take place at ! the nth harmonic (including relativity and Doppler effects) is: ! ! γ(ω - k∥⋅v∥) = n ω_ce ! ! where ω_ce = eB/m, γ = 1/√(1 - β²). ! Defining the adimensional velocity u̅ = p̅/mc = γβ̅, the Lorentz ! factor can be written as ! ! γ² = 1 + γ²β² = 1 + u∥² + u⊥² ! ! Substituting this expression for γ, the condition can be rewritten ! as the equation for an ellipse in the (u∥, u⊥) plane: ! ! (u∥ - u∥₀)²/α∥² + u⊥²/α⊥² = 1 ! ! where u∥₀ = nYN∥/(1 - N∥²) ! α∥ = √(N∥² + n²Y² - 1)/(1 - N∥²) ! α⊥ = √(N∥² + n²Y² - 1)/√(1 - N∥²) ! Y = ω_ce/ω ! N̅ = ck̅/ω ! Assuming thermal equilibrium, the electrons are distributed ! according to a Boltzmann factor ! ! exp(-E/kT) = exp(-mc²(γ-1)/kT) ! = exp(-μ⋅√(1+u²)), (μ=mc²/kT) ! ! which is constant along circles centered at the origin: ! ! u∥² + u⊥² = u² ! ! A harmonic number is possible when its ellipse intersects a region ! containing a significant fraction of the electron population. ! ! subroutine arguments ! CMA Y variable: Y=ω_c/ω real(wp_), intent(in) :: Y ! squared parallel refractive index: N∥²=N²cos²θ real(wp_), intent(in) :: Npl2 ! reciprocal of adimensional temperature: μ=mc²/kT real(wp_), intent(in) :: mu ! whether to use the weakly-relativistic limit logical, intent(in) :: weakly ! min, max harmonic numbers integer, intent(out) :: nhmin, nhmax ! local variables integer :: nh, nhc real(wp_) :: Yc, Yn, gamma, rdu2, argexp, uu2 ! local constants ! Maximum value for μ⋅(γ-1) above which the ! distribution function is considered to be 0. real(wp_), parameter :: expcr = 16.0_wp_ nhmin = 0 nhmax = 0 ! Compute the critical Y value for which the first ! harmonic ellipse collapses into a single point. ! ! The condition is α∥⋅α⊥=0 (with n=1), from which: ! ! Yc = √(1 - N∥²) (fully relativistic) ! Yc = 1 -½ N∥² (weakly relativistic) ! if (weakly) then Yc = max(one - npl2/2, zero) else Yc = sqrt(max(1 - npl2, zero)) end if ! At the given Y and Yc, the "critical" harmonic ! numeber, i.e. the first non-collapsed, is ⌈Yc/Y⌉. ! Note: to see why consider the axes can be written ! ! α∥ = √(n²Y² - Yc²)/Yc² ! α⊥ = √(n²Y² - Yc²)/Yc ! nhc = ceiling(Yc/Y) ! Check a few numbers starting with nhc do nh = nhc, nhc + 10 Yn = Y*nh ! Find the closest point to the origin ! ! The solutions are simply: ! ! u⊥ = 0 ! u∥ = u∥₀ - sgn(N∥)⋅α∥ ! ! Note that N∥ can be negative, in which case u∥₀ is ! in the second quadrant and we must swap the sign. ! The value of γ in this point is ! ! γ² = 1 + u∥² ! = 1 + u∥₀² + α∥² - 2⋅sgn(N∥)⋅u∥₀α∥ ! ! The sign can be simplified by noting that ! ! sgn(N∥)⋅u∥₀ = sgn(N∥)⋅N∥ nY/(1 - N∥²) ! = |N∥|⋅Yn/Yc² ! ! Putting everything back together: ! ! γ = √[1 + u∥₀² + α∥² - 2⋅sgn(N∥)⋅u∥₀α∥] ! = √[Yc⁴ + N∥²Yn² + (Yn²-Yc²) - 2⋅Yn √(N∥²(Yn²-Yc²))]/Yc² ! = √[Yn² + N∥²(Yn² - Yc²) - 2⋅Yn √(N∥²(Yn²-Yc²))]/Yc² ! = |Yn - √(N∥²(Yn²-Yc²))|/Yc² ! if (weakly) then rdu2 = 2*(Yn - Yc) uu2 = npl2 + rdu2 - 2*sqrt(Npl2*rdu2) argexp = mu*uu2/2 else rdu2 = Yn**2 - Yc**2 gamma = (Yn - sqrt(Npl2*rdu2))/(1 - Npl2) argexp = mu*(gamma - one) end if if (argexp <= expcr) then ! The are enough electrons with this energy, ! update the max harmonic number nhmax = nh ! Also update the min if this is the first check if (nhmin == 0) nhmin = nh else ! Too few electrons, exit now as higher harmonics ! will have even fewer particles than the current if (nhmin > 0) exit end if end do end subroutine harmnumber subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & error, Npr, e, & model, nlarmor, max_iters, fallback) ! Solves the warm plasma dispersion relation using an iterative ! method starting from the cold plasma solution ! ! Reference: https://doi.org/10.13182/FST08-A1660 use, intrinsic :: ieee_arithmetic, only : ieee_is_finite use gray_errors, only : gray_error, warmdisp_convergence, warmdisp_result, & raise_error ! subroutine arguments ! Inputs ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω real(wp_), intent(in) :: X, Y ! reciprocal of adimensional temperature: μ=mc²/kT real(wp_), intent(in) :: mu ! solution of the cold plasma dispersion: N⊥ real(wp_), intent(in) :: Npr_cold ! parallel refractive index: N∥ real(wp_), intent(in) :: Npl ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X integer, intent(in) :: sox ! Outputs ! GRAY error code integer(kind=gray_error), intent(out) :: error ! orthogonal refractive index: N⊥ complex(wp_), intent(out) :: Npr ! polarization vector: ε (local cartesian frame, z∥B̅) complex(wp_), intent(out), optional :: e(3) ! Parameters ! switch for the dielectric tensor model: ! 1: weakly relativistic ! 2: fully relativistic (faster variant) ! 3: fully relativistic (slower variant) integer, intent(in) :: model ! order of the electron Larmor radius expansion integer, intent(in) :: nlarmor ! max number of iterations integer, intent(in) :: max_iters ! whether to fall back to the result of the first ! iteration in case the method fails to converge logical, intent(in) :: fallback ! local variables ! iteration counter integer :: i ! squares of the refractive index complex(wp_) :: Npr2 real(wp_) :: Npl2 ! full dielectric tensor: ε̅ complex(wp_) :: eps(3, 3) ! N⊥-independent parts: ε₃₃⁰, ε̅^(l) complex(wp_) :: e330, epsl(3, 3, nlarmor) ! N⊥, N⊥² from the previous iteration complex(wp_) :: Npr_prev, Npr2_prev ! N⊥² from the first iteration, needed when `fallback` is true complex(wp_) :: Npr2_fallback ! "error" estimated by the change in successive iterations real(wp_) :: Npr2_err ! change in N⊥ between successive iterations complex(wp_) :: delta_Npr ! change in complex argument of N⊥ between successive iterations real(wp_) :: angle, angle_norm ! factor for controlling the speed of the iterations ! Note: the factor is used as follow: ! ! N⊥²(i) ← N⊥²(i-1)⋅(1 - speed) + N⊥²(i)⋅speed ! real(wp_) :: speed ! Use the cold N⊥ as the initial "previous" N⊥, ! set the ΔN⊥ to a small value Npr2_prev = Npr_cold**2 delta_Npr = 1e-10_wp_ Npl2 = Npl**2 Npr2_err = 1 ! Compute the N⊥-independent part of the dielectric tensor call dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl, error) iterate_to_solve_disperion: do i = 1, max_iters Npr_prev = sqrt(Npr2_prev) ! Compute the "new" dielectric tensor ! ! ε̅ = Σ_l N⊥^(2l - 2) ε̅^(l) (eq. 10) compute_eps: block integer :: l eps = 0 do l = 1, nlarmor eps = eps + Npr2_prev**(l-1) * epsl(:,:,l) end do end block compute_eps ! Compute the speed factor if (i < 10) then speed = 1 elseif (i < 50) then ! Linear function that is 1 at i=10 and 0.1 at i=50 speed = 1 - 0.9_wp_*real(i - 10, wp_)/40_wp_ else speed = 0.1_wp_ end if ! Check if convergence is reached ! ! Make sure we do at least two steps in the iteration ! and the relative error is within tolerance if (i > 2 .and. Npr2_err < 1.0e-3_wp_*speed) exit solve_quadratic: block ! Solve the quadratic equation for N⊥² ! ! AN⊥⁴ + Bn⊥² + C = 0 (eq. 29) ! ! with complex coefficients A,B,C. integer :: s complex(wp_) :: A, B, C, delta A = (eps(1,1) - Npl2)*(1 - eps(3,3)) + (eps(1,3) + Npl)**2 B = -(1 - eps(3,3)) * ( (eps(1,1) - Npl2)*(eps(2,2) - Npl2) + eps(1,2)**2 ) & -e330*(eps(1,1) - Npl2) + 2*eps(1,2)*eps(2,3) * (eps(1,3) + Npl) & +eps(2,3)**2 * (eps(1,1) - Npl2) - (eps(1,3) + Npl)**2 * (eps(2,2) - Npl2) C = e330 * ((eps(1,1) - Npl2)*(eps(2,2) - Npl2) + eps(1,2)**2) ! polynomial discriminant delta = B**2 - 4*A*C ! Choose the sign of the solution if (Y > 1) then s = sox if (delta%im <= 0) s = -s else s = -sox if (delta%re <= 0 .and. delta%im >= 0) s = -s end if ! Solution Npr2 = (-B + s*sqrt(delta)) / (2*A) end block solve_quadratic ! A new value for Npr2 has been established: estimate the error ! and store the new value as the latest prediction, Npr2_prev Npr2_err = abs(1 - abs(Npr2)/abs(Npr2_prev)) ! Apply a sequence of modification to N⊥ to improve the stability modify_fixed_point: block ! fixed-point prediction of N⊥ complex(wp_) :: Npr_fix, change ! Save the original fixed-point iteration Npr_fix = sqrt(Npr2) ! If real and imaginary parts of the refractive index have opposite ! sign, then the wave is a backward wave, typically a Bernstein ! branch (or we have an instability). We then have to find our way ! back to the EM branch, or at least make sure this wave number is ! as emission. if (Npr_fix%re * Npr_fix%im < 0) Npr2 = conjg(Npr_fix)**2 ! Change the speed when the current/previous step go out of phase if (i > 1) then change = (Npr_prev - sqrt(Npr2)) / delta_Npr angle = atan2(change%im, change%re) angle_norm = modulo(angle/pi + 1, two) - 1 if (abs(angle_norm) > 0.4_wp_) then ! Change in speed is linear from 1.0 at abs(angle)=0.4*angle ! to 0.3 at angle=π speed = speed * (1 - 0.7_wp_ * (abs(angle_norm) - 0.4_wp_) / (1 - 0.4_wp_)) end if end if ! In the iteration process we do not want huge steps; do not allow steps > 0.05 Npr2 = Npr2_prev*(1 - speed) + Npr2*speed if (abs(Npr2 - Npr2_prev) > 0.05_wp_) then Npr2 = Npr2_prev + 0.05_wp_ * (Npr2 - Npr2_prev) / abs(Npr2 - Npr2_prev) ! Again, make sure that we have a damped EM wave and not a ! Bernstein-like wave (see above) if (real(sqrt(Npr2)) * aimag(sqrt(Npr2)) < zero) & Npr2 = conjg(sqrt(Npr2))**2 end if end block modify_fixed_point ! Store the value of N⊥² and its change between iterations Npr2_prev = Npr2 delta_Npr = Npr_prev - sqrt(Npr2_prev) ! Store the result of the first iteration as a fallback ! in case the method fails to converge if (i == 1 .and. fallback) Npr2_fallback = Npr2 end do iterate_to_solve_disperion ! Handle convergence failures if (i > max_iters) then if (fallback) then ! first try failed, use fallback error = raise_error(error, warmdisp_convergence, 0) Npr2 = Npr2_fallback else ! first try failed, no retry error = raise_error(error, warmdisp_convergence, 1) end if end if ! Catch NaN and ±Infinity values if (.not. ieee_is_finite(Npr2%re) .or. .not. ieee_is_finite(Npr2%im)) then Npr2 = 0 error = raise_error(error, warmdisp_result, 0) end if if (Npr2%re < 0 .and. Npr2%im < 0) then Npr2 = sqrt(0.5_wp_) ! Lorenzo, what should we have here?? error = raise_error(error, warmdisp_result, 1) end if ! Final result Npr = sqrt(Npr2) if (.not. present(e)) return ! Compute the polarization vector ! TODO: find a reference for the formula polarization: block complex(wp_) :: den e = 0 if (abs(Npl) > 1.0e-6_wp_) then den = Npr*eps(1,2)*eps(2,3) - Npr*(eps(1,3) + Npl)*(eps(2,2) - Npr2 - Npl2) e(2) = -(eps(1,2)*Npr*(eps(1,3) + Npl) + (eps(1,1) - Npl2)* Npr*eps(2,3))/den e(3) = (eps(1,2)**2 + (eps(2,2) - Npr2 - Npl2)*(eps(1,1) - Npl2))/den e(1) = sqrt(1/(1 + abs(e(3))**2 + abs(e(2))**2)) e(2) = e(2)*e(1) e(3) = e(3)*e(1) else if (sox < 0) then e(3) = 1 else e(1) = sqrt(one/(1 + abs(-eps(1,1)/eps(1,2))**2)) e(2) = -e(1)*eps(1,1)/eps(1,2) end if end if end block polarization end subroutine warmdisp subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl, error) ! Returns the N⊥-independent parts of the dielectric tensor, namely ! ε₃₃⁰ and the 3x3 tensors ε̅^(l) for l=1,nlarmor. ! ! Reference: https://doi.org/10.13182/FST08-A1660 use gray_errors, only : gray_error ! subroutine arguments ! Inputs ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω real(wp_), intent(in) :: X, Y ! reciprocal of adimensional temperature: μ=mc²/kT real(wp_), intent(in) :: mu ! parallel refractive index: N∥ real(wp_), intent(in) :: Npl ! Parameters ! switch for the dielectric tensor model: ! 1: weakly relativistic ! 2: fully relativistic (faster variant) ! 3: fully relativistic (slower variant) integer, intent(in) :: model ! order of the electron Larmor radius expansion integer, intent(in) :: nlarmor ! Outputs ! element ε₃₃⁰, defined in eq. 9 complex(wp_), intent(out) :: e330 ! tensors ε̅_ij^(l), defined in eq. 10 complex(wp_), intent(out) :: epsl(3,3,nlarmor) ! GRAY error code integer(kind=gray_error), intent(out) :: error ! local variables integer :: l, n real(wp_) :: Npl2, alpha_nl, beta_nl, A_l complex(wp_) :: Qp_0l, Qm_0l, Qp_1l, Qm_1l, Qp_2l ! weakly relativistic functions complex(wp_), dimension(0:nlarmor,0:2) :: Fp, Fm ! fully relativistic functions real(wp_) :: rr(-nlarmor:nlarmor, 0:2, 0:nlarmor) real(wp_) :: ri(nlarmor, 0:2, nlarmor) Npl2 = Npl**2 if (model == 1) then ! Compute the (combinations of) Shkarofsky functions F_q(n) ! in which the functions Q⁺_hl(n), Q⁻_hl(n) can be expressed call fsup(nlarmor, Y, Npl, mu, Fp, Fm, error) else ! Compute the real and imaginary parts of the integrals ! defining the functions Q_hl(n), Q_hl(-n) block real(wp_) :: cr, ci, cmxw cmxw = 1 + 15/(8*mu) + 105/(128 * mu**2) cr = -mu**2/(sqrt_pi * cmxw) ci = sqrt(2*pi*mu) * mu**2/cmxw ! real part if (model == 2) then ! fast variant call hermitian(rr, Y, mu, Npl, cr, model, nlarmor) else ! slow variant call hermitian_2(rr, Y, mu, Npl, cr, model, nlarmor, error) if (error /= 0) error = 2 end if ! imaginary part call antihermitian(ri, Y, mu, Npl, ci, nlarmor) end block end if block ! Compute all factorials incrementally... integer :: fact_l ! l! integer :: fact_2l ! (2l)! integer :: fact_2l_l ! (2l)! / l! integer :: fact_l2 ! (l!)² integer :: fact_lpn ! (l+n)! integer :: fact_lmn ! (l-n)! ! ...starting from 1 fact_l = 1 fact_2l = 1 fact_2l_l = 1 fact_l2 = 1 epsl = 0 do l = 1, nlarmor fact_l = fact_l * l ! Coefficient A_l, eq. 11 if (model == 1) then fact_2l_l = fact_2l_l * (4*l - 2) ! see http://oeis.org/A001813 A_l = (1/(2 * Y**2 * mu))**(l - 1)/2 * fact_2l_l else fact_2l = fact_2l * 2*l * (2*l - 1) fact_l2 = fact_l2 * l**2 A_l = - fact_2l/(4**l * Y**(2*(l - 1)) * fact_l2) end if ! When n=0, (l+n)! = (l-n)! = l! fact_lpn = fact_l fact_lmn = fact_l do n = 0, l ! Coefficients α_nl and β_nl, eq. 12 alpha_nl = one * (-1)**(l - n)/(fact_lpn * fact_lmn) beta_nl = alpha_nl*(n**2 + one * (2*(l - n)*(l + n)*(l - 1)) / (2*l - 1)) ! Compute the next factorial values fact_lpn = fact_lpn * (l + n+1) if (n0) then Qp_0l = Qp_0l + im*ri(n,0,l) + rr(-n,0,l) Qm_0l = Qm_0l + im*ri(n,0,l) - rr(-n,0,l) Qp_1l = Qp_1l + im*ri(n,1,l) + rr(-n,1,l) Qm_1l = Qm_1l + im*ri(n,1,l) - rr(-n,1,l) Qp_2l = Qp_2l + im*ri(n,2,l) + rr(-n,2,l) end if end if ! Components of the ε̅^(l) tensors, eq. 11 epsl(1,1,l) = epsl(1,1,l) - n**2 * alpha_nl * Qp_0l epsl(1,2,l) = epsl(1,2,l) + im * l * n*alpha_nl * Qm_0l epsl(2,2,l) = epsl(2,2,l) - beta_nl * Qp_0l epsl(1,3,l) = epsl(1,3,l) - n*alpha_nl * Qm_1l / Y epsl(2,3,l) = epsl(2,3,l) - im * l*alpha_nl * Qp_1l / Y epsl(3,3,l) = epsl(3,3,l) - alpha_nl * Qp_2l / Y**2 end do ! Add the common factor A_l epsl(:,:,l) = epsl(:,:,l) * A_l end do end block ! Multiply by X and add the Kroneker δ_l,1 to diagonal elements epsl = epsl * X epsl(1,1,1) = epsl(1,1,1) + 1 epsl(2,2,1) = epsl(2,2,1) + 1 ! The special ε₃₃⁰ element if (model == 1) then Qp_2l = mu*Fp(0,1) + mu**2 * Npl2*(Fp(0,2) + Fp(0,0) - 2*Fp(0,1)) else Qp_2l = -rr(0,2,0) end if e330 = 1 - X * Qp_2l ! Other components given by symmetry epsl(2,1,:) = - epsl(1,2,:) epsl(3,1,:) = + epsl(1,3,:) epsl(3,2,:) = - epsl(2,3,:) end subroutine dielectric_tensor subroutine hermitian(rr,yg,mu,npl,cr,fast,lrm) use eierf, only : calcei3 ! subroutine arguments integer :: lrm,fast real(wp_) :: yg,mu,npl,cr real(wp_), dimension(-lrm:lrm,0:2,0:lrm) :: rr ! local variables integer :: i,k,n,n1,nn,m,llm real(wp_) :: mu2,mu4,mu6,npl2,npl4,bth,bth2,bth4,bth6,t,rxt,upl2, & upl,gx,exdx,x,gr,s,zm,zm2,zm3,fe0m,ffe,sy1,sy2,sy3 ! do n=-lrm,lrm do k=0,2 do m=0,lrm rr(n,k,m)=zero end do end do end do ! llm=min(3,lrm) ! bth2=2.0_wp_/mu bth=sqrt(bth2) mu2=mu*mu mu4=mu2*mu2 mu6=mu4*mu2 ! do i = 1, npts+1 t = ttv(i) rxt=sqrt(one+t*t/(2.0_wp_*mu)) x = t*rxt upl2=bth2*x**2 upl=bth*x gx=one+t*t/mu exdx=cr*extv(i)*gx/rxt*dtex ! n1=1 if(fast.gt.2) n1=-llm ! do n=n1,llm nn=abs(n) gr=npl*upl+n*yg zm=-mu*(gx-gr) s=mu*(gx+gr) zm2=zm*zm zm3=zm2*zm call calcei3(zm,fe0m) ! do m=nn,llm if(n.eq.0.and.m.eq.0) then rr(0,2,0) = rr(0,2,0) - exdx*fe0m*upl2 else if (m.eq.1) then ffe=(one+s*(one-zm*fe0m))/mu2 else if (m.eq.2) then ffe=(6.0_wp_-2.0_wp_*zm+4.0_wp_*s+s*s*(one+zm-zm2*fe0m))/mu4 else ffe=(18.0_wp_*s*(s+4.0_wp_-zm)+6.0_wp_* & (20.0_wp_-8.0_wp_*zm+zm2)+s**3*(2.0_wp_+zm+zm2-zm3*fe0m))/mu6 end if ! rr(n,0,m) = rr(n,0,m) + exdx*ffe rr(n,1,m) = rr(n,1,m) + exdx*ffe*upl rr(n,2,m) = rr(n,2,m) + exdx*ffe*upl2 ! end if ! end do end do end do ! if(fast.gt.2) return ! sy1=one+yg sy2=one+yg*2.0_wp_ sy3=one+yg*3.0_wp_ ! bth4=bth2*bth2 bth6=bth4*bth2 ! npl2=npl*npl npl4=npl2*npl2 ! rr(0,2,0) = -(one+bth2*(-1.25_wp_+1.5_wp_*npl2) & +bth4*(1.71875_wp_-6.0_wp_*npl2+3.75_wp_*npl2*npl2) & +bth6*3.0_wp_*(-65.0_wp_+456.0_wp_*npl2-660.0_wp_*npl4 & +280.0_wp_*npl2*npl4)/64.0_wp_+bth6*bth2*15.0_wp_ & *(252.853e3_wp_-2850.816e3_wp_*npl2+6942.720e3_wp_*npl4 & -6422.528e3_wp_*npl4*npl2+2064.384e3_wp_*npl4*npl4) & /524.288e3_wp_) rr(0,1,1) = -npl*bth2*(one+bth2*(-2.25_wp_+1.5_wp_*npl2) & +bth4*9.375e-2_wp_*(6.1e1_wp_-9.6e1_wp_*npl2+4.e1_wp_*npl4 & +bth2*(-184.5_wp_+4.92e2_wp_*npl2-4.5e2_wp_*npl4 & +1.4e2_wp_*npl2*npl4))) rr(0,2,1) = -bth2*(one+bth2*(-0.5_wp_+1.5_wp_*npl2)+0.375_wp_*bth4 & *(3.0_wp_-15.0_wp_*npl2+10.0_wp_*npl4)+3.0_wp_*bth6 & *(-61.0_wp_+471.0_wp_*npl2-680*npl4+280.0_wp_*npl2*npl4) & /64.0_wp_) rr(-1,0,1) = -2.0_wp_/sy1*(one+bth2/sy1*(-1.25_wp_+0.5_wp_*npl2/sy1) & +bth4/sy1*(-0.46875_wp_+(2.1875_wp_+0.625_wp_*npl2)/sy1 & -2.625_wp_*npl2/sy1**2+0.75_wp_*npl4/sy1**3)+bth6/sy1 & *(0.234375_wp_+(1.640625_wp_+0.234375_wp_*npl2)/sy1 & +(-4.921875_wp_-4.921875_wp_*npl2)/sy1**2 & +2.25_wp_*npl2*(5.25_wp_+npl2)/sy1**3 - 8.4375_wp_*npl4/sy1**4 & +1.875_wp_*npl2*npl4/sy1**5)+bth6*bth2/sy1*(0.019826889038*sy1 & -0.06591796875_wp_+(-0.7177734375_wp_ - 0.1171875_wp_*npl2)/sy1 & +(-5.537109375_wp_ - 2.4609375_wp_*npl2)/sy1**2 & +(13.53515625_wp_ + 29.53125_wp_*npl2 + 2.8125_wp_*npl4)/sy1**3 & +(-54.140625_wp_*npl2 - 32.6953125_wp_*npl4)/sy1**4 & +(69.609375_wp_*npl4 + 9.84375_wp_*npl2*npl4)/sy1**5 & -36.09375_wp_*npl2*npl4/sy1**6 + 6.5625_wp_*npl4**2/sy1**7)) rr(-1,1,1) = -npl*bth2/sy1**2*(one+bth2*(1.25_wp_-3.5_wp_/sy1 & +1.5_wp_*npl2/sy1**2)+bth4*9.375e-2_wp_*((5.0_wp_-7.e1_wp_/sy1 & +(126.0_wp_+48.0_wp_*npl2)/sy1**2-144.0_wp_*npl2/sy1**3 & +4.e1_wp_*npl4/sy1**4)+bth2*(-2.5_wp_-3.5e1_wp_/sy1+(3.15e2_wp_ & +6.e1_wp_*npl2)/sy1**2+(-4.62e2_wp_-5.58e2_wp_*npl2)/sy1**3 & +(9.9e2_wp_*npl2+2.1e2_wp_*npl4)/sy1**4-6.6e2_wp_*npl4/sy1**5+ & 1.4e2_wp_*npl4*npl2/sy1**6))) rr(-1,2,1) = -bth2/sy1*(one+bth2*(1.25_wp_-1.75_wp_/sy1+1.5_wp_*npl2/sy1**2) & +bth4*3.0_wp_/32.0_wp_*(5.0_wp_-35.0_wp_/sy1 & +(42.0_wp_+48.0_wp_*npl2)/sy1**2-108.0_wp_*npl2/sy1**3 & +40.0_wp_*npl4/sy1**4+0.5_wp_*bth2*(-5.0_wp_-35.0_wp_/sy1 & +(21.e1_wp_+12.e1_wp_*npl2)/sy1**2-(231.0_wp_+837.0_wp_*npl2) & /sy1**3+12.0_wp_*npl2*(99.0_wp_+35.0_wp_*npl2)/sy1**4 & -1100.0_wp_*npl4/sy1**5 + 280.0_wp_*npl2*npl4/sy1**6))) ! if(llm.gt.1) then ! rr(0,0,2) = -4.0_wp_*bth2*(one+bth2*(-0.5_wp_+0.5_wp_*npl2)+bth4 & *(1.125_wp_-1.875_wp_*npl2+0.75_wp_*npl4)+bth6*3.0_wp_ & *(-61.0_wp_+157.0_wp_*npl2-136.0_wp_*npl4+40.0_wp_*npl2*npl4) & /64.0_wp_) rr(0,1,2) = -2.0_wp_*npl*bth4*(one+bth2*(-1.5_wp_+1.5_wp_*npl2)+bth4 & *(39.0_wp_-69.0_wp_*npl2+30.0_wp_*npl4)/8.0_wp_) rr(0,2,2) = -2.0_wp_*bth4*(one+bth2*(0.75_wp_+1.5_wp_*npl2)+bth4* & (13.0_wp_-48.0_wp_*npl2 +40.0_wp_*npl4)*3.0_wp_/32.0_wp_) rr(-1,0,2) = -4.0_wp_*bth2/sy1*(one+bth2* & (1.25_wp_-1.75_wp_/sy1+0.5_wp_*npl2/sy1**2)+bth4* & (0.46875_wp_-3.28125_wp_/sy1+(3.9375_wp_+1.5_wp_*npl2)/sy1**2 & -3.375_wp_*npl2/sy1**3+0.75_wp_*npl4/sy1**4) & +bth4*bth2*3.0_wp_/64.0_wp_*(-5.0_wp_-35.0_wp_/sy1 & +(210.0_wp_+40.0_wp_*npl2)/sy1**2-3.0_wp_* & (77.0_wp_+93.0_wp_*npl2)/sy1**3+(396.0*npl2+84.0_wp_*npl4) & /sy1**4-220.0_wp_*npl4/sy1**5+40.0_wp_*npl4*npl2/sy1**6)) rr(-1,1,2) = -2.0_wp_*bth4*npl/sy1**2*(one+bth2 & *(3.0_wp_-4.5_wp_/sy1+1.5_wp_*npl2/sy1**2)+bth4 & *(20.0_wp_-93.0_wp_/sy1+(99.0_wp_+42.0_wp_*npl2)/sy1**2 & -88.0_wp_*npl2/sy1**3+20.0_wp_*npl4/sy1**4)*3.0_wp_/16.0_wp_) rr(-1,2,2) = -2.0_wp_*bth4/sy1*(one+bth2 & *(3.0_wp_-2.25_wp_/sy1+1.5_wp_*npl2/sy1**2)+bth4 & *(40.0_wp_*npl4/sy1**4-132.0_wp_*npl2/sy1**3 & +(66.0_wp_+84.0_wp_*npl2)/sy1**2-93.0_wp_/sy1+40.0_wp_) & *3.0_wp_/32.0_wp_) rr(-2,0,2) = -4.0_wp_*bth2/sy2*(one+bth2 & *(1.25_wp_-1.75_wp_/sy2+0.5_wp_*npl2/sy2**2)+bth4 & *(0.46875_wp_-3.28125_wp_/sy2+(3.9375_wp_+1.5_wp_*npl2) & /sy2**2-3.375_wp_*npl2/sy2**3+0.75_wp_*npl4/sy2**4)+bth4*bth2 & *3.0_wp_/64.0_wp_*(-5.0_wp_-35.0_wp_/sy2 & +(210.0_wp_+40.0_wp_*npl2)/sy2**2-3.0_wp_ & *(77.0_wp_+93.0_wp_*npl2)/sy2**3 & +(396.0*npl2+84.0_wp_*npl4)/sy2**4-220.0_wp_*npl4/sy2**5 & +40.0_wp_*npl4*npl2/sy2**6)) rr(-2,1,2) =-2.0_wp_*bth4*npl/sy2**2*(one+bth2 & *(3.0_wp_-4.5_wp_/sy2+1.5_wp_*npl2/sy2**2)+bth4 & *(20.0_wp_-93.0_wp_/sy2+(99.0_wp_+42.0_wp_*npl2)/sy2**2 & -88.0_wp_*npl2/sy2**3+20.0_wp_*npl4/sy2**4)*3.0_wp_/16.0_wp_) rr(-2,2,2) = -2.0_wp_*bth4/sy2*(one+bth2 & *(3.0_wp_-2.25_wp_/sy2+1.5_wp_*npl2/sy2**2)+bth4 & *(40.0_wp_*npl4/sy2**4-132.0_wp_*npl2/sy2**3 & +(66.0_wp_+84.0_wp_*npl2)/sy2**2-93.0_wp_/sy2+40.0_wp_) & *3.0_wp_/32.0_wp_) ! if(llm.gt.2) then ! rr(0,0,3) = -12.0_wp_*bth4*(one+bth2*(0.75_wp_+0.5_wp_*npl2)+bth4 & *(1.21875_wp_-1.5_wp_*npl2+0.75_wp_*npl2*npl2)) rr(0,1,3) = -6.0_wp_*npl*bth6*(1+bth2*(-0.25_wp_+1.5_wp_*npl2)) rr(0,2,3) = -6.0_wp_*bth6*(one+bth2*(2.5_wp_+1.5_wp_*npl2)) rr(-1,0,3) = -12.0_wp_*bth4/sy1*(one+bth2 & *(3.0_wp_-2.25_wp_/sy1+0.5_wp_*npl2/sy1**2)+bth4 & *(3.75_wp_-8.71875_wp_/sy1+(6.1875_wp_+2.625_wp_*npl2) & /sy1**2-4.125_wp_*npl2/sy1**3+0.75*npl2*npl2/sy1**4)) rr(-1,1,3) = -6.0_wp_*npl*bth6/sy1**2* & (one+bth2*(5.25_wp_-5.5_wp_/sy1+1.5_wp_*npl2/sy1**2)) rr(-1,2,3) = -6.0_wp_*bth6/sy1* & (one+bth2*(5.25_wp_-2.75_wp_/sy1+1.5_wp_*npl2/sy1**2)) ! rr(-2,0,3) = -12.0_wp_*bth4/sy2 & *(one+bth2*(3.0_wp_-2.25_wp_/sy2+0.5_wp_*npl2/sy2**2) & +bth4*(3.75_wp_-8.71875_wp_/sy2+(6.1875_wp_+2.625_wp_*npl2) & /sy2**2-4.125_wp_*npl2/sy2**3+0.75*npl2*npl2/sy2**4)) rr(-2,1,3) = -6.0_wp_*npl*bth6/sy2**2 & *(one+bth2*(5.25_wp_-5.5_wp_/sy2+1.5_wp_*npl2/sy2**2)) rr(-2,2,3) = -6.0_wp_*bth6/sy2 & *(one+bth2*(5.25_wp_-2.75_wp_/sy2+1.5_wp_*npl2/sy2**2)) ! rr(-3,0,3) = -12.0_wp_*bth4/sy3 & *(one+bth2*(3.0_wp_-2.25_wp_/sy3+0.5_wp_*npl2/sy3**2) & +bth4*(3.75_wp_-8.71875_wp_/sy3+(6.1875_wp_+2.625_wp_*npl2) & /sy3**2-4.125_wp_*npl2/sy3**3+0.75*npl2*npl2/sy3**4)) rr(-3,1,3) = -6.0_wp_*npl*bth6/sy3**2 & *(one+bth2*(5.25_wp_-5.5_wp_/sy3+1.5_wp_*npl2/sy3**2)) rr(-3,2,3) = -6.0_wp_*bth6/sy3 & *(one+bth2*(5.25_wp_-2.75_wp_/sy3+1.5_wp_*npl2/sy3**2)) ! end if end if ! end subroutine hermitian ! ! ! subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm,error) use gray_errors, only : gray_error, dielectric_tensor, raise_error use quadpack, only : dqagsmv ! local constants integer,parameter :: lw=5000,liw=lw/4,npar=7 real(wp_), parameter :: epsa=zero,epsr=1.0e-4_wp_ ! arguments integer :: lrm,fast real(wp_) :: yg,mu,npl,cr real(wp_), dimension(-lrm:lrm,0:2,0:lrm) :: rr integer(kind=gray_error), intent(out) :: error ! local variables integer :: n,m,ih,k,n1,nn,llm,neval,ier,last,ihmin integer, dimension(liw) :: iw real(wp_) :: mu2,mu4,mu6,npl2,bth,bth2,bth4,bth6 real(wp_) :: sy1,sy2,sy3,resfh,epp real(wp_), dimension(lw) :: w real(wp_), dimension(npar) :: apar ! do n=-lrm,lrm do k=0,2 do m=0,lrm rr(n,k,m)=zero end do end do end do ! llm=min(3,lrm) ! bth2=2.0_wp_/mu bth=sqrt(bth2) mu2=mu*mu mu4=mu2*mu2 mu6=mu4*mu2 ! n1=1 if(fast.gt.10) n1=-llm ! apar(1) = yg apar(2) = mu apar(3) = npl apar(4) = cr ! do n=n1,llm nn=abs(n) apar(5) = real(n,wp_) do m=nn,llm apar(6) = real(m,wp_) ihmin=0 if(n.eq.0.and.m.eq.0) ihmin=2 do ih=ihmin,2 apar(7) = real(ih,wp_) call dqagsmv(fhermit,-tmax,tmax,apar,npar,epsa,epsr,resfh, & epp,neval,ier,liw,lw,last,iw,w) if (ier /= 0) error = raise_error(error, dielectric_tensor, 1) rr(n,ih,m) = resfh end do end do end do if(fast.gt.10) return ! sy1=one+yg sy2=one+yg*2.0_wp_ sy3=one+yg*3.0_wp_ ! bth4=bth2*bth2 bth6=bth4*bth2 ! npl2=npl*npl ! rr(0,2,0) = -(one+bth2*(-1.25_wp_+1.5_wp_*npl2) & +bth4*(1.71875_wp_-6.0_wp_*npl2+3.75_wp_*npl2*npl2)) rr(0,1,1) = -npl*bth2*(one+bth2*(-2.25_wp_+1.5_wp_*npl2)) rr(0,2,1) = -bth2*(one+bth2*(-0.5_wp_+1.5_wp_*npl2)) rr(-1,0,1) = -2.0_wp_/sy1*(one+bth2/sy1*(-1.25_wp_+0.5_wp_*npl2 & /sy1)+bth4/sy1*(-0.46875_wp_+(2.1875_wp_+0.625_wp_*npl2) & /sy1-2.625_wp_*npl2/sy1**2+0.75_wp_*npl2*npl2/sy1**3)) rr(-1,1,1) = -npl*bth2/sy1**2*(one+bth2*(1.25_wp_-3.5_wp_/sy1 & +1.5_wp_*npl2/sy1**2)) rr(-1,2,1) = -bth2/sy1*(one+bth2*(1.25_wp_-1.75_wp_/sy1+1.5_wp_ & *npl2/sy1**2)) ! if(llm.gt.1) then ! rr(0,0,2) = -4.0_wp_*bth2*(one+bth2*(-0.5_wp_+0.5_wp_*npl2) & +bth4*(1.125_wp_-1.875_wp_*npl2+0.75_wp_*npl2*npl2)) rr(0,1,2) = -2.0_wp_*npl*bth4*(one+bth2*(-1.5_wp_+1.5_wp_*npl2)) rr(0,2,2) = -2.0_wp_*bth4*(one+bth2*(0.75_wp_+1.5_wp_*npl2)) rr(-1,0,2) = -4.0_wp_*bth2/sy1*(one+bth2 & *(1.25_wp_-1.75_wp_/sy1+0.5_wp_*npl2/sy1**2)+bth4 & *(0.46875_wp_-3.28125_wp_/sy1+(3.9375_wp_+1.5_wp_*npl2) & /sy1**2-3.375_wp_*npl2/sy1**3+0.75_wp_*npl2*npl2/sy1**4)) rr(-1,1,2) = -2.0_wp_*bth4*npl/sy1**2*(one+bth2 & *(3.0_wp_-4.5_wp_/sy1+1.5_wp_*npl2/sy1**2)) rr(-1,2,2) = -2.0_wp_*bth4/sy1*(one+bth2 & *(3.0_wp_-2.25_wp_/sy1+1.5_wp_*npl2/sy1**2)) rr(-2,0,2) = -4.0_wp_*bth2/sy2*(one+bth2 & *(1.25_wp_-1.75_wp_/sy2+0.5_wp_*npl2/sy2**2)+bth4 & *(0.46875_wp_-3.28125_wp_/sy2+(3.9375_wp_+1.5_wp_*npl2) & /sy2**2-3.375_wp_*npl2/sy2**3+0.75_wp_*npl2*npl2/sy2**4)) rr(-2,1,2) = -2.0_wp_*bth4*npl/sy2**2*(one+bth2 & *(3.0_wp_-4.5_wp_/sy2+1.5_wp_*npl2/sy2**2)) rr(-2,2,2) = -2.0_wp_*bth4/sy2*(one+bth2 & *(3.0_wp_-2.25_wp_/sy2+1.5_wp_*npl2/sy2**2)) ! if(llm.gt.2) then ! rr(0,0,3) = -12.0_wp_*bth4*(1+bth2*(0.75_wp_+0.5_wp_*npl2)+bth4 & *(1.21875_wp_-1.5_wp_*npl2+0.75_wp_*npl2*npl2)) rr(0,1,3) = -6.0_wp_*npl*bth6*(1+bth2*(-0.25_wp_+1.5_wp_*npl2)) rr(0,2,3) = -6.0_wp_*bth6*(1+bth2*(2.5_wp_+1.5_wp_*npl2)) rr(-1,0,3) = -12.0_wp_*bth4/sy1 & *(one+bth2*(3.0_wp_-2.25_wp_/sy1+0.5_wp_*npl2/sy1**2) & +bth4*(3.75_wp_-8.71875_wp_/sy1 & +(6.1875_wp_+2.625_wp_*npl2)/sy1**2 & -4.125_wp_*npl2/sy1**3+0.75_wp_*npl2*npl2/sy1**4)) rr(-1,1,3) = -6.0_wp_*npl*bth6/sy1**2 & *(one+bth2*(5.25_wp_-5.5_wp_/sy1+1.5_wp_*npl2/sy1**2)) rr(-1,2,3) = -6.0_wp_*bth6/sy1 & *(one+bth2*(5.25_wp_-2.75_wp_/sy1+1.5_wp_*npl2/sy1**2)) ! rr(-2,0,3) = -12.0_wp_*bth4/sy2 & *(one+bth2*(3.0_wp_-2.25_wp_/sy2+0.5_wp_*npl2/sy2**2) & +bth4*(3.75_wp_-8.71875_wp_/sy2 & +(6.1875_wp_+2.625_wp_*npl2)/sy2**2 & -4.125_wp_*npl2/sy2**3+0.75_wp_*npl2*npl2/sy2**4)) rr(-2,1,3) = -6.0_wp_*npl*bth6/sy2**2 & *(one+bth2*(5.25_wp_-5.5_wp_/sy2+1.5_wp_*npl2/sy2**2)) rr(-2,2,3) = -6.0_wp_*bth6/sy2 & *(one+bth2*(5.25_wp_-2.75_wp_/sy2+1.5_wp_*npl2/sy2**2)) ! rr(-3,0,3) = -12.0_wp_*bth4/sy3 & *(one+bth2*(3.0_wp_-2.25_wp_/sy3+0.5_wp_*npl2/sy3**2) & +bth4*(3.75_wp_-8.71875_wp_/sy3 & +(6.1875_wp_+2.625_wp_*npl2)/sy3**2 & -4.125_wp_*npl2/sy3**3+0.75_wp_*npl2*npl2/sy3**4)) rr(-3,1,3) = -6.0_wp_*npl*bth6/sy3**2 & *(one+bth2*(5.25_wp_-5.5_wp_/sy3+1.5_wp_*npl2/sy3**2)) rr(-3,2,3) = -6.0_wp_*bth6/sy3 & *(one+bth2*(5.25_wp_-2.75_wp_/sy3+1.5_wp_*npl2/sy3**2)) ! end if end if ! end subroutine hermitian_2 ! ! ! function fhermit(t,apar,npar) use eierf, only : calcei3 ! arguments integer, intent(in) :: npar real(wp_), intent(in) :: t real(wp_), dimension(npar), intent(in) :: apar ! local variables integer :: n,m,ih real(wp_) :: yg,mu,npl,cr,mu2,mu4,mu6,bth,bth2,rxt,x,upl2, & upl,gx,exdxdt,gr,zm,s,zm2,zm3,fe0m,ffe,uplh ! external functions/subroutines real(wp_) :: fhermit ! yg = apar(1) mu = apar(2) npl = apar(3) cr = apar(4) n = int(apar(5)) m = int(apar(6)) ih = int(apar(7)) ! bth2=2.0_wp_/mu bth=sqrt(bth2) mu2=mu*mu mu4=mu2*mu2 mu6=mu4*mu2 ! rxt=sqrt(one+t*t/(2.0_wp_*mu)) x = t*rxt upl2=bth2*x**2 upl=bth*x gx=one+t*t/mu exdxdt=cr*exp(-t*t)*gx/rxt gr=npl*upl+n*yg zm=-mu*(gx-gr) s=mu*(gx+gr) zm2=zm*zm zm3=zm2*zm call calcei3(zm,fe0m) ffe=zero uplh=upl**ih if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2 if(m.eq.1) ffe=(one+s*(one-zm*fe0m))*uplh/mu2 if(m.eq.2) ffe=(6.0_wp_-2.0_wp_*zm+4.0_wp_*s+s*s*(one+zm-zm2*fe0m))*uplh/mu4 if(m.eq.3) ffe=(18.0_wp_*s*(s+4.0_wp_-zm)+6.0_wp_*(20.0_wp_-8.0_wp_*zm+zm2) & +s**3*(2.0_wp_+zm+zm2-zm3*fe0m))*uplh/mu6 fhermit= exdxdt*ffe ! end function fhermit ! ! ! subroutine antihermitian(ri,yg,mu,npl,ci,lrm) ! local constants integer, parameter :: lmx=20,nmx=lmx+2 ! arguments integer :: lrm real(wp_) :: yg,mu,npl,ci real(wp_) :: ri(lrm,0:2,lrm) ! local variables integer :: n,k,m,mm real(wp_) :: cmu,npl2,dnl,ygn,rdu2,rdu,du,ub,aa,up,um,gp,gm,xp,xm, & eep,eem,ee,cm,cim,fi0p0,fi1p0,fi2p0,fi0m0,fi1m0,fi2m0, & fi0p1,fi1p1,fi2p1,fi0m1,fi1m1,fi2m1,fi0m,fi1m,fi2m real(wp_), dimension(nmx) :: fsbi ! do n=1,lrm do k=0,2 do m=1,lrm ri(n,k,m)=zero end do end do end do ! npl2=npl*npl dnl=one-npl2 cmu=npl*mu ! do n=1,lrm ygn=n*yg rdu2=ygn**2-dnl if(rdu2.gt.zero) then rdu=sqrt(rdu2) du=rdu/dnl ub=npl*ygn/dnl aa=mu*npl*du if (abs(aa).gt.5.0_wp_) then up=ub+du um=ub-du gp=npl*up+ygn gm=npl*um+ygn xp=up+one/cmu xm=um+one/cmu eem=exp(-mu*(gm-one)) eep=exp(-mu*(gp-one)) fi0p0=-one/cmu fi1p0=-xp/cmu fi2p0=-(one/cmu**2+xp*xp)/cmu fi0m0=-one/cmu fi1m0=-xm/cmu fi2m0=-(one/cmu**2+xm*xm)/cmu do m=1,lrm fi0p1=-2.0_wp_*m*(fi1p0-ub*fi0p0)/cmu fi0m1=-2.0_wp_*m*(fi1m0-ub*fi0m0)/cmu fi1p1=-((one+2.0_wp_*m)*fi2p0-2.0_wp_*(m+one)*ub*fi1p0 & +up*um*fi0p0)/cmu fi1m1=-((one+2.0_wp_*m)*fi2m0-2.0_wp_*(m+one)*ub*fi1m0 & +up*um*fi0m0)/cmu fi2p1=(2.0_wp_*(one+m)*fi1p1-2.0_wp_*m* & (ub*fi2p0-up*um*fi1p0))/cmu fi2m1=(2.0_wp_*(one+m)*fi1m1-2.0_wp_*m* & (ub*fi2m0-up*um*fi1m0))/cmu if(m.ge.n) then ri(n,0,m)=0.5_wp_*ci*dnl**m*(fi0p1*eep-fi0m1*eem) ri(n,1,m)=0.5_wp_*ci*dnl**m*(fi1p1*eep-fi1m1*eem) ri(n,2,m)=0.5_wp_*ci*dnl**m*(fi2p1*eep-fi2m1*eem) end if fi0p0=fi0p1 fi1p0=fi1p1 fi2p0=fi2p1 fi0m0=fi0m1 fi1m0=fi1m1 fi2m0=fi2m1 end do else ee=exp(-mu*(ygn-one+npl*ub)) call ssbi(aa,n,lrm,fsbi) block ! incrementally compute m! integer :: fact_m fact_m = 1 do m=n,lrm fact_m = fact_m * m cm=sqrt_pi*fact_m*du**(2*m+1) cim=0.5_wp_*ci*dnl**m mm=m-n+1 fi0m=cm*fsbi(mm) fi1m=-0.5_wp_*aa*cm*fsbi(mm+1) fi2m=0.5_wp_*cm*(fsbi(mm+1)+0.5_wp_*aa*aa*fsbi(mm+2)) ri(n,0,m)=cim*ee*fi0m ri(n,1,m)=cim*ee*(du*fi1m+ub*fi0m) ri(n,2,m)=cim*ee*(du*du*fi2m+2.0_wp_*du*ub*fi1m+ub*ub*fi0m) end do end block end if end if end do ! end subroutine antihermitian ! ! ! subroutine ssbi(zz,n,l,fsbi) ! local constants integer, parameter :: lmx=20,nmx=lmx+2 real(wp_), parameter :: eps=1.0e-10_wp_ ! arguments integer :: n,l real(wp_) :: zz real(wp_), dimension(nmx) :: fsbi ! local variables integer :: k,m,mm real(wp_) :: c0,c1,sbi ! do m=n,l+2 c0=one/gamma(dble(m)+1.5_wp_) sbi=c0 do k=1,50 c1=c0*0.25_wp_*zz**2/(dble(m+k)+0.5_wp_)/dble(k) sbi=sbi+c1 if(c1/sbi.lt.eps) exit c0=c1 end do mm=m-n+1 fsbi(mm)=sbi end do ! end subroutine ssbi ! ! ! subroutine expinit ! local variables integer :: i ! do i = 1, npts+1 ttv(i) = -tmax+dble(i-1)*dtex extv(i)=exp(-ttv(i)*ttv(i)) end do ! end subroutine expinit pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm, error) use gray_errors, only : gray_error, dielectric_tensor, raise_error ! subroutine arguments integer, intent(in) :: lrm real(wp_), intent(in) :: yg, npl, mu complex(wp_), intent(out), dimension(0:lrm,0:2) :: cefp, cefm integer(kind=gray_error), intent(out) :: error ! local constants real(wp_), parameter :: apsicr=0.7_wp_ ! local variables integer :: is,l,iq,ir,iflag real(wp_) :: psi,apsi,alpha,phi2,phim,xp,yp,xm,ym,x0,y0, & zrp,zip,zrm,zim,zr0,zi0 complex(wp_) :: czp,czm,cf12,cf32,cphi,cz0,cdz0,cf0,cf1,cf2 ! psi=sqrt(0.5_wp_*mu)*npl apsi=abs(psi) ! do is=0,lrm alpha=npl*npl/2.0_wp_+is*yg-one phi2=mu*alpha phim=sqrt(abs(phi2)) if (alpha.ge.0) then xp=psi-phim yp=zero xm=-psi-phim ym=zero x0=-phim y0=zero else xp=psi yp=phim xm=-psi ym=phim x0=zero y0=phim end if call zetac (xp,yp,zrp,zip,iflag) if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) call zetac (xm,ym,zrm,zim,iflag) if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) ! czp=dcmplx(zrp,zip) czm=dcmplx(zrm,zim) cf12=czero if (alpha.ge.0) then if (alpha.ne.0) cf12=-(czp+czm)/(2.0_wp_*phim) else cf12=-im*(czp+czm)/(2.0_wp_*phim) end if ! if(apsi.gt.apsicr) then cf32=-(czp-czm)/(2.0_wp_*psi) else cphi=phim if(alpha.lt.0) cphi=-im*phim call zetac (x0,y0,zr0,zi0,iflag) if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) cz0=dcmplx(zr0,zi0) cdz0=2.0_wp_*(one-cphi*cz0) cf32=cdz0 end if ! cf0=cf12 cf1=cf32 cefp(is,0)=cf32 cefm(is,0)=cf32 do l=1,is+2 iq=l-1 if(apsi.gt.apsicr) then cf2=(one+phi2*cf0-(iq+0.5_wp_)*cf1)/psi**2 else cf2=(one+phi2*cf1)/dble(iq+1.5_wp_) end if ir=l-is if(ir.ge.0) then cefp(is,ir)=cf2 cefm(is,ir)=cf2 end if cf0=cf1 cf1=cf2 end do ! if(is.ne.0) then ! alpha=npl*npl/2.0_wp_-is*yg-one phi2=mu*alpha phim=sqrt(abs(phi2)) if (alpha.ge.zero) then xp=psi-phim yp=zero xm=-psi-phim ym=zero x0=-phim y0=zero else xp=psi yp=phim xm=-psi ym=phim x0=zero y0=phim end if call zetac (xp,yp,zrp,zip,iflag) if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) call zetac (xm,ym,zrm,zim,iflag) if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) ! czp=dcmplx(zrp,zip) czm=dcmplx(zrm,zim) ! cf12=czero if (alpha.ge.0) then if (alpha.ne.zero) cf12=-(czp+czm)/(2.0_wp_*phim) else cf12=-im*(czp+czm)/(2.0_wp_*phim) end if if(apsi.gt.apsicr) then cf32=-(czp-czm)/(2.0_wp_*psi) else cphi=phim if(alpha.lt.0) cphi=-im*phim call zetac (x0,y0,zr0,zi0,iflag) if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) cz0=dcmplx(zr0,zi0) cdz0=2.0_wp_*(one-cphi*cz0) cf32=cdz0 end if ! cf0=cf12 cf1=cf32 do l=1,is+2 iq=l-1 if(apsi.gt.apsicr) then cf2=(one+phi2*cf0-(iq+0.5_wp_)*cf1)/psi**2 else cf2=(one+phi2*cf1)/dble(iq+1.5_wp_) end if ir=l-is if(ir.ge.0) then cefp(is,ir)=cefp(is,ir)+cf2 cefm(is,ir)=cefm(is,ir)-cf2 end if cf0=cf1 cf1=cf2 end do ! end if ! end do ! end subroutine fsup ! ! PLASMA DISPERSION FUNCTION Z of complex argument ! Z(z) = i sqrt(pi) w(z) ! Function w(z) from: ! algorithm 680, collected algorithms from acm. ! this work published in transactions on mathematical software, ! vol. 16, no. 1, pp. 47. ! pure subroutine zetac (xi, yi, zr, zi, iflag) ! ! given a complex number z = (xi,yi), this subroutine computes ! the value of the faddeeva-function w(z) = exp(-z**2)*erfc(-i*z), ! where erfc is the complex complementary error-function and i ! means sqrt(-1). ! the accuracy of the algorithm for z in the 1st and 2nd quadrant ! is 14 significant digits; in the 3rd and 4th it is 13 significant ! digits outside a circular region with radius 0.126 around a zero ! of the function. ! all real variables in the program are real(8). ! ! ! the code contains a few compiler-dependent parameters : ! rmaxreal = the maximum value of rmaxreal equals the root of ! rmax = the largest number which can still be ! implemented on the computer in real(8) ! floating-point arithmetic ! rmaxexp = ln(rmax) - ln(2) ! rmaxgoni = the largest possible argument of a real(8) ! goniometric function (cos, sin, ...) ! the reason why these parameters are needed as they are defined will ! be explained in the code by means of comments ! ! ! parameter list ! xi = real part of z ! yi = imaginary part of z ! u = real part of w(z) ! v = imaginary part of w(z) ! iflag = an error flag indicating whether overflow will ! occur or not; type integer; ! the values of this variable have the following ! meaning : ! iflag=0 : no error condition ! iflag=1 : overflow will occur, the routine ! becomes inactive ! xi, yi are the input-parameters ! u, v, iflag are the output-parameters ! ! furthermore the parameter factor equals 2/sqrt(pi) ! ! the routine is not underflow-protected but any variable can be ! put to 0 upon underflow; ! ! reference - gpm poppe, cmj wijers; more efficient computation of ! the complex error-function, acm trans. math. software. ! real(wp_), intent(in) :: xi, yi real(wp_), intent(out) :: zr, zi integer, intent(out) :: iflag real(wp_) :: xabs,yabs,x,y,qrho,xabsq,xquad,yquad,xsum,ysum,xaux,daux, & u,u1,u2,v,v1,v2,h,h2,qlambda,c,rx,ry,sx,sy,tx,ty,w1 integer :: i,j,n,nu,kapn,np1 real(wp_), parameter :: factor = 1.12837916709551257388_wp_, & rpi = 2.0_wp_/factor, & rmaxreal = 0.5e+154_wp_, & rmaxexp = 708.503061461606_wp_, & rmaxgoni = 3.53711887601422e+15_wp_ iflag=0 xabs = abs(xi) yabs = abs(yi) x = xabs/6.3_wp_ y = yabs/4.4_wp_ ! ! the following if-statement protects ! qrho = (x**2 + y**2) against overflow ! if ((xabs>rmaxreal).or.(yabs>rmaxreal)) then iflag=1 return end if qrho = x**2 + y**2 xabsq = xabs**2 xquad = xabsq - yabs**2 yquad = 2*xabs*yabs if (qrho<0.085264_wp_) then ! ! if (qrho<0.085264_wp_) then the faddeeva-function is evaluated ! using a power-series (abramowitz/stegun, equation (7.1.5), p.297) ! n is the minimum number of terms needed to obtain the required ! accuracy ! qrho = (1-0.85_wp_*y)*sqrt(qrho) n = nint(6 + 72*qrho) j = 2*n+1 xsum = 1.0_wp_/j ysum = 0.0_wp_ do i=n, 1, -1 j = j - 2 xaux = (xsum*xquad - ysum*yquad)/i ysum = (xsum*yquad + ysum*xquad)/i xsum = xaux + 1.0_wp_/j end do u1 = -factor*(xsum*yabs + ysum*xabs) + 1.0_wp_ v1 = factor*(xsum*xabs - ysum*yabs) daux = exp(-xquad) u2 = daux*cos(yquad) v2 = -daux*sin(yquad) u = u1*u2 - v1*v2 v = u1*v2 + v1*u2 else ! ! if (qrho>1.o) then w(z) is evaluated using the laplace ! continued fraction ! nu is the minimum number of terms needed to obtain the required ! accuracy ! ! if ((qrho>0.085264_wp_).and.(qrho<1.0)) then w(z) is evaluated ! by a truncated taylor expansion, where the laplace continued fraction ! is used to calculate the derivatives of w(z) ! kapn is the minimum number of terms in the taylor expansion needed ! to obtain the required accuracy ! nu is the minimum number of terms of the continued fraction needed ! to calculate the derivatives with the required accuracy ! if (qrho>1.0_wp_) then h = 0.0_wp_ kapn = 0 qrho = sqrt(qrho) nu = int(3 + (1442/(26*qrho+77))) else qrho = (1-y)*sqrt(1-qrho) h = 1.88_wp_*qrho h2 = 2*h kapn = nint(7 + 34*qrho) nu = nint(16 + 26*qrho) endif if (h>0.0_wp_) qlambda = h2**kapn rx = 0.0_wp_ ry = 0.0_wp_ sx = 0.0_wp_ sy = 0.0_wp_ do n=nu, 0, -1 np1 = n + 1 tx = yabs + h + np1*rx ty = xabs - np1*ry c = 0.5_wp_/(tx**2 + ty**2) rx = c*tx ry = c*ty if ((h>0.0_wp_).and.(n<=kapn)) then tx = qlambda + sx sx = rx*tx - ry*sy sy = ry*tx + rx*sy qlambda = qlambda/h2 endif end do if (h==0.0_wp_) then u = factor*rx v = factor*ry else u = factor*sx v = factor*sy end if if (yabs==0.0_wp_) u = exp(-xabs**2) end if ! ! evaluation of w(z) in the other quadrants ! if (yi<0.0_wp_) then if (qrho<0.085264_wp_) then u2 = 2*u2 v2 = 2*v2 else xquad = -xquad ! ! the following if-statement protects 2*exp(-z**2) ! against overflow ! if ((yquad>rmaxgoni).or. (xquad>rmaxexp)) then iflag=1 return end if w1 = 2.0_wp_*exp(xquad) u2 = w1*cos(yquad) v2 = -w1*sin(yquad) end if u = u2 - u v = v2 - v if (xi>0.0_wp_) v = -v else if (xi<0.0_wp_) v = -v end if zr = -v*rpi zi = u*rpi end subroutine zetac ! end module dispersion