From d5b7177321b5cae84ef55405d7e63423297f27c0 Mon Sep 17 00:00:00 2001 From: kamischi <karl-michael.schindler@web.de> Date: Fri, 23 Aug 2024 18:04:29 +0200 Subject: [PATCH] test case for WFW version call works, but calculation gives rubbish. The selection of the term has changed and needs to be reflected in the input parameters. --- tests/SetEpsTestCasesFromScratch/WFW/Makefile | 74 ++++ tests/SetEpsTestCasesFromScratch/WFW/README | 13 + .../WFW/calltree.txt | 16 + .../WFW/change_working_dir.f90 | 25 ++ tests/SetEpsTestCasesFromScratch/WFW/doRun.sh | 15 + .../SetEpsTestCasesFromScratch/WFW/doeels.f90 | 261 ++++++++++++++ tests/SetEpsTestCasesFromScratch/WFW/eels | Bin 0 -> 143688 bytes tests/SetEpsTestCasesFromScratch/WFW/eels.f90 | 160 +++++++++ .../WFW/epsLog/seteps001.log | 331 +++++++++++++++++ .../WFW/epsLog/seteps004.log | 331 +++++++++++++++++ .../WFW/epsLog/seteps006.log | 333 ++++++++++++++++++ .../WFW/extend3.f90 | 22 ++ .../SetEpsTestCasesFromScratch/WFW/fint1.f90 | 64 ++++ .../SetEpsTestCasesFromScratch/WFW/fint2.f90 | 66 ++++ .../SetEpsTestCasesFromScratch/WFW/fint3.f90 | 59 ++++ tests/SetEpsTestCasesFromScratch/WFW/fun.f90 | 23 ++ .../WFW/get_commandline_options.f90 | 93 +++++ .../WFW/myEels20-seteps.f90 | 71 ++++ .../SetEpsTestCasesFromScratch/WFW/param.f90 | 5 + .../SetEpsTestCasesFromScratch/WFW/phint.f90 | 74 ++++ tests/SetEpsTestCasesFromScratch/WFW/qrat.f90 | 14 + .../SetEpsTestCasesFromScratch/WFW/quanc8.f90 | 240 +++++++++++++ .../SetEpsTestCasesFromScratch/WFW/queels.f90 | 92 +++++ .../SetEpsTestCasesFromScratch/WFW/surlos.f90 | 174 +++++++++ .../SetEpsTestCasesFromScratch/WFW/usurlo.f90 | 23 ++ .../WFW/version.inc | 2 + 26 files changed, 2581 insertions(+) create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/Makefile create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/README create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/calltree.txt create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/change_working_dir.f90 create mode 100755 tests/SetEpsTestCasesFromScratch/WFW/doRun.sh create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/doeels.f90 create mode 100755 tests/SetEpsTestCasesFromScratch/WFW/eels create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/eels.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps001.log create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps004.log create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps006.log create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/extend3.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/fint1.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/fint2.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/fint3.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/fun.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/get_commandline_options.f90 create mode 100755 tests/SetEpsTestCasesFromScratch/WFW/myEels20-seteps.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/param.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/phint.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/qrat.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/quanc8.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/queels.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/surlos.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/usurlo.f90 create mode 100644 tests/SetEpsTestCasesFromScratch/WFW/version.inc diff --git a/tests/SetEpsTestCasesFromScratch/WFW/Makefile b/tests/SetEpsTestCasesFromScratch/WFW/Makefile new file mode 100644 index 0000000..8c01114 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/Makefile @@ -0,0 +1,74 @@ +.SUFFIXES: .f90 .F90 .o + +# fortran compiler +FC = gfortran + +# fortran compiler options +FFLAGS = -g -gdwarf-2 -fbounds-check -fcheck=all -ffpe-trap=invalid -O0 -Wall +# gfortran version 4.8 does not know -fdiagnostics-color +# safeguard for major version >= 5 +GFORTAN_VERSION_GE_5 := $(shell echo `gcc -dumpversion | cut -f1 -d. ` \>= 5 | bc) +ifeq "$(GFORTAN_VERSION_GE_5)" "1" + FFLAGS += -fdiagnostics-color=auto +endif + +# implicit rules for compiling fortran files +%.o: %.f90 ; $(FC) $(FFLAGS) -c $< +%_mod.mod: %.f90 ; $(FC) $(FFLAGS) -c $< + +all: eels + +eelsmods = quanc8_mod.mod queels_mod.mod seteps_mod.mod param_mod.mod +doeels.o: doeels.f90 $(eelsmods) + +queels_mod.mod: quanc8_mod.mod fint1_mod.mod fint2_mod.mod fint3_mod.mod +quanc8_mod.mod: fun_mod.mod +fint1_mod.mod: surlos_mod.mod param_mod.mod +fint2_mod.mod: surlos_mod.mod param_mod.mod +fint3_mod.mod: surlos_mod.mod param_mod.mod +fun_mod.mod: param_mod.mod + +seteps_mod.mod: myEels20-seteps.f90 + $(FC) $(FFLAGS) -c -o seteps.o myEels20-seteps.f90 + +getoptsubs = getopt.o text.o system.o constants.o date_and_time.o kinds.o dummy_variables.o + +getoptmods = sufr_getopt.mod sufr_text.mod sufr_system.mod sufr_constants.mod sufr_date_and_time.mod sufr_kinds.mod sufr_dummy.mod + +sufr_getopt.mod: sufr_text.mod + $(FC) $(FFLAGS) -ffree-line-length-180 -c -o getopt.o ../getopt-libs/libsufr-0.7.7/src/getopt.f90 + +sufr_text.mod: sufr_system.mod + $(FC) $(FFLAGS) -ffree-line-length-160 -c -o text.o ../getopt-libs/libsufr-0.7.7/src/text.f90 + +sufr_system.mod: sufr_constants.mod sufr_dummy.mod + $(FC) $(FFLAGS) -c -o system.o ../getopt-libs/libsufr-0.7.7/src/system.f90 + +sufr_constants.mod: sufr_date_and_time.mod + $(FC) $(FFLAGS) -ffree-line-length-256 -c -o constants.o ../getopt-libs/libsufr-0.7.7/src/constants.f90 + +sufr_date_and_time.mod: sufr_kinds.mod + $(FC) $(FFLAGS) -ffree-line-length-150 -c -o date_and_time.o ../getopt-libs/libsufr-0.7.7/src/date_and_time.f90 + +sufr_kinds.mod: + $(FC) $(FFLAGS) -c -o kinds.o ../getopt-libs/libsufr-0.7.7/src/kinds.f90 + +sufr_dummy.mod: + $(FC) $(FFLAGS) -c -o dummy_variables.o ../getopt-libs/libsufr-0.7.7/src/dummy_variables.f90 + +eelssubs = doeels.o usurlo.o quanc8.o fun.o queels.o fint1.o fint2.o fint3.o surlos.o seteps.o phint.o qrat.o param.o +eels: eels.f90 change_working_dir.o $(eelssubs) $(eelsmods) + $(FC) $(FFLAGS) -o eels eels.f90 change_working_dir.o $(eelssubs) + +get_commandline_options.o: get_commandline_options.f90 $(getoptmods) $(getoptsubs) + $(FC) $(FFLAGS) -c -o get_commandline_options.o get_commandline_options.f90 + +clean: + rm -f *.o + rm -rf *.dSYM + rm -f *.mod + rm -f *.so + rm -f eelsou + rm -f eels eels.exe + +.PHONY: all clean diff --git a/tests/SetEpsTestCasesFromScratch/WFW/README b/tests/SetEpsTestCasesFromScratch/WFW/README new file mode 100644 index 0000000..d483945 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/README @@ -0,0 +1,13 @@ +Command for compilation: + +> make + +Command for execution: + +> doRun.sh + +Cleanup: + +> make clean + +KMS \ No newline at end of file diff --git a/tests/SetEpsTestCasesFromScratch/WFW/calltree.txt b/tests/SetEpsTestCasesFromScratch/WFW/calltree.txt new file mode 100644 index 0000000..78ec920 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/calltree.txt @@ -0,0 +1,16 @@ +Calltrees of EELS + +EELS + change_working_dir + doeels + quanc8 + fun + queels + quanc8 + fun + fint1, fint2, fint3 + usurlo + surlos + phint + seteps + (extend3) diff --git a/tests/SetEpsTestCasesFromScratch/WFW/change_working_dir.f90 b/tests/SetEpsTestCasesFromScratch/WFW/change_working_dir.f90 new file mode 100644 index 0000000..4685606 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/change_working_dir.f90 @@ -0,0 +1,25 @@ +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 + + 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 + + return +end subroutine change_working_dir diff --git a/tests/SetEpsTestCasesFromScratch/WFW/doRun.sh b/tests/SetEpsTestCasesFromScratch/WFW/doRun.sh new file mode 100755 index 0000000..f5a8e51 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/doRun.sh @@ -0,0 +1,15 @@ +#!/bin/sh -v +cp ../inputFIles/eelsin001 eelsin +./eels +mv seteps.log epsLog/seteps001.log + +cp ../inputFIles/eelsin004 eelsin +./eels +mv seteps.log epsLog/seteps004.log + +cp ../inputFIles/eelsin006 eelsin +./eels +mv seteps.log epsLog/seteps006.log + +rm eelsin +rm EELSOU diff --git a/tests/SetEpsTestCasesFromScratch/WFW/doeels.f90 b/tests/SetEpsTestCasesFromScratch/WFW/doeels.f90 new file mode 100644 index 0000000..972e962 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/doeels.f90 @@ -0,0 +1,261 @@ +module doeels_mod +contains +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, debug, 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. * +! * * +! * e0: impact energy (eV) * +! * theta: incidence angle (°) * +! * phia, phib: angular apertures of the elliptic detector (°) * +! * wmin, wmax, dw: energy-loss interval and step size (cm**-1) * +! * comment: comment lines * +! * comment_size: number of comment lines * +! * layers: number of layers * +! * neps: number of epsilon values * +! * nper: number of periodic layers * +! * name: layer names * +! * name_size: number of layers * +! * thick: layer thickness * +! * epsinf: epsilon at infinite frequency * +! * nos: number of oscillators * +! * osc: oscillator parameters, wTO, Q, lambda * +! * osc_size: number of oscillators * +! * contrl: 'image' for image-charge screening * +! * mode: 'kurosawa' for kurosawa model * +! * wn_array: frequencies * +! * debug: logical * +! * f_array: spectrum * +! * wn_array_size: number of frequencies * +! * * +! ****************************************************************** + + use quanc8_mod + use queels_mod + use seteps_mod + use param_mod + + 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, intent(in), optional :: debug + + integer :: i, iw, neps, nofu, nout, nw, lmax + double precision :: a, aerr, alpha, b, beta, & + c1, c2, dx, ener, facru, f, f0, & + f1, fpic, fun, pi, prefac, psia, psii, qrat, rerr, t, & + width, wpic, x, xmin, xmax, z, z1, z2 + double precision :: table(nt) + double complex, allocatable :: eps(:) + +! **** log modification start + integer :: j +! **** log modification end + + external fun, qrat + + data aerr / 0.0d0 /, rerr / 1.0d-06 /, f / 0.0d0 /, f1 / 0.0d0 / + +! **** log modification start + open(unit = 99, file = 'seteps.log') +! **** log modification end + + 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 + + pi = 4 * datan(1.0d0) + + dlimf = 0.0d0 + rational = .false. + +! *** read target specifications + + user = layers == 0 + 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 + + lmax = size(thick) + nw = size(wn_array) + if (debug) write (*,*) 'lmax: ', lmax + allocate(eps(lmax)) + if (debug) write (*,*) 'eps: ', size(eps) + nout = 1 + nw / 20 + ener = 8065 * e0 + psia = phia / 180 * pi + psii = theta / 180 * pi + cospsi = dcos(psii) + sinpsi = dsin(psii) + tanpsi = dtan(psii) + prefac = dsqrt(255500 / e0)/(137 * cospsi) + facru = psia / cospsi * dsqrt(0.2624664d0 * e0) + elleps = (1.0d0 - phia / phib) * (1.0d0 + phia / phib) + acoef = sinpsi**2 + elleps * cospsi**2 + bcoef = sinpsi * cospsi + 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 + if (debug) write (*,*) 'rational: ', rational + 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 + + if (debug) write(*, 110) + +! **** log modification start + write (99, '(i5, i5)') neps, layers + write (99, '(a)') mode + do i = 1, neps + write (99, '(a, g15.7, i5)') name(i), epsinf(i), nos(i) + do j = 1, nos(i) + write (99, '(3g15.7)') osc(1,j), osc(2,j), osc(3,j) + enddo + enddo + write (99, *) +! **** log modification end + + do iw = 1, nw + f0 = f1 + f1 = f + f = 0.0d0 + wn = wmin + (iw - 1) * dw +! if (debug) write (*,*) 'wn: ', wn + if (wn >= 0.0d0) then + if (wn /= 0.0d0) then + if (.not. user) call seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode) + + x = wn / (2 * ener * psia) + if (rational) then + f = qrat(x, alpha, beta, c1, c2) * dimag(-2 / (1.0d0 + eps(1))) + else + call queels(x, f, aerr, rerr, facru, eps, thick, layers, nper) + endif + f = prefac * f / wn + endif ! wn /= 0.0d0 + + wn_array(iw) = wn + f_array(iw) = f + +! *** localize a peak using a parabolic interpolation + + if ((iw >= 3) .and. (f1 - f0 > aerr) .and. (f1 - f > aerr)) then + a = (f1 - f0) + (f1 - f) + if (a > 4 * rerr * f1) then + b = 0.5d0 * (f1 - f0 + 3 * (f1 - f)) + t = b / a + wpic = wn - t * dw + fpic = f + 0.5d0 * b * t + width = dsqrt(8 * fpic / a) * dw + if (debug) write(*, 120) wpic, fpic, width + endif ! a > 4 * rerr * f1 + endif ! iw >= 3 ... + endif ! wn >= 0.0d0 + if (mod(iw, nout) == 0) then + if (debug) write(*, 130) 100 * iw / nw, wn, f + endif + enddo + +! **** log modification start + close (unit = 99) +! **** log modification end + + 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 +end module doeels_mod diff --git a/tests/SetEpsTestCasesFromScratch/WFW/eels b/tests/SetEpsTestCasesFromScratch/WFW/eels new file mode 100755 index 0000000000000000000000000000000000000000..27712386d9484900bda8e1e215ae9c17c5fa86b8 GIT binary patch literal 143688 zcmeFaeSB2aoj-nNCOifN2^i!dnGh8#Dk4-s#7q(>rM9I7cImflO+rAVJOl_XVvvkN z%ePH`C0DxA6_q;)NQvFbmbNs-WX7_s->vORp<BDPLJ+N1S4C;7QyKDmzt1`M&Y8)~ zFai92_mAZDn%r~m`JB)BoX`0@orinw^M8Bm+yEgw0)JlohU0gnw-8?uex!uB62DL4 zS5fh~(pfXVICJhCK27|Sg%cljxg(z60kEQC{>(4TPt2C(gGpHQeJ_AIKbESfShwWf zb@B4a_<lK6HD+CZg^%}am;09^f^Hg#ii+BG_bjS)A0^`p%<xF&e4xqWy@{`{#s?nC zzlw_L+I4G}ty=2dC$@LWCltPg*DC<t`!p;sCT`vXU(Kr8c!QJa`}SIeFVKh3Lh#<i zr=G_DS5z#yyW+03RV$aMXUX_}w?W~%7Bav;-e<zce*HFn6%{i|=a*KL&-v^;)w}3F z$+7@4`KS9+1mvB>NxXy^j9*2?(wcZH)8qS`?g!#C=_0w;;phC3t|@?3bU9CPjU>^9 zd-)eo7vgk&d{R-dxN2RM%2W4jXZ+s5N4=4Mz-Q9sJa(=X6^m-)($<x};hMfLI0bMu z@o8$^T*}AFl{KrJ2g&VSt?`+>&2fNnO~yyQDSiL+?9bjYbH;4dn|SI~umL3HpOVy1 zPEdEmW8#{G9%sL@oF5fV)c5`zuecXQS(-eH^!AE;6bxSJ6{m*?k&9=a!;gQ@>t6w` zZ{cV5FK!e5ZW`_tE~!S&Wk~z+qkfLRr*`f570VWlU%X+(V%+m@Bz~mh&u3n&{o;>q z+PUSyTYuB>pBuh^GoJbJQ*9y^m6m<c@1T#gsP@S$aV?raJh!b@_s)OTI|V=HDUMN6 zdyBx@)m7{6CN)d%s#&{kZPlt9KYk<VTDAz;`0HAZlz5DCYwoF9yb?HWyeU&T@8VLJ zPfxt<sA1lQ+I34--dJXq8wL8&2F0hE)$7KusHv`6A-P2wtE)*A{s`yCmGlKx)Ko2A zvbKUP!L9sbKIKR9I}s0%Oqett|HoHVt*yTMmPt2MtzNxi$qkFE)~;W+>gI{#CrrBW z<{NK{<4Hr0#G{~uxCROK!?z(vg#Vbnr(E%Gw7y?`y;uCs6nvT6^RI^kJsjxaKo19c zIMBm^9uD+ypoarJ9O&Ud4+nZU(8GZq4)k!KhXXwv=;1&Q2YNWr!+{<S^l+eu13et* z;Xn@udN|O-fgTR@aG-|+JsjxaKo19cIMBm^9uD+ypoarJ9O&Ud4+nZU(8GZq4)k!K zhXeojIMC@E@#|P!?!MSzuXSeLfPH^iIB?(l)tBvyE+4cnBE)`CDJ;CT(^@FR8r;th z;+dyoxW%imueA1*iElY@MItyCZydFTE4)$LH`ppXJ7Ayi1drNFy-`2%?0I=pD|uaz zH#J=1!1PrxvmBTq8fM7UaGr+Q=-U?{o{<{n%;2dJ;(34Y)bKkVUDxgVf&wt_dsJP= z22GW9zA$KNIPAboQFT4*i3TJ+m+87Ln;PD(VKz=wFgrD#GXoUN`vazimlCEg+BHzY z)Mz~0a}~^vT(rx9aoVNv%<He~>Ob{!1I4!_|NJ8?d!mR|=AbRhMXN1Dqp0j}*>gk` z?K^7UA)?_1Pix2KzBk(oMQf!8a20^7_E`T5d^d5}8iyDRm?YTG12%Lj&tf?*aIC>S z)1OVk(UAd1fH*b}ba7N_9Mh9<_<O0gq@rUyaU720n4@uApM)cj0Y?FGG{<q=u5pY^ z!V%1XV<2&y0-rTMZqPXTC*iO&;5b_f9KKu^$0&^>isxp(L^9xbQ{ouv;>guFUQfay zdK>Z)1l<cn)TVqa6j31_?9$)=iSknR5BF+4aCl~`y)e)!as=L##BY?kQlG2(i1{}1 z3;T*S9WU$}!Dk*k3*i~mdq{~8%b#&{zH)@p`Hm6RcYxnM{1w`Scdz$|`Oy2LC;NzK zhfhS=c2DI<7xpp00{^==0IoC{E;!1C+YUHSj`(%c=Al-H&(k7K-+v76AU5z?Dnx_q ztA4?w`(@prS>!E`IJDWLT|5uye6;s+Tem;b*_mJH_x&f{?i<Q0acu*f*-rt@V@GZy zv=3U>wgR>pWtx28Ids%x_Xmx>mPXcv=MN^~vd#YTLq^bF9XW~pW#X^V_<zbiCH_dR z#D6cIn|k|I4cA0lNVpCKw~BCvZf?|YiyS!r0Eurg;SBv-n+#W};J!#WL&xscaAgjC zwu1W{p0oXBxCUvvv_0$4?Z>mzA~Ha>-=CrV<2C+8-r!M|@l+0!^iIKZQ%^tPz<~}4 z7g2CG63)=Y+cex5v|qw?D7dR6TwLF7)^IO7`23g2dPm?n+cf}Jly}=D`Y1h2B_rno zpaHVWa$d-3KaJx<JU4OpGvIg&ur_db{2mw2`+NqRf5J6<s5JILab>~Dt-D5q52|zr z($Ld2qHo~7jy_6<BhU*k<hM%mQv@!v{=o9}Wp4TQ1En$IJ29>vuSy>|8tMP~dj(?Y z<EM)^6?7a-kF&_ZX<ueN3!G-#f_)SDWBRN8fDIkET(;{$(8jW%d4t^h!@S?D^2^X~ zKf*oRy@%JD88ZIp$lo=>4)<!c!?aIgF6~udoM%4s?ZTW^`-Sq>!W=Q7$OpTIXVM1X zxu3fHg0I!2wX$C#t)@OUpib($x3`ubF-XeAw=_*()iilF7S5A(f{u=76M6QsPh1xQ zJ*0znl=ihEO+SU88=wukGRRsB+$No(zoJv-R}PlC6d9D>miz<Cn|xl)dxOs_NITm0 zcesY11uq$=Wr(lrLhx~;-Ct}BYzM#J^}~*P#Oa5M21Fh~yWU*^S!3LX*QAGfaHWmq z^^{Mvy!@Ky`(q}adaqu1-&_W`F*={JSt>laO~6f_ybRv)S_zmu=pv4#7z+Nibojh2 z;?JAl2fPj9gZg7<@XXFDi;@o1y#_p)uXt5p3ARh;kw*F4_M_jVPipi83h5)D%YN&F z66qh0)_cXfl|zD$+#$;%*3{)s0x#;KFM}=it8#%N>lIzj9tT*-m%w&=yo^bY+GD*@ zA%YLt`QE6ERf_>5d7D0-ycphk#o33t%UAYI9{j1jy)-2cpF;bbJeP+ZMJai>H;IRT zlVuV6YTW5~xZUC5lYn)2_>_#1kJ@`ZiibOJrSEEHES-5A@&Am66{+R$PVzzLxjYnO zQu6TclX$pKmWABt^67Y3=J0U6A%}HZ4r{d>R^v(@)@V6An|^$BaXkFs!sHNoZNHOR z4hNDCI?v^yKQ#|~CGqeJvMl6AmruvTcW+lV!0u!4Fjw=?r+IkB0~-Jyc6y>ehOA0I zQk&EIqgwbI{0e=r>w9(F{4U~VJWGn38$Bv+ZWMtr*6&d#VH}_E%l?l2QsY27!PXCv zwm!%{OcUq72z;c2wmCG8^5}W<I@ms@jd&b3Y>cvD9mQ5FX&A3y8Iv=e09$#L3;R95 zHo+!_^2?$)0jEL6<<cG_E{8oX)p0rf^)F0Xo@BpjT98h_AH3Sd|1fEJO2JB6kWPTL zuW?};2n*ZNw4)R@0{#GYtjHM`U`+co>xW<aXIzz?S3K>|JauBWN@EO7oG$;?8~j5$ zfsV*nO~<v?GSX3J@DJ$(*vgN&uy+y`y4+M<mXhXeX8nDP`jJk+8NANLIh#0hjrx&J zfVIcDu(uNydeC%6^9+6G@Qg9oU92BG`!ud5&qS^fbJ$7_?Ibxg<E0XfC%~&CkJ+>* za=kSSWwL4yVEYr+XpjHi${!;ZSo7U0l<qlt_31&0dKEC+)1lf^iP+B6E3-Yhtmkk% zRxLE!m5rVSZm>>F2aas?46wfeY_@uK2Cy!Vv(Yo)f10#pqi2BqSHfn{vscafv(Xvg ze3Uq|(HX#6gw3EcKV$un**%0aZ0X9wiG5ya_IU(zUl@an3}eu560ak>*|aNigLM~a z$fjL@T}apr?b^or(XJX?U+*P`y0m1YM}Zrya^la%E&z54VGTW+3O$m(?YUm92cN@y z5q^c@N;zhsepUS#`h{nvZ>yB!^br=X0yv9Fe^~kR&|r+eWZag1zLo9EgWs~zKE&bD zZ!z70c_hZ(5<l#o#2>lQ$|ZiZkFfL=Oees??m4i6uxMWhcCX3{Uv+Plv=cZFY<}`{ z`cB5h0)C$33dS4m_y+S<9N%!<!to8qErjLuB*r2|@PWnfK_@X6BoL2dOd`iQ<fX<r z7rrXj;q@Y~x1s3XML!1Inv;6G2%6^)Mcp_fY$)+ef?o@PkAQg>^sEUYU$)C%U_E<V zoX4b7J`aqSu_e<S<H%>=R~WFM<LEAoae^4*43YhVxJdUA<Dx%1e8w0@ma80ZJ*mst zFL=-gu#qq7*!-Bjo>H-u$y@T!eoe!lkk^Oo!#ak4-V<Hqg*;QHa?tKO&<C1^XLQ>< z;wDY^O_qiFLob#_?H5Whk1w7RIl>ZW?XTO1IPO>S)umiug0-3DUPHN4N#zc5jGKAR z-aZ7`#u#@f+ChGiw_xW(V!&@6K_8?Y<91JmCs<zujoHcY1ZxgqwT@D5ypS9F;LO$l zj=Q5=NF7L%pGu~2Qr<&l7_;y4df!|LzNeCDj>)%d5<6e^<w$w1NGZ>r+`suIwS>U4 z641f1JLS6>x>oJ=H5A`l{o@emtO73PQ_g{BjUV$<C;O;563E*8VXn;Qw0xX0zcjVX z_a`5aSEzMU%4y){6f!@-@&cAJA4VVAkomupWg)M+JZ1hbU1a`d>+NZJEKm>FLtkmN zH)~n1)7J((H|&WA_Jn-_`OKBFzs(c<!B$W7C(E_WulGcsfUNF448A~)_duSTAeXx# zpWgv3`<6paQBUaLS6b^uuQ%TWk6Hj%FZ#dv#@<~c>OEr3j4K28y#anzAYDm&hkWYA z$z04O!$x9$McGI(QRz-#lJ!%x!^H<!i4XHm7mDwDCO*_3&O{&T3jxNV@6{CaJ!0Ym zyu_D|zAW&475MfD-*bClBleu#^Ysdmqj(8kS7f3U`Lq=dt(a?4yo^k8d07M8E<V6Y zeCcS-0^dRtAK-Uoq7U_j0OQd2N(%aBnfL%N@uj0L3w#rSZ;S9gw*|c1a(4IEsc&xr zUrRmO-b=juMhe_P$|dqco8z=n|2XGHJ@trNTz-xOeitX;A=mG|0l80>&F<jrZQ=wS zg&BBC8k+#~GW5q8FD29X-mQr=N_;8WO?}J;-*MpU58o9AU4Hr$q+iwRWU}tc$u576 z0k#A@C!hB;KsR#4yRRZm{og%2aNn!21!34x$sgp0wx!{DaEgnw-NXralz;bC@cd$M z?lW<MjuP0fSNC?sNgmm#vk5TdIc*kVCM6dJjsIlg1H8nSqTSizdjR;ZJ{$V_F3?4u z`;mUd<ax(bmp}gq*u_|LoQWU&;n51rS$qRW!hTBhYG^)D8ZE!yxAQIGYiQnz^hV## z>pk#qfUAdJ901&37C8-1dipnfSqx}6aqsf8g&VFrTMhmB0PBrqxc)b;PRvpSc=ow4 z<25l&(PB@t=7AQ%HsLDwD<EBn{tCip@%f|9^HT$@P;ss`O+=>OmJ^viy>Ddt?E@my zO9n@#7hM~OjugS^qjG}NNB0d*FBlM<UN|^7-E%GW5peIqvxsxJE{b`7NsC^u3k(d( z{RgS$NXk$zbb79{mYQxYIfS_^xt2U2zLs2yXUS{HO78Wz!GD`I4ERIO!so)?m4L6D zzhXK;UhLZxp8($v@HXrjVPy=2bOLN-x(j=5Dqyir-2}f>h&2c~H%1?oAP={@@b3Y> z34VvLay|-a2NwO%=n=w&eT}f_m*&Ip>nBR1MH_Q=BIar6vSy9>Cbm)XMz2{*9$?K{ z^}%h{ZzP?%Z}gh=6DiiLUnU&JO&&dNa`~mlO>zCw^F>%kY-1UWm0W)5agz({%ohRn z`x-xEnkG5!)I3Al;e({znU6x5Jt)%zoG#ChPJlJ=KSJ1dRIFG8y(&c<H3n;Fj5~`q zdh&`On@t!eKo6In%-=P_fyHUo4FS&eVEn*!G_+++4s4>c-f8D(-O}<L)bflzJ0b#D zlf~HI-z)um&o;IjWA8>>RiEhjSQC#Jrq<pv%!A#p@oeH%k9QupF>!9{8zU2G^<M_v z!rmWAYZh};i-Fsmo4S{EIrA3dj5$Q46Z*lK1GYWZ6yiuTHzoVmnF9vw4S>ydZtBw} zEl*O9y3PRu{}B)7fF+-Se;8?bO2JB<L)xJGhzE1P4s3tIW|*70-K;;&9FD9X^SKVr zzf1wnG;=r-7W26d>>Gs5Fo#os`lYRW9{X;dpP6>wjkJ}>FHE+Tc$PHBa)j##PW``N z{h4hg_@!+nc-s#647T!!$FY@wJwVv!Gup}{o(^p*0lydU>1-wGb!{bJA0=#-w(?h! zj`;o_*$0lT1g@{?eCfwDpHtb&#{k#YBZjiPvoEDGu50^67IT8KE=RYrkte^`szSY4 z$rE7zp0F9@>Fdb<$m9w2-@^K{l_$Vo2l%Yy39yBP%~qZOp9lD?<O#5S37e%nO_g+{ zktg8#%Vg52eQySN>IXdQFh^&~Q^b^~fbtK1?h!NA>mrPu^gc4y^^x_b?X`Xi{8{%W z;D5xnr;nu!`{1>nlzht8K6tHNfX}KwLC3cUn`M9gOwy62KY{B3oqs<2@K1!Jj~4lI zAE9Aq5Hs9p)#!YQN6(Aj`Zu+YkpA|BSEML6#;RgD+6>0odad}L4^%l*b}9!czwgA{ z5kpRv69?jv6}XyZg1St2J@z^s!}?UY+@pd$+DU80IS+Q(6Z$#e3YDw<X<V<pQ`hks zUf1jWRUB&ortlTdoRhx09)}H&o-FZR<{2zDm~;fV#~3k7IBB2p1mJ3~E{%mNN@FLc z)fb#AuPg{}E{mNQSzplMySG5{<bX3Dn|^<^eE|A+fqRZJ?u!B53dD}d`;FXvF(r&c z;(Lo!n+vQz-2xuzJdP>x{Isxr&K%?sci!FavrYom&oQFjTQrIM<(_a|o^z0TeFOTX z=iNJe)~|FqdjVh%U_7)?&ppn=oFexXmFv01Imnl5D-Qqkp7?q3ee?hC_{aFOEC2pk z@=xWt{0ru!<e!zqzgJ{g?9nJu<<s%+ewTlMmHg9tIp*oT9Obw&cAw*oh6h0g^<J-M z&>!qikDa69tY9x|HS0{y7rov5yOGcM*+z`*H`gfEds}V9r{sZ)Hxe)RMdF3OPsq3X zxOnF>ALFBh{9s=<|5MC|olMBD?C0i}b<IEGuC1V5w<15~8c8ka*7E7n9q6y{V*m2> z#H;!5=G(b${#7zRiLT%PH-8lCU052m7nCOQs|NLg&-=i$7VzaswEeqi-@kyC@1efE z*bDoFZ>Y5gdwX|dzkW0JBsXFIGWR-v2k$e?L4Uq8PekPv?@Y)m-kBJL{kQMT649p) z^t0Z9jeNUAM87w2DCha3Z%r1_f5AS}xA5M@zrfbLH6GW){j9gfiRgci711|~Mf8nA z5&hFh5#4vJpY^AF%w6SSzdiO)b_~L`cCgif`79fISvl_zM8Gcg50kRo;p^M7!jojP zFITow!bK1xfKJ%@F7I|=pN|_Sy!`P*9JYz0L&X_^zKL;$S^h<oH}{Qyz<cnFKE;cA z{>*#u(YZg(`{wxlf8(Bczr+=~;l|Gqla4g+A((Wu5-<8t)`LFzH;wm4xc&gP7<O{5 z?TnQ;9uh;{HRqE9FmBE?9>N~+uO`tah8VmS`D)ESX}oUIXOSMns3FLuJ61*77=NDj zVDC71a9Usw@RtEkIVP+et0J8M8yw`qRuC3)+*Do)``vZF(ss9Q+5<gKUXV0qw+FlU z?vQmWSXFltY~*qmHX!T9zN7N8ly$FT-RR%>l7=*~6!?3i&KHbjbU!_;=cT}t{m!^5 zp!v>tz~t`@D9^QW><2yyA9i#a?qhZB`(i#}*@dDt1U(An2Cxq_(0V2Epd*bQf4<d$ zG9z(M`M3ht1@Cq4K>4|`x}Sb0;;Z;=XZb+9vnCWC(y>}3yrV{$^<XZASJ9DYot=dC zYPcdvk7#`z&p4mt!CoP|yr1Ru2aj@1pEL+N-+Br8O?V#i^|S0#eXUd2i%gh-pvNx8 z{Ci<(^goDqju(5+mF?rYD)S}~W=<L88T;}f$C3^kdH_4(_d-U;iB=i+jTEgB&_iC> z)xvrU_?mHj8Dj|2;fa73*c0$q)`RDy*~Asdfa~`f*BG^5AN!w&OI)Y$T(;|QU(0Te zx9bHw4~>Mp=EF{-uK&WljWm6YY!~uA<TK>`99NT$$_#Y;7_e*;-<^?dvK1ZA;CW~w zc#FIT!2{xC`l%$`!3?-V8uu8re;%@PrEFu1#`P$!@_irO_h?Hi?_+g=eUKqb$x=RK z$%pqJ1>S$m!TXQB@cv_Oy#Lr|--qxW@FOvM2jq0Sh&F8#(dGsb-Mv*r_iPr?C+bCX zZykIA@drg$8N{5S<P&W6Lec~{58x=*e$ZkcqkabD{!RN3;+8{Y(HyK1bzuG7hOY|4 zZ|%;XwU75FeZ5+|XYW7uF|3u(7X4cK;>vY;pU3;;%g7t!@jgOXIY~J^d*87eUb*ks zXT-pk0l3<)8TA}QJqL|?;0MF-iM#VZVbrtf*f{JJyhG%+^v9L;l<RsHCDgO=Sn+>v zJk|_f+dL}}Efc<$CZq+Pacq)M?*Y_%0Cb>ze}kT3o#JSq^S)zMuibYnv>yCj9{?{r zEjH2tJhN+)%0GqjPpR@Bp!}66FHpXS<!eyBMwdrgu>5kA|Ildnliq0PN$eToyJgt7 zBj@Qs%U^*jAGjjt?mKq!kHEDQxRz>MNJsDta+%ncy`W{Uq6K?sTDX@Xc&_eP=necT zP`*N!M>>dSoDZ<`%Hw_AjPlK@{ArZGoNYn<MXY}T$}iC6kru3fUQ+p8D8EaUe;efo zk^UIye--8Dp!^(N9_a|4$@fT0P^Nqc-d!!(Ctwq7<9*v8ZP0=tQ|-d?XyH}j-S1+I z6M)~8aSom_&XIB~ZHP*HDz8wo4I8z4B4oQbMOi#^%T^ATF%t2ZW&a6fCGOyG%YFxa z&-Zdpc{m^5C?Z!{b}rtrr4R4piH3bX^)Bo&Piy4>z`fvUwafch53~MU#b4Uo$R?p| zuErVBIHk=+S!r{1*+!IY29B56hqSq3gtQ_5Oq_V$kVIcc23+@RTw_45#D(<;>f~36 z3%X5P`ZCrpLP(b&?d*-Y3Tqk8-k2cQzmQMf!%iF}|K!!D`5qar_zv8k0A2?tuc`1= z0v~+dQSK#Ry+*u(@n;;@L<1M$7?YgjIxP5U(w8m1Lg4ct<_kTKda?hi<m{Av^8Q&7 zr4J6x>pk@uUKjM9O8P=CK*nEGZMBcVCWRsU9QPhZ`4<D&4;DxqJLKoIhUVdYDD)Bi z$1UVB(yZ5%C;9#q<SA9(Uk+II{ZE~~$6P!6YXot@MmoF=AzhM&zp}rOpSHhjIe(8h zWTf@`o4eUxz?;xt2CiofT%<3n{>p&w0PvMUZq5FZ_m`@_^3Y%9@&5V+>SKSgUbDaa z+^?UezxDt&L05ySzrIIY@GWkCA)T?mke{}{Dn}{(!kWV>;B|2FdS3kk-h}=#a4j@& zg(R(6=vM}Ovw+W&)L-)6=&!1a(XT=1uZi*gx|2Lcf3aS(zXBQhYbs#bUl?=h*iZFW z8F9gvJG>1cow2`=pSHiSFILf0XnBFx!I@Qm0dGQo8Mxk=5XVLOveK^%_+A9QkF>vX z(O<>!{(2AfvA<Za*<Y0z`s-JKWq&m}_AsLQ>$k*(*va87<EV`N#W*Xq|Jv_Sd+d$~ z;kO<JUI%AZ{RO;6e<k4BY~Uh&S@l;2d^NzAuD?pP|4P;`+#e0WhZXS(A11G-;MZQG z??w!BLishMrC%F^XCC0n$GTGSsK9pm9u=GV<oiv{hkIL#@2>nHWnMGp6_FPdYL7Je z_bH4&;WLkxZ4^8G7_)muiM$ewac_uq&d0m7*V%m_BYANdS%x;Uzv(}Lr_s})`nwv> zFHwIt>i*u({vKiU(FOE(0<Jy=t_$q%1bk<1`sn*R^cv!rQ%P}*4}F69eCjIuA&0z2 z`Yq4{y~Vgv#l^u4dV2z}j^3ioX^$8w<Cxd+oZ}Qn$3jS#r14`ijzNA}y$#B>iv1!| zWPJyC9h_P9Bk(5Zt%2*G3|yoyE4|Ht?*ZU@fcj^|o$}uBi_pspk2^bkt)VkGw?OyT zx5;CqS+Ch&c830{1uXlk41Q7emtXbQM&iP_(BW+e>5Tn_{Ivar_e~`|cCmFE@H#lN z>M!6;=r04;#|&JgFRT8_fNvP^rL(iLpD<3Sx|sg@5dC#3-e1>~$KXBdHTx@)p}+b7 zmi@(Xo$N21Z$n!>n7A+wcX%5@I%9t!KW%^cuTuSmy-`QTceB5MH=(}_Tx|v}(w9|# zWx%%&_&(D9dLR9DEZ$!c)W`l}z2tr2;u*2<PH)WS{6GNvOt3~R(?OZW86axT2hZ1V z524KS`z#OV+gJu^&bcw4bF(fy{A7spw}b&M${x#9qE4SZQpQA~&N4X{F;DSyjfZ<) zo%&E-*7u^WpJyte{u0!aW9(I|#9l=w_Wd-+`iyP+B#iMJJ!;PxeC9RQXV`wmpV$v> z!UtgE^gjCp_)@~ZQx+`{ff%oaIf2-T-1>sf!&L>cZHi~HA&5D$rXzTbHBRF#B;Gjh z-8Abf9K5Q4>;5sU`#Ui1bDT31#y$Y&IrtyDL}Rx79yQJazs6b~)i-jU1$A_H%$9)b z?Hl8Ic!6WK1bn{-KB+q&A`YxC*eBCx^nA#HKFZD=?4#yGd|2nloR=T{Y<B|wO+6oi zIkHe6%+dRj)>4;ZPl&X~dX5Qc$X?3!OJ9%tw6@5x_rbB&zZg8oil6ChJ#al@;JRqv z0eoNmDET=adjyO*MS(MgDvfwd?)^$LPf-WD=r4YXzC)SQo`~`nTkt%zo%?oFtoJqC zGrd+`&oei>^uo{K@rrZl&*%W&3fQQWdtx|`UqYXeu%`nyw7|NGH0wM+?3m)sweX8_ zo)YU#@?BKu<h9l)z*3H#Jux-1Easzi`IE4_dQS{uavf*d*IIpbIeP(MW&81->B1EI zIyj$=eH~_hMbw^fH+~k}JC&xt2BDsW__;&%*KkRvYX38!k!|z@<bGhZ>pJU?H=sWk zmSJB*(ms$g@L`AJK8*O;&mREZ9zV17<S)N5_nBXQcdch*(axf?mB)%U`troa!aTna zcL=Vvj)Uy&2!QW_nD8Gt+L5o$=IYGdR1n?{-?vTY@!YNfI&b8rg7%%*YXw^$-UgWM z=YeSeOhZ{L2AT>3+Ftvm50%gRfoAu4B+c^<{oJ8-z|W6%6e)VR=gajuz)L?M#J?T& zgWl0NXY6?H{RN$%IlL8gNWV;6fmr*tve<F@35^pxKRx>H8L>{n!-vWCO;bFOKCVHv zll2vt{s%lV^NCaXC;B1qNY+`Y>1X}2je*!XCl54}o|tnlVFNLi<(X-+e8j`}L(z}& zvOc3Nq*cOo7;piF_XC{c_G!q}1o~l?2})WLY>q$4=5%o1T^bp_6a6S<81Fe$(&mU! zGDclWnS<<M4sI^uo5mv4zfM1!&#OVN+g42@@f55e?>EW*VEfvUX1@iIr{k}fsdt)Y zv`39wE4gPR4b2lY&3YVk#1m08->T{3c=!R-Qv}|W>Stqdm38Bp6USA~MI2?=>LA|X z0_`pha$iUq8U||`9>jSdD09RkMoYe3sq4$dwQz8y(k;kP_<hf%id!0W+yXkUvyNRK zk2SLLvx#%V3AiE#t_zG?67YQ&_%1YVvEPJ!G-8)#?7eF8;XID=JS)^Z$l86PEc!|w z)}IkazcL6{$m%PYJKTNha*OZn?LIZc;`^_!VZ8N6tZn1HV(!s6jq?~RUVS$AOL$sj z-*t?({OtG0&DO9!c#jVDumF8W+d#bzFW2^Qxi`kNd<K}&fN28^@rg>MpQB()c1=?@ z%*~6aJZS@nuNru&%VMm%VD*d`&rYOHkD^XzkVl*QWzxvw3h<hGTn3*e+awC4Jl=!n z)Z?YNN;)C`kdxNm;2j62Ero-XeQXC!qd|ZBG^MLJk2%I?;R@Y;&~0y3_A$I2J^=PC z47)4aMHzn<`I|H3Z$`VHRcA@gMtiw`7G*HUDUOU@HzU@LJ-^K7Sl}GS0)>r>WFK@c zzrP?{kFtlc2Q%ztY&+N9g?D9k&1|*tD;yesckm?MlM3M(T&{Z#Z2<kF--l@wVw}~R zJONCBSOmQA8|)L-%l71<4RW1D$HPdMWQvE8pC%qQbYYzJr*W4`PsdrmHE><Lo{qDA z3VfGJPn#fXA?Toe_;Twzn+IFX(9=`U!@bZ`rH7st=%v&}oLwSy(F<LqY<eOBGTXQ1 z_~q=MuW+p0<(y#1N~1?y&p8jolQ&q~$%~=toZt{*MPojr(PQ7>!hH>Jh>?CR{Q~MK zVO+uUBTfyrO78ZbZHHcm^OP(&{+@bHx&+|q)9IJv_afE~y@=ljlpH}{r2HO+ZE*61 zDo;H*Qe0o~Z-_C?a%TLI31^r(HVAq!Na@hO!(WkR`rme3`J6Pi<0^F+KDYg_mg~dN zWdo*@aABPn_QrTl4fi-vC*|HJ>Xm-)G{C-5?@~$qxdptW{>Ye9>tTW3-xmhX*Ab&U zz_U~{*s|cT^tNm&@Tan6WvmzSH3?VAKW)nppA2QavaRINQt~LXov7-29y@`vaL!(L z9y^hZjN-lRH-P87GRpV60Y@2i<Z}Gom1pV4n&J1FisAR*`)ugca^gem-8ohJC-(a| zB^%Nnrs_Ae=l|&TnEmz>r{93%9O7ih<~nj#Ngs7@_CBiPiuC<@81<*>*I$w@^lLk= zX1|U@znU`Cp~hy2d-R^4d{q}|6qN<oCy6~j12Lw{=erMhmnCV>PBZXJ9-=PBwGrNf zw@#e<H1EN0=l;jM2j89hAMn0ft+}xLe%v#!kyk_3E0f~Lz;J_3zhBYWKRulrG~WAg zMcdr=+9sZ0k;Ee~w`{i8+~2fR<6%4;!u(>Lc@KNMUf)+H)e|x6u~j{loKG_NYS!}^ zjfe4g$X6cyg1-M0uc&tz=1u*{K3fcP_iTr*KA$d*x2{a0QCw;86l+WD`$z^}Tt!;a z=xadw9=Yz(C;~TG=RSshKFIOA-mi_c<5wgs-gB1y0@yQv^+0En&$tYoD#JQ1&Yy>T zkS}@XtzW)dEbBkc`Z;&&)Q@yx{c<ir!X5)G^jX;q%-cG3Qa<~huTFJNn#A)A>(uLL za-X=dRwmo5&PkK7KXU5C9?#I3c%789myyRdKgVk(emiwXZnB!PtsAgA084(`sP_t- z10dyUN3Yf&>_9AtU!l0dm8+u5Ay;^oHWmcF&A=!7(y=L~T!}H~)v<4~ojX#;DvEFN zy=w4nf^|3BsrML29vk}zB&>R`TEc!ww=-eiDdUy*<NRhHZ|g?=pX^5c(*etU?M2x0 z>W3VQzFNkT>aXcW{nuu$ALIDx7WqJ4jZUuLk@2MZ=dgar+jvO_`+ul>zY=(kFlLeX z(9ivJK6K@BTg%Hj?fa5^rcrLZ_2IRY<HT6!$T{_1ro^q~#@2G<>=!icV=ql(Z(M(X zc|0kzJ28IQ37N&ue1D~=8ZwJ#Niutc;}_MI3D$2xQ;<HdL&ueM-k7}CacVqI7|*fZ zrN?fp6MXpmXpDdCdNmKm*pxAZEF192`Kfq1AkY8jK9PULnSADR#(EvoWITJsStAq9 z8T(@p14^5s$5KdV97`cT)mSP(`;4`JUIor@dl>EK*noRtc|Gs`x@22s;M#29a_rxE z_pfHaR|9+xz*mNza`vxI+9&V#>h+p=y|_+;eFXA)G;|RCa!7DYd(=LFH0+waORw29 z<6hbe%|E0w@(=lG_y^sZ<nr$mMw_$ZU!t8#z%|jp<=D3K@-G9vLf}g^7A+;uQml1V zk=I!70#A>Iu-+xFgRlmMc%YIz#XfKwHiL6XxPORy#V^GkT`7~c)=9|ZyJNd4lLwQ0 zx`FGsfy;?Q&dbvb_<jL=rHsv<cVu$PK6!8KZFxq@WN*r5mp#=sVi`NEW%4lir`KPe zfjy9Xz`U8Gn@^#xX5=@W>V@%2AG`yVW6eO_oSTzZ&dni4<lLN^6Yybf&fBsD_G}C6 zPBUT{&dbsFFm}AkuWcRjIFFoi-U)LyLE1x%2c5ZHyB}=N0>n^Q|Nacu%aJzM%T@nP z)HZdJ^#JI$u_i-(kn810cYzOl+=agf@TtaTIe0#j>muMms2si;=Xs66dYD&u8@>ek zi;(6!BwSzTxW8@Orh<Z3#Q08JIp-C|@9g2qf|u}XEWW>>W5KQw!ua-50Bwn39$l_$ z=(24X=X1W2V}8zA_QkmAC0w~yLYv+`?*4+&sKd$cEUzl4hrO3PlRiVl7ICbFxh%F{ zzH5MVMq7jYG`8k6=clDSRNi8Z13w&`yq?!q18-ui8^<-=z;z+>Ub4>(+SpgAYz=fE zrL8g7)syF^yD{%|M9+I2(eqxgr)u8IKTyj$?xn16Rr67@pCJEo4hsEh_M5#kL+pun zP-Q=K++zLvHHba2#_o<OkWSz;-a&O>+X3rAyk*vF_rZHQ_ZFSaSRc+MlsIrs-@mZF zT}FLKC)6jVy0A@VeL9A*2RZd2KEj?D+CNE8hvrQkcr#7&!d7t*=hl+X`GOs}&MYU> zx{hqhLuS_mFVpdjQ8(xKV$juexS1b@?uFND+awg-g@YGCFE!o3=gcd%uP=*{ey#(2 z0lH}74V!I6nbEK#W;?=~7S7{4ah#k#AG|`Xd)8t;4EWn?F;6rSG}UODAXAhx#(=~b zt|?P<`J_D`^M!yf1P;>6JeTf>rkm}X47kaJTdd$v59{Duk`r5=VEK`lS2yNJu--3y zhclNW3f%Ghe&JW+t^m)(G3Fd1&`Vpte|3|T=lERAXXL%wBf)XeY~CM^=g-7F#|c;S ziugA)Zf2`@z%wIwE#2rD37dQe`}%=iiSw2M_W{P8PMmMjX)8K|8R+ayI#6E_@0}HG z6eGEg?;j8uVZVm8{_XL#{-K&a(*6K^SrPhXIG%C7jL+EL%1-!NhG^LLu2yX^=4p+w z%LMC}I$v<iqsJSsob5Vx2^1>chm+oqdk*F0TrhkYZ8vC`Prud&_47NNQU-&R1tH&$ zeP<Tdz3a;NotZZS@59c-`>>zD`>?a{K5QV`S)Q-noAvdJzk`i+JHE>;@d%IWyCY6q z;lEYc2%LBM&!EMDqkZ;pUmo58r~Z6UjBlI4KG=ie_sMpEearRTfb7d@mOT>jAm0(5 zCe%A#c(+>0*LdJ<7S^{>Uo);vct^OB`#10n6z<)`Tov}(3Tr!g1OFW^^u)r2rAi;s z2J)f<c0s-?jCa8KuJC{_KX|mW`ugW$b>4lk&3s?j2ma)MKfU%vCnBDj1wNGHJ>e3( z4?7L<(PX?QjOR~`M~pO1^Jj1h{zRrJ{^0Dk&pZ49oXekGv!i)K13dfw-BR#$BI98^ zLyQ!)*Z0FXV|LUt3TZwY9>5qxy%S*9BCmFKv=V1}zkxU`7v%!64aaIuZa6kb3~Ko- zu6Fb6=uqI^rEyOpZakB?Yfv6(iJQ+%+{=--e0DSl+^++7Ka>-A2lWQLZ#t0g9pmb! zo<BdNHR8iHH&6A;b3<Ci8qDcFI|T3l&%?I{=C&4L-wUsW_lh+u7~?UG@3=ijIYpXs zKL)r0{>TWfZwREhzQHT>Jy)zAJumd^T#Uzo$CuZNer_ym#6HoAeq5VvEc_PI3w8S2 zNH5UoN0FYV(>Bs``o-*>v#BFddm?2J`|h@*o#S<S2h!tY8t<Lz{Us$zRwJ<aVxeb# zcu*ebKzsY-Nf^p?xHczdYq;=%oLIOXc+fAvuVOT$OOQ?#Q^~!0$WIfa9r5_(TnpyS zZ?}GWRX1Za;5B@C0<NbFTp>wo7VoQNz_%0lJXq5XE#%(E{`+Klp0oGSJM{tZEOgMT z#x2$0|56b>P=V|E*-^f~+~Lc6j%TOX<z8zeaPhru`;I{L1ilBsb3pN~%a9=STlfX- z&s6Z!-FIuo-Q24dpbf@ZhP;w*>w(WHi$0l3xpw;o>5P4Y{Iq?eVpLhjLW7^govN=Z zK^t^ly{ikn34LSW`lNv?Bx%i}Z!+MU1bi}nu4CU^ruL=QI(_5Q`%-<XZ#JV3wu<Nv z>T#_>-z;a}41Lb7Vc+<yVZenm_FL@bfoM3nZ?5Wo-wXym**7RF`$o#*rPZ517hgKP zIcea!6ngV3;FG$tlYKKp>CJYhZ*sKW<fy)BMjz~j-s}QwgVvdC@JV|Eu;X#RAaL%1 zUEygB)%L<VWU=+ctU&ZtjMaBPIXl|CcQ(fD!!T~Yq9r0qers=oP1Jf8p*@EFG4FBE z)q(3@adrAy<@=C7tmR>uKdj|#=EL^M_!Ksi_K-2PifiwZ<J|Ba{bKD0V1H|8$n{T^ zr}piHx97y#x2MV5tn*Oz1Y+2Pa?QxwoTglz&I8^q<r*~3<`j85k%#BuRP{ZX6ARGp z=jFAw=FMyU0O=p*ZEuAR$h{Bn3G!Wjet#lB`%_t>$4$bT4<3cf)ZPc4&A|Bz4R_ju zSOoBhTZQ#0^77Kkzr;HUcxjWeo`+3s=OwNp1J^~$zr;5X_%4;-czcZa*8AFTbSB%T zYrAEeByE2Jz94-L%1YZLZLzjZNSCCs^U|jy|2%#}*6}OU(XGA#-V}aA;`)h!>mvPz z#P<a7T`Iq^6}Y5Nhu=7n+&9;Ezi%D{KG`=YEBodn(3`sqezK3w>(evn%{&9wrO=xa z;QL7Y#zLIA0=vfcUhl*!uveG|P<9PxmO1vSfNcsWyGB|P<PztVIdDS(7wU~~zc&|K zMOYi1g?ZIsTpO<LrG1kyeA5;1P5JmXa3ATP<R10+u%6l{U=;=X$UR{6SJdx_>#y=# z-YcTN%5V7;>Y~4*jH!J3t1DXGkoi}%lrbOnyAa<$V?Joy@2)@kz1E)q$9$}#UhdfX z4&U(WhW8Pl2T#@je*S}cU6EsjB8>St#%-@l<BQf~Y;dV%)Am|n(8)f!XqygP-!*Vu zv}^+3w}J0c+4K#-m3Ay}Q#b7fXj=z-LE3bb1^+JfSm3lk9bYi`$v!%--AHHCforCL z>!NKs@J#@|OJ&oC02ktzRAYfgF)6FQ!F-J#r-8PCz=yt(Wo6%78ok-?j-gS<yI7~X z2)#+bb=<&pfqIjG?-#)Lk=XPF7z@n9SOE4n(^!DEHhiXU%JIM>ppD~!_i#-$9(YLR z4{Lb?*Ge5H?}i;cEYq;3*vH_-cSkC;Uvc7Pew$0$U==SP;ogM*4ICAL#CTc59ntZ! zgxf?~Xsa=Q#BZjQB;m%ILKoa?P^Kgaw}wl=y;R~(>j(Gi_+HZS4}eYU2UUD8;pPIa z41TW+K5!5EqW$hgXNz!8pGe<Wad*{OXS{+o{hIHkaE@N?Geti4KP%sr<M^&d?YoZf zyy{ye&2isF{$nomx@z#gQ2VaO^j?}uV=VK#%wrj&tvOZ^>PRz|G5D8b<&&06t-t;p z%LBY<%Z2*AxV{0`+l6sl7p=d*_j};GRDSOeaHSc`G>R$Rw5_1+N5Cg-E6Rd@m&)%! zHsrc5<U7ZD#Na3UsGEK-0oPUo*9H2$1bnr?cd7i|eBeqmmT45Xb-!=s0H5p|l$Cw+ z5$MfCgP&RT4QMdp`2@Z38@Mio-V6Y~kHqiI%Q%)n+%Ly6Qg+`gNNdCQd(>D)!kq+M zDAib|{L+nO!l;8|nSAgi)mY}B%)g>#5c6T<jj@d351{ja>&{pP^XjdzuXsN>5S3#Z z%!A1}yX5)2S(wioG&>gFHrt(F3s=mJaoxG?%qOunjXccP#jw5>UWhcFAIE*$`^f8K z<Sj7rK180+$eU;6bs}$$x_A6zc-CzB-fFlU^Ou8E3?L8#P_|Q!Tm8MP0M`Qoy{*;Y zOW{TpU!266up#GIyY215-32(=&haVQnMYa;dd+s`<M~)4&ur&NJTEr#%yt&yc~UzW z2e7SNn^$qb)hZ6?PMbUWSQE2pGvGc3IMx}jtMY7nxoUs71pQtD8tZWnnFvqA{SMqu z)A>94$#m!8$^s!y9PPL~{_c8sGWuznnjf|&t9=YCdmONQ=lH~O?0>prcB~C;ZI}C| z@b37;ez6eRX+J-=wM2;JZ_dG<Ti#<`k@vY=OT>B}=boK+!SwoGtQPN0QjRIpA)kVO z1MA8*(&FC2d0a>2+9TI39UPp;misGeA@^wCS;(*RjO&4gT<gR0T_fuYD2Gm6UqD^$ zlfj3Hv)#QR^H_f$T^H91b${8;`eA-vYjZx@f?uIO{$AuW!_Y=N!w}VcFS1eZE5=$= zUn?>idQBMUYQ~t_v2PuHtpg~NCZ>+tH^A)!TxigC@;6~#+l`U&JqoEufPEaWCHE>l z!Wyug*LLCX{Qw8<Taw<ivGfrSzBS>%ehskT$GhBjPk)`Ln*&_9dcZN3X6!7_dcI2h zx|DG==84LJePk>xWq9+)l)hxcFJ}yb{a8zZKMlX)aW3G#1h~xnx=hnN(3%U_^!$?h zx{e6o{3PJ){7n4PehaYE0h^v*at==7oCvs-{L<qN#j9Q}uL^}vjxqdLo6uv7{gk~T zNgHI3dsVay&vD*OgPp+Icg)YR1@fA}^JMw1;=3M6@-6F7`{wrxoa^5gb#PvdWs!GB zhOr&|o#X{*%d`GAD(%>=QTuQScN%cT8h<6v<bXZM;5!oPIEk`p#(D{LyZ|_k?PyCW zqs@@ZYRF|h#(0v~#(ZYhZB_CK9_3l#^lja*$H)@yhuv+fQO6VM+nP|vqseWBOu6lJ z#;xc}d_%OWOrgKOf%4LxX_=zjtkp6#5AQK>uLf;=7`{i^UA1>3#(f*|Ic%?duJ&=n zOy5J_)CqsmzFqsM?eIge%N~4}6tacyyL<B3{~<i1PYKYcB<0J!9?E`yp7#5+rcK&+ z(yZ-!p%7~Z&O~|GPoz`DA&8@890EBUYJHacO&y2q*Kvr1n*q3FpG-aDIwkj>gb|;F zv5xELUABHDQSXLY*P)&!&|Z;b`{Q#o()NR<VOAkvE0S!#3#V+qgc}C95NsKEo64W; z*S256_6KZ<w*3OOpX;RV98N^qehG*D9!g%a*4ts$nUUG(?Jx`91^O1}XNs^7ya+zU z(cdQc5zfyf>uzR#N!t%S9A+H@-t_#EwqL^iJK!?&OWQl(YzJ(5eo5OeasB{sQoqi} zFIV@6S$k0?jh!=ef0*?c;8OFeSDar7HXwyx=@@D?I5ej4D_CH;WU*9&%;OseE1h~X z$5L@0(1GtFpiG+Bpwaa+fctyEWrz(hKcV#ruy+79O>D5=<LDFMN+rD+?2C$}(C=4R zw*WRxY_Q+s=o8?s2VBb7;KE}m#H(`u1N1fD8VdYr_!XBW!1V`QW`1dX0&M4qO#IUN z1h~HdE<-E@`{Tq9`PLhNP0uf_Pr%s$xRm_rRxEXe^&gV9k203Zw|<B^vW_oYUJSSP zqFhNJaSmSF;<$`g4!0fy+=azb!>w<kY#JMrP{$Uw;UkNshFdGrw>2)G9mB0f-EFI* z11LK$eOnXin4R3#ZpBjh)>M>F7fUVB-zjzDDIF8(y|6N#;uu%HE6Q_0WSoUF(=bL? zdt|Y<3*+8!9bzb*$FY8U-Rv0m$3DZ@NycD-SSQEFIv#UiQ{j1A$5MFLW@w>tt_kAu z(*k>>5l6Mg>dLTJy&QY7XJ9Y(Ozg!DMD<?n2JYLA)wQW_g$<Ht6P_0SE_<+HYp>)u zY(KwgcKS;BF1dumeoVDTeUKX8b5D4g<e`X#^VRuql9oz!HWc3lJB~R3^BL#2$d?$( z%I6VvRvgQ7-o`mMhHy?k=Na0Ouigp6Sx$IIjOVcN47_k2>dp(qBrMuO9;m(Qht>Dd zHuL*v@o!q$<FN-2`_G%kVqbbO_N5nMU;0SwOV7u?^gPXzn^N#Zz00^?1V>s=Iy?cK z%af6Tu6xJvOzs`uf%i3$mV3wf%-lP^9eMcXBlnJ{-Ityp=(^_}&m`_`D37$n&1WX= z2IMuA%02Jt_ND7R+<ITPdhcNN|AuT9>V1k@M%g!#$5<Bzk2SwlKc7RMBg-c+-y`9K z0f#il9`bny`$+A33~$AI23t#G)Qj*^kIGL4I}UnOtzbVWR%ZqN1Uf*zffluI9O)9I zQ>~-OJ!Z&Hv&XDauX*nmVw5!ree2-lH5I;uy~)7q%)9tG{#Sd@8;PsPz!j3TX0bOp z1HOU4C+GfWVNYvkAGQCrMDKsi)BDEr^}g{r*#C+><v+lFaeGH;)W-W4a&Pi9+(Rau zHinSS$T#Gt;+t%%TF;Vg{p0X%@(p+s_-5cbX5b=ye82C!?+RzY7Y07gkMVnE;m%yO zr@O-8o8Hr%ulID<fIrypeE@s8?QP)O0`d*}?H9scl;`4`tGPG}tiNR!<2$T7OVxRO zc>#4^U$a-8*VpV(=k-<RCY{%ZeZ|FrX!F<r_7#hk0QXb*uf`sI?Bh6HVvRe~xr5(} z);byEdmZ5#PgH(KdTFUUU)pwwWIsV(gQKi_48C8K>;u;#1J^~%KJa~(eaLqv{-?=4 zXd7#^H7mZQlU3jvYT&wP*#|yBzFh*@-vimVA^SOar=%J3-{gZH6yq6mKz<V#Iv~FZ z3>~1ou%QE-ld!9EExw1sJ&~{B+l+t2_Y+^mJ7<m_K=!e}3~?pCYf8D7b5#6hI&6M| z9@O&e1FT6!l>ac}Bl*12E8ijIS(T4XyT2fY_r&FM?Qa-QIPao1VT^3%@vJ#{uYFIF zFFV3>Fi972KaTg&&tZ;C!fRT>dX0neh{B`J8gbIHP52hPQ-9k>p}#n*^(*<^w83e8 z&1h?_f$O667qs08eE%c$7qm?_+Ims?3tU$jxGq|Mfp0MRb_w)%ce&DEUjXlaLVrWh z->h^O@Je@cpu6f^3rBYsK_Bs*Q~Iy{qH&^C_ilO}PS6GI%VuDIaO9)Vh2UuGr&nAe zT>!493|tqj3&6J%`2I)g0%%)jwDqEN0l2CRTo<hi!1o#Q?W53zZ2edUUEujHRnUhb zLmv*FPSOXAg%f<*dQViYL%iT^J@7&qzt_@QI3T{heR2fW|L_dz<*sl0QHNgNM%-Os z^#xB@F5{Y9uz=(FMqw9PXNN&A2B^8JD!yL{AFbb4X!PJ5M+fdbz;Qlb&f&v{m*Ct^ z&UFNb%J%>oMF8g$A|{sK;yR??H-A#UZ~heDH}}fAnsF^dp8@ZFPvC0U5?@O%&>Wfv z`-8j(F{jWU_k>#ndL$ga`$M=7K||JU!`YEP%cgCBJEYqdLc7?$BA#i8`gtb7Y2=ZX zF%HePqUCql(1Nq|B>zxHQ$|{fUAS*Mv>fhA%da#o;F01>hoa?&qy;qY<rRIOan0B1 z`@mJ!{eUmgeP5~H!E@&+oxTTLt;U;`OaXo=;Lj^lfcp~QI5#Qx@PhYua*ZfS#ySr# zDd2t?%7J-CYajRtFJz?%dncI3?*V<8^<o~3JgD~yXT9`@UW0F>egvGhaIHk^*)`6Z zB;rlZMaa62=(GLhGa<)b#9hgl6Z1T<g_v7Y=RkCB*56Ct+zWF{>imOuMv9m{7IRCO zQ)<GTQZwe1cIRVmNzEzk7lFh%B|qm1B3D_w3dq<Mpc8dTICoA7F*Wz(upfASn!OA6 zcwaR>m$Tv64LG+zo-x4l3rzT(fZwU$|A?_>AM|$w-$H-$mHUq6;#=mwz}3b$)S1J4 z0KPc99^>Qnf!HUwhJ-vlXR}{)sCgdF0fl$4ovK_X@3#wi=2LiwewS4~({rEOaGo@t zJNu0|PZZv+^BR!XfOYC^`hF|s@sO5tdMFo0o=lr?kc;qU)`$6`&8ogwIp$Q6mixcQ zAs=<{oEDZRz3s%$dh3A`aHQcU^#RT`f!FHX8PddZPj3z2toTI*k74X^2<4wf*^Rjy z3l8D=#%UW1jw8+UrryEaBlKr(0Cm?7g$(lC;ww?ddZF+Jy?pP!6*991^?<K|*t?Kl zS*Jh#%`HFoSebmb4L;X`&zpS~e9C^FoATP*ojVHn&H-@A_Mom889SvbUxV_H7to;Q z$B-^TI@NrZoF7Aent7fhoYV4CPW{(dE7)F*lh;)E66VK%*YPpHk$~%e8Ms1{)-2|E zGT@sDd@`=WoLIQCx0)ZDtmk=-dDQ#iFL>1Z;h0C@SqgsWg<Xy`=4CiHEO?&YM5Q}U zJ44W;jJ!jBTHaNTa(Oq}Xm3`$1KtGQ8MyizxJX}Cb7UFt{WVX=R+twHcj8P6##sv- z-g(uz2QTQk#iig6=Eo%Os*x@P?=ased51a2rfRPhs`glRwa?m}kNK{Vf#?pO_*Sy~ zN__|@ng0N>1<#4PCXGxVhklyyQU~gw({<8*Hx~Q`_8aR0eD9slxR0d+>7n40tQXIS zBf%DFU1~>O$7{atN6x(_zsK_u^q=bxa-9KqsMpjxReyYq;hm8@8>nnw=jLAe9H6*A zu;Bwj@BvNm0nP9My9dD+U|vMx_*@EKU|*x;1mE)7j=pf<Twjn2TkiORt>F1qW!uyG z0uz2S;5RGyw7vj+nBWUA4;ros#O_mbCGY`i-lY>fq%C#U2PB>TxZPsDd{(LEOq_3T zoRc;aHWzrfo+Dv|D(CDUmCsV;kHbDv;5~tB7=13`A|7>~9Quy3!1<ss?0Nuk?&#mC zZzMVMEjS!U*<RpweHzM>4$4Xx<=aq}bkOH8{-7K=^^h-I*T6Xvv0?|8Q9tF1bvrPW zC!R4zy4v7J3I>0<fctizMtU@!JN{)<<BEbh_!Q17Z5+3;;Ay1Eug7q<p(%qV?@U_Q z2fU{rm37l^!e@~elyj&3jGcH-UpP8P%xU|uW7Ubmwj0BM$M>^hw?N+H`)u0hGy$*l zI~|zM9<1}=Pb5ErNxls-lGe8YkIQd6E+gUPs2As>IKJ(sx3q5qF7hmlx}<MQSAHqV zL;hX+jC4l(jQq6r+3{`oHro$}bkny1Z-RX`a5WjYLXy@jd|L*5-vmDC+u+Z_oxPNA zgTJDE2yNHCZHM-4@PG7a(mqc^8upj5skG1Ilr3<12fL7wcgRo6JI5b^cXt`>CGPXa z1;CrYI|J7|16K%eS^2gM_)5q-?a#uU{giKmzf!!@zHNv0ZSa3Ke46ClLZsopc$S&u z9pZ(iV)(XV__kvBHmrj;bKcvdeT*T$Qda^gd|L{cPU+jU&!Ig_@NMuzj&J)C+U5Av zF1}6E(hm5Hz74uYJ1^;lUqavM^>g?&Usu0&-f>Q43f~sIR{1u3%k>ZV4vwq~aISAN z$1M%udBb^pn+b3D!nD3^Y?5!Qg$~sQVqbQAo3`b&k&b;lR~Wld(&<-y==M>ix~H9{ zzlbII3HlOQPWY@*u+5Cm>1Vhm!Sii-&*u(|%!9p8_@*DmI5NLO`#58)!}6|81>Mpg z&}Py<!Jc1;{r?gAMB0DIoZalRnox&h{}&I6`-UzwjZ^LaC}^@An&9uBLSNa36Y08L z)5Y(_gyFX&ZNVhpgmP(plgUTh@J-MM#;uNT`r4b?HvyOQO{gnr+=Y2Hqda{Z<d-&3 zo;i$kNg8`7$6d%zV_)~{an}*i@iFT%@<Zd~^}O~Kcn$lSfa^p0bc4RE>}v*muK}O* zP4G40dOhyK7)sjL1KKwo)V>LRiS|_5SNf@9*jM-@X<y;bE|Fi~&v$M8oVx+<!bW?u z;@#FH-Wj;QYv8&-zn(zbx5+#DCioh879i6YLn+>A-;{I~pyVC>6nqKCu#$K1XHDhs zP37=S<?(UXRH1!@A-~X-M)kc)wcjVjxXYDkZTDPwsRNu#LA(q8o=)&hc;@(~AE8~+ zH{n?q-=t}2M_xwX1YM(@mh{43pzow_a>recf7+bu&ch|2L;7|}pOvZ3B9%B-ckxw6 zJeAiexe1K3CZG)toa?K~Q|!gUGmOn3*Vu1`wA_ouXXak4a^|VMSbxTTlC<Ntl0bAA z&L9nY)%m6$=lleoN&F?ifwaWWXD0q>$eUIgwXu)vufAc=rSpXtyWxA;eXYW?H|^tl z13bSF@-)|;mov4J*9AFK#lx>1EyQ<e!;`V+YBKg`!OqLMxQQM)&(@A-a=fAQIUmP+ zj<=`ies{({;qmxJDBA7X8>CCp*e_`hk$+x$7#U~%g7&{V_7Hee+C$(vXy6J-TC<4N zGvI3mK4}jfo#wt3L#H2L%$}sv2cXYk=(2ORXJ{$psXBnZQRmkRo>REM+j_pf!#4pi zUu^ex=I;G`H&C2@sAxds5%PTn=2$4pJUiBdt5;v~J;N5SIQvk&@V>bWFk^Hc=iEzi zu51(VR$xBU{5EzzWRvCdv4*0~#g2a)8+14F+t^<MFPuEmg6G)7alKqqCJeugEuRIh zxBd>Wv~8V`{f;5QN4_J=LSE~YKRpS$aLxmH(#UUP%W@spTcx_3y#ufd9>hDU+P`cU z$~JDsm9oAK`O~03m`CAzwVY1@4)vYsd0wO;|2w>`;bw0u;fBk6JlpGSJ-D~572kbn zEfmAuZ!vjBVSNtIFhy&w&oy#=PSO^<!SaK)(B21I$$P$Q&GtF(>d1CrtxLj+8?9k% zFLW5NZ}Oc^rW0V*9MKU0*#3kC-liSojrYy#03Slykntm);EgWsbC=<LrJV|vxH)H- z0BetTVb2hDRv>yO;?$yz!*)L7^)ztJhHIRml2Z64Z$k<4Lg<&K&E%uEfqtRnUjJFP z)dn4G18HK~XhWPPe81?wP!HEA-1Z`!0PD~M*#96ba2IMmLWbS7OYjfxCQMZH1NJcK zmwN1#e56hrZN)m|$6eST6Bg}gMmw-}+fZ?D)!Fv(YOhX|->KM`w=?L)x1#Wk3#^s3 z*t1H}UQffIK3Gfhi8Z2c;69$Y3%b?_p0l0}jrbPSBdFK0HL^d_(FoXWfK5jS<Nxp} z5#zTv!Y2fNm=8LI-vZd~g&z25$ZxZ^^%m$UX`_rV-peR2kqOpX(kJCbe=jbbya0A7 zc?dlS&7+O*EI+w-*9eEFLh`ioZ2LsTQ~B)=d`k&^?(&X0r*tddS;K4O8lIy=bAX@u zf@^qM7e_m5c+J#5F_-)gJ_Y_yraWTpB}c3o!2Xz~?j3uWeBMEV?BJU>lh|%4JG$Lz z`{icqV}QNTe)()K^v8}e%JUjA0c%7fJyD#yjyYHOfm+NRgO{Y;hqi>!X7V@uJo1hS zoNLQ+KGzPF4nL0_!nbVB;X52z*`dHh>kMGCvO|D<ov>NS{6y;|!0wj3&B?5nz<U(1 zrtNvfwLLiv(;@Sj_+9{AI0Sn%3g3tM1LjoadjtF)4d!4sVI8jx@l5y_;!<9phu+}Y zJ`(zZbMM=aL5H3fv8;F%nP`0{TV4V7v21w-*oU&^6=1hsB(FZieDNRfTV1@d%eQ%K z*zmJh1L*`GxxUhl>qzkO__#WGeKfaD?zcS7`3L+4z|Q;uSFV+w0AKkXpvS@Y<Cx=U z!!z~?=RzF#r!n``h8X|&IO4y*pq=wki2Id{d%M~L$nqVMk8yiI-I(Us1B~6j(@VyN zYbpmq|EWjNABhh(;Twn{F#jBVYIz{mKpj+e<@1Ee(C2ag5}0I-WBrh0*f8a%kv9Am zVjAgR_6xxJC9N6#^nUH9CHzppXYkXYL;7h6+Yhis@J|oWKf^ys{|x^`|9m5KWIXi1 z^wA~vR-R>yekaOtzULhJ{CCCp)?L!~nAqp+H|}$#&hR@+?3<VRVGUzkO~E-_Uz&!! zutuK~o_$JKu3b>4ViF$j0O1;&y2*uid{aC87Jb%bIo#(Cp9MKM3z}ZSFV%OVIzX5C zov3!`@YA5L9Xk9dWfNEGI-j>s+f?uv=;65JGRS1Yx1#_j%K=W7BRtBr;kj9^y}Ys@ zg5OE}&JIxH*%I*0p+CH!RQlfVLR_`){h%1|;3wO|G<i#3moq#bCv<Fp48lkAn^pAP z^w~0<!hipp4>CG0qu-9p559-D1$1Wc+j02;>;}RbF`OeG?Rbte(PLTH3wggQ+j;@J zFipMRXT5Vc#>CueG3@bJ_(SXkhi{~9#``R3<3Gd;GXC2yA|JQPiC_Aels-JJpMaeL z*mS;x^Q?@k+UAtT{s5V#Y{zhw_78UN5T18x``7lOh|%7KUqFAp06BUQzdX>Lj~~9* z7v}f-vY|h4i<Orx{ebP8fquq@^o8BhDYS{YKs!yF%q!O=$ph1FC+igTU+NU#9)~VI zjcY1<@7Q<QuRibw-#Mg>_As_XymN()cYerthkgGf_{l5%@?OM02k|?E-?Jhr<>Wa$ zM_i=*Mf^Ltr8jqti4Z5A)?+is$t_lEF3x8p?oQwUz2Oh(ce2#Az+`J5%4el(fPI{> zDRk`=o=aWJw%*8Ot0CKZ0UJzH?+dIKx|Sh!i|bTiinWS(ve7BPRs;4@=oIy-Z5R5c zSx9|iU;QDMItICxx`@3Q&@<_)AhU<?oG|r(c@#0mOVr0*(8FfLgRt>Y>MU&h?|H@j zi@4Hf(8t|A?Y@F~z*BbVyBytZhh9^s+n-dpW&L}N`tjVYAJ3ioA3^;u0?*UBe)<E( zIp5B$FZeIm5bgosn7Lru23g;osF(G(E1f`nsCTc5gD~#-e7ZfJ7{A2(DKf>%VZUY* zn*cW2AF%0T_XNCe12)5$+OVlptW$tZ7nj@yn=~CZDV=;cHtTPPD+^x4?{md<UE-5U z!~q}5*rZOzq>rN?e?JbfBd*60-+xK_{r>R#FCu*$&nZjH|2=FP;pyicIQsdQ@H>ZJ z_(Q}7um!tN&KyIGfV{Qk-Y5HvV+biDT%*PM6!!*k&qrG)Y%TiLiGdE`o^{h-IDSIX z0U0&Fk?F{I_*vk24!<+N`96M-O&QB5UkDo*h=mWU=ixA}2a$gWzc-P7N5tA8W1KH< z&-KLMZ}1%vJbwy5ErV@`6@6{tK#VckA>5O`L%1h>uFmN890Dy;XF!+K8Qe>q!F__x zfS2KEUB^&vy$5u?Pw9GbFYCp9%6iczSuffr>&3mS7xyXaO&Bu_lVb*H>(v~BwDq@3 zS&xqyq<=yE@FzGM2-icnf)DU3;4!Z7E6bg;>oKO-?p0$7KlZ;%TW}0*VVkM*&tq*6 zX{R6geb}uU{(NJQFZ0Z%?8YE*eGE*s`m;UR=s93J`zbxoFjp9lApw5|@Y#()rdqE9 zHl40#8iPy+{r{83AaB7|IP%)9c<DtwmUtnc;vQUOyoBeO;-z*8hp`Ia=yRUNmGRQk zcy5++#v0C8DFQla6T<L2QvM)&VaOomk9Lgv%TmRvBk?>z*J@xVY6G$E+F3P@uPbQd ze&~_#ulS|1vy2@%H(<iGJH8zDC=8$c2l!pSdr3dI6aMlnY_u~DBwZW_elErLg1d|Z z!MmRT7kNrQ`V(A#0$bgN=U+1B8j@pH!pm}iljUR|+$YPM<u(E5F3{$T1vfzExwa|C zg4zchLYtmOn=Xic^rxf)eJSa{Rnm{=CLLYqmvCylE6V|%KI$iUZkF2w`^<T*U62vZ zZMl9$&v~`2hb+{Ac9v%!puF@`kSC|j6Vc{D9Qy*6dnVfmOT2{Dvgg2F3s~y*=_2** zSMEP_${oj;mt);_%{M1LA`QgLy^n1e6SID2Ol{6*JM%g6y${1~s`wvszP(|ePjUa{ zssb7Rolx-^=dMiqT!=keoO^{n&F31YvxdpNiKDZ$1%agb2*!bOuCbfnt!AB^k9rOD zo@e~gsOO3rMP#b=wZ0k0AMh7)u1LaeB5cz5L!Mppm(KZF*BQ9gT9Iv?fUQbVX9r}e z9nbC8)L09B(B<qu<-6s09<<+T%^{vNV-E?7IR!u08v#25u<2q;?z`f6xDE0VevLXt z8N;4f^c%<ZoReX!_95(6E1nak1TfU4ddMSXlXBUnWtV%U{y*%mpd1B_dBuCskLl(T z6WSJ-W(_AV()jU&wgGk^VbjFx$#}hhO&6~x<9#1%S?PTE?ciBShIpO(*G3`#kmP5b zjCFUCf4F}%cSFIIQYSVPES!e;QtEafR>v#&OMX5LdAJ|H-{;<6@FDyk=Z(1+?jHaz z^;g@b|EqpVubnwQ_m7AJoi*qn`kCnu!87iq2yi_c`$!-It(~Qq??YSOE}6A5)j3R^ zn{(xR0_GV^vV5@1J`}8PPg61-yv^E%x*RyJ-Fr}%t<PYh4xh}|-w8%PNS(%gi81Gh z_gY#*;{&Y@4$|7;^IR9J(`O|4E=#+ACe_}VJrm>K@s&6mQ_wHuJ4214W4cx0&<{A5 z{&K8s@jJfI|5EU@80)#*J0thG@Ju75Q|_a2&NSL6cAoM1o_iZ>vV%}I6018F{Zrks zE5(&9zr$7TnE`&ZRpRHmd#HqKyn$$Ft0&6z7<?yKKEpm5>~VOmC`VYK%_^_RkF{;U z%(zxv%dQGU9~9o_$^g?;?};|mmqxGJ=-r8PT%IdK+5g7576Vbo7JT>jAinpzL*%zS zjO*@Mf#?>@X>Gwfe7oyB(cN{W(Jg0pfBg->^Bcr7@D1X7@b2C}iD51CaZP~7`@y@{ zEBMy|&+iK70RIHwwZ!EugNgrOX>`}M;$0iQD1bc<+&g3A8P_8OWz+V@_gVyRS2Be@ z`CD(S#T`P~tHitCHOu0eE_+(w%sRds3_O@OyTUE|Cd!6Oe66JQYiQe5kjdSXeXV;Z zm$&Z5Z<pw|{FN7Wjo1YlJoybVe-EDJV2?yt<Sc&%eNDRbJRIbs#Cn0WOx}!giF+ho zQFR~j1m&Kr{rv9ZF9Ev;=cnuvz0}$Z(xlggwa1Z83HuCSw`}z7T<pzh*z)I1$6f^w zcjHW-EuvQg`{XM)*LJta$=!l;ZFj%sZI$zh+c9>=d0jHyjPzd2-9CxmQ|Jq<E$@Q8 z>!jbP@x)qQD^qQs=4<`VwDQ*QJY7!N4}}Bs=s3G7iu0|ay8*LH^mlnvRSVwWSsLE# z=h<wH<V{I=>lVz9?k3MmGVlrC%$9A(zTgMRC&;MFC!`bj6p{M5pXaHpC+z0ZXc+we zG-Sa>ny}OvxBrE3aW9vA_3WG>avE~*+k(9zlDAtGZ^<LBeQ;l474}%~0sXtdE5+-; z_O$)%(2qUVbF-x%u(P^DKj5bk9(~ykI??wT+rjmt@N!SA+ieg&t0>zx0QQRRv;pw_ z0e?Q)DW@jgUx5u`ADjD-UUhWL;XQDl!TRZL$P4sU>6n%s7uJy-z`jmc$jh$2iLz6f zN_GJM2f#m}<>B9??qwklDEA*M2Re3%-U&3hdT-Eln6RMfgh-?*r91=oLBO|TE$=MH zkn>7qJa`;2NE@Cp9&AHAa|q8letw*DU5EkqKyE$(I|81x?D?Sn*dClA@H5=+fiC<J z=N#>UF6<G#8=gLoj6H_@$3`L!T3#A^1hJ2^-f{@{hfrS|u1*Z~4;XXbkC<sB#x`|G zJ8S9bVyG`6hH8T?l4GfPutR7A+u5#cTr4YFV*9Kb@Fy!<0@%fb&DNFx{tJLVpB$fs z9A``GjBe06)uA;T*#)gP5FU1;4SmLISoiHr*g(i(%NEGN7RUkYTo;@8IAA&6pgeId zQ;rjn?u5+UiGHF?JqFn-z%!1g&moQs<JX3{;<J=RtTD?uTinaJQ`}P?AH&$B4L(Np zG4^|LExH{r?aFt>{sQ}(fp_~oO79y*AjkSkPR4%T?@>A?Vc#HZ=6*hcI90_Az<nI> z=j&(s44w(pf$>vTG<W1!|0i3T0o&FcngRc>S<uY+t3B|^Y3yp8zrgnx@MW;8asC4K zA;QwG<|Nv%K-{jz`3v~35<Y!wkU&QrU^Cd#1Ugm|HnS~Fpkp!Mllz%xzU`K=LGK2} z2In&Ca{{fO>jtee9a_`ta{{eX0q@xHUeMiC_M1B1271Fa*z*p(q~1F=JzKq(`=2o9 zW$v?j1^Ys_j1og&FNPo<TEqDL?8uL9pH&;q+-k=hSUb*H<$03QUQ{zS#dGF8iS$#R z*q@+#X3W5G3+PLT8MM!h_fxQ!h4&3I`G~ko0ro}0W@|eD|2*K6Whyz|q#iy2`KQk} z;~UV$@k;`&hk%cIcp@pjabX=hhr0I@mU@_tOaXp3+ot0i#=99}t^}GM)-<Jwxf1x- z&`p}^37?H^0R8I#`=pGsa#O`w33MzYEO=AjFVRLBv0q%@fcr~?hyK)i0<oVXe!2tk z5#Lk!Idt#m@SAsedN=$W@za~W-YrkSUi@6-Hv9ql=jiD9VlBq=jL*)&r*>CQX_t0E zhPd`~2)fR_kg0U^Y2-bMF#%yX&s9H;>%<!i>Y!VkyLt)!coWAA_;D{J;m?lTDCf3l zLkRz(goocwf|oi9`duA`y^-f0C+VmfNA`8rIUz@KuXH2Fi#hCTj49=Kahe*JE;`$p z&-rpwE;*LmgE4nnzZDs)e7C>1^$ylDIS=8=Cej=q>hZvS58ms=xD)mFwq7A@rm?)* zrzK(k3$Rbr88JQ531d`eUKe=(gRnR&Y7h1ZwcmT!S?;rlp|2S4bmF=nXR5yno4N-! z6|p$x8gTXw&evf+V<`5A8B58zN5sOMb9@?c@y2mD=Lgph(Wl(kR}bDlinzQJzai*H z*2fr3;=xtMWOyd=$a27gauN@&jLnG8#KZoTbD-!??x&2ER~GET4}JCP7a@C2?9Z_W z=cG9gdOzZP+Uh584)<5E$CEy)4d-v~!TlE8GuA|m5L5FmeG+5>cJDOjSL$V-XqnLS zG!K8MWFjEvO5}N^d>00D=)EnT<r3)2_XdP|)^Ux+Gm-s#hh{zP>4S6a?V`1E0N<Cv zTqNE*d;oRu$}`kL(DBE4kNKgYm}kMbYpxA>I{>?AZx3)EKUj(%`ij>=^i@$G%sU{z zPz-hVwPXHfehAM%qWgXg;304F{Tj$nAM2B>2eu73QmrE#(fa}=>~z3}AY*+0LB4Z= zwDVpF@E+0ob|vi1fCZhu=6fHQJH@;q&WO|P`+}^KwB{x7+b?yXk5#PmQT}qL-9Od) zlf957xpoCylL5zcfPH|uQM93z_bKO&eV6IE<L}Ay``^jM`bnXPI%mvC9CnI%*hn9x zn|5C-hBYq-&Y6cL&4r8cZKQhQ<9sRJ$%IVa$~Itci}`F8pJDBTXCe?z&V$wBeq=zt zlkPk}0a@l59=tCA9Ou$_776#oa?P3dlxdD}Wc{zqkndJ={|ez!?U@ZiKFqS=@dDp+ z!n!5O$T^yE>i!(-#yTS70q%wAL`=oCJ;pvSp}vkY;`=tv`6<BrjdHIr*B9U~o(Amb z!Bui!)~H6D0R}tH@mBjd#($d%$Oo<;9RhByWp2W=$AR-j;3Pb0qwnE!r!J@L(||9) zxnhr}$Rl3zh-EsFA08`WevH|8{t(xa8FR`u0Y7;r-^0XuA$dC*<xDv4Q?~K@nLdxg zeFt-?kbH9au6g6)tB&I%KOHy+-@hi~s{lT%SB9%W6L3kI7UJ1bJgZQ2Qik@SoV}py z^B?m0JoMLQ?BfJKKPav&_!P!{{h%A_?3}oKaBmg+kmtS)q8;Y_VI?;Oug{2`F5V>P z8T%qn@&foI&rm;MrEH+AQ_snR*qfaKN5>1h<eG{3K9_$m{0G{|wi1tQ=N(3S&^K)7 z1DZeKI{-(1a9{VVu5e2XIN8_bfSCoJZALltE$2DH<w|zl=T86QIrWYHhOqQQtot7{ zPPea15^<_N#JNa-lk4|ReZcGFO+_Bp>*ZMr(@?$yG|l54E%?l7NH4$-X_v1HH9y$S zIe;gO1JAM!Jm{437$*Eyz(dB^ALJGLjCiT9@;UpAX~?_O3HJTu)MX^@$%@YK99)6N z(Tng_O*3fVxlN8O3Im?_&yjY(B$NRTr=2`YV>IOrwDUdn<KQ>j#j>)0(QYU21oGO^ ze&&VIKDRF7=QH4C8%YQ2WZf>ENp-fP4)W6Fb%Vm&hUct<JUK@@u?=!EPSGUCr^u7E z<5>)HVd`g?v`<uY(QdbsF8M58k1Ic8LE~6;j+THwp--Zmy_3tc&W_4&s>JFJsc*;_ z{>M4v4F07N9@;(sFrwCbPJ8_Q@Qm?DKWib@%pEw#54k=>`yleou|u6mJ9d(MqaC*M z0@0DY0-kp8e%Na225mNFihTM#=<ZSu>!#-`N8PL-e@tGWO%5+&bw87Q!2KH5<=Gwl zBKZ0qb*Ny}rh=l4o_Sd3`go=I;1+rQ8DuyDTKN4Pw#C!g=N7i_uZRh50>1`(15x2E z{H;v)Hq!l#^kqhRh>^a+NRKem1x9*|k-o-AUuUFmG}1R4>06BStw#EGBVA^sKVhUl zWu!lAq~{vxzc<ofG}3n(=|x8RE+hSABfZK<uQk&DXr%8m(wmI*gGTzRM*8bU`kO}j zAtU{;k$%KTKW?O(jPw&KJ$Fg<J!@;1{o@jUb=8U${&h9}RW+;Dx(`>ct*KtJ_@1>( z{9~6aSy4Myl&<ryShi}3e}cc}F8^K2RxI)J`Hgo?h0yRX@lJ?O*VL@`FI%@{ZPmIZ zwf<$R?pn4AcPln1XyM|UG$l2{*{c>Wxz|7T>J^K}5@7MNl}lFDE~{Ck2&h`SwrYca zY<0~_JYF}}zi7#dn)QDCmaO$Js<~&?VwOQ!A-adBYEcb1a?k43UGT{E{A1PHWmSvN z_PeWUQRk9X{*^V0m)*5&$zs&AZpl)VLix3Kk(l*&BOAP}zGp?%x@D`D`q$mP#9y@% z#H}N7E0?USS-XLJ2m5QP!6n6+(sk>WtX!=M#Pa~Oc1hLZ4KDCvze~&5tI<6!Kd0QB znxCtxR>~S?u3cNR)<HwEu3mAEf6)eXlB$39D%HEP+u6NiulHB2l`NJ0KV_nVTX*-e z>g7l-SzU|H6=GG5Lg2{Bl!>Y7$f`d`Y^F{{y>2}pzvz0>cRK3HieAZ#%mjA_PbPx0 z!lODRBS9|zCSL^qoO&`4oE5zgo{Yq1+Yb^?24YOPnVP!Y4qv9;be{Y)hR`jaM*iLW zG;&*0TMa5HcbWBY>MasXlq%~Xty0z6rT5St`o{vxSpTwG*z|RduER$0&#ga=-F5TN z!=DK^rDo8&m8&b*hV#)s;id_ao~-0UHmByEb+yYDFJZSVxp(!F>UB#N``1>jTDru4 zEiC=5_<hW1AMmH_z_j&AyQ%$3H|k4^Pgbw-X{NxFWz`1RPtwCFU$$~t?aHck)pz^v zg6U4}3EdjVi}9Kl#FLgjv%PNVOG;nL8j|VL{*I1t#>~=Be|k=N>HL}gdsaDKaWP^5 z)4xr)DYb)h<aWYMw_Ieur;jDF=BJdO%yuYkznF3(@mwH1UGQY4C(ggr4#CZzklL`t z^HUqPcz$Zb<K|CHjeq^B3aL!&1OA=YzD%4X;Sie0zXUs3z3!f>6$on_bE)EujP`QE z#58<KX@Aor%t(ihBT^BPv7F<FQxwctPWsCfg}||l<)mC*Ksl!qQ~0KM|EG};lh4Wh zcUKC8srudR@C<}F<(%=(q}1}5mHx)7Nvp(eyC+TRp4~F_J6=y(CC`dp@G&i+nlIVn zNlS>vll}O!`ffV;iuE;Xmt(wAv3S{940qf(MzTSUN4vtsN2fxh%gvsS#WiXqtQz3- zTU<I+zEfwZlF{0AYdOeTv1HZKb$2r&tG&BsEk?x{4OXp^L&C*Ns#jF4MJ!W;(Y18& z_pDrrfu1X?aDOtwQC{7tH*RGU>#eR?wQlVSAuuqN<GsY(4Do;|Pwu$78}WeFm&x(C zF`??L<VKmztSWW*FeM&_cE3JtZ?mb;tuHQc-LKEF)7ez0^~9|=9`+>E`~SChE$~eh z>HZ{5dA6k@J{FY?4`~Y}ZPGXLQj6;Yc?c?tEa5aw+AL|3CMhk5SOpan1*Mgj)cRU{ zAfVy`Dr;9!QBh%61ym40ab*!vQ4mq?_st{c<mBYRb-BNLe|HYQ$;p4dnQy-N&o^`C z%uGXA9>Won;~n3>k~1n?o-$lT*ppOw%6OGrp{o8~8pH9ED%VBl5XsgRn;=|YWm*|g zUoKbAJJkskx$<1om0R@-2<Nvd^AG*}`W+8kj(AddzNBQ<tn65AxRe+1+z*dIpe4p5 zHf*0wzaxg_h58K_mKP_;U##B&l4|!5z8=SawkdY#ccfD(C*HF!Tuw=X_KY!*hVz%{ zw>OvDRPX0XlPmsgQ@tPS%k}wI^vCySTP%N$H;;?$X0hL6)z2-GKmH6RY<~=YO{$$q zIq_Oq#j`B|k2ZZjNXm}zcow_PN`3X)>(t7LUFYQIQFXz(fX9b(IJ~%h%43zv`8+BT z!4GYQw<r?rTh!LcUPr069Hn<tIWH$E9!R?<_s9y@OVtn6C0^`>DIRWO;Bvn&Aj4ij zJ!~pCV3{8lwvwEuHN+p&BKWl`eyFx1sd_^>v0I)rqtcF6r6twng!!;3&jEe8(jWBY z%5zdzUZ7vkaQOxL_9ewbZNObs9e2A);h|7@o~i9TrE=ns9<EPeg8CLK&k401;ru11 zXeTL-t!ui1^3(*+)tL3w(G#b-$<N#UAx@gz1TAkHX(6`jfie4gVqjn4Wm{#%u?5 zy?={6I(p*SqwtqtJSa)f-&OM{UH?|ill0^^)q1kNT>pu|?YHGCFG9NcRqLht@<jRL z^JBBd@834nJV&>Do4)@h|D1Kzc~$$*iZ+Pn4_A0ZB(JZwF8(?Tm2=g+9g39cO+*4N zFILm6I;_CuCF)^RPp)=D9l2`XJ)-=i`kl)oS0#9qb30Z0>i^wx;&n?YLaAzgRvLr9 zr0VCyLkICIcq)40RUpihy4$GQD|tVJ6DuJrpFW14>&j?A<lGIz@*>4gA`+ZhxwZu5 z+M>%<{2$8=Cpfio3lfxD5M8dS!$Y~@1j(06UI%cCk6gU{lu?b!dBI#>sPG(-AS_p= zkr8rOS0&|y*FzMkW<8{^6cciNMauY~<WS2M$B2+8<=0?T{N9n|`vB1suXm{Gqx9zp zd&Bir`j5UmMg#Hu#WDEja@BaNtG^@$pW*xo>{pGGy86o$ztxu~Xz#KF{FEiAziL{d zTfTaqTb+Q{-<q#akRsP7NRjIs7?!Kv({#%(PEfw;fvcZipTKbWs`nCg{^ZA%>COPG zGH`wRY@EuxB<Kv``#%{yRQ)KtYgaREIGl*I<qx0<*?^_H4+INt^lT12)`2GjPjMS4 zT(1IU+}AaL*S|pNuW{3<wjtxEQ`{$hx+Lu@Owzt0<#jb~xvKFxemcdC<EK-eH*wPy zE8~i8I^9UsbEG&yIm!FKzZ#y$++D-|qL@eI{IPl}^K&@TtMn^HUoFIk)sgeSZm%R| z%Oy5Emjwrks{(QjIW=s~9_#eOx>p74gTnh)KFYJ@T(hg#VuRDB@TJIN<zK=v_$5a& zJM+zh%@W+N%f<cEkW8TT%9>m=)QBO8;|eN8428gj<y2X%mO`RG`WxDUm)Yrc2CAmP znc`p&I;$(-tD*5BAGzn6gD@O$?U3G4)34tEtJJSwKXZM+2YbG5ufHy2229dpe1%-8 ziIpVm`@y;j*e8Z{KinXe&9c`W^o0UGzuN%^q1<wV8xE_$Na2LT#~#$<_D+>O?n-kA zp7+j>IS7a0;1DvLA%^qOK~`e6yCBz`W%kI^oPq3MPOjODXPO6_VUpvv<>h6W&1MgL zw>3a2t{ZR(U~?AP7xg-WL5L~B6fp3?DPKT^i0a|6R<60;&!x32V^)=1lLx{qATe9E z=I5F#twok1LTR=Ya+cV#vj$}4koHE5HlhV2Z55;i?6T%#<G8FgdZVb!Qb=1Sy8tuB zPLy3p46gUYGZB`bsWz~aq_>rF!?Fvhe!bF-(31En;e#Jo%+5kbgo;zA&Emd@F~obb zm^;!4a9$Dt{>UPj0ogVu)Vn07JX?}Z4kfN4!S9xY1N@B#xiAhdY=3%|0Zsx2;1ndD zB6rDfAit8dyb=i~51f?)%_JAtP0oEg;8cJ)JLq%{1c!p->D7=thor%C18#2}9#(ME zg9$-!uaHs*)59Ka9v_`XOLXsOpu~<0gv-ONd<4Xm;0yJjfo3>5MAD=%eQ*ScV7?_x z0La+Z28WWnMy|pG4W1gg0V6>lk%3DOF*`j#1avyH9D;sCPKH*%5el!F<U+onqtLzz z-3)q(?BR+@ey9{2`iHK^rDoX=16BiM$JBT@5&J2Or%o8e;B<v@nFkLZJOs}9d-14k zr4Rn(fy%Op8_FJkH4aJG=HSAt|2VACQ8sOqIje#WZY*yeRp@{}+mId8XliGfg$%Ax zWH$EcK%EzkisYEh^>C*H5&*x$lU?u4!I&KL3{Z!z=#;0x*ha?;GyGkG-{TI|fpfrB z@2(70^A3{jE`okql#>l(FJ5%~1dk=LfD5qSmgHm?=CghYu3o8h-mrbx&uPynAPKG5 z*(<GuY<xI(K>4_zdv+Y=HJ|zE&@EeszOx0d6GK=1;j`)8R(&({#%}nniQPj(tDUfs zm}{<f!k?(XsEN%}7laWG$7?ezJvp%9z}^ExzrShbC)<Xe<=nXFzUJdYJ)?&YX?SJp z%a(avdw=*$`(clKx9sD_n~%IqDhD6t{{8!#aXkk5BVIY?!Q`H6rag!vnhpinADs=b zDMcIw-o4;-5ABI0o7D4-UKM$VcMg5jKe*HON45?<@x-_PxVmEW(CmVIN*qll(fZ9E zdim;?mVUQ=`jBoc-Z1?d$<WZlkeNnD0{rp9^C7(d@~<cH^0RZ`lF0<I8I9@WzWfD4 z3;dTR7?pOx*3PC(__A=lp=q6=`LBlipEoR8Ygqi8Vac<GrSP}4<igXE=(-1gpnA(` zw8Zq9mz!Erh`#g-O!#8Fdlxk*mYstj4cfVf=^KHACtlRft?j6te{OMet%D%NuHCGf zQY^jL`AGebQaej8B3i?GJ5#-D4K2^{`nPsAoePN;K5J-t#?ZXRaQ|w<qE&{)PaBr3 zG%S6}u<S`g>Zm6SrcBWE@y?I7@A%-uf4;wKce_!<lvsA~`pveHW;l6-nZfGUW}<dn z-tq@tzt15>;`&X5R9HWA0oA*H08?y&zc7|Ew)UlZ_nTV2phcuzu%<6<-_K1=`%KMy zP4|ChTC~Ts_*2u8-6m5L<lo8b1T~6PpF&pd<jJ($!>m|VpZw+RT2j?DN!0%6KUfb* zruru0FD>9jlw|B-?eA)Fk_e~;jagXP<S<<hiLnjh*vs{w2vGf*?WlUL-$;4NZpmW$ z+h|PcsJ^BXsn`;<fn)K(+jnIH)xRTBKhw6nPfKs6ztw2k3~6`^$I{E~&b^uHZB0V2 z396Pls+VaV(I4OLN=S>mdv2rp*Yt|yg*M{y@h#VLUQ$hH>%!}ru5BK7|JX&>EFQCD zbn3{FsiUqo&4bKg=V>>RbqS&n{be?UwJT^j3rH+{5VaW*{(|We`O7#+gLWTb`c|F= z057Vqu|2sj_q0qQU(zlB4<e)5pPSsx)%RB|axHc)sWi=lY+;MB1IKCy-oE)8sr|Df z+o9DMZdw@dcmsqW(e7P@hFE%$7p%v`SE=3(eBwsyF?xBYmX72rs67Gcw7jHf#(UI- zq1z^^uN=NdUzv@xn9N2x%w{&y`Bi2knl4G1jo|X+POF&-KW9NPjw&ec$nNleDG=P5 z-QhQc2u^2r_&)#$hTlRUcWeX9e+IkLZl2kACc~7c%tqRrnT=<&_;c7DM`q0bJ9h8R z?&q;P)*a>WxA5fNgWc)VFtgFja8Gup?w;A$o8dm}ej&T}Wq0~woY_bp^O=ooI^Iac z;zKYX7!V8y1_T3w0l|P^KrkQ}5DW+g1OtKr!GK^uFd!HZ3<w4U1A+m;fM7r{AQ%t~ z2nGZLf&syRU_dY+7!V8y1_T3w0l|P^KrkQ}5DW+g1OtKr!GK^uFd!HZ3<w4U1A+m; zfM7r{AQ%t~2nGZLr^vwFUq9b^^z}7E4D8~+qjQRCL~_A^U_dY+7!V8y1_T3w0l|P^ zKrkQ}5DW+g1OtKr!GK^uFd!HZ3<w4U1A+m;fM7r{AQ%t~2nGZLf&syRU_dY+7!V8y z1_T3w0l|P^KrkQ}5DW+g1OtKr!GK^uFd!HZ3<w4U1A+m;fM7r{AQ%t~2nK#517^rS z8p1IM{U8`1@GF&nf&01I@59<jQlpX@zZr}WK#5@zggFqNgzyQ3t_I043_<|HlMwbp z$S_KVLI@QQ?u4)g!iNyrrAP(~gvk)@f$$Q9V-SX!B*Sb7FF}w}B||xcsi{)RS_o|r z-i7cHgwG)yhVUzdPVJ<WZV>uFus|3Hp&abC+pAo@Kqw%4?T+e7cfjs0jEM8q_~n4p z9t;KC-YT`E${CUaRhr~(h;pb?xxFE0z$<&~e%UcaS5W8GihM!4GZ63v>~ct*BT(lJ zxoezM7E=%-Q@tGULIpI^pfh9#i-T%y80YdkwNXBQ$nEn+rVQB~9$zprX~^#LJG~LP zL-v4EuGERLSNgo#-1Px>NH<O=vCps5F0`fW&|Wp_BC&;oE@uE*Sq{h!XqEU<dz2xH z2Ji`wvms6jZ1t)bB`{f>+#X+*+W{^U+XCRzqia42cs|+TgXWS$vezcr9dZQ#5(6LA z@HJQ)0QYhQoWM4eDA{XdkH-hj8+1eS!3QLQyyT-AI<8$hl46p!d;z^XwbJKwdV=;E zU#0!(QRUDw2&Hx|sGD>Z^cT4%ER&iHwRO@5#-KDGvY+@gMUqmiVJYUYO2cgt5jN>o zLqSALfi%ZplI9&ukxbI(KjG)$#`cob399IHN%t5k$o(UOL)ves^tfwW(k}+5$K&(| zrCG*+G~Wo-@dw?~V@9YEetF7RC#^Hqd!;vxCg~XHGD(fMfJP%xw8!8JVydH1=p2YP zNejs5V)EHa63W;epcUqVrUyaE<3zTGe6A;-Z6sZd^o%hGuJ9_P*-YZMlFwtiV1zWv zUPo!TrnPklrT2^x<kr>^a%Ta1tuJELH@#zmoX!x`E}RuQEmn4JN?ibYWH?%S8FJ0O z1JehIMB0!t-5Kyn+aP}J7L2z`X*Jb$XM@vG7joLGp`(Jahf}~fNosGGwx>9xr%kRp zuU&fIR3n3Hd}Q)hL#x@PJ*L`#4EOz}I-=Yz9W()rpMirIt7UJM(_ZfjOu<11hgE5X z$>$9Pd>*^>E+pGSNLaH0EA-MH1AHDnW`fTnx23}8>`zGfV8WK*VNdxM)B34L4C`+= zZ1|+lg0vm!t?6rqFEFg^KRbO6Pj?GEpfaH3EoI=UkZ&&F9`NPFrSN<K+=Ozd@8&v+ zA4UufF11noB*PsT{u{#`8U6>uofzKFaA$^(Gn~e720U6&e>%gx7|vk0nBguAU&(ML z!=n**fM+ek6%3!xaEM_u!*?*8&G3AN^BHbpxRhb8Zy3W%nY>&hzntNfOg@6)zcD<T z;jIjh)8O$8?`86d3?E^5GQ;hGE39uN!`!}V4fZgMZCZ&R?CK8rD>;CBjRw~Pc1h6I z_iuJ#GR$g$5AH8(LEcA&QQilFzZgE?LOfC<$NYo*RQ3|f&oR$mPI4>xx%?GK-&f_2 zTM?_;8}rAj-(E<suJ3WcSiioLBJ0=4?2|z^){kRu-w4)$5M%wgd@7_zd$_-$Jf1ys zAg$V-CcwIU2PhZcfSmBnFt?}lB1#yyXA6^e=m_>e9i=wfUo#Pd4l~2O8Rq`qpJDF* z`3!UaPgI{`NUyHXC4jL$Rn3w7npV;JEC)H(hhtu!BeXh@lk&^ue_{3vWA-PqXT8Rr zP0SvxKObZE{Js;lhhuKf`T>#ra(OAFM}9jndqyzK$G^!8^YJglFdzTsFwDol2N>q# z-x`Mb__v8+KK|`M4DHFspG15dfb{D2`5rK~kCu-$_|5@UzN$0jgJaIeVNLsR`4~u# z_VDp9kv$GbueQer812c@^vBj+)SkDPJsfj;CU%W%A1>d=?BU~YB71(&*dwLs+jHU@ zYR?3aqyKWu?a7338ZX2+K5%(0q(^_@&#QXGz!!fWVETCD#au|Q=Bo*?9$#rc(y|YM z9PQ&6_4k3WLt`J8r>0YT`12@{J?B7rwLQH6>)G@932M&_kfS{ub9-7f_Hg+^W)FWJ zCbH)-jXlpWd$i;G>=Zg4re#ojIOg`O=>;Ugg&5mUEyw3qVPBg6LYALno__(!t>owO zL69EXn?KJ-GCaAB>cjl;+FyqBYQCoeM!vQEcN4R38_1Dwj=6m+&!Pq*#`<yj0cH=M zPy7PdOp}sNOdrp__FW?F>k1g{(~kfCLA1UekfVJZqki<y*&2V~@(`p){`mYRkv+F- z>{$@Qo+Fn~d%j}!aLn!LfM4Lk8}h~FhnYRe&aaG_k^H0sMt-<|Rsl9+LF}5T!Oa@H zQiC^X@W&eby#}YB5m{bO4KCB*F&gaD;2Sj<`BMAKe2}aC_ih-s)p!wL<l~L=o_8Sy zm>>Hae;?Tga`YFDQI7sn4goL3*gk5xx_q3rbpbh_ADzcAAD?Xu^YQQshWU6j88OTu z{B6`<>KSfe_#TEE8Gel6ISg-Pcmc!jGrXALgABJa+zxET@*ZcnJHu-jwlTb(;VT(# zV|Wt7n;8x<yp`cOh#~O#%tH+G=jRIy^YP)I4D<2uAj5oo!x!K#25>UCKqQ?*urfc+ zMGW=X!^$saaz1|;#4w-lU&1i&|5q~1``;Lb`TTGKU?>-hCYM)(@x@myU#`KcH24h- z-mAesX>eCqr%~&NW!MNjLW5z6HbNfM;5#+AMT4Kw;7uC*i3T6m;Etxq^357tqQTc{ z@H7oxrornpc&7&ctij!2JwVN00mWc7y~;Frm<C^_!B=VU7!4k;!IL!Dp}|u$*ssB` zoR0Q~?Y&WhzlLfYg76K5!w|lO@EwHjA>cdD4-jDLCLM)<`w2fmI1b@w2rz|`eu3~S zgp&{mQBar>3xr^TkO~3EjP?-lnbi?OCkUM(q(MlBkO2Yvd?ti5Ai&mwbS8vu5YB>d zHiZ8xcA8;1m#w3QSMS)eDhS9jX*ixPnzCr^(kYA9uAf45<hm?fI*p2^E2&Y@4zey9 z9ZlCsqvGirY*aK`l8uUoWkkA|8ZD_qULq}LYq?SJY-u+tz7F>RG&<^VNj54L)`KCA zF2jb|W~;J9ptmp!0=oJd&cPOAQ4qNt%hq6{1!PUuZjWAw#odail~h>9XDg|8vM>uP z_qb%wmSyd9Ro2cIWyJ=U*x>TQc9?eii*FMt`=3gh|A!k~_Ny<xVx-+p1_fWeyV3=l zW^B*XQew5_OR$L<2v{K24?DpjJKIpTSo19<QE@g)aa3H1r6?+{%mSOuk@{=}mV&6b zqUdr93M^J_x`HCK)9=PlcMuGA)CGd>sZL8irq#q~(^!kq7#S-G2DmVXi_|+JU#X>3 zqp38iw8E&;3Zv4LMx}xBqfEkldDp^N9n7n7!%*n-<dr+8Iz6z(9LO6!e8gCrxyl_H z5S;3zC@<s;hJtxxouT3W;5hhxsT_0$FAMl;#ya33<EYNN{<7=yVDk*ND=n@vP(D=- z<hgtv*z^kKIV$q}r8PxGm4!iT;q;;gm&@TQ4p~d`OC9-eReAHq<>fgXu0LE>U1GEN zq_E4-@i)?2B3MKIjnP9nHQ_&AKb?ZW*yvbrc8!(!l95h2PLpL!Iam{Ovcr2>gk`aO z(?}n3r^#BT9IOpLOCzkx2k$h}>A`8TfGH<@#h8{vjvoK}MxzMEaP%<3oaHppg_?mY zIw^`Mm^GEs@BibSj&9Rq@i>FjohU|Fgu|wh%;ZlO8F;jU+N1mWp!1l(ft%Q5JB>Dx zk(0jCM7_I2Foc6Yv8LcOT2r7f6!pp!Z4I5kpC)sda^gX!_Tms>XpASTv1ZzEc<R~m z2hi7((^Uy_5QUN&<FzQlS}U0vgkNZk<R#`b*+z8(zx+)p!V)?}z%1c3m5#4v$P0Tt zi!PhI<({U(sczz0kuTF?8Gy|ABaW=WGC&?o=JMb)oa`Yp#2VS*heM*lyox%vr!o&; z4D1yZL~e%}q@6_LX(u>eND|T*_T2-o@b)T)!&;a}KQSd!CCdRvb#dVU*$+nw2UN;| zdbihBkmqq%K+?Rl{{cIooB_}aZI*mqszgh-SsZE}Ad;Dnmc_6;LXkzz1_sExl*N~$ zezStLz};zWu+q1UuKwX7pwtcGU=`Z5QG%m#k`ew%(w`uVM<HCjC1bANY{Z@HBrZcE zPh`M{u~g{G@_#V)iz)xQK4?B@3*41EaL;tZ!!sv*y7~CDUQ0Ls`1H*ypB?Xe{VP|S zY5TTqQ;eg=r=Fa#cE<GEKH7NK<#TgaB$7|R_mf38^v-y==?+)PxT1~~-|n2@?!SG< z%(?IUu*#LTv8wyVKRw_7wHcYVM>kJ6w(bm1kNK;IRlPBL>&AskT07r1_L}!TKTvwk zpq01ZSN3tojGeD7z191`x-pO5wNr9?Gg9B~_j2ojg5%%+{b%#=^>;4W{)OSGgD38I z;KN?`4!ice?%8(R5X-dQ%j(}P^Bmf>W9;P*KA!Gf`Ohmx|Le<p&fovzn5SCirayGf zrgih~bASIt>vJ=$mYZh1yX)S=YxB*WkA86N3p+;{S{6L8<%*)*ksrOdegDZ7D`s>` z&)Ib(<L!rE9rr;;Yu65s96Y*WKt}p6pH7;Q^X_wPw=R5duj}?>?H~DUmHmptJtp>; z@W;<43>^RBeJTF53H4w7@KeaTth%Vy`geKI+Ps@*b~*Cw?hRG28a844=gr<lhdNG~ zKjPdSkCb#D^oILn&*sh(I~7{`{Bh%u1C62irk848ZF{`-FJ0bSUcd65-lobt?~Sv@ zoRQyW%v(bqIrPE7#~%LW&HKLF`t{wr&)W6KL*L)sZo&A!Uwu`@TaRzfoV09mw~^<6 zH?I7O1Lk$h=Zv{w)8of4?P{4@@N3qM_AaeYHP3&0;LI`0hZmT#Q?C2czv<1Dw#)u} zu=)2__srQgXX>Ez{u@e+1#9+}`VW0O$@JiL*Z#C*{DW1SOgoPMbmEr@ukV`j{_om8 zyU5_ED(tbe=j*RLQ8MCBz4x`Oy?K4pUmt$7d-*Ek#(TfM>Y-I5-nsPMRbOsc+qeHw z@2%JE`*B|JxV8yfkB%{YomSp$?>#?W?z`!(X5ZN#_j&p4)s;`r&HCM;<?`72tDZ0W z=D=%>DQ!J(nDJo0&t^TjYluAEHqU;CJn6pG?e2X2?Ki)CY`weBv8k&acWlXd^H7tq zN5joeEW2WV%7&qTysN=r-Su<ris`NYYBt~S(&V=c9Tq?F#f*DrocZh<^Li}0$g=9+ z@6~Ty>$&)r<Ev&kuD*6+&#snzYhV29Z8_7jwl&?->F1(Z*Ik^*zkjQI*~0f%&z!qF zz3i<^{<i4e?sdlpzqj<>hh}Yf_8z(KpyH`}eyjdO`G;Tsq4fB^bzQG~@|OG33@a)& zeDm?>Z-*bOzUiXpUiom#SJ#{vapl>i7x>1$w)(~5FXv3&m{s{fm;Qb27t}86eqrO= t%R1Ti4H$RwF-KWWb85RzKi}fG@X5?~zuap)^6v#*Z+u~2@$GE1`7f*GPa*&S literal 0 HcmV?d00001 diff --git a/tests/SetEpsTestCasesFromScratch/WFW/eels.f90 b/tests/SetEpsTestCasesFromScratch/WFW/eels.f90 new file mode 100644 index 0000000..3eb1277 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/eels.f90 @@ -0,0 +1,160 @@ +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. * +! * * +! ****************************************************************** + + use doeels_mod + + 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(:, :) + logical :: debug + integer :: ioStatus + +! integer, parameter :: name_length = 10 +! *** read spectrometer parameters + + call change_working_dir() + open(unit = 11, file = 'eelsin') +! impact energy (eV) + read(11, *) e0 +! incidence angle (°) + read(11, *) theta +! angular apertures of the elliptic detector (°) + read(11, *) phia + read(11, *) phib +! energy-loss interval and step size (cm**-1) + read(11, *) wmin + read(11, *) wmax + read(11, *) dw +! comment lines + read(11, '(a72)') (comment(k), k = 1, 2) + + write(*,*) 'program eels (September 2020)' + write(*,'(a, f6.2, a, f5.1, a, f5.2, a, f5.2, a)') & + ' e0 = ', e0, ' eV, theta = ', theta, '°, phia = ', & + phia, '°, phib = ', phib, '°' + write(*,'(a, g11.4, a, g11.4, a, g11.4, a)') & + ' energy losses from', wmin, ' to', wmax, ', step = ', dw, ' cm**-1' + write(*,*) comment(1) + write(*,*) comment(2) + + if (phia <= 0.0d0 .or. phib <= 0.0d0) then + stop '*** wrong input ***' + endif + if (e0 <= 0.0d0 .or. theta + phia >= 90.0d0) then + stop '*** bad input ***' + endif + +! *** read target specifications + + read(11, *) layers, nper, mode + user = layers == 0 + if (.not. user) then + neps = layers + if (nper == -1) then + neps = layers + 1 + nper = 1 + write(*,*) 'the substrate is a anisotropic uniaxial material' + endif + if (layers < 0 .or. nper < 1 .or. nper > layers) then + stop '*** invalid target specifications ***' + endif + write(6, 101) layers, nper + jos = 0 + allocate (name(neps), thick(neps), epsinf(neps), nos(neps)) + do l = 1, neps + if (l <= layers) then + read(11, 102) name(l), thick(l) + endif + read(11, *) epsinf(l), nos(l) + write(6, 103) + 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) + else + do j = 1, nos(l) + if (.not. allocated(osc)) then + allocate(osc(3, nos(l))) + else if (j == 1) then + ! call extend3(osc, nos(j)) + old_size_1 = size(osc, 1) + old_size_2 = size(osc, 2) + allocate(tmp_osc(old_size_1, old_size_2 + nos(l))) + tmp_osc(1:old_size_1, 1:old_size_2) = osc + deallocate(osc) + call move_alloc(tmp_osc, osc) + endif + jos = jos + 1 + read(11, *) (osc(k, jos), k = 1, 3) + if ((j == nos(l) / 2 + 1) .and. (nos(l) > 1)) then + write(6, 106) + write(6, 107) + endif + if (j == 1) then + if (l <= layers) then + 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 + enddo + endif + enddo + write(*,*) + read(11, 102, IOSTAT = ioStatus) contrl + if (ioStatus /= 0) then + contrl = '' + endif + endif + + close (unit = 11) + + nw = 1 + int((wmax - wmin) / dw) + allocate (wn_array(nw), f(nw)) + + debug = .false. + 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, debug, f, size(wn_array)) + + open (unit = 12, file = 'eelsou') + write (12, 207) e0, theta, phia, phib + write (12, '(a72)') comment(1) + do i = 1, nw + write (12, '(2e15.7)') wn_array(i), f(i) + enddo + close (unit = 12) + + 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) +end program eels diff --git a/tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps001.log b/tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps001.log new file mode 100644 index 0000000..a1dd64d --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps001.log @@ -0,0 +1,331 @@ + 1 1 +No-layers +MnO 4.950000 1 + 269.0000 16.00000 0.5000000E-01 + + 50.00000 4.950000 0.000000 + 52.00000 4.950000 0.000000 + 54.00000 4.950000 0.000000 + 56.00000 4.950000 0.000000 + 58.00000 4.950000 0.000000 + 60.00000 4.950000 0.000000 + 62.00000 4.950000 0.000000 + 64.00000 4.950000 0.000000 + 66.00000 4.950000 0.000000 + 68.00000 4.950000 0.000000 + 70.00000 4.950000 0.000000 + 72.00000 4.950000 0.000000 + 74.00000 4.950000 0.000000 + 76.00000 4.950000 0.000000 + 78.00000 4.950000 0.000000 + 80.00000 4.950000 0.000000 + 82.00000 4.950000 0.000000 + 84.00000 4.950000 0.000000 + 86.00000 4.950000 0.000000 + 88.00000 4.950000 0.000000 + 90.00000 4.950000 0.000000 + 92.00000 4.950000 0.000000 + 94.00000 4.950000 0.000000 + 96.00000 4.950000 0.000000 + 98.00000 4.950000 0.000000 + 100.0000 4.950000 0.000000 + 102.0000 4.950000 0.000000 + 104.0000 4.950000 0.000000 + 106.0000 4.950000 0.000000 + 108.0000 4.950000 0.000000 + 110.0000 4.950000 0.000000 + 112.0000 4.950000 0.000000 + 114.0000 4.950000 0.000000 + 116.0000 4.950000 0.000000 + 118.0000 4.950000 0.000000 + 120.0000 4.950000 0.000000 + 122.0000 4.950000 0.000000 + 124.0000 4.950000 0.000000 + 126.0000 4.950000 0.000000 + 128.0000 4.950000 0.000000 + 130.0000 4.950000 0.000000 + 132.0000 4.950000 0.000000 + 134.0000 4.950000 0.000000 + 136.0000 4.950000 0.000000 + 138.0000 4.950000 0.000000 + 140.0000 4.950000 0.000000 + 142.0000 4.950000 0.000000 + 144.0000 4.950000 0.000000 + 146.0000 4.950000 0.000000 + 148.0000 4.950000 0.000000 + 150.0000 4.950000 0.000000 + 152.0000 4.950000 0.000000 + 154.0000 4.950000 0.000000 + 156.0000 4.950000 0.000000 + 158.0000 4.950000 0.000000 + 160.0000 4.950000 0.000000 + 162.0000 4.950000 0.000000 + 164.0000 4.950000 0.000000 + 166.0000 4.950000 0.000000 + 168.0000 4.950000 0.000000 + 170.0000 4.950000 0.000000 + 172.0000 4.950000 0.000000 + 174.0000 4.950000 0.000000 + 176.0000 4.950000 0.000000 + 178.0000 4.950000 0.000000 + 180.0000 4.950000 0.000000 + 182.0000 4.950000 0.000000 + 184.0000 4.950000 0.000000 + 186.0000 4.950000 0.000000 + 188.0000 4.950000 0.000000 + 190.0000 4.950000 0.000000 + 192.0000 4.950000 0.000000 + 194.0000 4.950000 0.000000 + 196.0000 4.950000 0.000000 + 198.0000 4.950000 0.000000 + 200.0000 4.950000 0.000000 + 202.0000 4.950000 0.000000 + 204.0000 4.950000 0.000000 + 206.0000 4.950000 0.000000 + 208.0000 4.950000 0.000000 + 210.0000 4.950000 0.000000 + 212.0000 4.950000 0.000000 + 214.0000 4.950000 0.000000 + 216.0000 4.950000 0.000000 + 218.0000 4.950000 0.000000 + 220.0000 4.950000 0.000000 + 222.0000 4.950000 0.000000 + 224.0000 4.950000 0.000000 + 226.0000 4.950000 0.000000 + 228.0000 4.950000 0.000000 + 230.0000 4.950000 0.000000 + 232.0000 4.950000 0.000000 + 234.0000 4.950000 0.000000 + 236.0000 4.950000 0.000000 + 238.0000 4.950000 0.000000 + 240.0000 4.950000 0.000000 + 242.0000 4.950000 0.000000 + 244.0000 4.950000 0.000000 + 246.0000 4.950000 0.000000 + 248.0000 4.950000 0.000000 + 250.0000 4.950000 0.000000 + 252.0000 4.950000 0.000000 + 254.0000 4.950000 0.000000 + 256.0000 4.950000 0.000000 + 258.0000 4.950000 0.000000 + 260.0000 4.950000 0.000000 + 262.0000 4.950000 0.000000 + 264.0000 4.950000 0.000000 + 266.0000 4.950000 0.000000 + 268.0000 4.950000 0.000000 + 270.0000 4.950000 0.000000 + 272.0000 4.950000 0.000000 + 274.0000 4.950000 0.000000 + 276.0000 4.950000 0.000000 + 278.0000 4.950000 0.000000 + 280.0000 4.950000 0.000000 + 282.0000 4.950000 0.000000 + 284.0000 4.950000 0.000000 + 286.0000 4.950000 0.000000 + 288.0000 4.950000 0.000000 + 290.0000 4.950000 0.000000 + 292.0000 4.950000 0.000000 + 294.0000 4.950000 0.000000 + 296.0000 4.950000 0.000000 + 298.0000 4.950000 0.000000 + 300.0000 4.950000 0.000000 + 302.0000 4.950000 0.000000 + 304.0000 4.950000 0.000000 + 306.0000 4.950000 0.000000 + 308.0000 4.950000 0.000000 + 310.0000 4.950000 0.000000 + 312.0000 4.950000 0.000000 + 314.0000 4.950000 0.000000 + 316.0000 4.950000 0.000000 + 318.0000 4.950000 0.000000 + 320.0000 4.950000 0.000000 + 322.0000 4.950000 0.000000 + 324.0000 4.950000 0.000000 + 326.0000 4.950000 0.000000 + 328.0000 4.950000 0.000000 + 330.0000 4.950000 0.000000 + 332.0000 4.950000 0.000000 + 334.0000 4.950000 0.000000 + 336.0000 4.950000 0.000000 + 338.0000 4.950000 0.000000 + 340.0000 4.950000 0.000000 + 342.0000 4.950000 0.000000 + 344.0000 4.950000 0.000000 + 346.0000 4.950000 0.000000 + 348.0000 4.950000 0.000000 + 350.0000 4.950000 0.000000 + 352.0000 4.950000 0.000000 + 354.0000 4.950000 0.000000 + 356.0000 4.950000 0.000000 + 358.0000 4.950000 0.000000 + 360.0000 4.950000 0.000000 + 362.0000 4.950000 0.000000 + 364.0000 4.950000 0.000000 + 366.0000 4.950000 0.000000 + 368.0000 4.950000 0.000000 + 370.0000 4.950000 0.000000 + 372.0000 4.950000 0.000000 + 374.0000 4.950000 0.000000 + 376.0000 4.950000 0.000000 + 378.0000 4.950000 0.000000 + 380.0000 4.950000 0.000000 + 382.0000 4.950000 0.000000 + 384.0000 4.950000 0.000000 + 386.0000 4.950000 0.000000 + 388.0000 4.950000 0.000000 + 390.0000 4.950000 0.000000 + 392.0000 4.950000 0.000000 + 394.0000 4.950000 0.000000 + 396.0000 4.950000 0.000000 + 398.0000 4.950000 0.000000 + 400.0000 4.950000 0.000000 + 402.0000 4.950000 0.000000 + 404.0000 4.950000 0.000000 + 406.0000 4.950000 0.000000 + 408.0000 4.950000 0.000000 + 410.0000 4.950000 0.000000 + 412.0000 4.950000 0.000000 + 414.0000 4.950000 0.000000 + 416.0000 4.950000 0.000000 + 418.0000 4.950000 0.000000 + 420.0000 4.950000 0.000000 + 422.0000 4.950000 0.000000 + 424.0000 4.950000 0.000000 + 426.0000 4.950000 0.000000 + 428.0000 4.950000 0.000000 + 430.0000 4.950000 0.000000 + 432.0000 4.950000 0.000000 + 434.0000 4.950000 0.000000 + 436.0000 4.950000 0.000000 + 438.0000 4.950000 0.000000 + 440.0000 4.950000 0.000000 + 442.0000 4.950000 0.000000 + 444.0000 4.950000 0.000000 + 446.0000 4.950000 0.000000 + 448.0000 4.950000 0.000000 + 450.0000 4.950000 0.000000 + 452.0000 4.950000 0.000000 + 454.0000 4.950000 0.000000 + 456.0000 4.950000 0.000000 + 458.0000 4.950000 0.000000 + 460.0000 4.950000 0.000000 + 462.0000 4.950000 0.000000 + 464.0000 4.950000 0.000000 + 466.0000 4.950000 0.000000 + 468.0000 4.950000 0.000000 + 470.0000 4.950000 0.000000 + 472.0000 4.950000 0.000000 + 474.0000 4.950000 0.000000 + 476.0000 4.950000 0.000000 + 478.0000 4.950000 0.000000 + 480.0000 4.950000 0.000000 + 482.0000 4.950000 0.000000 + 484.0000 4.950000 0.000000 + 486.0000 4.950000 0.000000 + 488.0000 4.950000 0.000000 + 490.0000 4.950000 0.000000 + 492.0000 4.950000 0.000000 + 494.0000 4.950000 0.000000 + 496.0000 4.950000 0.000000 + 498.0000 4.950000 0.000000 + 500.0000 4.950000 0.000000 + 502.0000 4.950000 0.000000 + 504.0000 4.950000 0.000000 + 506.0000 4.950000 0.000000 + 508.0000 4.950000 0.000000 + 510.0000 4.950000 0.000000 + 512.0000 4.950000 0.000000 + 514.0000 4.950000 0.000000 + 516.0000 4.950000 0.000000 + 518.0000 4.950000 0.000000 + 520.0000 4.950000 0.000000 + 522.0000 4.950000 0.000000 + 524.0000 4.950000 0.000000 + 526.0000 4.950000 0.000000 + 528.0000 4.950000 0.000000 + 530.0000 4.950000 0.000000 + 532.0000 4.950000 0.000000 + 534.0000 4.950000 0.000000 + 536.0000 4.950000 0.000000 + 538.0000 4.950000 0.000000 + 540.0000 4.950000 0.000000 + 542.0000 4.950000 0.000000 + 544.0000 4.950000 0.000000 + 546.0000 4.950000 0.000000 + 548.0000 4.950000 0.000000 + 550.0000 4.950000 0.000000 + 552.0000 4.950000 0.000000 + 554.0000 4.950000 0.000000 + 556.0000 4.950000 0.000000 + 558.0000 4.950000 0.000000 + 560.0000 4.950000 0.000000 + 562.0000 4.950000 0.000000 + 564.0000 4.950000 0.000000 + 566.0000 4.950000 0.000000 + 568.0000 4.950000 0.000000 + 570.0000 4.950000 0.000000 + 572.0000 4.950000 0.000000 + 574.0000 4.950000 0.000000 + 576.0000 4.950000 0.000000 + 578.0000 4.950000 0.000000 + 580.0000 4.950000 0.000000 + 582.0000 4.950000 0.000000 + 584.0000 4.950000 0.000000 + 586.0000 4.950000 0.000000 + 588.0000 4.950000 0.000000 + 590.0000 4.950000 0.000000 + 592.0000 4.950000 0.000000 + 594.0000 4.950000 0.000000 + 596.0000 4.950000 0.000000 + 598.0000 4.950000 0.000000 + 600.0000 4.950000 0.000000 + 602.0000 4.950000 0.000000 + 604.0000 4.950000 0.000000 + 606.0000 4.950000 0.000000 + 608.0000 4.950000 0.000000 + 610.0000 4.950000 0.000000 + 612.0000 4.950000 0.000000 + 614.0000 4.950000 0.000000 + 616.0000 4.950000 0.000000 + 618.0000 4.950000 0.000000 + 620.0000 4.950000 0.000000 + 622.0000 4.950000 0.000000 + 624.0000 4.950000 0.000000 + 626.0000 4.950000 0.000000 + 628.0000 4.950000 0.000000 + 630.0000 4.950000 0.000000 + 632.0000 4.950000 0.000000 + 634.0000 4.950000 0.000000 + 636.0000 4.950000 0.000000 + 638.0000 4.950000 0.000000 + 640.0000 4.950000 0.000000 + 642.0000 4.950000 0.000000 + 644.0000 4.950000 0.000000 + 646.0000 4.950000 0.000000 + 648.0000 4.950000 0.000000 + 650.0000 4.950000 0.000000 + 652.0000 4.950000 0.000000 + 654.0000 4.950000 0.000000 + 656.0000 4.950000 0.000000 + 658.0000 4.950000 0.000000 + 660.0000 4.950000 0.000000 + 662.0000 4.950000 0.000000 + 664.0000 4.950000 0.000000 + 666.0000 4.950000 0.000000 + 668.0000 4.950000 0.000000 + 670.0000 4.950000 0.000000 + 672.0000 4.950000 0.000000 + 674.0000 4.950000 0.000000 + 676.0000 4.950000 0.000000 + 678.0000 4.950000 0.000000 + 680.0000 4.950000 0.000000 + 682.0000 4.950000 0.000000 + 684.0000 4.950000 0.000000 + 686.0000 4.950000 0.000000 + 688.0000 4.950000 0.000000 + 690.0000 4.950000 0.000000 + 692.0000 4.950000 0.000000 + 694.0000 4.950000 0.000000 + 696.0000 4.950000 0.000000 + 698.0000 4.950000 0.000000 + 700.0000 4.950000 0.000000 diff --git a/tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps004.log b/tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps004.log new file mode 100644 index 0000000..d8ef70b --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps004.log @@ -0,0 +1,331 @@ + 1 1 +No-layers +Pt 8.900000 1 + 185541.0 -1.000000 0.1200000 + + 50.00000 8.900000 0.000000 + 52.00000 8.900000 0.000000 + 54.00000 8.900000 0.000000 + 56.00000 8.900000 0.000000 + 58.00000 8.900000 0.000000 + 60.00000 8.900000 0.000000 + 62.00000 8.900000 0.000000 + 64.00000 8.900000 0.000000 + 66.00000 8.900000 0.000000 + 68.00000 8.900000 0.000000 + 70.00000 8.900000 0.000000 + 72.00000 8.900000 0.000000 + 74.00000 8.900000 0.000000 + 76.00000 8.900000 0.000000 + 78.00000 8.900000 0.000000 + 80.00000 8.900000 0.000000 + 82.00000 8.900000 0.000000 + 84.00000 8.900000 0.000000 + 86.00000 8.900000 0.000000 + 88.00000 8.900000 0.000000 + 90.00000 8.900000 0.000000 + 92.00000 8.900000 0.000000 + 94.00000 8.900000 0.000000 + 96.00000 8.900000 0.000000 + 98.00000 8.900000 0.000000 + 100.0000 8.900000 0.000000 + 102.0000 8.900000 0.000000 + 104.0000 8.900000 0.000000 + 106.0000 8.900000 0.000000 + 108.0000 8.900000 0.000000 + 110.0000 8.900000 0.000000 + 112.0000 8.900000 0.000000 + 114.0000 8.900000 0.000000 + 116.0000 8.900000 0.000000 + 118.0000 8.900000 0.000000 + 120.0000 8.900000 0.000000 + 122.0000 8.900000 0.000000 + 124.0000 8.900000 0.000000 + 126.0000 8.900000 0.000000 + 128.0000 8.900000 0.000000 + 130.0000 8.900000 0.000000 + 132.0000 8.900000 0.000000 + 134.0000 8.900000 0.000000 + 136.0000 8.900000 0.000000 + 138.0000 8.900000 0.000000 + 140.0000 8.900000 0.000000 + 142.0000 8.900000 0.000000 + 144.0000 8.900000 0.000000 + 146.0000 8.900000 0.000000 + 148.0000 8.900000 0.000000 + 150.0000 8.900000 0.000000 + 152.0000 8.900000 0.000000 + 154.0000 8.900000 0.000000 + 156.0000 8.900000 0.000000 + 158.0000 8.900000 0.000000 + 160.0000 8.900000 0.000000 + 162.0000 8.900000 0.000000 + 164.0000 8.900000 0.000000 + 166.0000 8.900000 0.000000 + 168.0000 8.900000 0.000000 + 170.0000 8.900000 0.000000 + 172.0000 8.900000 0.000000 + 174.0000 8.900000 0.000000 + 176.0000 8.900000 0.000000 + 178.0000 8.900000 0.000000 + 180.0000 8.900000 0.000000 + 182.0000 8.900000 0.000000 + 184.0000 8.900000 0.000000 + 186.0000 8.900000 0.000000 + 188.0000 8.900000 0.000000 + 190.0000 8.900000 0.000000 + 192.0000 8.900000 0.000000 + 194.0000 8.900000 0.000000 + 196.0000 8.900000 0.000000 + 198.0000 8.900000 0.000000 + 200.0000 8.900000 0.000000 + 202.0000 8.900000 0.000000 + 204.0000 8.900000 0.000000 + 206.0000 8.900000 0.000000 + 208.0000 8.900000 0.000000 + 210.0000 8.900000 0.000000 + 212.0000 8.900000 0.000000 + 214.0000 8.900000 0.000000 + 216.0000 8.900000 0.000000 + 218.0000 8.900000 0.000000 + 220.0000 8.900000 0.000000 + 222.0000 8.900000 0.000000 + 224.0000 8.900000 0.000000 + 226.0000 8.900000 0.000000 + 228.0000 8.900000 0.000000 + 230.0000 8.900000 0.000000 + 232.0000 8.900000 0.000000 + 234.0000 8.900000 0.000000 + 236.0000 8.900000 0.000000 + 238.0000 8.900000 0.000000 + 240.0000 8.900000 0.000000 + 242.0000 8.900000 0.000000 + 244.0000 8.900000 0.000000 + 246.0000 8.900000 0.000000 + 248.0000 8.900000 0.000000 + 250.0000 8.900000 0.000000 + 252.0000 8.900000 0.000000 + 254.0000 8.900000 0.000000 + 256.0000 8.900000 0.000000 + 258.0000 8.900000 0.000000 + 260.0000 8.900000 0.000000 + 262.0000 8.900000 0.000000 + 264.0000 8.900000 0.000000 + 266.0000 8.900000 0.000000 + 268.0000 8.900000 0.000000 + 270.0000 8.900000 0.000000 + 272.0000 8.900000 0.000000 + 274.0000 8.900000 0.000000 + 276.0000 8.900000 0.000000 + 278.0000 8.900000 0.000000 + 280.0000 8.900000 0.000000 + 282.0000 8.900000 0.000000 + 284.0000 8.900000 0.000000 + 286.0000 8.900000 0.000000 + 288.0000 8.900000 0.000000 + 290.0000 8.900000 0.000000 + 292.0000 8.900000 0.000000 + 294.0000 8.900000 0.000000 + 296.0000 8.900000 0.000000 + 298.0000 8.900000 0.000000 + 300.0000 8.900000 0.000000 + 302.0000 8.900000 0.000000 + 304.0000 8.900000 0.000000 + 306.0000 8.900000 0.000000 + 308.0000 8.900000 0.000000 + 310.0000 8.900000 0.000000 + 312.0000 8.900000 0.000000 + 314.0000 8.900000 0.000000 + 316.0000 8.900000 0.000000 + 318.0000 8.900000 0.000000 + 320.0000 8.900000 0.000000 + 322.0000 8.900000 0.000000 + 324.0000 8.900000 0.000000 + 326.0000 8.900000 0.000000 + 328.0000 8.900000 0.000000 + 330.0000 8.900000 0.000000 + 332.0000 8.900000 0.000000 + 334.0000 8.900000 0.000000 + 336.0000 8.900000 0.000000 + 338.0000 8.900000 0.000000 + 340.0000 8.900000 0.000000 + 342.0000 8.900000 0.000000 + 344.0000 8.900000 0.000000 + 346.0000 8.900000 0.000000 + 348.0000 8.900000 0.000000 + 350.0000 8.900000 0.000000 + 352.0000 8.900000 0.000000 + 354.0000 8.900000 0.000000 + 356.0000 8.900000 0.000000 + 358.0000 8.900000 0.000000 + 360.0000 8.900000 0.000000 + 362.0000 8.900000 0.000000 + 364.0000 8.900000 0.000000 + 366.0000 8.900000 0.000000 + 368.0000 8.900000 0.000000 + 370.0000 8.900000 0.000000 + 372.0000 8.900000 0.000000 + 374.0000 8.900000 0.000000 + 376.0000 8.900000 0.000000 + 378.0000 8.900000 0.000000 + 380.0000 8.900000 0.000000 + 382.0000 8.900000 0.000000 + 384.0000 8.900000 0.000000 + 386.0000 8.900000 0.000000 + 388.0000 8.900000 0.000000 + 390.0000 8.900000 0.000000 + 392.0000 8.900000 0.000000 + 394.0000 8.900000 0.000000 + 396.0000 8.900000 0.000000 + 398.0000 8.900000 0.000000 + 400.0000 8.900000 0.000000 + 402.0000 8.900000 0.000000 + 404.0000 8.900000 0.000000 + 406.0000 8.900000 0.000000 + 408.0000 8.900000 0.000000 + 410.0000 8.900000 0.000000 + 412.0000 8.900000 0.000000 + 414.0000 8.900000 0.000000 + 416.0000 8.900000 0.000000 + 418.0000 8.900000 0.000000 + 420.0000 8.900000 0.000000 + 422.0000 8.900000 0.000000 + 424.0000 8.900000 0.000000 + 426.0000 8.900000 0.000000 + 428.0000 8.900000 0.000000 + 430.0000 8.900000 0.000000 + 432.0000 8.900000 0.000000 + 434.0000 8.900000 0.000000 + 436.0000 8.900000 0.000000 + 438.0000 8.900000 0.000000 + 440.0000 8.900000 0.000000 + 442.0000 8.900000 0.000000 + 444.0000 8.900000 0.000000 + 446.0000 8.900000 0.000000 + 448.0000 8.900000 0.000000 + 450.0000 8.900000 0.000000 + 452.0000 8.900000 0.000000 + 454.0000 8.900000 0.000000 + 456.0000 8.900000 0.000000 + 458.0000 8.900000 0.000000 + 460.0000 8.900000 0.000000 + 462.0000 8.900000 0.000000 + 464.0000 8.900000 0.000000 + 466.0000 8.900000 0.000000 + 468.0000 8.900000 0.000000 + 470.0000 8.900000 0.000000 + 472.0000 8.900000 0.000000 + 474.0000 8.900000 0.000000 + 476.0000 8.900000 0.000000 + 478.0000 8.900000 0.000000 + 480.0000 8.900000 0.000000 + 482.0000 8.900000 0.000000 + 484.0000 8.900000 0.000000 + 486.0000 8.900000 0.000000 + 488.0000 8.900000 0.000000 + 490.0000 8.900000 0.000000 + 492.0000 8.900000 0.000000 + 494.0000 8.900000 0.000000 + 496.0000 8.900000 0.000000 + 498.0000 8.900000 0.000000 + 500.0000 8.900000 0.000000 + 502.0000 8.900000 0.000000 + 504.0000 8.900000 0.000000 + 506.0000 8.900000 0.000000 + 508.0000 8.900000 0.000000 + 510.0000 8.900000 0.000000 + 512.0000 8.900000 0.000000 + 514.0000 8.900000 0.000000 + 516.0000 8.900000 0.000000 + 518.0000 8.900000 0.000000 + 520.0000 8.900000 0.000000 + 522.0000 8.900000 0.000000 + 524.0000 8.900000 0.000000 + 526.0000 8.900000 0.000000 + 528.0000 8.900000 0.000000 + 530.0000 8.900000 0.000000 + 532.0000 8.900000 0.000000 + 534.0000 8.900000 0.000000 + 536.0000 8.900000 0.000000 + 538.0000 8.900000 0.000000 + 540.0000 8.900000 0.000000 + 542.0000 8.900000 0.000000 + 544.0000 8.900000 0.000000 + 546.0000 8.900000 0.000000 + 548.0000 8.900000 0.000000 + 550.0000 8.900000 0.000000 + 552.0000 8.900000 0.000000 + 554.0000 8.900000 0.000000 + 556.0000 8.900000 0.000000 + 558.0000 8.900000 0.000000 + 560.0000 8.900000 0.000000 + 562.0000 8.900000 0.000000 + 564.0000 8.900000 0.000000 + 566.0000 8.900000 0.000000 + 568.0000 8.900000 0.000000 + 570.0000 8.900000 0.000000 + 572.0000 8.900000 0.000000 + 574.0000 8.900000 0.000000 + 576.0000 8.900000 0.000000 + 578.0000 8.900000 0.000000 + 580.0000 8.900000 0.000000 + 582.0000 8.900000 0.000000 + 584.0000 8.900000 0.000000 + 586.0000 8.900000 0.000000 + 588.0000 8.900000 0.000000 + 590.0000 8.900000 0.000000 + 592.0000 8.900000 0.000000 + 594.0000 8.900000 0.000000 + 596.0000 8.900000 0.000000 + 598.0000 8.900000 0.000000 + 600.0000 8.900000 0.000000 + 602.0000 8.900000 0.000000 + 604.0000 8.900000 0.000000 + 606.0000 8.900000 0.000000 + 608.0000 8.900000 0.000000 + 610.0000 8.900000 0.000000 + 612.0000 8.900000 0.000000 + 614.0000 8.900000 0.000000 + 616.0000 8.900000 0.000000 + 618.0000 8.900000 0.000000 + 620.0000 8.900000 0.000000 + 622.0000 8.900000 0.000000 + 624.0000 8.900000 0.000000 + 626.0000 8.900000 0.000000 + 628.0000 8.900000 0.000000 + 630.0000 8.900000 0.000000 + 632.0000 8.900000 0.000000 + 634.0000 8.900000 0.000000 + 636.0000 8.900000 0.000000 + 638.0000 8.900000 0.000000 + 640.0000 8.900000 0.000000 + 642.0000 8.900000 0.000000 + 644.0000 8.900000 0.000000 + 646.0000 8.900000 0.000000 + 648.0000 8.900000 0.000000 + 650.0000 8.900000 0.000000 + 652.0000 8.900000 0.000000 + 654.0000 8.900000 0.000000 + 656.0000 8.900000 0.000000 + 658.0000 8.900000 0.000000 + 660.0000 8.900000 0.000000 + 662.0000 8.900000 0.000000 + 664.0000 8.900000 0.000000 + 666.0000 8.900000 0.000000 + 668.0000 8.900000 0.000000 + 670.0000 8.900000 0.000000 + 672.0000 8.900000 0.000000 + 674.0000 8.900000 0.000000 + 676.0000 8.900000 0.000000 + 678.0000 8.900000 0.000000 + 680.0000 8.900000 0.000000 + 682.0000 8.900000 0.000000 + 684.0000 8.900000 0.000000 + 686.0000 8.900000 0.000000 + 688.0000 8.900000 0.000000 + 690.0000 8.900000 0.000000 + 692.0000 8.900000 0.000000 + 694.0000 8.900000 0.000000 + 696.0000 8.900000 0.000000 + 698.0000 8.900000 0.000000 + 700.0000 8.900000 0.000000 diff --git a/tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps006.log b/tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps006.log new file mode 100644 index 0000000..cd282a8 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/epsLog/seteps006.log @@ -0,0 +1,333 @@ + 2 2 +No-layers +MnO 4.950000 1 + 269.0000 16.00000 0.5000000E-01 +Platinum 8.900000 1 + 269.0000 16.00000 0.5000000E-01 + + 50.00000 4.950000 0.000000 8.900000 0.000000 + 52.00000 4.950000 0.000000 8.900000 0.000000 + 54.00000 4.950000 0.000000 8.900000 0.000000 + 56.00000 4.950000 0.000000 8.900000 0.000000 + 58.00000 4.950000 0.000000 8.900000 0.000000 + 60.00000 4.950000 0.000000 8.900000 0.000000 + 62.00000 4.950000 0.000000 8.900000 0.000000 + 64.00000 4.950000 0.000000 8.900000 0.000000 + 66.00000 4.950000 0.000000 8.900000 0.000000 + 68.00000 4.950000 0.000000 8.900000 0.000000 + 70.00000 4.950000 0.000000 8.900000 0.000000 + 72.00000 4.950000 0.000000 8.900000 0.000000 + 74.00000 4.950000 0.000000 8.900000 0.000000 + 76.00000 4.950000 0.000000 8.900000 0.000000 + 78.00000 4.950000 0.000000 8.900000 0.000000 + 80.00000 4.950000 0.000000 8.900000 0.000000 + 82.00000 4.950000 0.000000 8.900000 0.000000 + 84.00000 4.950000 0.000000 8.900000 0.000000 + 86.00000 4.950000 0.000000 8.900000 0.000000 + 88.00000 4.950000 0.000000 8.900000 0.000000 + 90.00000 4.950000 0.000000 8.900000 0.000000 + 92.00000 4.950000 0.000000 8.900000 0.000000 + 94.00000 4.950000 0.000000 8.900000 0.000000 + 96.00000 4.950000 0.000000 8.900000 0.000000 + 98.00000 4.950000 0.000000 8.900000 0.000000 + 100.0000 4.950000 0.000000 8.900000 0.000000 + 102.0000 4.950000 0.000000 8.900000 0.000000 + 104.0000 4.950000 0.000000 8.900000 0.000000 + 106.0000 4.950000 0.000000 8.900000 0.000000 + 108.0000 4.950000 0.000000 8.900000 0.000000 + 110.0000 4.950000 0.000000 8.900000 0.000000 + 112.0000 4.950000 0.000000 8.900000 0.000000 + 114.0000 4.950000 0.000000 8.900000 0.000000 + 116.0000 4.950000 0.000000 8.900000 0.000000 + 118.0000 4.950000 0.000000 8.900000 0.000000 + 120.0000 4.950000 0.000000 8.900000 0.000000 + 122.0000 4.950000 0.000000 8.900000 0.000000 + 124.0000 4.950000 0.000000 8.900000 0.000000 + 126.0000 4.950000 0.000000 8.900000 0.000000 + 128.0000 4.950000 0.000000 8.900000 0.000000 + 130.0000 4.950000 0.000000 8.900000 0.000000 + 132.0000 4.950000 0.000000 8.900000 0.000000 + 134.0000 4.950000 0.000000 8.900000 0.000000 + 136.0000 4.950000 0.000000 8.900000 0.000000 + 138.0000 4.950000 0.000000 8.900000 0.000000 + 140.0000 4.950000 0.000000 8.900000 0.000000 + 142.0000 4.950000 0.000000 8.900000 0.000000 + 144.0000 4.950000 0.000000 8.900000 0.000000 + 146.0000 4.950000 0.000000 8.900000 0.000000 + 148.0000 4.950000 0.000000 8.900000 0.000000 + 150.0000 4.950000 0.000000 8.900000 0.000000 + 152.0000 4.950000 0.000000 8.900000 0.000000 + 154.0000 4.950000 0.000000 8.900000 0.000000 + 156.0000 4.950000 0.000000 8.900000 0.000000 + 158.0000 4.950000 0.000000 8.900000 0.000000 + 160.0000 4.950000 0.000000 8.900000 0.000000 + 162.0000 4.950000 0.000000 8.900000 0.000000 + 164.0000 4.950000 0.000000 8.900000 0.000000 + 166.0000 4.950000 0.000000 8.900000 0.000000 + 168.0000 4.950000 0.000000 8.900000 0.000000 + 170.0000 4.950000 0.000000 8.900000 0.000000 + 172.0000 4.950000 0.000000 8.900000 0.000000 + 174.0000 4.950000 0.000000 8.900000 0.000000 + 176.0000 4.950000 0.000000 8.900000 0.000000 + 178.0000 4.950000 0.000000 8.900000 0.000000 + 180.0000 4.950000 0.000000 8.900000 0.000000 + 182.0000 4.950000 0.000000 8.900000 0.000000 + 184.0000 4.950000 0.000000 8.900000 0.000000 + 186.0000 4.950000 0.000000 8.900000 0.000000 + 188.0000 4.950000 0.000000 8.900000 0.000000 + 190.0000 4.950000 0.000000 8.900000 0.000000 + 192.0000 4.950000 0.000000 8.900000 0.000000 + 194.0000 4.950000 0.000000 8.900000 0.000000 + 196.0000 4.950000 0.000000 8.900000 0.000000 + 198.0000 4.950000 0.000000 8.900000 0.000000 + 200.0000 4.950000 0.000000 8.900000 0.000000 + 202.0000 4.950000 0.000000 8.900000 0.000000 + 204.0000 4.950000 0.000000 8.900000 0.000000 + 206.0000 4.950000 0.000000 8.900000 0.000000 + 208.0000 4.950000 0.000000 8.900000 0.000000 + 210.0000 4.950000 0.000000 8.900000 0.000000 + 212.0000 4.950000 0.000000 8.900000 0.000000 + 214.0000 4.950000 0.000000 8.900000 0.000000 + 216.0000 4.950000 0.000000 8.900000 0.000000 + 218.0000 4.950000 0.000000 8.900000 0.000000 + 220.0000 4.950000 0.000000 8.900000 0.000000 + 222.0000 4.950000 0.000000 8.900000 0.000000 + 224.0000 4.950000 0.000000 8.900000 0.000000 + 226.0000 4.950000 0.000000 8.900000 0.000000 + 228.0000 4.950000 0.000000 8.900000 0.000000 + 230.0000 4.950000 0.000000 8.900000 0.000000 + 232.0000 4.950000 0.000000 8.900000 0.000000 + 234.0000 4.950000 0.000000 8.900000 0.000000 + 236.0000 4.950000 0.000000 8.900000 0.000000 + 238.0000 4.950000 0.000000 8.900000 0.000000 + 240.0000 4.950000 0.000000 8.900000 0.000000 + 242.0000 4.950000 0.000000 8.900000 0.000000 + 244.0000 4.950000 0.000000 8.900000 0.000000 + 246.0000 4.950000 0.000000 8.900000 0.000000 + 248.0000 4.950000 0.000000 8.900000 0.000000 + 250.0000 4.950000 0.000000 8.900000 0.000000 + 252.0000 4.950000 0.000000 8.900000 0.000000 + 254.0000 4.950000 0.000000 8.900000 0.000000 + 256.0000 4.950000 0.000000 8.900000 0.000000 + 258.0000 4.950000 0.000000 8.900000 0.000000 + 260.0000 4.950000 0.000000 8.900000 0.000000 + 262.0000 4.950000 0.000000 8.900000 0.000000 + 264.0000 4.950000 0.000000 8.900000 0.000000 + 266.0000 4.950000 0.000000 8.900000 0.000000 + 268.0000 4.950000 0.000000 8.900000 0.000000 + 270.0000 4.950000 0.000000 8.900000 0.000000 + 272.0000 4.950000 0.000000 8.900000 0.000000 + 274.0000 4.950000 0.000000 8.900000 0.000000 + 276.0000 4.950000 0.000000 8.900000 0.000000 + 278.0000 4.950000 0.000000 8.900000 0.000000 + 280.0000 4.950000 0.000000 8.900000 0.000000 + 282.0000 4.950000 0.000000 8.900000 0.000000 + 284.0000 4.950000 0.000000 8.900000 0.000000 + 286.0000 4.950000 0.000000 8.900000 0.000000 + 288.0000 4.950000 0.000000 8.900000 0.000000 + 290.0000 4.950000 0.000000 8.900000 0.000000 + 292.0000 4.950000 0.000000 8.900000 0.000000 + 294.0000 4.950000 0.000000 8.900000 0.000000 + 296.0000 4.950000 0.000000 8.900000 0.000000 + 298.0000 4.950000 0.000000 8.900000 0.000000 + 300.0000 4.950000 0.000000 8.900000 0.000000 + 302.0000 4.950000 0.000000 8.900000 0.000000 + 304.0000 4.950000 0.000000 8.900000 0.000000 + 306.0000 4.950000 0.000000 8.900000 0.000000 + 308.0000 4.950000 0.000000 8.900000 0.000000 + 310.0000 4.950000 0.000000 8.900000 0.000000 + 312.0000 4.950000 0.000000 8.900000 0.000000 + 314.0000 4.950000 0.000000 8.900000 0.000000 + 316.0000 4.950000 0.000000 8.900000 0.000000 + 318.0000 4.950000 0.000000 8.900000 0.000000 + 320.0000 4.950000 0.000000 8.900000 0.000000 + 322.0000 4.950000 0.000000 8.900000 0.000000 + 324.0000 4.950000 0.000000 8.900000 0.000000 + 326.0000 4.950000 0.000000 8.900000 0.000000 + 328.0000 4.950000 0.000000 8.900000 0.000000 + 330.0000 4.950000 0.000000 8.900000 0.000000 + 332.0000 4.950000 0.000000 8.900000 0.000000 + 334.0000 4.950000 0.000000 8.900000 0.000000 + 336.0000 4.950000 0.000000 8.900000 0.000000 + 338.0000 4.950000 0.000000 8.900000 0.000000 + 340.0000 4.950000 0.000000 8.900000 0.000000 + 342.0000 4.950000 0.000000 8.900000 0.000000 + 344.0000 4.950000 0.000000 8.900000 0.000000 + 346.0000 4.950000 0.000000 8.900000 0.000000 + 348.0000 4.950000 0.000000 8.900000 0.000000 + 350.0000 4.950000 0.000000 8.900000 0.000000 + 352.0000 4.950000 0.000000 8.900000 0.000000 + 354.0000 4.950000 0.000000 8.900000 0.000000 + 356.0000 4.950000 0.000000 8.900000 0.000000 + 358.0000 4.950000 0.000000 8.900000 0.000000 + 360.0000 4.950000 0.000000 8.900000 0.000000 + 362.0000 4.950000 0.000000 8.900000 0.000000 + 364.0000 4.950000 0.000000 8.900000 0.000000 + 366.0000 4.950000 0.000000 8.900000 0.000000 + 368.0000 4.950000 0.000000 8.900000 0.000000 + 370.0000 4.950000 0.000000 8.900000 0.000000 + 372.0000 4.950000 0.000000 8.900000 0.000000 + 374.0000 4.950000 0.000000 8.900000 0.000000 + 376.0000 4.950000 0.000000 8.900000 0.000000 + 378.0000 4.950000 0.000000 8.900000 0.000000 + 380.0000 4.950000 0.000000 8.900000 0.000000 + 382.0000 4.950000 0.000000 8.900000 0.000000 + 384.0000 4.950000 0.000000 8.900000 0.000000 + 386.0000 4.950000 0.000000 8.900000 0.000000 + 388.0000 4.950000 0.000000 8.900000 0.000000 + 390.0000 4.950000 0.000000 8.900000 0.000000 + 392.0000 4.950000 0.000000 8.900000 0.000000 + 394.0000 4.950000 0.000000 8.900000 0.000000 + 396.0000 4.950000 0.000000 8.900000 0.000000 + 398.0000 4.950000 0.000000 8.900000 0.000000 + 400.0000 4.950000 0.000000 8.900000 0.000000 + 402.0000 4.950000 0.000000 8.900000 0.000000 + 404.0000 4.950000 0.000000 8.900000 0.000000 + 406.0000 4.950000 0.000000 8.900000 0.000000 + 408.0000 4.950000 0.000000 8.900000 0.000000 + 410.0000 4.950000 0.000000 8.900000 0.000000 + 412.0000 4.950000 0.000000 8.900000 0.000000 + 414.0000 4.950000 0.000000 8.900000 0.000000 + 416.0000 4.950000 0.000000 8.900000 0.000000 + 418.0000 4.950000 0.000000 8.900000 0.000000 + 420.0000 4.950000 0.000000 8.900000 0.000000 + 422.0000 4.950000 0.000000 8.900000 0.000000 + 424.0000 4.950000 0.000000 8.900000 0.000000 + 426.0000 4.950000 0.000000 8.900000 0.000000 + 428.0000 4.950000 0.000000 8.900000 0.000000 + 430.0000 4.950000 0.000000 8.900000 0.000000 + 432.0000 4.950000 0.000000 8.900000 0.000000 + 434.0000 4.950000 0.000000 8.900000 0.000000 + 436.0000 4.950000 0.000000 8.900000 0.000000 + 438.0000 4.950000 0.000000 8.900000 0.000000 + 440.0000 4.950000 0.000000 8.900000 0.000000 + 442.0000 4.950000 0.000000 8.900000 0.000000 + 444.0000 4.950000 0.000000 8.900000 0.000000 + 446.0000 4.950000 0.000000 8.900000 0.000000 + 448.0000 4.950000 0.000000 8.900000 0.000000 + 450.0000 4.950000 0.000000 8.900000 0.000000 + 452.0000 4.950000 0.000000 8.900000 0.000000 + 454.0000 4.950000 0.000000 8.900000 0.000000 + 456.0000 4.950000 0.000000 8.900000 0.000000 + 458.0000 4.950000 0.000000 8.900000 0.000000 + 460.0000 4.950000 0.000000 8.900000 0.000000 + 462.0000 4.950000 0.000000 8.900000 0.000000 + 464.0000 4.950000 0.000000 8.900000 0.000000 + 466.0000 4.950000 0.000000 8.900000 0.000000 + 468.0000 4.950000 0.000000 8.900000 0.000000 + 470.0000 4.950000 0.000000 8.900000 0.000000 + 472.0000 4.950000 0.000000 8.900000 0.000000 + 474.0000 4.950000 0.000000 8.900000 0.000000 + 476.0000 4.950000 0.000000 8.900000 0.000000 + 478.0000 4.950000 0.000000 8.900000 0.000000 + 480.0000 4.950000 0.000000 8.900000 0.000000 + 482.0000 4.950000 0.000000 8.900000 0.000000 + 484.0000 4.950000 0.000000 8.900000 0.000000 + 486.0000 4.950000 0.000000 8.900000 0.000000 + 488.0000 4.950000 0.000000 8.900000 0.000000 + 490.0000 4.950000 0.000000 8.900000 0.000000 + 492.0000 4.950000 0.000000 8.900000 0.000000 + 494.0000 4.950000 0.000000 8.900000 0.000000 + 496.0000 4.950000 0.000000 8.900000 0.000000 + 498.0000 4.950000 0.000000 8.900000 0.000000 + 500.0000 4.950000 0.000000 8.900000 0.000000 + 502.0000 4.950000 0.000000 8.900000 0.000000 + 504.0000 4.950000 0.000000 8.900000 0.000000 + 506.0000 4.950000 0.000000 8.900000 0.000000 + 508.0000 4.950000 0.000000 8.900000 0.000000 + 510.0000 4.950000 0.000000 8.900000 0.000000 + 512.0000 4.950000 0.000000 8.900000 0.000000 + 514.0000 4.950000 0.000000 8.900000 0.000000 + 516.0000 4.950000 0.000000 8.900000 0.000000 + 518.0000 4.950000 0.000000 8.900000 0.000000 + 520.0000 4.950000 0.000000 8.900000 0.000000 + 522.0000 4.950000 0.000000 8.900000 0.000000 + 524.0000 4.950000 0.000000 8.900000 0.000000 + 526.0000 4.950000 0.000000 8.900000 0.000000 + 528.0000 4.950000 0.000000 8.900000 0.000000 + 530.0000 4.950000 0.000000 8.900000 0.000000 + 532.0000 4.950000 0.000000 8.900000 0.000000 + 534.0000 4.950000 0.000000 8.900000 0.000000 + 536.0000 4.950000 0.000000 8.900000 0.000000 + 538.0000 4.950000 0.000000 8.900000 0.000000 + 540.0000 4.950000 0.000000 8.900000 0.000000 + 542.0000 4.950000 0.000000 8.900000 0.000000 + 544.0000 4.950000 0.000000 8.900000 0.000000 + 546.0000 4.950000 0.000000 8.900000 0.000000 + 548.0000 4.950000 0.000000 8.900000 0.000000 + 550.0000 4.950000 0.000000 8.900000 0.000000 + 552.0000 4.950000 0.000000 8.900000 0.000000 + 554.0000 4.950000 0.000000 8.900000 0.000000 + 556.0000 4.950000 0.000000 8.900000 0.000000 + 558.0000 4.950000 0.000000 8.900000 0.000000 + 560.0000 4.950000 0.000000 8.900000 0.000000 + 562.0000 4.950000 0.000000 8.900000 0.000000 + 564.0000 4.950000 0.000000 8.900000 0.000000 + 566.0000 4.950000 0.000000 8.900000 0.000000 + 568.0000 4.950000 0.000000 8.900000 0.000000 + 570.0000 4.950000 0.000000 8.900000 0.000000 + 572.0000 4.950000 0.000000 8.900000 0.000000 + 574.0000 4.950000 0.000000 8.900000 0.000000 + 576.0000 4.950000 0.000000 8.900000 0.000000 + 578.0000 4.950000 0.000000 8.900000 0.000000 + 580.0000 4.950000 0.000000 8.900000 0.000000 + 582.0000 4.950000 0.000000 8.900000 0.000000 + 584.0000 4.950000 0.000000 8.900000 0.000000 + 586.0000 4.950000 0.000000 8.900000 0.000000 + 588.0000 4.950000 0.000000 8.900000 0.000000 + 590.0000 4.950000 0.000000 8.900000 0.000000 + 592.0000 4.950000 0.000000 8.900000 0.000000 + 594.0000 4.950000 0.000000 8.900000 0.000000 + 596.0000 4.950000 0.000000 8.900000 0.000000 + 598.0000 4.950000 0.000000 8.900000 0.000000 + 600.0000 4.950000 0.000000 8.900000 0.000000 + 602.0000 4.950000 0.000000 8.900000 0.000000 + 604.0000 4.950000 0.000000 8.900000 0.000000 + 606.0000 4.950000 0.000000 8.900000 0.000000 + 608.0000 4.950000 0.000000 8.900000 0.000000 + 610.0000 4.950000 0.000000 8.900000 0.000000 + 612.0000 4.950000 0.000000 8.900000 0.000000 + 614.0000 4.950000 0.000000 8.900000 0.000000 + 616.0000 4.950000 0.000000 8.900000 0.000000 + 618.0000 4.950000 0.000000 8.900000 0.000000 + 620.0000 4.950000 0.000000 8.900000 0.000000 + 622.0000 4.950000 0.000000 8.900000 0.000000 + 624.0000 4.950000 0.000000 8.900000 0.000000 + 626.0000 4.950000 0.000000 8.900000 0.000000 + 628.0000 4.950000 0.000000 8.900000 0.000000 + 630.0000 4.950000 0.000000 8.900000 0.000000 + 632.0000 4.950000 0.000000 8.900000 0.000000 + 634.0000 4.950000 0.000000 8.900000 0.000000 + 636.0000 4.950000 0.000000 8.900000 0.000000 + 638.0000 4.950000 0.000000 8.900000 0.000000 + 640.0000 4.950000 0.000000 8.900000 0.000000 + 642.0000 4.950000 0.000000 8.900000 0.000000 + 644.0000 4.950000 0.000000 8.900000 0.000000 + 646.0000 4.950000 0.000000 8.900000 0.000000 + 648.0000 4.950000 0.000000 8.900000 0.000000 + 650.0000 4.950000 0.000000 8.900000 0.000000 + 652.0000 4.950000 0.000000 8.900000 0.000000 + 654.0000 4.950000 0.000000 8.900000 0.000000 + 656.0000 4.950000 0.000000 8.900000 0.000000 + 658.0000 4.950000 0.000000 8.900000 0.000000 + 660.0000 4.950000 0.000000 8.900000 0.000000 + 662.0000 4.950000 0.000000 8.900000 0.000000 + 664.0000 4.950000 0.000000 8.900000 0.000000 + 666.0000 4.950000 0.000000 8.900000 0.000000 + 668.0000 4.950000 0.000000 8.900000 0.000000 + 670.0000 4.950000 0.000000 8.900000 0.000000 + 672.0000 4.950000 0.000000 8.900000 0.000000 + 674.0000 4.950000 0.000000 8.900000 0.000000 + 676.0000 4.950000 0.000000 8.900000 0.000000 + 678.0000 4.950000 0.000000 8.900000 0.000000 + 680.0000 4.950000 0.000000 8.900000 0.000000 + 682.0000 4.950000 0.000000 8.900000 0.000000 + 684.0000 4.950000 0.000000 8.900000 0.000000 + 686.0000 4.950000 0.000000 8.900000 0.000000 + 688.0000 4.950000 0.000000 8.900000 0.000000 + 690.0000 4.950000 0.000000 8.900000 0.000000 + 692.0000 4.950000 0.000000 8.900000 0.000000 + 694.0000 4.950000 0.000000 8.900000 0.000000 + 696.0000 4.950000 0.000000 8.900000 0.000000 + 698.0000 4.950000 0.000000 8.900000 0.000000 + 700.0000 4.950000 0.000000 8.900000 0.000000 diff --git a/tests/SetEpsTestCasesFromScratch/WFW/extend3.f90 b/tests/SetEpsTestCasesFromScratch/WFW/extend3.f90 new file mode 100644 index 0000000..1ebde11 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/extend3.f90 @@ -0,0 +1,22 @@ +subroutine extend3(array, extension) + +! increase the size of a allocatable array by extension. +! if the array is not allocated, create one of size extension. + + implicit none + + double precision, intent(in out), allocatable :: array(:, :) + integer, intent(in) :: extension + + integer :: old_size_1, old_size_2 + double precision, allocatable :: tmp_array(:, :) + + old_size_1 = size(array, 1) + old_size_2 = size(array, 2) + allocate(tmp_array(old_size_1, old_size_2 + extension)) + tmp_array(1:old_size_1, 1:old_size_2) = array + deallocate(array) + call move_alloc(tmp_array, array) + + return +end subroutine extend3 diff --git a/tests/SetEpsTestCasesFromScratch/WFW/fint1.f90 b/tests/SetEpsTestCasesFromScratch/WFW/fint1.f90 new file mode 100644 index 0000000..62361a0 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/fint1.f90 @@ -0,0 +1,64 @@ +module fint1_mod +contains +double precision function fint1(u, eps, thick, layers, nper, eps_size) + +! ****************************************************************** +! * * +! * integration over the azimutal angle from 0.0 to pi * +! * * +! * eps: epsilon * +! * thick: thickness * +! * layers: number of layers * +! * nper: number of periodic layers * +! * eps_size: number of layers * +! * * +! ****************************************************************** + + use surlos_mod + use param_mod + + implicit none + + double precision, intent(in) :: u + double complex, intent(in) :: eps(eps_size) + double precision, intent(in) :: thick(eps_size) + integer, intent(in) :: layers, nper, eps_size + + double precision :: den, dif, e, pi, rom, rop, sum, t, u2 + double precision :: usurlo + + data pi / 3.141592653589793238d0 / + +! write (*,*) 'fint1:' +! write (*,*) 'thick: ', size(thick) +! write (*,*) 'eps: ', size(eps) + + if (u == 0.0d0) then + fint1 = 0.0d0 + return + endif + e = tanpsi * u + u2 = u**2 + rom = (1.0d0 - e)**2 + u2 + rop = (1.0d0 + e)**2 + u2 + sum = rop + rom + rom = dsqrt(rom) + rop = dsqrt(rop) + dif = rop - rom + den = dsqrt((2.0d0 - dif) * (2.0d0 + dif)) * rop * rom + fint1 = pi * u2 * (4 * sum - dif**2 * (sum - rop * rom)) / den**3 + if (rational) then + return + endif + if (user) then + fint1 = fint1 * usurlo(ru * u, wn) + else + fint1 = fint1 * surlos(ru * u, eps, thick, layers, nper) + 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 + endif + return +end function fint1 +end module fint1_mod diff --git a/tests/SetEpsTestCasesFromScratch/WFW/fint2.f90 b/tests/SetEpsTestCasesFromScratch/WFW/fint2.f90 new file mode 100644 index 0000000..12baaf6 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/fint2.f90 @@ -0,0 +1,66 @@ +module fint2_mod +contains +double precision function fint2(u, eps, thick, layers, nper, eps_size) + +! ****************************************************************** +! * * +! * integration over the azimutal angle from 0.0 to phi < pi * +! * * +! * eps: epsilon * +! * thick: thickness * +! * layers: number of layers * +! * nper: number of periodic layers * +! * eps_size: number of layers * +! * * +! ****************************************************************** + + use surlos_mod + use param_mod + + implicit none + + double precision, intent(in) :: u + double complex, intent(in) :: eps(eps_size) + double precision, intent(in) :: thick(eps_size) + integer, intent(in) :: layers, nper, eps_size + + double precision :: arg, b2, c, phi, t, x + double precision :: phint, usurlo + +! write (*,*) 'fint2:' +! write (*,*) 'thick: ', size(thick) +! write (*,*) 'eps: ', size(eps) + + if (u == 0.0d0) then + fint2 = 0.0d0 + return + endif + b2 = bcoef**2 + c = (1.0d0 - elleps) * (cospsi * u)**2 + (bcoef - ccoef) * (bcoef + ccoef) + if (dabs(acoef * c) > 1.0d-03 * b2) then + x = (bcoef - dsqrt(b2 - acoef * c)) / acoef + else + x = acoef * c / b2 + x = 0.5d0 * c * (1.d0 + 0.25d0 * x * (1.d0 + 0.5d0 * x * (1.d0 + 0.625d0 * x))) / bcoef + endif + arg = x / u + if (dabs(arg) > 1.0d0) then + arg = dsign(1.0d0, arg) + endif + phi = dacos(arg) + fint2 = phint(phi, tanpsi, u) + if (rational) then + return + endif + if (user) then + fint2 = fint2 * usurlo(ru * u, wn) + else + fint2 = fint2 * surlos(ru * u, eps, thick, layers, nper) + 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 + endif + return +end function fint2 +end module fint2_mod diff --git a/tests/SetEpsTestCasesFromScratch/WFW/fint3.f90 b/tests/SetEpsTestCasesFromScratch/WFW/fint3.f90 new file mode 100644 index 0000000..7135e6f --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/fint3.f90 @@ -0,0 +1,59 @@ +module fint3_mod +contains +double precision function fint3(u, eps, thick, layers, nper, eps_size) + +! ****************************************************************** +! * * +! * integration over the azimutal angle from phi1 > 0 to phi2 < pi * +! * * +! * eps: epsilon * +! * thick: thickness * +! * layers: number of layers * +! * nper: number of periodic layers * +! * eps_size: number of layers * +! * * +! ****************************************************************** + + use surlos_mod + use param_mod + + implicit none + + double precision, intent(in) :: u + double complex, intent(in) :: eps(eps_size) + double precision, intent(in) :: thick(eps_size) + integer, intent(in) :: layers, nper, eps_size + + double precision :: arg, phi1, phi2, rac, t + double precision :: phint, usurlo + +! write (*,*) 'fint3:' +! write (*,*) 'thick: ', size(thick) +! write (*,*) 'eps: ', size(eps) + + if (u == 0.0d0) then + fint3 = 0.0d0 + return + endif + rac = dsign(1.0d0, acoef) * cospsi * dsqrt((1.0d0 - elleps) * acoef * (um - u) * (um + u)) + arg = (bcoef - rac) / (u * acoef) + if (dabs(arg) > 1.0d0) arg = dsign(1.0d0, arg) + phi2 = dacos(arg) + fint3 = phint(phi2, tanpsi, u) + arg = (bcoef + rac) / (u * acoef) + if (dabs(arg) > 1.0d0) arg = dsign(1.0d0, arg) + phi1 = dacos(arg) + fint3 = fint3 - phint(phi1, tanpsi, u) + if (rational) return + if (user) then + fint3 = fint3 * usurlo(ru * u, wn) + else + fint3 = fint3 * surlos(ru * u, eps, thick, layers, nper) + 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 + endif + return +end function fint3 +end module fint3_mod diff --git a/tests/SetEpsTestCasesFromScratch/WFW/fun.f90 b/tests/SetEpsTestCasesFromScratch/WFW/fun.f90 new file mode 100644 index 0000000..b077a2a --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/fun.f90 @@ -0,0 +1,23 @@ +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. * +! * * +! ****************************************************************** + + use param_mod + + implicit none + + double precision, intent(in) :: phi + + double precision :: sinphi + + sinphi = dsin(phi) + fun = dsqrt((1.0d0 - elleps + elleps * sinphi**2) * & + (1.0d0 - sinpsi * sinphi) * & + (1.0d0 + sinpsi * sinphi)) + return +end function fun diff --git a/tests/SetEpsTestCasesFromScratch/WFW/get_commandline_options.f90 b/tests/SetEpsTestCasesFromScratch/WFW/get_commandline_options.f90 new file mode 100644 index 0000000..e9977be --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/get_commandline_options.f90 @@ -0,0 +1,93 @@ +module get_commandline_options_mod +contains +subroutine get_commandline_options(eelsin_name, eelsou_name, bosin_name, bosou_name) + +! This routine defines the commandline options, parses the commandline and sets the +! filenames of I/O files. +! +! It uses sufr_getopt from libSUFR and has been derived from getopt_long_example.f90 +! +! KMS + + use sufr_getopt, only: getopt_t, getopt_long, longOption, optArg, getopt_long_help + + implicit none + + character (len = :), allocatable, intent(in out) :: eelsin_name, eelsou_name + character (len = :), allocatable, intent(in out) :: bosin_name, bosou_name + + include 'version.inc' + + character :: option + integer :: status + + ! Set up the longOpts struct to define the valid options: + ! short option, long option, argument (no = 0 / yes = 1), short description + type(getopt_t) :: longOpts(8) = & + [ & + getopt_t('v', 'version', 0, 'Print version'), & + getopt_t('V', 'version', 0, 'Print version'), & + getopt_t('h', 'help', 0, 'Print options'), & + getopt_t('d', 'dir', 1, 'Change I/O directory'), & + getopt_t('e', 'eelsin', 1, 'Name of EELS input file (eelsin)'), & + getopt_t('f', 'eelsou', 1, 'Name of EELS output file (optional)'), & + getopt_t('b', 'bosin', 1, 'Name of BOSON input file (bosin)'), & + getopt_t('c', 'bosou', 1, 'Name of BOSON output file (bosou)') & + ] + + do ! scan all the command-line parameters + + ! getopt_long() returns a single character" ">","!",".", or the short-option character (e.g. "a" for -a). + ! It also sets two 'global' variables through the SUFR_getopt module: + ! - longOption: the full option (e.g. "-a" or "--all") including the dashes + ! - optArg: the argument following the option (if required and present) + option = getopt_long(longOpts) + + ! Do different things depending on the option returned: + select case(option) + case('d') ! Change I/O directory + status = chdir(trim(optArg)) + if (status /= 0) then + write (*,*) 'WARNING: change directory failed!' + write (*,*) 'Directory tried: ', trim(optarg) + write (*,*) 'Error code (see: man chdir): ', status + write (*,*) 'Continuing in the start directory.' + write (*,*) '' + end if + case('e') + eelsin_name = optArg + case('f') + eelsou_name = optArg + case('b') + bosin_name = optArg + case('c') + bosou_name = optArg + case('v') + write (*,*) 'eels-boson version: ' // version + stop + case('V') + write (*,*) 'eels-boson version: ' // version + stop + case('h') + call getopt_long_help(longOpts) + stop + case('!') ! Unknown option (starting with "-" or "--") + write (*,*) 'WARNING: unknown option: ' // trim(optArg) + call getopt_long_help(longOpts) + stop + case('.') ! Parameter is not an option (i.e., it doesn't start with "-" or "--") + write (*,*) 'WARNING: parameter without option: ' // trim(optArg) + call getopt_long_help(longOpts) + stop + case default + write (*,*) 'WARNING: unknown option: ' // trim(longOption) + call getopt_long_help(longOpts) + stop + case('>') ! Last parameter. Exit case statement + exit + end select + end do + + return +end subroutine get_commandline_options +end module get_commandline_options_mod diff --git a/tests/SetEpsTestCasesFromScratch/WFW/myEels20-seteps.f90 b/tests/SetEpsTestCasesFromScratch/WFW/myEels20-seteps.f90 new file mode 100755 index 0000000..27e82c5 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/myEels20-seteps.f90 @@ -0,0 +1,71 @@ +module seteps_mod +contains +!subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps) +! Angepasste Version +subroutine seteps(nLayer, 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) :: nLayer +! integer, dimension(nLayer),intent(in) :: nos +! integer, intent(in) :: osc_size +! double precision, dimension(3,osc_size),intent(in) :: osc + integer, intent(in) :: nos(:) + double precision, intent(in) :: osc(:,:) + character (len=10), intent(in) :: name(:) + integer, intent(in) :: layers + character (len=10), intent(in) :: mode + + !f2py depend(osc_size) osc + double precision, dimension(nLayer),intent(in) :: epsinf + double precision, intent(in) :: wn + double complex, dimension(nLayer), intent(out) :: eps + !f2py depend(nLayer) nos, epsinf, eps + + double complex :: nomi, deno, addeps + double precision :: wn2, b + integer j, k, l, m + logical debugFirstRun + + common /control/ debugFirstRun + + j = 0 + do l = 1, nLayer ! loop over different thin film layers + m = nos(l)/2 ! m number of TO modes = offset to reach the LO mode in the joint TO-LO list + nomi = dcmplx(1.0d0, 0.0d0) + deno = dcmplx(1.0d0, 0.0d0) + addeps = dcmplx(0.0d0, 0.0d0) + wn2 = wn**2 + do k = 1, m ! loop over all TO modes + j = j + 1 + + if (osc(1,j) > 0.) then ! positive TO mode: 'Kurosawa' form: _Multiplicative_ phonon mode + b = wn/osc(1, j+m) + nomi = nomi * osc(1, j+m)**2 * (1.0 - b * dcmplx(b, osc(3,j+m)/osc(1, j+m)) ) + deno =deno * (osc(1,j)**2 - wn * dcmplx( wn, osc(3,j) ) ) + + else if (osc(1,j) < 0.) then! Negative TO mode means: _Additive_ Lorentz oscillator with Q + addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 - wn2, -1*wn*osc(3,j)) ! Sign of imaginary part changed (WFW) + + else ! osc(1,j) = 0 -> it is a Drude term + addeps = addeps - dcmplx(osc(1,j+m)**2, wn*(osc(3,j)-osc(3,j+m))) /dcmplx(wn2, wn*osc(3,j)) + end if + + enddo + j = j+m ! we have already looped over the LO modes, therefore increase the index + eps(l) = epsinf(l) * (nomi / deno + addeps) ! brackets changed by HHe 230915 + enddo + debugFirstRun = .false. + +! **** log modification start + write (99, '(30g15.7)') wn, (eps(j), j = 1, nLayer) +! **** log modification end + + return +end subroutine seteps +end module seteps_mod diff --git a/tests/SetEpsTestCasesFromScratch/WFW/param.f90 b/tests/SetEpsTestCasesFromScratch/WFW/param.f90 new file mode 100644 index 0000000..a688b3c --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/param.f90 @@ -0,0 +1,5 @@ +module param_mod + double precision :: acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi + double precision :: ru, um, dlimf, wn + logical :: user, rational +end module param_mod diff --git a/tests/SetEpsTestCasesFromScratch/WFW/phint.f90 b/tests/SetEpsTestCasesFromScratch/WFW/phint.f90 new file mode 100644 index 0000000..57f15b5 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/phint.f90 @@ -0,0 +1,74 @@ +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 * +! * * +! * Reference: * +! * Ph. Lambin, J. P. Vigneron, and A. A. Lucas, * +! * Phys. Rev. B 32 (1985) 8203-8215. * +! * * +! ******************************************************************* + + 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 + + pi = 3.141592653589793238d0 + c = dcos(phi) + s = dsin(phi) + u2 = u**2 + e = a*u + 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) + else + rm = dsqrt((1.0d0 - e)**2 + u2) + tm = 0.5d0 * datan2(u, 1.0d0 - e) + rp = dsqrt((1.0d0 + e)**2 + u2) + tp = 0.5d0 * datan2(u, 1.0d0 + e) + root = dsqrt(rm * rp) + cpr = dcos(tm + tp) + spr = dsin(tm + tp) + x = 0.0d0 ! ensure initialization + if (c >= 0.0d0) then + x = s / (1.0d0 + c) + elseif (dabs(s) > 1.0d-07) then + x = (1.0d0 - c) / s + endif + if ((c >= 0.0d0) .or. (dabs(s) > 1.0d-07)) then + zeta = dsqrt(rm / rp) + zetar = -zeta * dsin(tm - tp) + zetai = zeta * dcos(tm - tp) + br = 0.5d0 * dlog(((zetar + x)**2 + zetai**2) / ((zetar - x)**2 + zetai**2)) + bi = datan2(zetai, zetar + x) - datan2(zetai, zetar - x) + rr = -(br * spr - bi * cpr) / root + ri = -(bi * spr + br * cpr) / root + d = e * s / ((1.0d0 - e * c)**2 + u2) + ar = d * (1.0d0 - e * c) - rr + u * ri + ai = -d * u - ri - u * rr + else + rr = -pi / root * cpr + ri = pi / root * spr + ar = -rr + u * ri + ai = -ri - u * rr + endif + qr = (ar * (cpr - spr) * (cpr + spr) + 2 * ai * cpr * spr) / (rm * rp) + phint = 0.5d0 * (ri / u - qr) + endif + return +end function phint diff --git a/tests/SetEpsTestCasesFromScratch/WFW/qrat.f90 b/tests/SetEpsTestCasesFromScratch/WFW/qrat.f90 new file mode 100644 index 0000000..7bbb966 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/qrat.f90 @@ -0,0 +1,14 @@ +double precision function qrat(x, alpha, beta, c1, c2) + +! 1 + x * (beta + c1 * x) +! qrat = --------------------------------------------- +! (1 + x * (beta + c2 * x)) * (1 + alpha * x)^2 + + 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 diff --git a/tests/SetEpsTestCasesFromScratch/WFW/quanc8.f90 b/tests/SetEpsTestCasesFromScratch/WFW/quanc8.f90 new file mode 100644 index 0000000..86ee4f3 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/quanc8.f90 @@ -0,0 +1,240 @@ +module quanc8_mod +contains +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. +! +! additional input .. +! +! eps epsilon +! thick thickness +! layers number of layers +! nper number of periodic layers + + 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) + + levmin = 1 + levmax = 30 + levout = 6 + nomax = 5000 + nofin = nomax - 8 * (levmax - levout + 2**(levout + 1)) + +! trouble when nofun reaches nofin + + w0 = 3956.0d0 / 14175.0d0 + w1 = 23552.0d0 / 14175.0d0 + w2 = -3712.0d0 / 14175.0d0 + w3 = 41984.0d0 / 14175.0d0 + w4 = -18160.0d0 / 14175.0d0 + +! initialize running sums to zero. + + flag = 0.0d0 + result = 0.0d0 + cor11 = 0.0d0 + errest = 0.0d0 + area = 0.0d0 + nofun = 0 + if (a == b) return + +! *** stage 2 *** initialization for first interval + + lev = 0 + nim = 1 + x0 = a + x(16) = b + qprev = 0.0d0 + f0 = fun(x0, eps, thick, layers, nper, size(eps)) + stone = (b - a) / 16 + x(8) = (x0 + x(16)) / 2 + x(4) = (x0 + x(8)) / 2 + x(12) = (x(8) + x(16)) / 2 + x(2) = (x0 + x(4)) / 2 + x(6) = (x(4) + x(8)) / 2 + x(10) = (x(8) + x(12)) / 2 + x(14) = (x(12) + x(16)) / 2 +!*** !$OMP PARALLEL shared (f, x, eps, thick, layers, nper) private (j) +!*** !$OMP DO + do j = 2, 16, 2 + f(j) = fun(x(j), eps, thick, layers, nper, size(eps)) + enddo +!*** !$OMP END DO NOWAIT +!*** !$OMP END PARALLEL + nofun = 9 + + 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. + + x(1) = (x0 + x(2)) / 2 + f(1) = fun(x(1), eps, thick, layers, nper, size(eps)) + do j = 3, 15, 2 + x(j) = (x(j - 1) + x(j + 1)) / 2 + enddo +!*** !$OMP PARALLEL shared (f, x, eps, thick, layers, nper) private (j) +!*** !$OMP DO SCHEDULE(auto) + do j = 3, 15, 2 + f(j) = fun(x(j), eps, thick, layers, nper, size(eps)) + enddo +!*** !$OMP END DO NOWAIT +!*** !$OMP END PARALLEL + nofun = nofun + 8 + step = (x(16) - x0) / 16.0d0 + qleft = (w0 * (f0 + f(8)) + w1 * (f(1) + f(7)) + w2 * (f(2) + f(6)) & + + w3 * (f(3) + f(5)) + w4 * f(4)) * step + 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 + qnow = qleft + qright(lev + 1) + qdiff = qnow - qprev + area = area + qdiff + +! *** stage 4 *** interval convergence test + + esterr = dabs(qdiff) / 1023 + tolerr = dmax1(abserr, relerr * dabs(area)) * (step / stone) + + if (lev >= levmin) then + if (lev >= levmax) then +! current level is levmax. + flag = flag + 1.0d0 + else + 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) + else + if (esterr > tolerr) then +! *** stage 5 *** no convergence +! locate next interval. + nim = 2 * nim + lev = lev + 1 +! store right hand elements for future use. + do i = 1, 8 + fsave(i, lev) = f(i + 8) + xsave(i, lev) = x(i + 8) + enddo +! assemble left hand elements for immediate use. + qprev = qleft + do i = 1, 8 + f(18 - 2 * i) = f(9 - i) + x(18 - 2 * i) = x(9 - i) + enddo + cycle + endif + endif + endif + +! *** stage 7 *** interval converged +! add contributions into running sums. + result = result + qnow + errest = errest + esterr + cor11 = cor11 + qdiff / 1023 +! locate next interval. + do while (nim /= 2 * (nim / 2)) + nim = nim / 2 + lev = lev - 1 + enddo + nim = nim + 1 + + if (lev <= 0) exit + +! assemble elements required for the next interval. + qprev = qright(lev) + x0 = x(16) + f0 = f(16) + do i = 1, 8 + f(2*i) = fsave(i, lev) + x(2*i) = xsave(i, lev) + enddo + cycle + else +! *** stage 5 *** no convergence +! locate next interval. + nim = 2 * nim + lev = lev + 1 +! store right hand elements for future use. + do i = 1, 8 + fsave(i, lev) = f(i + 8) + xsave(i, lev) = x(i + 8) + enddo +! assemble left hand elements for immediate use. + qprev = qleft + do i = 1, 8 + f(18 - 2 * i) = f(9 - i) + x(18 - 2 * i) = x(9 - i) + enddo + endif + + enddo + + ! *** stage 8 *** finalize and return + result = result + cor11 + +! make sure errest not less than roundoff level. + if (errest /= 0.0d0) then + temp = dabs(result) + errest + do while (temp == dabs(result)) + errest = 2 * errest + temp = dabs(result) + errest + enddo + endif + return +end subroutine quanc8 +end module quanc8_mod diff --git a/tests/SetEpsTestCasesFromScratch/WFW/queels.f90 b/tests/SetEpsTestCasesFromScratch/WFW/queels.f90 new file mode 100644 index 0000000..d2128df --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/queels.f90 @@ -0,0 +1,92 @@ +module queels_mod +contains +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 * +! * * +! * eps: epsilon * +! * thick: thickness * +! * layers: number of layers * +! * nper: number of periodic layers * +! * * +! ****************************************************************** + + use quanc8_mod + use fint1_mod + use fint2_mod + use fint3_mod + use param_mod + + 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 + + double precision :: error, flag + double precision :: u1, u2, ut, y + integer :: ie, nofu + dimension error(3), flag(3) + +! write (*,*) 'queels:' +! write (*,*) 'eps: ', size(eps) +! write (*,*) 'thick: ', size(thick) + + f = 0.0d0 + if (x <= 0.0d0) then + return + endif + ru = facru * x + ccoef = cospsi**2 / x + ut = ccoef - bcoef + u1 = dabs(ut) + u2 = ccoef + bcoef + if (ut > 0.0d0) then + call quanc8(fint1, 0.0d0, u1, aerr, rerr, y, error(1), nofu, flag(1), eps, thick, layers, nper) + f = y + else + flag(1) = 0.0d0 + endif + if (u2 > u1) then + call quanc8(fint2, u1, u2, aerr, rerr, y, error(2), nofu, flag(2), eps, thick, layers, nper) + f = f + y + else + flag(2) = 0.0d0 + endif + if (dabs(acoef) > x * (1.0d0 - elleps) * bcoef) then + um = dsqrt(ccoef / x / (1.0d0 - elleps) + bcoef**2 / acoef) + if (um > u2) then + call quanc8(fint3, u2, um, aerr, rerr, y, error(3), nofu, flag(3), eps, thick, layers, nper) + f = f + y + endif + 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 + do ie = 1, 3 + 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 + f = (2 / 3.141592653589793238d0)**2 * f + return +end subroutine queels +end module queels_mod diff --git a/tests/SetEpsTestCasesFromScratch/WFW/surlos.f90 b/tests/SetEpsTestCasesFromScratch/WFW/surlos.f90 new file mode 100644 index 0000000..66552e9 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/surlos.f90 @@ -0,0 +1,174 @@ +module surlos_mod +contains + +logical function zero(z) + double complex, intent(in) :: z + zero = (dble(z) == 0.0d0) .and. (dimag(z) == 0.0d0) +end function zero + +double precision function surlos(dk, eps, thick, layers, nper) + +! ****************************************************************** +! * * +! * eels surface loss function for an arbitrary multilayered target* +! * * +! * dk: * +! * eps: epsilon * +! * thick: thickness * +! * layers: number of layers * +! * nper: number of periodic layers * +! * * +! ****************************************************************** + + 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, 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 + + epsmac = epsilon(1.0d0) + argmin = dsqrt(2 * epsmac) + epsmac = epsmac / 2 + argmax = 0.5d0 * dlog(2 / epsmac) + +! write (*,*) 'surlos:' +! write (*,*) 'thick: ', size(thick) +! write (*,*) 'eps: ', size(eps) + + lmax = size(eps) + allocate (arg(lmax)) + lstart = layers - nper + 1 + static = .true. + skip = .false. + + do n = 1, layers + arg(n) = dk * thick(n) + if (arg(n) > argmax .or. zero(eps(n))) then + csi = eps(n) + skip = .true. + exit + endif + static = .not. (n >= lstart .and. arg(n) > argmin) + enddo + + if (.not. skip) then +! *** periodic continued fraction, period = nper + + do ! Do loop is only for the option to skip out of it. + if (nper <= 1) then + 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 + + n = lstart + + exit + enddo + endif ! .not. skip + +! *** backward algorithm + + do + n = n - 1 + if (n <= 0) then + a = csi + 1.0d0 + if (zero(a)) then + surlos = 2 / epsmac + else + surlos = dimag(-2 / a) + endif + return + endif ! n <= 0 + + if (arg(n) /= 0.0d0) then + tn = dtanh(arg(n)) + b = eps(n) + csi * tn + if (zero(b)) then + surlos = 0.0d0 + return + endif + csi = eps(n) * (csi + tn * eps(n)) / b + endif + enddo +end function surlos +end module surlos_mod diff --git a/tests/SetEpsTestCasesFromScratch/WFW/usurlo.f90 b/tests/SetEpsTestCasesFromScratch/WFW/usurlo.f90 new file mode 100644 index 0000000..222946c --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/usurlo.f90 @@ -0,0 +1,23 @@ +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 out) :: dq + double precision, intent(in out) :: wn + + if ((wn .GT. 0) .AND. (dq .GT. 0)) then + write(*,*) 'hello, here is the user loss function' + endif + usurlo = 1.0d0 + return +end function usurlo diff --git a/tests/SetEpsTestCasesFromScratch/WFW/version.inc b/tests/SetEpsTestCasesFromScratch/WFW/version.inc new file mode 100644 index 0000000..453b451 --- /dev/null +++ b/tests/SetEpsTestCasesFromScratch/WFW/version.inc @@ -0,0 +1,2 @@ +! version of eels-boson + character (len = *), parameter :: version = '1.0.0' -- GitLab