From 3a0e45e95b9744c14790791f31c4a3d06c5ece2c Mon Sep 17 00:00:00 2001
From: kamischi <karl-michael.schindler@web.de>
Date: Fri, 13 Jan 2023 17:52:54 +0100
Subject: [PATCH] commented source code added.
---
.../boson_all.f90.fcat-analysis.txt | 972 ++++++++++++
.../eels_all.f90.fcat-analysis.txt | 1353 +++++++++++++++++
2 files changed, 2325 insertions(+)
create mode 100644 source/f90/fcat-analysis/boson_all.f90.fcat-analysis.txt
create mode 100644 source/f90/fcat-analysis/eels_all.f90.fcat-analysis.txt
diff --git a/source/f90/fcat-analysis/boson_all.f90.fcat-analysis.txt b/source/f90/fcat-analysis/boson_all.f90.fcat-analysis.txt
new file mode 100644
index 0000000..8d97d00
--- /dev/null
+++ b/source/f90/fcat-analysis/boson_all.f90.fcat-analysis.txt
@@ -0,0 +1,972 @@
+program boson
+
+! *******************************************************************
+! * *
+! * perform the quantum-mechanical complement to the classical step *
+! * of the dielectric theory of eels in specular geometry using a *
+! * suitable thermodynamical average of the quantized surface *
+! * harmonic oscillators *
+! * *
+! *******************************************************************
+
+ implicit none
+
+ integer, parameter :: mmax=14, nmax=2**mmax
+ double precision :: asym, emax, emin, gauss, t, width, wmin, wmax, y
+ double precision :: p2(nmax)
+ integer :: i, np, ioStatus
+ character (len = 72) :: comment(2)
+ double precision, allocatable :: xout(:), yout(:)
+ integer :: nout
+
+! *** read input parameters
+
+1 call change_working_dir()
+1 open(unit = 13, file = 'bosin')
+! target temperature (Kelvin)
+1 read(13, *) t
+! width of the instrumental response (cm**-1)
+1 read(13, *) width
+! fraction of gaussian for the instrumental response
+1 read(13, *) gauss
+! asymmetry of the instrumental response
+1 read(13, *) asym
+! lower and upper energy losses for this computation (cm**-1)
+1 read(13, *) emin
+1 read(13, *) emax
+1 close(unit=13)
+
+1 write(*,*) 'program boson (September 2020)'
+1 write(*,'(a, f6.1, a, f7.2, a)') ' t =', t, ' K, width =', width, ' cm**-1'
+1 write(*,'(a, f5.2, a, f5.2)') ' gauss =', gauss, ', asym =', asym
+1 write(*,'(a, g11.4, a, g11.4, a)') ' energy losses from', emin, ' to', emax, ' cm**-1'
+
+! *** read the table of values of the classical loss spectrum
+
+1 open(unit = 12, file = 'eelsou')
+1 read(12, '(a48)') comment(1)
+1 read(12, '(a72)') comment(2)
+1 np = 0
+1 do
+327 read(12, *, IOSTAT = ioStatus) wmax, y
+327 if (ioStatus /= 0) exit
+326 if (wmax < 0.0d0) cycle
+326 np = np + 1
+326 p2(np) = y
+326 if (np == 1) wmin = wmax
+326 if (np < nmax) cycle
+*> enddo
+1 close(unit = 12)
+
+1 if (np <= 0) stop '*** no input values for pcl ***'
+
+1 write(*,*) comment(2)
+1 write(*,'(a, i6, a, g15.7, a, g15.7)') ' read', np, ' values of pcl from', wmin, ' to', wmax
+
+! length calculation for xout, yout.
+!
+1 nout = 2**14
+1 allocate (xout(nout))
+1 allocate (yout(nout))
+
+1 call doboson(t, width, gauss, asym, emin, emax, wmin, wmax, np, p2, xout, yout, nout)
+
+1 open(unit = 14, file = 'bosou')
+1 write(14, '(a, a, f6.1, a, f5.2)') comment(1), 'T =', t, ' GAUSS =', gauss
+1 write(14, *) comment(2)
+1 do i = 1, nout
+851 write(14, '(2e15.7)') xout(i), yout(i)
+851 end do
+1 close(unit = 14)
+1 write(*,*) nout, ' values written on disk'
+
+1 deallocate (xout, yout)
+1 stop
+end program boson
+subroutine change_working_dir()
+
+! This routine gets the first argument of the commandline and takes it
+! as the path to change the working directory
+! used intrinsic routines:
+! iarg returns the number of commandline arguments without the program cname.
+! chdir changes the directory and returns 0 on success.
+! trim removes blanks from strings.
+
+ character (len = 256) :: argument
+ integer :: status
+
+1 if (iargc() == 1) then
+*> call getarg(1, argument)
+*> status = chdir(trim(argument))
+*> if (status /= 0) then
+*> write (*,*) '*** change directory failed ***'
+*> write (*,*) 'directory tried: ', trim(argument)
+*> write (*,*) 'error code (see: man chdir): ', status
+*> write (*,*) 'continuing in the start directory'
+*> end if
+*> end if
+
+1 return
+end subroutine change_working_dir
+double precision function respon(w, width)
+
+!*******************************************************************
+! *
+! instrumental response of the spectrometer for the frequency w *
+! width is the full width at half maximum of the response function *
+! *
+!*******************************************************************
+
+ implicit none
+ double precision, intent(in) :: w, width
+
+ double precision :: u
+
+*> u = 1.25d0 * w / width
+*> if (u <= -1.0d0 .or. u >= 2.0d0) then
+*> respon = 0.0d0
+*> elseif (u > 1.0d0) then
+*> respon = 4.0d0 + 2 * dsqrt((u - 1.0d0)**3) - 3 * u
+*> elseif (u > 0.0d0) then
+*> respon = 2.0d0 - 4 * dsqrt(u**3) + 3 * u
+*> else
+*> respon = 2 * dsqrt((u + 1.0d0)**3)
+*> endif
+*> return
+end function respon
+module sicot_mod
+contains
+subroutine sicot(f, m, h, x0)
+
+! *******************************************************************
+! * *
+! * integral transform g(y, x0) of a real function f(x) defined by *
+! * *
+! * g(x0, y) = integral of f(x)/tanh(x/x0)*(1-cos(x*y)) dx *
+! * *
+! * f(x) is tabulated on the mesh of points xj = (j-1/2)*h , *
+! * j = 1, 2, ..., n , with n = 2**m *
+! * g(x0, y) is computed on the mesh of points yk = k*2pi/(n*h) , *
+! * k = 0, 1, 2, ..., n-1 *
+! * *
+! * input ... *
+! * f(j) contains the value of f(xj) , j = 1, 2, ..., n *
+! * m is such that n = 2**iabs(m) *
+! * h is the step size (must be positive) *
+! * *
+! * output ... *
+! * f(k+1) contains g(x0, yk), k = 0, 1, ..., n-1 *
+! * *
+! * computing remarks ... *
+! * f is a real array with dimension n or more *
+! * it is supposed that f(x) is zero for x > n*h *
+! * external reference : sintr (sine transform of a real function). *
+! * *
+! *******************************************************************
+
+ use sintr_mod
+
+ implicit none
+
+ double precision, intent(in out) :: f(*)
+ integer, intent(in) :: m
+ double precision, intent(in) :: h
+ double precision, intent(in) :: x0
+
+ double precision :: a, b, c, d, e, s, sa
+ integer :: i, j, msign, n
+
+ logical :: debug
+
+1 debug = .false.
+1 if (m == 0) then
+*> f(1) = 0.0d0
+*> return
+*> endif
+
+1 if (h <= 0.0d0) then
+*> if (debug) write(*,*) ' *** incorrect step size in <sicot>, h =', h, ' ***'
+*> stop
+*> endif
+
+1 if (x0 < 0.0d0) then
+*> if (debug) write(*,*) ' *** incorrect input in <sicot>, x0 =', x0, ' ***'
+*> stop
+*> endif
+
+1 n = 2**iabs(m)
+
+! *** evaluate the integral from xj to infinity of f(x)/tanh(x/x0) dx
+! and store the result in f(j)
+
+1 c = 0.0d0
+1 if (x0 > h / 16) then
+1 c = dexp(-h / x0)
+1 endif
+1 s = 1.0d0 - c
+1 e = c**2
+1 a = 0.25d0 * h
+1 b = 0.25d0 * x0
+1 do i = 1, n
+4096 if (i < n) then
+4095 f(i) = f(i) + f(i+1)
+4095 endif
+4096 d = a
+4096 if (s /= 1.0d0) then
+3902 sa = s
+3902 c = c * e
+3902 s = 1.0d0 - c
+3902 d = d + b * dlog(s / sa)
+3902 endif
+4096 f(i) = f(i) * d
+4096 enddo
+
+1 j = n
+1 do i = 2, n
+4095 j = j - 1
+4095 f(j) = f(j) + f(j + 1)
+4095 enddo
+! alternative, but not yet tested:
+! do j = n - 1, 1, -1
+! f(j) = f(j) + f(j + 1)
+! enddo
+1 msign = -iabs(m)
+1 call sintr(f, msign, h)
+1 return
+end subroutine sicot
+end module sicot_mod
+module sintr_mod
+contains
+subroutine sintr(f, msign, h)
+
+! *******************************************************************
+! * *
+! * integral sine transform of a real function f(x) *
+! * *
+! * g(y) = integral from zero to infinity of f(x)*sin(x*y) dx *
+! * if msign >= 0 or *
+! * g(y) = y * integral from zero to infinity of f(x)*sin(x*y) dx *
+! * if msign < 0 *
+! * *
+! * f(x) is tabulated on the mesh of points xj = (j-1/2)*h , *
+! * j = 1, 2, ..., n , with n = 2**iabs(msig) *
+! * g(y) is computed on the mesh of points yk = k*2pi/(n*h) , *
+! * k = 0, 1, 2, ..., n-1 *
+! * *
+! * input ... *
+! * f(j) contains the value of f(xj) , j = 1, 2, ..., n *
+! * msign is such that n = 2**iabs(msign) *
+! * h is the step size (must be positive) *
+! * *
+! * output ... *
+! * f(k+1) contains g(yk), k = 0, 1, ..., n-1 *
+! * *
+! * computing remarks ... *
+! * f is a real array with dimension n or more *
+! * it is supposed that f(x) is zero for x > n*h *
+! * *
+! *******************************************************************
+
+ implicit none
+
+ double precision, intent(in out) :: f(*)
+ integer, intent(in) :: msign
+ double precision, intent(in) :: h
+
+ double precision :: a, c, ca, d, e, fnp1, s, sa, ti, ti1, ti2, tr, tr1, tr2
+ double precision :: ui, ur, wi, wr
+ integer :: i, i2, ip2, j, j2, k, l, le, le1, m2, n, n2, nm1, nv2
+
+ logical :: debug
+
+2 debug = .false.
+2 if (h <= 0.0d0) then
+*> if (debug) write(*,*) ' *** incorrect step size in <sintr>, h: ', h, ' ***'
+*> stop
+*> endif
+
+2 if (msign == 0) then
+*> f(1) = 0.0d0
+*> return
+*> endif
+
+! *** for j and k = 0, 1, 2 ... n/2-1, compute
+! s1 = sum ( f(2*j+1)*cos(4*pi*k*j/n) - f(2*j+2)*sin(4*pi*k*j/n) )
+! s2 = sum ( f(2*j+1)*sin(4*pi*k*j/n) + f(2*j+2)*cos(4*pi*k*j/n) )
+! store s1 in f(2*k+1) and s2 in f(2*k+2)
+! (adapted from cfft routine, cern library member d704)
+
+2 m2 = iabs(msign) - 1
+2 n = 2**m2
+2 if (n /= 1) then
+2 nv2 = n / 2
+2 nm1 = n - 1
+2 j = 1
+2 do i = 1, nm1
+4094 if (i < j) then
+1984 i2 = i + i
+1984 j2 = j + j
+1984 ti = f(j2)
+1984 f(j2) = f(i2)
+1984 f(i2) = ti
+1984 i2 = i2 - 1
+1984 j2 = j2 - 1
+1984 tr = f(j2)
+1984 f(j2) = f(i2)
+1984 f(i2) = tr
+1984 endif
+4094 k = nv2
+4094 do while (j > k)
+4072 j = j - k
+4072 k = k / 2
+4072 enddo
+4094 j = j + k
+4094 enddo
+2 do i = 1, n, 2
+2048 i2 = i + i
+2048 ti = f(i2 + 2)
+2048 f(i2 + 2) = f(i2) - ti
+2048 f(i2) = f(i2) + ti
+2048 i2 = i2 - 1
+2048 tr = f(i2 + 2)
+2048 f(i2 + 2) = f(i2) - tr
+2048 f(i2) = f(i2) + tr
+2048 enddo
+2 if (m2 /= 1) then
+2 c = 0.0d0
+2 s = 1.0d0
+2 le = 2
+2 do l = 2, m2
+20 wr = c
+20 wi = s
+20 ur = wr
+20 ui = wi
+20 c = dsqrt(c * 0.5d0 + 0.5d0 )
+20 s = s / (c + c)
+20 le1 = le
+20 le = le1 + le1
+20 do i = 1, n, le
+2046 i2 = i + i
+2046 ip2 = i2 + le
+2046 ti = f(ip2)
+2046 f(ip2) = f(i2) - ti
+2046 f(i2) = f(i2) + ti
+2046 i2 = i2 - 1
+2046 ip2 = ip2 - 1
+2046 tr = f(ip2)
+2046 f(ip2) = f(i2) - tr
+2046 f(i2) = f(i2) + tr
+2046 enddo
+20 do j = 2, le1
+4072 do i = j, n, le
+18434 i2 = i + i
+18434 ip2 = i2 + le
+18434 tr = f(ip2 - 1) * ur - f(ip2) * ui
+18434 ti = f(ip2) * ur + f(ip2 - 1) * ui
+18434 f(ip2) = f(i2) - ti
+18434 f(i2) = f(i2) + ti
+18434 i2 = i2 - 1
+18434 ip2 = ip2 - 1
+18434 f(ip2) = f(i2) - tr
+18434 f(i2) = f(i2) + tr
+18434 enddo
+4072 tr = ur * wr - ui * wi
+4072 ui = ui * wr + ur * wi
+4072 ur = tr
+4072 enddo
+20 enddo
+2 endif
+2 endif
+
+! *** for j and k = 0, 1 ... n-1, transform the array f so obtained into
+! sum f(j+1)*cos(2*pi*k*j/n) and sum f(j+1)*sin(2*pi*k*j/n)
+! and multiply the results by (1 - cos(2*pi*k/n)) and sin(2*pi*k/n)
+
+2 fnp1 = 4 * (f(1) - f(2))
+2 a = 3.141592653589793238d0 / n
+2 n2 = n + 1
+2 n = 2 * n
+2 if (n2 /= 2) then
+2 c = 1.0d0
+2 s = 0.0d0
+2 ca = dcos(a)
+2 sa = dsin(a)
+2 k = n + 1
+2 do j = 3, n2, 2
+2048 k = k - 2
+2048 d = c
+2048 c = d * ca - s * sa
+2048 s = d * sa + s * ca
+2048 tr1 = f(j) + f(k)
+2048 ti1 = f(j + 1) - f(k + 1)
+2048 d = f(j) - f(k)
+2048 e = f(j + 1) + f(k + 1)
+2048 tr2 = d * s + e * c
+2048 ti2 = e * s - d * c
+2048 f(j) = (1.0d0 - c) * (tr1 + tr2)
+2048 f(j + 1) = s * (ti1 + ti2)
+2048 f(k) = (1.0d0 + c) * (tr1 - tr2)
+2048 f(k + 1) = s * (ti2 - ti1)
+2048 enddo
+2 n2 = n2 - 1
+2 do j = 2, n2
+4094 j2 = j + j
+4094 f(j) = f(j2 - 1) + f(j2)
+4094 enddo
+2 do j = 2, n2
+4094 f(n + 2 - j) = f(j)
+4094 enddo
+2 n2 = n2 + 1
+2 endif
+2 f(n2) = fnp1
+2 if (msign >= 0) then
+
+! *** normalization
+
+1 c = 0.0d0
+1 a = (a + a) / h
+1 do j = 2, n
+4095 c = c + a
+4095 f(j) = f(j) / c
+4095 enddo
+1 endif
+2 f(1) = 0.0d0
+2 return
+end subroutine sintr
+end module sintr_mod
+module rcffi_mod
+contains
+subroutine rcffi(ar, ai, msign, h)
+
+! *******************************************************************
+! * *
+! * using a radix-two fast-fourier-transform technique, rcffi *
+! * computes the fourier integral transform of a real function *
+! * f(x) or the inverse fourier transform of a complex function *
+! * such that g(-y) = conj(g(y)). *
+! * (adapted from cfft routine, cern library member d704) *
+! * *
+! * g(y) = integral f(x) * cexp(-i * x * y) dx (msign < 0) *
+! * *
+! * f(x) = 1 / (2 * pi) integral g(y) * cexp(+i * y * x) dy (msign > 0) *
+! * *
+! * f(x) is tabulated on the meshes of points defined by *
+! * xj = j * hx and -xj = -j * hx, j = 0, +1, ..., n-2, n-1 *
+! * ar(j+1) contains f(+xj) and ai(j+1) contains f(-xj) *
+! * (input when msign < 0 , output when msign >0) *
+! * ar(1) may be different from ai(1) *
+! * *
+! * g(y) is tabulated on the mesh of points defined by *
+! * yk = k * hy , k = 0, 1, ..., n - 2, n - 1 *
+! * ar(j+1) contains real(g(yk)) and ai(j+1) contains aimag(g(yk)) *
+! * (output when msign > 0 , input when msign < 0) *
+! * remark : g(-yk) = conj(g(+yk)) *
+! * *
+! * n = 2**iabs(msign) (msign is an input) *
+! * *
+! * the step sizes satisfy *
+! * hy = (2 * pi) / (n * hx) or hx = (2 * pi) / (n * hy) *
+! * h is an input (h = hx if msign < 0 , h = hy if msign > 0) *
+! * h must be positive *
+! * *
+! * ar and ai are two real arrays with dimension n (input and *
+! * output) *
+! * *
+! *******************************************************************
+
+ implicit none
+
+ double precision, intent(in out) :: ar(*)
+ double precision, intent(in out) :: ai(*)
+ integer, intent(in) :: msign
+ double precision, intent(in) :: h
+
+ double precision :: a, as, c, s, t, ti, tr, ui, ur, wi, wr
+ integer :: i, ip, j, k, l, le, le1, m, n, nm1, nv2
+ logical :: debug
+
+1 debug = .false.
+1 as = 0.0d0
+1 if (msign == 0) then
+*> return
+*> endif
+
+1 if (h <= 0.0d0) then
+*> if (debug) write(*,*) ' *** negative step size in <rcffi>, h = ', h, ' *** '
+*> stop
+*> endif
+
+! *** initialization
+
+1 m = iabs(msign)
+1 n = 2**m
+1 if (msign > 0) then
+1 ar(1) = 0.5d0 * ar(1)
+*> else
+*> as = ar(1) - ai(1)
+*> ar(1) = 0.5d0 * (ar(1) + ai(1))
+*> do i = 2, n
+*> ar(i) = ar(i) + ai(n - i + 2)
+*> enddo
+*> do i = 1, n
+*> ai(i) = 0.0d0
+*> enddo
+*> endif
+! *** discrete fast-fourier transform
+1 nv2 = n / 2
+1 nm1 = n - 1
+1 j = 1
+1 do i = 1, nm1
+4095 if (i < j) then
+2016 tr = ar(j)
+2016 ar(j) = ar(i)
+2016 ar(i) = tr
+2016 ti = ai(j)
+2016 ai(j) = ai(i)
+2016 ai(i) = ti
+2016 endif
+4095 k = nv2
+4095 do while (j > k)
+4083 j = j - k
+4083 k = k / 2
+4083 enddo
+4095 j = j + k
+4095 enddo
+1 do i = 1, n, 2
+2048 tr = ar(i + 1)
+2048 ar(i + 1) = ar(i) - tr
+2048 ar(i) = ar(i) + tr
+2048 ti = ai(i + 1)
+2048 ai(i + 1) = ai(i) - ti
+2048 ai(i) = ai(i) + ti
+2048 enddo
+1 if (m == 1) return
+1 c = 0.0d0
+1 s = dble(isign(1, msign))
+1 le = 2
+1 do l = 2, m
+11 wr = c
+11 ur = wr
+11 wi = s
+11 ui = wi
+11 c = dsqrt(c * 0.5d0 + 0.5d0)
+11 s = wi / (c + c)
+11 le1 = le
+11 le = le1 + le1
+11 do i = 1, n, le
+2047 ip = i + le1
+2047 tr = ar(ip)
+2047 ar(ip) = ar(i) - tr
+2047 ar(i) = ar(i) + tr
+2047 ti = ai(ip)
+2047 ai(ip) = ai(i) - ti
+2047 ai(i) = ai(i) + ti
+2047 enddo
+11 do j = 2, le1
+4083 do i = j, n, le
+20481 ip = i + le1
+20481 tr = ar(ip) * ur - ai(ip) * ui
+20481 ti = ar(ip) * ui + ai(ip) * ur
+20481 ar(ip) = ar(i) - tr
+20481 ar(i) = ar(i) + tr
+20481 ai(ip) = ai(i) - ti
+20481 ai(i) = ai(i) + ti
+20481 enddo
+4083 tr = ur
+4083 ur = tr * wr - ui * wi
+4083 ui = tr * wi + ui * wr
+4083 enddo
+11 enddo
+
+! *** correction of the discrete fft, using a trapezoidal approximation
+! *** for the variations of the input function
+
+1 c = h
+1 if (msign >= 0) then
+1 c = c / 3.141592653589793238d0
+1 ai(1) = ar(1)
+1 do i = 2, n
+4095 ai(i) = ar(n - i + 2)
+4095 enddo
+1 endif
+1 ar(1) = c * ar(1)
+1 ai(1) = c * ai(1)
+1 a = 3.141592653589793238d0 / n
+1 wr = dcos(a)
+1 wi = dsin(a)
+1 t = dsqrt(c)
+1 a = a / t
+1 c = 0.0d0
+1 ur = 1.0d0
+1 ui = 0.0d0
+1 do i = 2, n
+4095 c = c + a
+4095 tr = ur
+4095 ur = wr * tr - wi * ui
+4095 ui = wi * tr + wr * ui
+4095 ti = ui / c
+4095 tr = ti**2
+4095 ar(i) = ar(i) * tr
+4095 ai(i) = ai(i) * tr
+4095 if (msign <= 0) then
+*> if (as /= 0.0d0) then
+*> ai(i) = ai(i) - as * (t - ur * ti) / (2 * c)
+*> endif
+*> endif
+4095 enddo
+1 return
+end subroutine rcffi
+end module rcffi_mod
+subroutine doboson(t, width, gauss, asym, emin, emax, wmin, wmax, np, p, xout, yout, nout)
+
+! *******************************************************************
+! * *
+! * perform the quantum-mechanical complement to the classical step *
+! * of the dielectric theory of eels in specular geometry using a *
+! * suitable thermodynamical average of the quantized surface *
+! * harmonic oscillators *
+! * *
+! * t: Target temperature *
+! * width: Width of instrumental response (cm**-1) *
+! * gauss: Fraction of gaussian for the instrumental response *
+! * asym: Asymmetry of the instrumental response *
+! * emin: Lower loss energy for this computation (cm**-1) *
+! * emax: Upper loss energy for this computation (cm**-1) *
+! * wmin: Lower loss energy of the eels computation (cm**-1) *
+! * wmax: Upper loss energy of the eels computation (cm**-1) *
+! * np: Number of points of the eels computation *
+! * p: Intensities of the eels computation *
+! * xout: Energies of this computation *
+! * yout: Intensities of this computation *
+! * nout: Number of points of this computation *
+! *******************************************************************
+
+ use rcffi_mod
+ use sicot_mod
+ use sintr_mod
+
+ implicit none
+
+ double precision, intent(in) :: t, width, gauss, asym, emin, emax, wmin
+ double precision, intent(in out) :: wmax
+ integer, intent(in) :: np
+ double precision, intent(in) :: p(np)
+ double precision, intent(in out) :: xout(nout), yout(nout)
+ integer, intent(in out) :: nout
+
+ integer, parameter :: mmax = 14, nmax = 2**mmax
+ double precision :: a, a1, a2, alfa, anorm, b, cp2, dwn, fac, fi
+ double precision :: fm, fm0, fm1, fp1, fmpic, fp, fp0, fppic, fr, g1, g2
+ double precision :: h, pi, sigma, sp2, test, u
+ double precision :: wmpic, wppic, wn, x, x1, x2, x3
+ double precision :: o1, o2, respon
+
+ external o1, o2
+
+ integer :: i, imax, imin, istep, j, jmin, jmax, m, n
+ logical :: picm, picp
+ double precision, allocatable :: p1(:), p2(:)
+ double precision :: r1(nmax), r2(nmax)
+
+ logical :: debug
+
+! remark : the two arrays r1 and r2 are used for fourier
+! transforming the user-supplied instrumental response of the
+! spectrometer that has to be coded in the external routine
+! respon called when the input parameter gauss is < 0 or > 1.
+! with 0 <= gauss <= 1, r1 and r2 are not used.
+
+! *** rational approximations for ei(u) * exp(-u) + e1(u) * exp(+u)
+! *** in the intervals (0, 1.3) and (1.3, infinity) <accuracy : 4.e-04>
+! *** used for fourier transforming half-lorentzian functions
+
+
+1 debug = .false.
+1 data fm1 / 0.0d0 /, fp1 / 0.0d0 /
+1 pi = 4 * datan(1.0d0)
+
+1 dwn = (wmax - wmin) / (np - 1)
+
+! *** redefine the frequency interval (0, wmax) to be used for the
+! *** fourier transforms
+
+1 if (debug) write(*,*) 'wmin: ', wmin, 'wmax: ', wmax, 'dwn: ', dwn
+1 jmin = int(0.5d0 + wmin / dwn)
+1 if ((jmin - 0.5d0) * dwn < wmin) then
+1 jmin = jmin + 1
+1 endif
+1 fac = ((jmin - 0.5d0) * dwn - wmin) / dwn
+1 jmax = int(0.5d0 + wmax / dwn)
+1 test = dmax1(1.5d0 * dabs(emin), 1.5d0 * dabs(emax), 6 * wmax)
+1 n = 2
+1 do m = 2, mmax
+11 n = 2 * n
+11 wmax = n * dwn
+11 if (wmax >= test) exit
+10 enddo
+1 if (wmax < test) then
+*> m = mmax
+*> n = nmax
+*> wmax = n * dwn
+*> if (debug) write(*,*) ' +++ n has been fixed at', nmax, ' = 2**', mmax, ' +++'
+*> endif
+1 if (debug) then
+*> write(*,*) ' classical spectrum redefined from 0.0 to', wmax
+*> write(*,*) ' step size =', dwn, ', ', n, ' (= 2**', m, ') points'
+*> endif
+
+1 allocate (p1(n))
+1 allocate (p2(n))
+
+1 p1 = 0
+
+! *** interpolate pcl on a suitable mesh in (0, wmax)
+1 i = 1
+1 do j = jmin, min(jmax, n)
+325 p1(j) = p(i) + fac * (p(i + 1) - p(i))
+325 i = i + 1
+325 enddo
+
+1 p2 = p1
+
+! *** characteristic function f(tau)
+
+1 call sicot(p1, m, dwn, 1.39d0 * t)
+1 call sintr(p2, m, dwn)
+1 h = (pi + pi) / wmax
+1 if (gauss < 0.0d0 .or. gauss > 1.0d0) then
+
+! *** broaden the spectrum by convoluting the characteristic function
+! *** with a user-supplied response function (arbitrary normalization)
+
+*> if (debug) write(*,*) '==> switch to a user-supplied instrumental response'
+*> alfa = dmax1(4 * dwn, width)
+*> if (alfa > width) then
+*> if (debug) write(*,*) ' +++ width has been enlarged to', alfa, ' +++'
+*> endif
+*> if (alfa < 10 * dwn) then
+*> if (debug) then
+*> write(*,*) ' ... poor representation of the response function ...'
+*> write(*,*) 'the step size (', dwn, ') should be reduced'
+*> endif
+*> endif
+! make a table of the response function
+*> r1(1) = respon(0.0d0, alfa)
+*> r2(1) = r1(1)
+*> do i = 2, n
+*> x = (i - 1) * dwn
+*> r1(i) = respon( x, alfa)
+*> r2(i) = respon(-x, alfa)
+*> enddo
+! fourier transform it
+*> call rcffi(r1, r2, -m, dwn)
+! normalization of the response function, and convolution
+*> anorm = r1(1)
+*> do i = 1, n
+*> fac = dexp(-p1(i)) / anorm
+*> cp2 = fac * dcos(p2(i))
+*> sp2 = fac * dsin(p2(i))
+*> p1(i) = r1(i) * cp2 + r2(i) * sp2
+*> p2(i) = r2(i) * cp2 - r1(i) * sp2
+*> enddo
+*> alfa = alfa / 2
+
+1 else
+
+! *** broaden the spectrum by convoluting the characteristic function by
+! *** a weighted sum of a lorentzian and a gaussian response functions
+
+1 alfa = dmax1(1.5d0 * dwn, width)
+1 if (alfa > width) then
+*> if (debug) write(*,*) ' +++ width has been enlarged to', alfa, ' +++'
+*> endif
+1 if (alfa > 0.5d0 * wmax / pi) then
+*> stop '*** width is too large, nothing done ***'
+*> endif
+1 alfa = 0.5d0 * alfa
+1 sigma = alfa / 1.66511d0
+1 a1 = (1.0d0 - asym) / 2 * (1.0d0 - gauss)
+1 a2 = (1.0d0 + asym) / 2 * (1.0d0 - gauss)
+1 if (a1 < 0.0d0 .or. a2 < 0.0d0) then
+*> stop '*** invalid input : asym should be in (-1, +1) ***'
+*> endif
+1 g1 = (1.0d0 - asym) * alfa
+1 g2 = (1.0d0 + asym) * alfa
+1 p1(1) = 1.0d0
+1 p2(1) = 0.0d0
+1 do i = 2, n
+4095 x = (i - 1) * h
+4095 fr = 0.0d0
+4095 fi = 0.0d0
+4095 if (a1 /= 0.0d0) then
+4095 x1 = g1 * x
+4095 if (x1 <= 100.0d0) then
+4095 fr = a1 * dexp(-x1)
+4095 endif
+4095 if (a1 /= a2) then
+4095 if (x1 <= 1.3d0) then
+193 fi = a1 * o1(x1)
+3902 else
+3902 fi = a1 * o2(x1)
+3902 endif
+4095 endif
+4095 endif
+4095 if (a2 /= 0.0d0) then
+4095 x2 = g2 * x
+4095 if (x2 <= 100.0d0) then
+4095 fr = fr + a2 * dexp(-x2)
+4095 endif
+4095 if (a2 /= a1) then
+4095 if (x2 <= 1.3d0) then
+104 fi = fi - a2 * o1(x2)
+3991 else
+3991 fi = fi - a2 * o2(x2)
+3991 endif
+4095 endif
+4095 endif
+4095 if (gauss /= 0.0d0) then
+4095 x3 = sigma * x
+4095 if (x3 <= 10.0d0) then
+1736 fr = fr + gauss * dexp(-x3**2)
+1736 endif
+4095 endif
+4095 fi = fi / pi
+4095 fac = dexp(-p1(i))
+4095 cp2 = fac * dcos(p2(i))
+4095 sp2 = fac * dsin(p2(i))
+4095 p1(i) = fr * cp2 + fi * sp2
+4095 p2(i) = fi * cp2 - fr * sp2
+4095 enddo
+1 if (dabs(fi) > 1.0d-03) then
+*> if (debug) then
+*> write(*,*) ' ... poor representation of the response function ...'
+*> write(*,*) 'the step size (', dwn, ') should be reduced'
+*> endif
+*> endif
+1 endif
+
+! *** full eels spectrum
+
+1 call rcffi(p1, p2, m, h)
+
+! *** output
+
+1 istep = max(int(alfa / dwn / 10), 1)
+1 nout = 0
+1 if (emin <= 0.0d0) then
+1 imin = n - int(dabs(emin) / dwn)
+1 if ((imin - n) * dwn < emin) then
+*> imin = imin + 1
+*> endif
+1 do i = imin, n
+251 j = n - i
+251 if (mod(j, istep) > 0) cycle
+251 x = -j * dwn
+251 if (x > emax) exit
+251 nout = nout + 1
+251 xout(nout)= x
+251 yout(nout)= p2(j + 1)
+251 enddo
+1 endif
+1 if (x <= emax) then
+1 if (emax >= dwn) then
+1 imax = int(emax / dwn) + 1
+1 if ((imax - 1) * dwn > emax) then
+*> imax = imax - 1
+*> endif
+1 if (imax >= 2) then
+1 do i = 2, imax
+600 if (mod(i - 1, istep) > 0) cycle
+600 x = (i - 1) * dwn
+600 if (x < emin) cycle
+600 nout = nout + 1
+600 xout(nout) = x
+600 yout(nout) = p1(i)
+600 enddo
+1 endif
+1 endif
+1 endif
+1 if (debug) write(*,*) nout, ' values, step size =', istep * dwn
+
+! *** analyze the spectrum
+
+1 if (debug) write(*,*) ' peak location amplitude'
+1 wn = 0.0d0
+1 fm = p1(1)
+1 fp = p2(1)
+1 if (p2(2) < fp .and. p1(2) < fp) then
+1 if (debug) write(*, '(f13.2, e13.4)') wn, fp
+1 endif
+1 fac = 5.0d-06 * fp * dwn
+1 imax = 2 + int(dmax1(dabs(emin), dabs(emax)) / dwn)
+1 do i = 2, imax
+601 fm0 = fm1
+601 fp0 = fp1
+601 fm1 = fm
+601 fp1 = fp
+601 fm = p2(i)
+601 fp = p1(i)
+601 if (i == 2) cycle
+600 picm = .false.
+600 picp = .false.
+600 wn = (i - 1) * dwn
+600 if ((fm1 >= fm0) .and. (fm1 >= fm)) then
+2 a = (fm1 - fm0) + (fm1 - fm)
+2 if (a >= fac) then
+2 b = 0.5d0 * ((fm1 - fm0) + 3 * (fm1 - fm))
+2 u = b / a
+2 wmpic = -wn + u * dwn
+2 fmpic = fm + 0.5d0 * b * u
+2 picm = .true.
+2 endif
+2 endif
+600 if ((fp1 >= fp0) .and. (fp1 >= fp)) then
+2 a = (fp1 - fp0) + (fp1 - fp)
+2 if (a >= fac) then
+2 b = 0.5d0 * ((fp1 - fp0) + 3 * (fp1 - fp))
+2 u = b / a
+2 wppic = wn - u * dwn
+2 fppic = fp + 0.5d0 * b * u
+2 picp = .true.
+2 if (picp) then
+2 if (picm) then
+*> if (debug) write(*, '(2(f13.2, e13.4, 5x))') wppic, fppic, wmpic, fmpic
+2 else
+2 if (debug) write(*, '(f13.2, e13.4)') wppic, fppic
+2 endif
+2 endif
+2 endif
+2 endif
+600 if (picm .and. (.not. picp)) then
+2 if (debug) write(*, '(33x, f13.2, e13.4)') wmpic, fmpic
+2 endif
+600 enddo
+1 deallocate (p1)
+1 deallocate (p2)
+1 return
+end subroutine doboson
+double precision function o1(u)
+
+ implicit none
+
+ double precision, intent(in) :: u
+ double precision :: u2
+
+297 u2 = u**2
+297 o1 = -dsinh(u) * dlog(u2) + u * ((0.03114d0 * u2 + 0.41666d0) * u2 + 0.84557d0)
+
+297 return
+end function o1
+double precision function o2(u)
+
+ implicit none
+
+ double precision, intent(in) :: u
+ double precision :: u2
+
+7893 u2 = 1 / u**2
+7893 o2 = (((202.91d0 * u2 + 932.21d0) * u2 + 41.740d0) * u2 + 2.0d0) / &
+ (((540.88d0 * u2 + 345.67d0) * u2 + 18.961d0) * u2 + 1.0d0) / u
+
+7893 return
+end function o2
diff --git a/source/f90/fcat-analysis/eels_all.f90.fcat-analysis.txt b/source/f90/fcat-analysis/eels_all.f90.fcat-analysis.txt
new file mode 100644
index 0000000..1feb9ef
--- /dev/null
+++ b/source/f90/fcat-analysis/eels_all.f90.fcat-analysis.txt
@@ -0,0 +1,1353 @@
+program eels
+
+! ******************************************************************
+! * *
+! * compute the classical eels spectrum of an arbitrary plane- *
+! * statified medium made from isotropic materials in specular *
+! * geometry using the dielectric theory of eels. *
+! * *
+! ******************************************************************
+
+ implicit none
+
+ double precision :: e0, theta, phia, phib, wmin, wmax, dw
+ integer :: i, j, jos, k, l, layers, neps, nper, nw
+ logical :: user
+ character (len = 72) :: comment(2)
+ character (len = 10) :: contrl, mode
+
+ double precision, allocatable :: thick(:), epsinf(:), osc(:, :)
+ integer, allocatable :: nos(:)
+ character (len = 10), allocatable :: name(:)
+ double precision, allocatable :: wn_array(:), f(:)
+
+ integer :: old_size_1, old_size_2
+ double precision, allocatable :: tmp_osc(:, :)
+ integer :: ioStatus
+
+! integer, parameter :: name_length = 10
+! *** read spectrometer parameters
+
+1 call change_working_dir()
+1 open(unit = 11, file = 'eelsin')
+! impact energy (ev)
+1 read(11, *) e0
+! incidence angle (%)
+1 read(11, *) theta
+! angular apertures of the elliptic detector (%)
+1 read(11, *) phia
+1 read(11, *) phib
+! energy-loss interval and step size (cm**-1)
+1 read(11, *) wmin
+1 read(11, *) wmax
+1 read(11, *) dw
+! comment lines
+1 read(11, '(a72)') (comment(k), k = 1, 2)
+
+1 write(*,*) 'program eels (September 2020)'
+1 write(*,'(a, f6.2, a, f5.1, a, f5.2, a, f5.2, a)') &
+ ' e0 = ', e0, ' eV, theta = ', theta, '°, phia = ', &
+ phia, '°, phib = ', phib, '°'
+1 write(*,'(a, g11.4, a, g11.4, a, g11.4, a)') &
+ ' energy losses from', wmin, ' to', wmax, ', step = ', dw, ' cm**-1'
+1 write(*,*) comment(1)
+1 write(*,*) comment(2)
+
+1 if (phia <= 0.0d0 .or. phib <= 0.0d0) then
+*> stop '*** wrong input ***'
+*> endif
+1 if (e0 <= 0.0d0 .or. theta + phia >= 90.0d0) then
+*> stop '*** bad input ***'
+*> endif
+
+! *** read target specifications
+
+1 read(11, *) layers, nper, mode
+1 user = layers == 0
+1 if (.not. user) then
+1 neps = layers
+1 if (nper == -1) then
+*> neps = layers + 1
+*> nper = 1
+*> write(*,*) 'the substrate is a anisotropic uniaxial material'
+*> endif
+1 if (layers < 0 .or. nper < 1 .or. nper > layers) then
+*> stop '*** invalid target specifications ***'
+*> endif
+1 write(6, 101) layers, nper
+1 jos = 0
+1 allocate (name(neps), thick(neps), epsinf(neps), nos(neps))
+1 do l = 1, neps
+2 if (l <= layers) then
+2 read(11, 102) name(l), thick(l)
+2 endif
+2 read(11, *) epsinf(l), nos(l)
+2 write(6, 103)
+2 if (nos(l) <= 0) then
+*> if (l <= layers) write(6, 104) l, name(l), thick(l), epsinf(l)
+*> if (l > layers) write(6, 105) epsinf(l)
+2 else
+2 do j = 1, nos(l)
+2 if (.not. allocated(osc)) then
+1 allocate(osc(3, nos(l)))
+1 else if (j == 1) then
+ ! call extend3(osc, nos(j))
+1 old_size_1 = size(osc, 1)
+1 old_size_2 = size(osc, 2)
+1 allocate(tmp_osc(old_size_1, old_size_2 + nos(l)))
+1 tmp_osc(1:old_size_1, 1:old_size_2) = osc
+1 deallocate(osc)
+1 call move_alloc(tmp_osc, osc)
+1 endif
+2 jos = jos + 1
+2 read(11, *) (osc(k, jos), k = 1, 3)
+2 if ((j == nos(l) / 2 + 1) .and. (nos(l) > 1)) then
+*> write(6, 106)
+*> write(6, 107)
+*> endif
+2 if (j == 1) then
+2 if (l <= layers) then
+2 write(6, 104) l, name(l), thick(l), epsinf(l), (osc(i, jos), i = 1, 3)
+*> else
+*> write(6, 105) epsinf(l), (osc(i, jos), i = 1, 3)
+*> endif
+*> else
+*> write(6, 108) (osc(i, jos), i = 1, 3)
+*> endif
+2 enddo
+2 endif
+2 enddo
+1 write(*,*)
+1 read(11, 102, IOSTAT = ioStatus) contrl
+1 if (ioStatus /= 0) then
+1 contrl = ''
+1 endif
+1 endif
+
+1 close (unit = 11)
+
+1 nw = 1 + int((wmax - wmin) / dw)
+1 allocate (wn_array(nw), f(nw))
+
+1 call doeels(e0, theta, phia, phib, wmin, wmax, dw, comment, size(comment), &
+ layers, neps, nper, name, size(name), thick, epsinf, nos, osc, size(osc, 2),&
+ contrl, mode, wn_array, f, size(wn_array))
+
+1 open (unit = 12, file = 'eelsou')
+1 write (12, 207) e0, theta, phia, phib, comment(1)
+1 do i = 1, nw
+326 write (12, 211) wn_array(i), f(i)
+326 enddo
+1 close (unit = 12)
+
+1 stop
+
+*> 101 format(i3, ' layer(s), nper = ', i2//' l', 2x, 'material', 7x, &
+ 'thickness', 5x, 'epsinf', 4x, 'wto , wp', 5x, 'q', 7x, 'gam/wto')
+*> 102 format(a10, d15.5)
+*> 103 format(1x, 72('-'))
+*> 104 format(1x, i3, 2x, a10, g15.3, f10.4, f12.4, f10.4, f9.4)
+*> 105 format(31x, f10.4, f12.4, f10.4, f9.4)
+*> 106 format(45x, 'wlo , wp', 5x, 'q', 7x, 'gam/wlo')
+*> 107 format(45x, 28('-'))
+*> 108 format(41x, f12.4, f10.4, f9.4)
+*> 207 format('e0 =', f6.2, ' theta =', f5.1, ' phia =', f5.2, ' phib =', f5.2 / a72)
+*> 211 format(2e15.7)
+end program eels
+
+double precision function qrat(x, alpha, beta, c1, c2)
+
+ implicit none
+
+ double precision, intent(in) :: x, alpha, beta, c1, c2
+
+*> qrat = (1.0d0 + x * (beta + c1 * x)) / ((1.0d0 + x * (beta + c2 * x)) * (1.0d0 + alpha * x)**2)
+
+*> return
+end function qrat
+
+double precision function usurlo(dq, wn)
+
+! ******************************************************************
+! * *
+! * user-supplied dielectric surface loss function aimag(g(dq, wn)) *
+! * input arguments : *
+! * dq : modulus of the two-dimensional surface wave vector *
+! * (angstroem**-1) *
+! * wn : frequency (cm**-1) *
+! * *
+! ******************************************************************
+
+ implicit none
+
+ double precision, intent(in) :: dq
+ double precision, intent(in) :: wn
+
+*> usurlo = 1.0d0
+*> return
+end function usurlo
+double precision function surlos(dk, eps, thick, layers, nper)
+
+! ******************************************************************
+! * *
+! * eels surface loss function for an arbitrary multilayered target*
+! * *
+! ******************************************************************
+
+ implicit none
+
+ double precision, intent(in) :: dk
+ double complex, intent(in) :: eps(:)
+ double precision, intent(in) :: thick(:)
+ integer, intent(in) :: layers, nper
+
+ integer :: lmax
+ logical :: static, zero, skip
+ integer :: lstart, n
+ double precision, allocatable :: arg(:)
+ double precision :: argmin, argmax, cn, cnm1, epsmac, sn, snm1, tn
+ double complex :: a, b, csi, pnm2, pnm1, pn, pp, qnm2, qnm1, qn, qp, z
+
+ common / mulayr / argmin, argmax, epsmac
+
+ zero(z) = (dble(z) == 0.0d0) .and. (dimag(z) == 0.0d0)
+
+! write (*,*) 'surlos:'
+! write (*,*) 'thick: ', size(thick)
+! write (*,*) 'eps: ', size(eps)
+
+40156 lmax = size(eps)
+40156 allocate (arg(lmax))
+40156 lstart = layers - nper + 1
+40156 static = .true.
+40156 skip = .false.
+
+40156 do n = 1, layers
+63246 arg(n) = dk * thick(n)
+63246 if (arg(n) > argmax .or. zero(eps(n))) then
+17169 csi = eps(n)
+17169 skip = .true.
+17169 exit
+*> endif
+46077 static = .not. (n >= lstart .and. arg(n) > argmin)
+46077 enddo
+
+40156 if (.not. skip) then
+! *** periodic continued fraction, period = nper
+
+22987 do ! Do loop is only for the option to skip out of it.
+22987 if (nper <= 1) then
+22987 csi = eps(layers)
+*> else
+*> if (static) then
+*> pn = 0.0d0
+*> qn = 0.0d0
+*> do n = lstart, layers
+*> pn = pn + thick(n) * eps(n)
+*> qn = qn + thick(n) / eps(n)
+*> enddo
+*> if (zero(qn)) then
+*> n = lstart
+*> if (n <= 1) then
+*> surlos = 0.0d0
+*> return
+*> endif
+*> n = n - 1
+*> csi = eps(n) / dtanh(arg(n))
+*> exit
+*> endif ! zero(qn)
+*> csi = cdsqrt(pn / qn)
+*> if ((dimag(csi) < 0.0d0) .and. (dble(qn) < 0.0d0)) then
+*> csi = -csi
+*> endif
+*> else ! static
+*> cn = dcosh(arg(lstart))
+*> sn = dsinh(arg(lstart))
+*> pnm1 = 1.0d0
+*> pn = cn
+*> pp = eps(lstart) * sn
+*> qnm1 = 0.0d0
+*> qn = sn / eps(lstart)
+*> qp = pn
+*> do n = lstart + 1, layers
+*> cnm1 = cn
+*> snm1 = sn
+*> cn = dcosh(arg(n))
+*> sn = dsinh(arg(n))
+*> a = eps(n) * sn
+*> pp = cn * pp + a * pn
+*> qp = cn * qp + a * qn
+*> b = (eps(n - 1) / eps(n)) * (sn / snm1)
+*> a = cnm1 * b + cn
+*> pnm2 = pnm1
+*> pnm1 = pn
+*> qnm2 = qnm1
+*> qnm1 = qn
+*> pn = a * pnm1 - b * pnm2
+*> qn = a * qnm1 - b * qnm2
+*> enddo
+*> if (zero(qn)) then
+*> a = qp - pn
+*> if (zero(a)) then
+*> n = lstart
+*> if (n <= 1) then
+*> surlos = 0.0d0
+*> return
+*> endif
+*> n = n - 1
+*> csi = eps(n) / dtanh(arg(n))
+*> exit
+*> endif
+*> csi = pp / a
+*> else
+*> a = 0.5d0 * (pn - qp) / qn
+*> b = cdsqrt(a**2 + pp / qn)
+*> pn = a - pn / qn
+*> if (cdabs(pn + b) > cdabs(pn - b)) then
+*> b = -b
+*> endif
+*> csi = a + b
+*> endif ! zero(qn)
+*> endif ! static
+! *** small-dk limit of the periodic tail
+
+*> endif ! nper <= 1
+
+22987 n = lstart
+
+22987 exit
+*> enddo
+22987 endif ! .not. skip
+
+! *** backward algorithm
+
+40156 do
+63246 n = n - 1
+63246 if (n <= 0) then
+40156 a = csi + 1.0d0
+40156 if (zero(a)) then
+*> surlos = 2 / epsmac
+40156 else
+40156 surlos = dimag(-2 / a)
+40156 endif
+40156 return
+*> endif ! n <= 0
+
+23090 if (arg(n) /= 0.0d0) then
+23090 tn = dtanh(arg(n))
+23090 b = eps(n) + csi * tn
+23090 if (zero(b)) then
+*> surlos = 0.0d0
+*> return
+*> endif
+23090 csi = eps(n) * (csi + tn * eps(n)) / b
+23090 endif
+23090 enddo
+end function surlos
+double precision function phint(phi, a, u)
+
+! ******************************************************************
+! * *
+! * evaluate the integral from zero to phi of *
+! * *
+! * u 2 *
+! * ( ----------------------------- ) dphi *
+! * 2 2 *
+! * (1 - a * u * cos(phi)) + u *
+! * *
+! * for 0 <= phi <= pi , u >= 0 and a >= 0 *
+! * *
+! ******************************************************************
+
+ implicit none
+
+ double precision, intent(in) :: phi
+ double precision, intent(in) :: a
+ double precision, intent(in) :: u
+
+ double precision :: ai, ar, bi, br, c, cpr, d, e, esr, pi, qr, ri, rm, root
+ double precision :: rp, rr, s, spr, tm, tp, u2, x, zeta, zetai, zetar, zr
+
+41010 pi = 3.141592653589793238d0
+41010 c = dcos(phi)
+41010 s = dsin(phi)
+41010 u2 = u**2
+41010 e = a*u
+41010 if (u < 1.0d0 .and. e < 1.0d-02 * (1.0d0 + u2)) then
+*> zr = 1.0d0 + u2
+*> esr = e / zr
+*> phint = u2 / zr**2 * ((( (4.0d0 / 3.0d0) * (2.0d0 + c**2) * s * (5.0d0 - 3 * u2) * &
+ esr + (phi + c * s) * (5.0d0 - u2)) * esr + 4 * s) * esr + phi)
+41010 else
+41010 rm = dsqrt((1.0d0 - e)**2 + u2)
+41010 tm = 0.5d0 * datan2(u, 1.0d0 - e)
+41010 rp = dsqrt((1.0d0 + e)**2 + u2)
+41010 tp = 0.5d0 * datan2(u, 1.0d0 + e)
+41010 root = dsqrt(rm * rp)
+41010 cpr = dcos(tm + tp)
+41010 spr = dsin(tm + tp)
+41010 x = 0.0d0 ! ensure initialization
+41010 if (c >= 0.0d0) then
+20833 x = s / (1.0d0 + c)
+19851 elseif (dabs(s) > 1.0d-07) then
+19851 x = (1.0d0 - c) / s
+19851 endif
+41010 if ((c >= 0.0d0) .or. (dabs(s) > 1.0d-07)) then
+40684 zeta = dsqrt(rm / rp)
+40684 zetar = -zeta * dsin(tm - tp)
+40684 zetai = zeta * dcos(tm - tp)
+40684 br = 0.5d0 * dlog(((zetar + x)**2 + zetai**2) / ((zetar - x)**2 + zetai**2))
+40684 bi = datan2(zetai, zetar + x) - datan2(zetai, zetar - x)
+40684 rr = -(br * spr - bi * cpr) / root
+40684 ri = -(bi * spr + br * cpr) / root
+40684 d = e * s / ((1.0d0 - e * c)**2 + u2)
+40684 ar = d * (1.0d0 - e * c) - rr + u * ri
+40684 ai = -d * u - ri - u * rr
+326 else
+326 rr = -pi / root * cpr
+326 ri = pi / root * spr
+326 ar = -rr + u * ri
+326 ai = -ri - u * rr
+326 endif
+41010 qr = (ar * (cpr - spr) * (cpr + spr) + 2 * ai * cpr * spr) / (rm * rp)
+41010 phint = 0.5d0 * (ri / u - qr)
+41010 endif
+41010 return
+end function phint
+double precision function fint1(u, eps, thick, layers, nper, eps_size)
+
+! ******************************************************************
+! * *
+! * integration over the azimutal angle from 0.0 to pi *
+! * *
+! ******************************************************************
+
+ implicit none
+
+ double precision, intent(in) :: u
+ integer, intent(in) :: layers, nper, eps_size
+ double complex, intent(in) :: eps(eps_size)
+ double precision, intent(in) :: thick(eps_size)
+
+ logical :: rational, user
+ double precision :: acoef, bcoef, ccoef, cospsi, den, dif, dlimf, e, elleps
+ double precision :: pi, rom, rop, ru, sinpsi, sum, um, t
+ double precision :: tanpsi, wn, u2
+
+ interface
+ double precision function usurlo(dq, wn)
+ double precision, intent(in) :: dq
+ double precision, intent(in) :: wn
+ end function usurlo
+ double precision function surlos(dk, eps, thick, layers, nper)
+ double precision, intent(in) :: dk
+ double complex, intent(in) :: eps(:)
+ double precision, intent(in) :: thick(:)
+ integer, intent(in) :: layers, nper
+ end function surlos
+ end interface
+
+ common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi, &
+ ru, um, dlimf, wn, user, rational
+
+ data pi / 3.141592653589793238d0 /
+
+! write (*,*) 'fint1:'
+! write (*,*) 'thick: ', size(thick)
+! write (*,*) 'eps: ', size(eps)
+
+14598 if (u == 0.0d0) then
+326 fint1 = 0.0d0
+326 return
+*> endif
+14272 e = tanpsi * u
+14272 u2 = u**2
+14272 rom = (1.0d0 - e)**2 + u2
+14272 rop = (1.0d0 + e)**2 + u2
+14272 sum = rop + rom
+14272 rom = dsqrt(rom)
+14272 rop = dsqrt(rop)
+14272 dif = rop - rom
+14272 den = dsqrt((2.0d0 - dif) * (2.0d0 + dif)) * rop * rom
+14272 fint1 = pi * u2 * (4 * sum - dif**2 * (sum - rop * rom)) / den**3
+14272 if (rational) then
+*> return
+*> endif
+14272 if (user) then
+*> fint1 = fint1 * usurlo(ru * u, wn)
+14272 else
+14272 fint1 = fint1 * surlos(ru * u, eps, thick, layers, nper)
+14272 if (dlimf > 0.0d0) then
+*> t = ru * u * dlimf
+*> fint1 = fint1 * (1.d0 + t * dlog(t / (t + 0.26d0)))**2 / (1.d0 + 1.40d0 * t)
+*> endif
+14272 endif
+14272 return
+end function fint1
+double precision function fint2(u, eps, thick, layers, nper, eps_size)
+
+! ******************************************************************
+! * *
+! * integration over the azimutal angle from 0.0 to phi < pi *
+! * *
+! ******************************************************************
+
+ implicit none
+
+ double precision, intent(in) :: u
+ integer, intent(in) :: layers, nper, eps_size
+ double complex, intent(in) :: eps(eps_size)
+ double precision, intent(in) :: thick(eps_size)
+
+ logical :: rational, user
+ double precision :: a, arg, b, b2, c, ccoef, cospsi, dlimf, elleps, phi
+ double precision :: phint, ru, sinpsi, um, t, tanpsi, wn, x
+
+ interface
+ double precision function usurlo(dq, wn)
+ double precision, intent(in) :: dq
+ double precision, intent(in) :: wn
+ end function usurlo
+ double precision function surlos(dk, eps, thick, layers, nper)
+ double precision, intent(in) :: dk
+ double complex, intent(in) :: eps(:)
+ double precision, intent(in) :: thick(:)
+ integer, intent(in) :: layers, nper
+ end function surlos
+ end interface
+
+ common / param / a, b, ccoef, elleps, cospsi, sinpsi, tanpsi, &
+ ru, um, dlimf, wn, user, rational
+
+! write (*,*) 'fint2:'
+! write (*,*) 'thick: ', size(thick)
+! write (*,*) 'eps: ', size(eps)
+
+10758 if (u == 0.0d0) then
+*> fint2 = 0.0d0
+*> return
+*> endif
+10758 b2 = b**2
+10758 c = (1.0d0 - elleps) * (cospsi * u)**2 + (b - ccoef) * (b + ccoef)
+10758 if (dabs(a * c) > 1.0d-03 * b2) then
+10758 x = (b - dsqrt(b2 - a * c)) / a
+*> else
+*> x = a * c / b2
+*> x = 0.5d0 * c * (1.d0 + 0.25d0 * x * (1.d0 + 0.5d0 * x * (1.d0 + 0.625d0 * x))) / b
+*> endif
+10758 arg = x / u
+10758 if (dabs(arg) > 1.0d0) then
+74 arg = dsign(1.0d0, arg)
+74 endif
+10758 phi = dacos(arg)
+10758 fint2 = phint(phi, tanpsi, u)
+10758 if (rational) then
+*> return
+*> endif
+10758 if (user) then
+*> fint2 = fint2 * usurlo(ru * u, wn)
+10758 else
+10758 fint2 = fint2 * surlos(ru * u, eps, thick, layers, nper)
+10758 if (dlimf > 0.0d0) then
+*> t = ru * u * dlimf
+*> fint2 = fint2 * (1.d0 + t * dlog(t / (t + 0.26d0)))**2 / (1.d0 + 1.40d0 * t)
+*> endif
+10758 endif
+10758 return
+end function fint2
+double precision function fint3(u, eps, thick, layers, nper, eps_size)
+
+! ******************************************************************
+! * *
+! * integration over the azimutal angle from phi1 > 0 to phi2 < pi *
+! * *
+! ******************************************************************
+
+ implicit none
+
+ double precision, intent(in) :: u
+ integer, intent(in) :: layers, nper, eps_size
+ double complex, intent(in) :: eps(eps_size)
+ double precision, intent(in) :: thick(eps_size)
+
+ logical :: rational, user
+ double precision :: a, arg, b, ccoef, cospsi, dlimf, elleps, phi1, phi2
+ double precision :: phint, sinpsi, rac, ru, um, t, tanpsi, wn
+
+ interface
+ double precision function usurlo(dq, wn)
+ double precision, intent(in) :: dq
+ double precision, intent(in) :: wn
+ end function usurlo
+ double precision function surlos(dk, eps, thick, layers, nper)
+ double precision, intent(in) :: dk
+ double complex, intent(in) :: eps(:)
+ double precision, intent(in) :: thick(:)
+ integer, intent(in) :: layers, nper
+ end function surlos
+ end interface
+
+ common / param / a, b, ccoef, elleps, cospsi, sinpsi, tanpsi, &
+ ru, um, dlimf, wn, user, rational
+
+! write (*,*) 'fint3:'
+! write (*,*) 'thick: ', size(thick)
+! write (*,*) 'eps: ', size(eps)
+
+15126 if (u == 0.0d0) then
+*> fint3 = 0.0d0
+*> return
+*> endif
+15126 rac = dsign(1.0d0, a) * cospsi * dsqrt((1.0d0 - elleps) * a * (um - u) * (um + u))
+15126 arg = (b - rac) / (u * a)
+15126 if (dabs(arg) > 1.0d0) arg = dsign(1.0d0, arg)
+15126 phi2 = dacos(arg)
+15126 fint3 = phint(phi2, tanpsi, u)
+15126 arg = (b + rac) / (u * a)
+15126 if (dabs(arg) > 1.0d0) arg = dsign(1.0d0, arg)
+15126 phi1 = dacos(arg)
+15126 fint3 = fint3 - phint(phi1, tanpsi, u)
+15126 if (rational) return
+15126 if (user) then
+*> fint3 = fint3 * usurlo(ru * u, wn)
+15126 else
+15126 fint3 = fint3 * surlos(ru * u, eps, thick, layers, nper)
+15126 if (dlimf > 0.0d0) then
+*> t = ru * u * dlimf
+*> fint3 = fint3 * (1.d0 + t * dlog(t / (t + 0.26d0)))**2 / (1.d0 + 1.40d0 * t)
+*> endif
+15126 endif
+15126 return
+end function fint3
+double precision function fun(phi)
+
+! ******************************************************************
+! * *
+! * integrand of the expression of the 1st order term in the *
+! * expansion of the eels integral for a homogeneous target. *
+! * *
+! ******************************************************************
+
+ implicit none
+
+ double precision, intent(in) :: phi
+
+ logical :: user, rational
+ double precision :: acoef, bcoef, ccoef, cospsi, dlimf, elleps, ru
+ double precision :: sinphi, sinpsi, tanpsi, um, wn
+
+ common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi, &
+ ru, um, dlimf, wn, user, rational
+
+*> sinphi = dsin(phi)
+*> fun = dsqrt((1.0d0 - elleps + elleps * sinphi**2) * &
+ (1.0d0 - sinpsi * sinphi) * &
+ (1.0d0 + sinpsi * sinphi))
+*> return
+end function fun
+subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, thick, layers, nper)
+
+! estimate the integral of fun(x) from a to b
+! to a user provided tolerance.
+! an automatic adaptive routine based on
+! the 8-panel newton-cotes rule (g. forsythe et al, 1977, p. 92)
+!
+! input ..
+!
+! fun the name of the integrand function subprogram fun(x).
+! a the lower limit of integration.
+! b the upper limit of integration.(b may be less than a.)
+! relerr a relative error tolerance. (should be non-negative)
+! abserr an absolute error tolerance. (should be non-negative)
+!
+! output ..
+!
+! result an approximation to the integral hopefully satisfying the
+! least stringent of the two error tolerances.
+! errest an estimate of the magnitude of the actual error.
+! nofun the number of function values used in calculation of result.
+! flag a reliability indicator. if flag is zero, then result
+! probably satisfies the error tolerance. if flag is
+! xxx.yyy , then xxx = the number of intervals which have
+! not converged and 0.yyy = the fraction of the interval
+! left to do when the limit on nofun was approached.
+
+ implicit none
+
+ double precision :: fun
+ double precision, intent(in) :: a
+ double precision, intent(in) :: b
+ double precision, intent(in out) :: abserr
+ double precision, intent(in) :: relerr
+ double precision, intent(out) :: result
+ double precision, intent(out) :: errest
+ integer, intent(out) :: nofun
+ double precision, intent(out) :: flag
+
+ external fun
+
+ double precision, intent(in) :: thick(:)
+ double complex, intent(in) :: eps(:)
+ integer, intent(in) :: layers, nper
+
+ double precision :: w0, w1, w2, w3, w4, area, x0, f0, stone, step, cor11, temp
+ double precision :: qprev, qnow, qdiff, qleft, esterr, tolerr
+ double precision :: qright(31), f(16), x(16), fsave(8, 30), xsave(8, 30)
+ double precision :: dabs, dmax1
+
+ integer :: levmin, levmax, levout, nomax, nofin, lev, nim, i, j
+
+! *** stage 1 *** general initialization
+! set constants.
+
+! write (*,*) 'quanc8:'
+! write (*,*) 'thick: ', size(thick)
+! write (*,*) 'eps: ', size(eps)
+
+978 levmin = 1
+978 levmax = 30
+978 levout = 6
+978 nomax = 5000
+978 nofin = nomax - 8 * (levmax - levout + 2**(levout + 1))
+
+! trouble when nofun reaches nofin
+
+978 w0 = 3956.0d0 / 14175.0d0
+978 w1 = 23552.0d0 / 14175.0d0
+978 w2 = -3712.0d0 / 14175.0d0
+978 w3 = 41984.0d0 / 14175.0d0
+978 w4 = -18160.0d0 / 14175.0d0
+
+! initialize running sums to zero.
+
+978 flag = 0.0d0
+978 result = 0.0d0
+978 cor11 = 0.0d0
+978 errest = 0.0d0
+978 area = 0.0d0
+978 nofun = 0
+978 if (a == b) return
+
+! *** stage 2 *** initialization for first interval
+
+978 lev = 0
+978 nim = 1
+978 x0 = a
+978 x(16) = b
+978 qprev = 0.0d0
+978 f0 = fun(x0, eps, thick, layers, nper, size(eps))
+978 stone = (b - a) / 16
+978 x(8) = (x0 + x(16)) / 2
+978 x(4) = (x0 + x(8)) / 2
+978 x(12) = (x(8) + x(16)) / 2
+978 x(2) = (x0 + x(4)) / 2
+978 x(6) = (x(4) + x(8)) / 2
+978 x(10) = (x(8) + x(12)) / 2
+978 x(14) = (x(12) + x(16)) / 2
+978 do j = 2, 16, 2
+7824 f(j) = fun(x(j), eps, thick, layers, nper, size(eps))
+7824 enddo
+978 nofun = 9
+
+978 do
+
+! *** stage 3 *** central calculation
+! requires qprev, x0, x2, x4, ..., x16, f0, f2, f4, ..., f16.
+! calculates x1, x3, ...x15, f1, f3, ...f15, qleft, qright, qnow, qdiff, area.
+
+3960 x(1) = (x0 + x(2)) / 2
+3960 f(1) = fun(x(1), eps, thick, layers, nper, size(eps))
+3960 do j = 3, 15, 2
+27720 x(j) = (x(j - 1) + x(j + 1)) / 2
+27720 f(j) = fun(x(j), eps, thick, layers, nper, size(eps))
+27720 enddo
+3960 nofun = nofun + 8
+3960 step = (x(16) - x0) / 16.0d0
+3960 qleft = (w0 * (f0 + f(8)) + w1 * (f(1) + f(7)) + w2 * (f(2) + f(6)) &
+ + w3 * (f(3) + f(5)) + w4 * f(4)) * step
+3960 qright(lev + 1) = (w0 * (f(8) + f(16)) + w1 * (f(9) + f(15)) + w2 * (f(10) + f(14)) &
+ + w3 * (f(11) + f(13)) + w4 * f(12)) * step
+3960 qnow = qleft + qright(lev + 1)
+3960 qdiff = qnow - qprev
+3960 area = area + qdiff
+
+! *** stage 4 *** interval convergence test
+
+3960 esterr = dabs(qdiff) / 1023
+3960 tolerr = dmax1(abserr, relerr * dabs(area)) * (step / stone)
+
+3960 if (lev >= levmin) then
+2982 if (lev >= levmax) then
+! current level is levmax.
+*> flag = flag + 1.0d0
+2982 else
+2982 if (nofun > nofin) then
+! *** stage 6 *** trouble section
+! number of function values is about to exceed limit.
+*> nofin = 2 * nofin
+*> levmax = levout
+*> flag = flag + (b - x0) / (b - a)
+2982 else
+2982 if (esterr > tolerr) then
+! *** stage 5 *** no convergence
+! locate next interval.
+513 nim = 2 * nim
+513 lev = lev + 1
+! store right hand elements for future use.
+513 do i = 1, 8
+4104 fsave(i, lev) = f(i + 8)
+4104 xsave(i, lev) = x(i + 8)
+4104 enddo
+! assemble left hand elements for immediate use.
+513 qprev = qleft
+513 do i = 1, 8
+4104 f(18 - 2 * i) = f(9 - i)
+4104 x(18 - 2 * i) = x(9 - i)
+4104 enddo
+513 cycle
+*> endif
+2469 endif
+2469 endif
+
+! *** stage 7 *** interval converged
+! add contributions into running sums.
+2469 result = result + qnow
+2469 errest = errest + esterr
+2469 cor11 = cor11 + qdiff / 1023
+! locate next interval.
+2469 do while (nim /= 2 * (nim / 2))
+2469 nim = nim / 2
+2469 lev = lev - 1
+2469 enddo
+2469 nim = nim + 1
+
+2469 if (lev <= 0) exit
+
+! assemble elements required for the next interval.
+1491 qprev = qright(lev)
+1491 x0 = x(16)
+1491 f0 = f(16)
+1491 do i = 1, 8
+11928 f(2*i) = fsave(i, lev)
+11928 x(2*i) = xsave(i, lev)
+11928 enddo
+1491 cycle
+978 else
+! *** stage 5 *** no convergence
+! locate next interval.
+978 nim = 2 * nim
+978 lev = lev + 1
+! store right hand elements for future use.
+978 do i = 1, 8
+7824 fsave(i, lev) = f(i + 8)
+7824 xsave(i, lev) = x(i + 8)
+7824 enddo
+! assemble left hand elements for immediate use.
+978 qprev = qleft
+978 do i = 1, 8
+7824 f(18 - 2 * i) = f(9 - i)
+7824 x(18 - 2 * i) = x(9 - i)
+7824 enddo
+978 endif
+
+978 enddo
+
+ ! *** stage 8 *** finalize and return
+978 result = result + cor11
+
+! make sure errest not less than roundoff level.
+978 if (errest /= 0.0d0) then
+978 temp = dabs(result) + errest
+978 do while (temp == dabs(result))
+*> errest = 2 * errest
+*> temp = dabs(result) + errest
+*> enddo
+978 endif
+978 return
+end subroutine quanc8
+subroutine queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
+
+! ******************************************************************
+! * *
+! * perform q-space integration for computing the eels spectrum of *
+! * a isotropic target using polar coordinates. *
+! * *
+! * x is the dimensionless energy loss hbar*omega/(2*e0*phia) *
+! * aerr and rerr are the desired absolute and relative accuracies *
+! * facru*x is the units of wavevectors omega/v_perpendicular *
+! * f is the q-integral multiplied by (2/pi)**2 *
+! * *
+! ******************************************************************
+
+ implicit none
+
+ double precision, intent(in) :: x
+ double precision, intent(out) :: f
+ double precision, intent(in out) :: aerr
+ double precision, intent(in out) :: rerr
+ double precision, intent(in) :: facru
+ double complex, intent(in) :: eps(:)
+ double precision, intent(in) :: thick(:)
+ integer, intent(in) :: layers, nper
+
+ logical :: rational, user
+ double precision :: acoef, bcoef, ccoef, cospsi, dlimf, elleps
+ double precision :: error, flag, ru, sinpsi
+ double precision :: u1, u2, um, ut, tanpsi, wn, y
+ integer :: ie, nofu
+ dimension error(3), flag(3)
+
+ interface
+ double precision function fint1(u, eps, thick, layers, nper, eps_size)
+ double precision, intent(in) :: u
+ integer, intent(in) :: layers, nper, eps_size
+ double complex, intent(in) :: eps(eps_size)
+ double precision, intent(in) :: thick(eps_size)
+ end function fint1
+ double precision function fint2(u, eps, thick, layers, nper, eps_size)
+ double precision, intent(in) :: u
+ integer, intent(in) :: layers, nper, eps_size
+ double complex, intent(in) :: eps(eps_size)
+ double precision, intent(in) :: thick(eps_size)
+ end function fint2
+ double precision function fint3(u, eps, thick, layers, nper, eps_size)
+ double precision, intent(in) :: u
+ integer, intent(in) :: layers, nper, eps_size
+ double complex, intent(in) :: eps(eps_size)
+ double precision, intent(in) :: thick(eps_size)
+ end function fint3
+ subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, thick, layers, nper)
+ double precision :: fun
+ double precision, intent(in) :: a
+ double precision, intent(in) :: b
+ double precision, intent(in out) :: abserr
+ double precision, intent(in) :: relerr
+ double precision, intent(out) :: result
+ double precision, intent(out) :: errest
+ integer, intent(out) :: nofun
+ double precision, intent(out) :: flag
+
+ external fun
+
+ double precision, intent(in) :: thick(:)
+ double complex, intent(in) :: eps(:)
+ integer, intent(in) :: layers, nper
+ end subroutine quanc8
+ end interface
+
+ common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi, &
+ ru, um, dlimf, wn, user, rational
+
+! write (*,*) 'queels:'
+! write (*,*) 'eps: ', size(eps)
+! write (*,*) 'thick: ', size(thick)
+
+326 f = 0.0d0
+326 if (x <= 0.0d0) then
+*> return
+*> endif
+326 ru = facru * x
+326 ccoef = cospsi**2 / x
+326 ut = ccoef - bcoef
+326 u1 = dabs(ut)
+326 u2 = ccoef + bcoef
+326 if (ut > 0.0d0) then
+326 call quanc8(fint1, 0.0d0, u1, aerr, rerr, y, error(1), nofu, flag(1), eps, thick, layers, nper)
+326 f = y
+*> else
+*> flag(1) = 0.0d0
+*> endif
+326 if (u2 > u1) then
+326 call quanc8(fint2, u1, u2, aerr, rerr, y, error(2), nofu, flag(2), eps, thick, layers, nper)
+326 f = f + y
+*> else
+*> flag(2) = 0.0d0
+*> endif
+326 if (dabs(acoef) > x * (1.0d0 - elleps) * bcoef) then
+326 um = dsqrt(ccoef / x / (1.0d0 - elleps) + bcoef**2 / acoef)
+326 if (um > u2) then
+326 call quanc8(fint3, u2, um, aerr, rerr, y, error(3), nofu, flag(3), eps, thick, layers, nper)
+326 f = f + y
+326 endif
+326 if (um < u1) then
+*> call quanc8(fint3, um, u1, aerr, rerr, y, error(3), nofu, flag(3), eps, thick, layers, nper)
+*> f = f - y
+*> endif
+*> else
+*> flag(3) = 0.0d0
+*> endif
+326 do ie = 1, 3
+978 if (flag(ie) == 0.0d0) cycle
+*> write(*,*) ' +++ flag(', ie, ') =', flag(ie), ', error =', error(ie), ' +++'
+*> if (flag(ie) - aint(flag(ie)) > 0.5d-02) then
+*> stop '*** execution aborted ***'
+*> endif
+*> enddo
+326 f = (2 / 3.141592653589793238d0)**2 * f
+326 return
+end subroutine queels
+subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
+
+! ******************************************************************
+! * *
+! * set up long-wavelength dielectric functions of the layers for *
+! * the present frequency wn (in cm**-1) *
+! * *
+! ******************************************************************
+
+ implicit none
+ integer, intent(in) :: neps
+ integer, intent(in) :: nos(:)
+ double precision, intent(in) :: osc(:, :)
+ double precision, intent(in) :: epsinf(:)
+ double precision, intent(in) :: wn
+ character (len=10), intent(in) :: name(:)
+ character (len=10), intent(in) :: mode
+
+ double complex, intent(in out) :: eps(:)
+ integer, intent(in) :: layers
+
+ double precision :: argmin, argmax, epsmac, x
+ double complex :: deno, nomi
+ integer :: j, k, l, m
+
+ common / mulayr / argmin, argmax, epsmac
+
+! write (*,*) 'seteps:'
+! write (*,*) 'nos: ', size(nos)
+! write (*,*) 'osc: ', size(osc)
+! write (*,*) 'epsinf: ', size(epsinf)
+! write (*,*) 'name: ', size(name)
+! write (*,*) 'eps: ', size(eps)
+! write (*,*) 'thick: ', size(thick)
+
+326 j = 0
+326 do l = 1, neps
+652 m = nos(l)/2
+652 nomi = dcmplx(1.0d0, 0.0d0)
+652 deno = dcmplx(1.0d0, 0.0d0)
+652 if (mode == 'kurosawa') then
+*> do k = 1, m
+*> j = j + 1
+! since osc(1,*) and wn are real, the following should be equivalent
+! Check required.
+! Furthermore the first term can be rewritten as
+! (osc(1, j + m) - wn) * (osc(1, j + m) + wn)
+! nomi = nomi * cmplx(osc(1, j + m)**2 - wn**2, - wn * osc(3, j + m))
+*> nomi = nomi * (osc(1, j + m)**2 - wn**2 - dcmplx(0.0d0, wn * osc(3, j + m)))
+*> deno = deno * (osc(1, j )**2 - wn**2 - dcmplx(0.0d0, wn * osc(3, j )))
+*> enddo
+*> eps(l) = epsinf(l) * nomi / deno
+*> elseif (name(l) == 'metal') then
+*> j = j + 1
+! suggestion for replacement
+! eps(l) = -osc(1, j)**2 / cmplx(wn**2, wn * osc(3, j))
+*> eps(l) = -osc(1, j)**2 / ( wn**2 + dcmplx(0.0d0, wn * osc(3, j)) )
+652 else
+652 eps(l) = epsinf(l)
+! The original version had this additional loop. It seems, it has been removed
+! because all cases of nos(l) > 1 are treated in case 1 above
+652 do k = 1, nos(l)
+652 j = j + 1
+652 x = wn / osc(1, j)
+652 deno = x * dcmplx(x, osc(3, j))
+652 if (osc(2, j) >= 0.0d0) then
+326 deno = 1.0d0 - deno
+326 endif
+652 if (cdabs(deno) == 0.0d0) then
+*> deno = epsmac
+*> endif
+652 eps(l) = eps(l) + osc(2, j) / deno
+652 enddo
+652 endif
+652 enddo
+326 if (neps == layers + 1) then
+! the substrate is a anisotropic uniaxial material
+*> eps(layers) = cdsqrt(eps(layers) * eps(layers + 1))
+*> if (dimag(eps(layers)) < 0.0d0) then
+*> eps(layers) = -eps(layers)
+*> endif
+*> endif
+326 return
+end subroutine seteps
+subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size, &
+ layers, neps, nper, name, name_size, thick, epsinf, nos, osc, osc_size,&
+ contrl, mode, wn_array, f_array, wn_array_size)
+
+! ******************************************************************
+! * *
+! * compute the classical eels spectrum of an arbitrary plane- *
+! * statified medium made from isotropic materials in specular *
+! * geometry using the dielectric theory of eels. *
+! * *
+! ******************************************************************
+
+ implicit none
+
+ integer, parameter :: nt = 5
+
+ double precision, intent(in) :: e0, theta, phia, phib, wmin, wmax, dw
+ character (len = 72) :: comment(comment_size)
+ character (len = 10) :: name(name_size)
+ double precision, intent(in) :: thick(name_size), epsinf(name_size), osc(3, osc_size)
+ character (len = 10) :: contrl, mode
+ integer, intent(in) :: comment_size, name_size, osc_size, wn_array_size
+ integer, intent(in out) :: layers, nper, nos(name_size)
+ double precision, intent(in out) :: wn_array(wn_array_size), f_array(wn_array_size)
+
+ logical :: rational, user, debug
+ integer :: i, iw, neps, nofu, nout, nw, lmax
+ double precision :: a, acoef, aerr, alpha, argmin, argmax, b, bcoef, beta, &
+ c1, c2, ccoef, cospsi, dlimf, dx, elleps, ener, epsmac, facru, f, f0, &
+ f1, fpic, fun, pi, prefac, psia, psii, qrat, rerr, ru, sinpsi, t, &
+ tanpsi, table, um, widt, wn, wpic, x, xmin, xmax, z, z1, z2
+ double complex, allocatable :: eps(:)
+ dimension table(nt)
+
+ external fun, qrat
+
+ interface
+ subroutine queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
+ double precision, intent(in) :: x
+ double precision, intent(out) :: f
+ double precision, intent(in out) :: aerr
+ double precision, intent(in out) :: rerr
+ double precision, intent(in) :: facru
+ double complex, intent(in) :: eps(:)
+ double precision, intent(in) :: thick(:)
+ integer, intent(in) :: layers, nper
+ end subroutine queels
+ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
+ integer, intent(in) :: neps
+ integer, intent(in) :: nos(:)
+ double precision, intent(in) :: osc(:, :)
+ double precision, intent(in) :: epsinf(:)
+ double precision, intent(in) :: wn
+ character (len=10), intent(in) :: name(:)
+ character (len=10), intent(in) :: mode
+ double complex, intent(in out) :: eps(:)
+ integer, intent(in) :: layers
+ end subroutine seteps
+ subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, thick, layers, nper)
+ double precision :: fun
+ double precision, intent(in) :: a
+ double precision, intent(in) :: b
+ double precision, intent(in out) :: abserr
+ double precision, intent(in) :: relerr
+ double precision, intent(out) :: result
+ double precision, intent(out) :: errest
+ integer, intent(out) :: nofun
+ double precision, intent(out) :: flag
+
+ external fun
+
+ double precision, intent(in) :: thick(:)
+ double complex, intent(in) :: eps(:)
+ integer, intent(in) :: layers, nper
+ end subroutine quanc8
+ end interface
+
+ common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi, &
+ ru, um, dlimf, wn, user, rational
+ common / mulayr / argmin, argmax, epsmac
+
+ data aerr / 0.0d0 /, rerr / 1.0d-06 /, f / 0.0d0 /, f1 / 0.0d0 /
+
+1 debug = .false.
+1 if (debug) then
+*> write (*,*) 'doeels:'
+*> write (*,*) 'comment: ', size(comment)
+*> write (*,*) 'name: ', size(name)
+*> write (*,*) 'thick: ', size(thick)
+*> write (*,*) 'epsinf: ', size(epsinf)
+*> write (*,*) 'osc: ', size(osc), size(osc, 1), size(osc, 2)
+*> write (*,*) 'nos: ', size(nos)
+*> write (*,*) 'wn_array: ', size(wn_array)
+*> write (*,*) 'f_array: ', size(f_array)
+*> endif
+
+! *** machine-dependent constants
+! *** epsmac + 1.0 = epsmac , cosh(argmin) = 1.0 , tanh(argmax) = 1.0
+
+1 pi = 4 * datan(1.0d0)
+1 epsmac = 1.0d0
+1 do while (1.0d0 + epsmac > 1.0d0)
+53 epsmac = epsmac / 2
+53 enddo
+1 argmin = dsqrt(2 * epsmac)
+1 argmax = 0.5d0 * dlog(2 / epsmac)
+
+1 dlimf = 0.0d0
+1 rational = .false.
+
+! *** read target specifications
+
+1 user = layers == 0
+1 if (user) then
+
+*> if (layers == 1) rational = .true.
+*> if (contrl == 'image') then
+! *** image-charge screening factor
+*> if (layers == 1 .and. neps == 2) then
+*> dlimf = dsqrt(epsinf(1) * epsinf(2))
+*> else
+*> dlimf = epsinf(1)
+*> endif
+*> dlimf = (dlimf - 1.0d0) / (dlimf + 1.0d0)
+*> endif
+*> endif
+
+! *** initialize constants
+
+1 lmax = size(thick)
+1 nw = size(wn_array)
+1 if (debug) write (*,*) 'lmax: ', lmax
+1 allocate(eps(lmax))
+1 if (debug) write (*,*) 'eps: ', size(eps)
+1 nout = 1 + nw / 20
+1 ener = 8065 * e0
+1 psia = phia / 180 * pi
+1 psii = theta / 180 * pi
+1 cospsi = dcos(psii)
+1 sinpsi = dsin(psii)
+1 tanpsi = dtan(psii)
+1 prefac = dsqrt(255500 / e0)/(137 * cospsi)
+1 facru = psia / cospsi * dsqrt(0.2624664d0 * e0)
+1 elleps = (1.0d0 - phia / phib) * (1.0d0 + phia / phib)
+1 acoef = sinpsi**2 + elleps * cospsi**2
+1 bcoef = sinpsi * cospsi
+1 if (dlimf > 0.0d0) then
+*> rational = .false.
+*> if (debug) write(*,*) ' = > electron attracted by an image charge = ', dlimf
+! *** dlimf : half the length unit imposed by the image force
+*> dlimf = 1.80d0 * dlimf/(e0 * cospsi**2)
+*> endif
+1 if (debug) write (*,*) 'rational: ', rational
+1 if (rational) then
+
+! *** set up coefficients for the rational approximation to the integral
+
+*> if (debug) write(*,*) ' = > set up a rational approximation to the integral'
+*> call quanc8(fun, 0.0d0, pi / 2, aerr, rerr, alpha, c1, nofu, c2, eps, thick, layers, nper)
+*> alpha = (2 / pi)**2 * alpha
+*> c1 = 2 / pi / dsqrt(1.0d0 - elleps) * sinpsi * alpha**2
+*> if (c1 > 0.99d0) then
+*> if (debug) write(*,*) ' ===> cannot do it'
+*> rational = .false.
+*> else
+*> c2 = 3 * alpha**2 / (1.0d0 - c1)
+*> c1 = c1 * c2
+*> xmin = wmin / (2 * ener * psia)
+*> xmax = wmax / (2 * ener * psia)
+*> if (xmin <= 0.0d0) xmin = 0.0d0
+*> dx = dmax1(0.02d0, (xmax - xmin) / nt)
+*> z1 = 0.0d0
+*> z2 = 0.0d0
+*> do i = 1, nt
+*> x = xmin + i * dx
+*> call queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
+*> table(i) = f
+*> f = f * (1.0d0 + alpha * x)**2
+*> if (dabs(c2 * f - c1) < c2 * rerr) cycle
+*> z = (1.0d0 - f) / (c2 * f - c1)
+*> if (z <= 0.0d0) cycle
+*> z1 = z1 + x * z * (x**2 - z)
+*> z2 = z2 + (x * z)**2
+*> enddo
+*> if (z2 == 0.0d0) then
+*> if (debug) write(*,*) ' ===> cannot do it'
+*> rational = .false.
+*> else
+*> beta = z1 / z2
+*> z = 0.0d0
+*> do i = 1, nt
+*> x = xmin + i * dx
+*> z = z + (table(i) - qrat(x, alpha, beta, c1, c2))**2
+*> enddo
+*> z = dsqrt(z) / nt
+*> if (z > 5.0d-03) then
+*> if (debug) write(*,*) ' ===> cannot do it'
+*> rational = .false.
+*> else
+*> if (debug) write(*, 100) alpha, c1, c2, beta, z
+*> endif ! z > 5.0d-03
+*> endif ! z2 == 0.0d0
+*> endif ! c1 > 0.99d0
+*> endif ! rational
+
+! *** loop over the energy losses
+
+1 if (debug) write(*, 110)
+1 do iw = 1, nw
+326 f0 = f1
+326 f1 = f
+326 f = 0.0d0
+326 wn = wmin + (iw - 1) * dw
+! if (debug) write (*,*) 'wn: ', wn
+326 if (wn >= 0.0d0) then
+326 if (wn /= 0.0d0) then
+326 if (.not. user) call seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
+
+326 x = wn / (2 * ener * psia)
+326 if (rational) then
+*> f = qrat(x, alpha, beta, c1, c2) * dimag(-2 / (1.0d0 + eps(1)))
+326 else
+326 call queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
+326 endif
+326 f = prefac * f / wn
+326 endif ! wn /= 0.0d0
+
+326 wn_array(iw) = wn
+326 f_array(iw) = f
+
+! *** localize a peak using a parabolic interpolation
+
+326 if ((iw >= 3) .and. (f1 - f0 > aerr) .and. (f1 - f > aerr)) then
+2 a = (f1 - f0) + (f1 - f)
+2 if (a > 4 * rerr * f1) then
+2 b = 0.5d0 * (f1 - f0 + 3 * (f1 - f))
+2 t = b / a
+2 wpic = wn - t * dw
+2 fpic = f + 0.5d0 * b * t
+2 widt = dsqrt(8 * fpic / a) * dw
+2 if (debug) write(*, 120) wpic, fpic, widt
+2 endif ! a > 4 * rerr * f1
+2 endif ! iw >= 3 ...
+326 endif ! wn >= 0.0d0
+326 if (mod(iw, nout) == 0) then
+19 if (debug) write(*, 130) 100 * iw / nw, wn, f
+19 endif
+326 enddo
+1 return
+*> 100 format(5x, 'alpha = ', f9.4, 4x, 'c1 = ', f9.4, 4x, 'c2 = ', f9.4, 4x, &
+ 'beta = ', f9.4/5x, 'accuracy = ', e9.2)
+*> 110 format(//' run (%) wn (cm**-1) pcl(wn) (cm) |', &
+ ' peak location amplitude width')
+*> 120 format(40x, f10.2, d12.4, f10.2)
+*> 130 format(2x, f5.1, 3x, f11.3, d14.5)
+end subroutine doeels
+subroutine change_working_dir()
+
+! This routine gets the first argument of the commandline and takes it
+! as the path to change the working directory
+! used intrinsic routines:
+! iarg returns the number of commandline arguments without the program cname.
+! chdir changes the directory and returns 0 on success.
+! trim removes blanks from strings.
+
+ character (len = 256) :: argument
+ integer :: status
+
+1 if (iargc() == 1) then
+*> call getarg(1, argument)
+*> status = chdir(trim(argument))
+*> if (status /= 0) then
+*> write (*,*) '*** change directory failed ***'
+*> write (*,*) 'directory tried: ', trim(argument)
+*> write (*,*) 'error code (see: man chdir): ', status
+*> write (*,*) 'continuing in the start directory'
+*> end if
+*> end if
+
+1 return
+end subroutine change_working_dir
--
GitLab