diff --git a/tests/SetEpsTestCasesFromScratch/f77/Makefile b/tests/SetEpsTestCasesFromScratch/f77/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..246d1222da1336f125659b20d1856a4905d1f98d
--- /dev/null
+++ b/tests/SetEpsTestCasesFromScratch/f77/Makefile
@@ -0,0 +1,20 @@
+# the fortran compiler
+FC = gfortran
+
+# the options for the fortran compiler
+FFLAGS = -g -gdwarf-2 -fbounds-check -fcheck=all -O0
+
+# the implicit rule for compiling for files
+%.o: %.for ; $(FC) $(FFLAGS) -c -o $@ $<
+
+all: eels
+
+eels: change_working_dir.o eels.for
+ $(FC) $(FFLAGS) -o eels eels.for change_working_dir.o
+
+clean:
+ rm -f *.o
+ rm -rf *.dSYM
+ rm -f eels
+
+.PHONY: all clean
diff --git a/tests/SetEpsTestCasesFromScratch/f77/README b/tests/SetEpsTestCasesFromScratch/f77/README
new file mode 100644
index 0000000000000000000000000000000000000000..d4839456790eb5b2cb8d560f56616f7d3432a629
--- /dev/null
+++ b/tests/SetEpsTestCasesFromScratch/f77/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/f77/change_working_dir.for b/tests/SetEpsTestCasesFromScratch/f77/change_working_dir.for
new file mode 100755
index 0000000000000000000000000000000000000000..6841ba89f000e319e3c37e61c774fbb097ad5784
--- /dev/null
+++ b/tests/SetEpsTestCasesFromScratch/f77/change_working_dir.for
@@ -0,0 +1,25 @@
+ subroutine change_working_dir()
+
+C This routine gets the first argument of the commandline and takes it
+C as the path to change the working directory
+C used intrinsic routines:
+C iarg returns the number of commandline arguments without the program cname.
+C chdir changes the directory and returns 0 on success.
+C trim removes blanks from strings.
+
+ character(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'
+ endif
+ endif
+
+ return
+ end
diff --git a/tests/SetEpsTestCasesFromScratch/f77/doRun.sh b/tests/SetEpsTestCasesFromScratch/f77/doRun.sh
new file mode 100755
index 0000000000000000000000000000000000000000..f5a8e515567c4f4d1e9bde28d4efd172018c2af5
--- /dev/null
+++ b/tests/SetEpsTestCasesFromScratch/f77/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/f77/eels.for b/tests/SetEpsTestCasesFromScratch/f77/eels.for
new file mode 100644
index 0000000000000000000000000000000000000000..768825cb831b12aad6368388adb7ebaf3656f791
--- /dev/null
+++ b/tests/SetEpsTestCasesFromScratch/f77/eels.for
@@ -0,0 +1,877 @@
+ PROGRAM EELS ABTI0030
+C ****************************************************************** ABTI0031
+C * * ABTI0032
+C * COMPUTE THE CLASSICAL EELS SPECTRUM OF AN ARBITRARY PLANE- * ABTI0033
+C * STATIFIED MEDIUM MADE FROM ISOTROPIC MATERIALS IN SPECULAR * ABTI0034
+C * GEOMETRY USING THE DIELECTRIC THEORY OF EELS. * ABTI0035
+C * * ABTI0036
+C ****************************************************************** ABTI0037
+ PARAMETER(LMAX=100,JMAX=300,NT=5) ABTI0038
+ LOGICAL RATION,USER ABTI0039
+ CHARACTER CONTRL*10,NAME*10,COMMEN*72 ABTI0040
+ DOUBLE PRECISION A,ACOEF,AERR,ALPHA,ARGMIN,ARGMAX,B,BCOEF,BETA, ABTI0041
+ ,C1,C2,CCOEF,COSPSI,DLIMF,DW,DX,E0,ELLEPS,ENER,EPSINF,EPSMAC, ABTI0042
+ ,FACRU,F,F0,F1,FPIC,FUN,OSC,PHIA,PHIB,PI,PREFAC,PSIA,PSII,QRAT, ABTI0043
+ ,RERR,RU,SINPSI,T,TANPSI,TABLE,THETA,THICK,UM,WIDT,WMIN,WMAX,WN, ABTI0044
+ ,WPIC,X,XMIN,XMAX,Z,Z1,Z2 ABTI0045
+ COMPLEX EPS ABTI0046
+ DIMENSION COMMEN(2),EPSINF(LMAX),NOS(LMAX),OSC(3,JMAX),TABLE(NT) ABTI0047
+ DIMENSION NAME(LMAX)
+ COMMON/PARAM/ACOEF,BCOEF,CCOEF,ELLEPS,COSPSI,SINPSI,TANPSI, ABTI0048
+ 1 RU,UM,DLIMF,WN,USER,RATION ABTI0049
+ COMMON/MULAYR/EPS(LMAX),THICK(LMAX),ARGMIN,ARGMAX,EPSMAC, ABTI0050
+ 2 LAYERS,NPER ABTI0051
+ EXTERNAL FUN ABTI0052
+ QRAT(X) = (1.0D0+X*(BETA+C1*X))/((1.0D0+X*(BETA+C2*X))* ABTI0053
+ * (1.0D0+ALPHA*X)**2) ABTI0054
+ DATA AERR/0.0D0/,RERR/1.0D-06/,F/0.0D0/,F1/0.0D0/ ABTI0055
+
+C **** log modification start
+ open(unit = 99, file = 'seteps.log')
+C **** log modification end
+
+C ABTI0056
+C *** MACHINE-DEPENDENT CONSTANTS ABTI0057
+C *** EPSMAC + 1.0 = EPSMAC , COSH(ARGMIN) = 1.0 , TANH(ARGMAX) = 1.0 ABTI0058
+C ABTI0059
+ PI = 4.0D0*ATAN(1.0D0) ABTI0060
+ EPSMAC = 1.0D0 ABTI0061
+ 1 EPSMAC = EPSMAC/2.0D0 ABTI0062
+ IF(1.0D0+EPSMAC.GT.1.0D0) GOTO 1 ABTI0063
+ ARGMIN = SQRT(2.0D0*EPSMAC) ABTI0064
+ ARGMAX = 0.5D0*LOG(2.0D0/EPSMAC) ABTI0065
+C ABTI0066
+C *** READ SPECTROMETER PARAMETERS ABTI0067
+C ABTI0068
+ call change_working_dir()
+ OPEN(UNIT=11,FILE='EELSIN') ABTI0069
+C IMPACT ENERGY (EV) ABTI0070
+ READ(11,*) E0 ABTI0071
+C INCIDENCE ANGLE (%) ABTI0072
+ READ(11,*) THETA ABTI0073
+C ANGULAR APERTURES OF THE ELLIPTIC DETECTOR (%) ABTI0074
+ READ(11,*) PHIA ABTI0075
+ READ(11,*) PHIB ABTI0076
+C ENERGY-LOSS INTERVAL AND STEP SIZE (CM**-1) ABTI0077
+ READ(11,*) WMIN ABTI0078
+ READ(11,*) WMAX ABTI0079
+ READ(11,*) DW ABTI0080
+C COMMENT LINES ABTI0081
+ READ(11,'(A72)') (COMMEN(K),K=1,2) ABTI0082
+C ABTI0083
+ WRITE(6,100) E0,THETA,PHIA,PHIB,WMIN,WMAX,DW,(COMMEN(K),K=1,2) ABTI0084
+ IF(PHIA.LE.0.0D0 .OR. PHIB.LE.0.0D0) STOP '*** WRONG INPUT ***' ABTI0085
+ IF(E0.LE.0.0D0 .OR. THETA+PHIA.GE.90.0D0) STOP '*** BAD INPUT ***'ABTI0086
+ DLIMF = 0.0D0 ABTI0087
+ RATION = .FALSE. ABTI0088
+C ABTI0089
+C *** READ TARGET SPECIFICATIONS ABTI0090
+C ABTI0091
+ READ(11,*) LAYERS,NPER ABTI0092
+ USER = LAYERS.EQ.0 ABTI0093
+ IF(USER) GOTO 15 ABTI0094
+ NEPS = LAYERS ABTI0095
+ IF(NPER.EQ.-1) THEN ABTI0096
+ NEPS = LAYERS+1 ABTI0097
+ NPER = 1 ABTI0098
+ WRITE(6,*) 'THE SUBSTRATE IS A ANISOTROPIC UNIAXIAL MATERIAL' ABTI0099
+ ENDIF ABTI0100
+ IF(LAYERS.LT.0 .OR. NPER.LT.1 .OR. NPER.GT.LAYERS) STOP ABTI0101
+ . '*** INVALID TARGET SPECIFICATIONS ***' ABTI0102
+ IF(NEPS.GT.LMAX) STOP '*** TOO MANY LAYERS ***' ABTI0103
+ WRITE(6,101) LAYERS,NPER ABTI0104
+ IF(LAYERS.EQ.1) RATION = .TRUE. ABTI0105
+ LSTART = LAYERS-NPER+1 ABTI0106
+ JOS = 0 ABTI0107
+ DO 10 L=1,NEPS ABTI0108
+ IF(L.LE.LAYERS) READ(11,102) NAME(L),THICK(L) ABTI0109
+ READ(11,*) EPSINF(L),NOS(L) ABTI0110
+ WRITE(6,103)
+C IF(L.EQ.LSTART) WRITE(6,103) ABTI0111
+ IF(NOS(L).LE.0) THEN ABTI0112
+ IF(L.LE.LAYERS) WRITE(6,104) L,NAME(L),THICK(L),EPSINF(L) ABTI0113
+ IF(L.GT.LAYERS) WRITE(6,105) EPSINF(L) ABTI0114
+ ELSE ABTI0115
+ DO 5 J=1,NOS(L) ABTI0116
+ JOS = JOS+1 ABTI0117
+ IF(JOS.GT.JMAX) STOP '*** TOO MANY OSCILLATORS ***' ABTI0118
+ READ(11,*) (OSC(K,JOS),K=1,3) ABTI0119
+ IF((J==NOS(L)/2+1).AND.(NOS(L).GT.1)) THEN
+ WRITE(6,114)
+ WRITE(6,115)
+ ENDIF
+ IF(J.EQ.1) THEN ABTI0120
+ IF(L.LE.LAYERS) WRITE(6,104) L,NAME(L),THICK(L), ABTI0121
+ , EPSINF(L),(OSC(I,JOS),I=1,3) ABTI0122
+ IF(L.GT.LAYERS) WRITE(6,105) ABTI0123
+ . EPSINF(L),(OSC(I,JOS),I=1,3) ABTI0124
+ ELSE ABTI0125
+ WRITE(6,106) (OSC(I,JOS),I=1,3) ABTI0126
+ ENDIF ABTI0127
+ 5 CONTINUE ABTI0128
+ ENDIF ABTI0129
+ 10 CONTINUE ABTI0130
+ WRITE(6,*) ABTI0131
+ READ(11,102,END=15) CONTRL ABTI0132
+ IF(CONTRL.EQ.'IMAGE') THEN ABTI0133
+C *** IMAGE-CHARGE SCREENING FACTOR ABTI0134
+ IF(LAYERS.EQ.1 .AND. NEPS.EQ.2) THEN ABTI0135
+ DLIMF = SQRT(EPSINF(1)*EPSINF(2)) ABTI0136
+ ELSE ABTI0137
+ DLIMF = EPSINF(1) ABTI0138
+ ENDIF ABTI0139
+ DLIMF = (DLIMF-1.0D0)/(DLIMF+1.0D0) ABTI0140
+ ENDIF ABTI0141
+C ABTI0142
+C *** INITIALIZE CONSTANTS ABTI0143
+C ABTI0144
+ 15 OPEN(UNIT=12,FILE='EELSOU') ABTI0145
+ WRITE(12,107) E0,THETA,PHIA,PHIB,COMMEN(1) ABTI0146
+ NW = 1+INT((WMAX-WMIN)/DW) ABTI0147
+ NOUT = 1+NW/20 ABTI0148
+ ENER = 8.065D+03*E0 ABTI0149
+ PSIA = PHIA/180.0D0*PI ABTI0150
+ PSII = THETA/180.0D0*PI ABTI0151
+ COSPSI = COS(PSII) ABTI0152
+ SINPSI = SIN(PSII) ABTI0153
+ TANPSI = TAN(PSII) ABTI0154
+ PREFAC = SQRT(2.555D+05/E0)/(1.37D+02*COSPSI) ABTI0155
+ FACRU = PSIA/COSPSI*SQRT(0.2624664D0*E0) ABTI0156
+ ELLEPS = (1.0D0-PHIA/PHIB)*(1.0D0+PHIA/PHIB) ABTI0157
+ ACOEF = SINPSI**2+ELLEPS*COSPSI**2 ABTI0158
+ BCOEF = SINPSI*COSPSI ABTI0159
+ IF(DLIMF.GT.0.0D0) THEN ABTI0160
+ RATION = .FALSE. ABTI0161
+ WRITE(6,108) DLIMF ABTI0162
+C *** DLIMF : HALF THE LENGTH UNIT IMPOSED BY THE IMAGE FORCE ABTI0163
+ DLIMF = 1.80D0*DLIMF/(E0*COSPSI**2) ABTI0164
+ ENDIF ABTI0165
+ IF(.NOT.RATION) GOTO 35 ABTI0166
+C ABTI0167
+C *** SET UP COEFFICIENTS FOR THE RATIONAL APPROXIMATION TO THE INTEGRALABTI0168
+C ABTI0169
+ WRITE(6,*) '==> SET UP A RATIONAL APPROXIMATION TO THE INTEGRAL' ABTI0170
+ CALL QUANC8(FUN,0.0D0,PI/2.0D0,AERR,RERR,ALPHA,C1,NOFU,C2) ABTI0171
+ ALPHA = (2.0D0/PI)**2*ALPHA ABTI0172
+ C1 = 2.0D0/PI/SQRT(1.0D0-ELLEPS)*SINPSI*ALPHA**2 ABTI0173
+ IF(C1.GT.0.99D0) GOTO 30 ABTI0174
+ C2 = 3.0D0*ALPHA**2/(1.0D0-C1) ABTI0175
+ C1 = C1*C2 ABTI0176
+ XMIN = WMIN/(2.0D0*ENER*PSIA) ABTI0177
+ XMAX = WMAX/(2.0D0*ENER*PSIA) ABTI0178
+ IF(XMIN.LE.0.0D0) XMIN = 0.0D0 ABTI0179
+ DX = MAX(0.02D0,(XMAX-XMIN)/NT) ABTI0180
+ Z1 = 0.0D0 ABTI0181
+ Z2 = 0.0D0 ABTI0182
+ DO 20 I=1,NT ABTI0183
+ X = XMIN+I*DX ABTI0184
+ CALL QUEELS(X,F,AERR,RERR,FACRU) ABTI0185
+ TABLE(I) = F ABTI0186
+ F = F*(1.0D0+ALPHA*X)**2 ABTI0187
+ IF(ABS(C2*F-C1).LT.C2*RERR) GOTO 20 ABTI0188
+ Z = (1.0D0-F)/(C2*F-C1) ABTI0189
+ IF(Z.LE.0.0D0) GOTO 20 ABTI0190
+ Z1 = Z1+X*Z*(X**2-Z) ABTI0191
+ Z2 = Z2+(X*Z)**2 ABTI0192
+ 20 CONTINUE ABTI0193
+ IF(Z2.EQ.0.0D0) GOTO 30 ABTI0194
+ BETA = Z1/Z2 ABTI0195
+ Z = 0.0D0 ABTI0196
+ DO 25 I=1,NT ABTI0197
+ X = XMIN+I*DX ABTI0198
+ Z = Z+(TABLE(I)-QRAT(X))**2 ABTI0199
+ 25 CONTINUE ABTI0200
+ Z = SQRT(Z)/NT ABTI0201
+ IF(Z.GT.5.0D-03) GOTO 30 ABTI0202
+ WRITE(6,109) ALPHA,C1,C2,BETA,Z ABTI0203
+ GOTO 35 ABTI0204
+ 30 WRITE(6,*) '===> CANNOT DO IT' ABTI0205
+ RATION = .FALSE. ABTI0206
+C ABTI0207
+C *** LOOP OVER THE ENERGY LOSSES ABTI0208
+C ABTI0209
+ 35 WRITE(6,110) ABTI0210
+
+C **** 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, *)
+C **** log modification end
+
+ DO 50 IW=1,NW ABTI0211
+ F0 = F1 ABTI0212
+ F1 = F ABTI0213
+ F = 0.0D0 ABTI0214
+ WN = WMIN+(IW-1)*DW ABTI0215
+ IF(WN.LT.0.0D0) GOTO 45 ABTI0216
+ IF(WN.EQ.0.0D0) GOTO 40 ABTI0217
+ IF(.NOT.USER) CALL SETEPS(NEPS,NOS,OSC,JOS,EPSINF,WN,NAME) ABTI0218
+ X = WN/(2.0D0*ENER*PSIA) ABTI0219
+ IF(RATION) THEN ABTI0220
+ F = QRAT(X)*AIMAG(-2.0/(1.0+EPS(1))) ABTI0221
+ ELSE ABTI0222
+ CALL QUEELS(X,F,AERR,RERR,FACRU) ABTI0223
+ ENDIF ABTI0224
+ F = PREFAC*F/WN ABTI0225
+ 40 WRITE(12,111) WN,F ABTI0226
+C *** LOCALIZE A PEAK USING A PARABOLIC INTERPOLATION ABTI0227
+ IF(IW.LT.3) GOTO 45 ABTI0228
+ IF(F1-F0.LE.AERR) GOTO 45 ABTI0229
+ IF(F1-F.LE.AERR) GOTO 45 ABTI0230
+ A = (F1-F0)+(F1-F) ABTI0231
+ IF(A.LE.4.0D0*RERR*F1) GOTO 45 ABTI0232
+ B = 0.5D0*(F1-F0+3.0D0*(F1-F)) ABTI0233
+ T = B/A ABTI0234
+ WPIC = WN-T*DW ABTI0235
+ FPIC = F+0.5D0*B*T ABTI0236
+ WIDT = SQRT(8.0D0*FPIC/A)*DW ABTI0237
+ WRITE(6,112) WPIC,FPIC,WIDT ABTI0238
+ 45 IF(MOD(IW,NOUT).EQ.0) WRITE(6,113) 100.0*IW/NW,WN,F ABTI0239
+ 50 CONTINUE ABTI0240
+ CLOSE(UNIT=11) ABTI0241
+ CLOSE(UNIT=12) ABTI0242
+
+C **** log modification start
+ close (unit = 99)
+C **** log modification end
+
+ STOP ABTI0243
+ 100 FORMAT(' PROGRAM EELS (MARCH 1990)'/' E0 =',F6.2,' EV , THETA =',ABTI0244
+ , F5.1,'% , PHIA =',F5.2,'% , PHIB =',F5.2,'%'/ ABTI0245
+ . ' ENERGY LOSSES FROM',G11.4,' TO',G11.4,' , STEP =',G11.4,ABTI0246
+ , ' CM**-1'//(1X,A72)) ABTI0247
+ 101 FORMAT(I3,' LAYER(S), NPER =',I2//' L',2X,'MATERIAL',7X, ABTI0248
+ , 'THICKNESS',5X,'EPSINF',4X,'WTO , WP',5X,'Q',7X,'GAM/WTO') ABTI0249
+ 102 FORMAT(A10,D15.5) ABTI0250
+ 103 FORMAT(1X,72('-')) ABTI0251
+ 104 FORMAT(1X,I3,2X,A10,G15.3,F10.4,F12.4,F10.4,F9.4) ABTI0252
+ 105 FORMAT(31X,F10.4,F12.4,F10.4,F9.4) ABTI0253
+ 106 FORMAT(41X,F12.4,F10.4,F9.4) ABTI0254
+ 107 FORMAT('E0 =',F6.2,' THETA =',F5.1,' PHIA =',F5.2,' PHIB =',F5.2/ ABTI0255
+ . A72) ABTI0256
+ 108 FORMAT(' ==> ELECTRON ATTRACTED BY AN IMAGE CHARGE =',F6.3) ABTI0257
+ 109 FORMAT(5X,'ALPHA =',F9.4,4X,'C1 =',F9.4,4X,'C2 =',F9.4,4X, ABTI0258
+ , 'BETA =',F9.4/5X,'ACCURACY =',E9.2) ABTI0259
+ 110 FORMAT(//' RUN (%) WN (CM**-1) PCL(WN) (CM) |', ABTI0260
+ , ' PEAK LOCATION AMPLITUDE WIDTH') ABTI0261
+ 111 FORMAT(2E15.7) ABTI0262
+ 112 FORMAT(40X,F10.2,D12.4,F10.2) ABTI0263
+ 113 FORMAT(2X,F5.1,3X,F11.3,D14.5) ABTI0264
+ 114 FORMAT(45X,'WLO , WP',5X,'Q',7X,'GAM/WLO')
+ 115 FORMAT(45X,28('-'))
+ END ABTI0265
+ SUBROUTINE QUEELS(X,F,AERR,RERR,FACRU) ABTI0266
+C ****************************************************************** ABTI0267
+C * * ABTI0268
+C * PERFORM Q-SPACE INTEGRATION FOR COMPUTING THE EELS SPECTRUM OF * ABTI0269
+C * A ISOTROPIC TARGET USING POLAR COORDINATES. * ABTI0270
+C * * ABTI0271
+C * X IS THE DIMENSIONLESS ENERGY LOSS HBAR*OMEGA/(2*E0*PHIA) * ABTI0272
+C * AERR AND RERR ARE THE DESIRED ABSOLUTE AND RELATIVE ACCURACIES * ABTI0273
+C * FACRU*X IS THE UNITS OF WAVEVECTORS OMEGA/V_PERPENDICULAR * ABTI0274
+C * F IS THE Q-INTEGRAL MULTIPLIED BY (2/PI)**2 * ABTI0275
+C * * ABTI0276
+C ****************************************************************** ABTI0277
+ LOGICAL RATION,USER ABTI0278
+ DOUBLE PRECISION ACOEF,AERR,BCOEF,CCOEF,COSPSI,DLIMF,ELLEPS, ABTI0279
+ ,ERROR,F,FACRU,FINT1,FINT2,FINT3,FLAG,RERR,RU,SINPSI,U1,U2,UM,UT, ABTI0280
+ ,TANPSI,WN,X,Y ABTI0281
+ DIMENSION ERROR(3),FLAG(3) ABTI0282
+ COMMON/PARAM/ACOEF,BCOEF,CCOEF,ELLEPS,COSPSI,SINPSI,TANPSI, ABTI0283
+ , RU,UM,DLIMF,WN,USER,RATION ABTI0284
+ EXTERNAL FINT1,FINT2,FINT3 ABTI0285
+ F = 0.0D0 ABTI0286
+ IF(X.LE.0.0D0) RETURN ABTI0287
+ RU = FACRU*X ABTI0288
+ CCOEF = COSPSI**2/X ABTI0289
+ UT = CCOEF-BCOEF ABTI0290
+ U1 = ABS(UT) ABTI0291
+ U2 = CCOEF+BCOEF ABTI0292
+ IF(UT.GT.0.0D0) THEN ABTI0293
+ CALL QUANC8(FINT1,0.0D0,U1,AERR,RERR,Y,ERROR(1),NOFU,FLAG(1)) ABTI0294
+ F = Y ABTI0295
+ ELSE ABTI0296
+ FLAG(1) = 0.0D0 ABTI0297
+ ENDIF ABTI0298
+ IF(U2.GT.U1) THEN ABTI0299
+ CALL QUANC8(FINT2,U1,U2,AERR,RERR,Y,ERROR(2),NOFU,FLAG(2)) ABTI0300
+ F = F+Y ABTI0301
+ ELSE ABTI0302
+ FLAG(2) = 0.0D0 ABTI0303
+ ENDIF ABTI0304
+ IF(ABS(ACOEF).GT.X*(1.0D0-ELLEPS)*BCOEF) THEN ABTI0305
+ UM = SQRT(CCOEF/X/(1.0D0-ELLEPS)+BCOEF**2/ACOEF) ABTI0306
+ IF(UM.GT.U2) THEN ABTI0307
+ CALL QUANC8(FINT3,U2,UM,AERR,RERR,Y,ERROR(3),NOFU,FLAG(3)) ABTI0308
+ F = F+Y ABTI0309
+ ENDIF ABTI0310
+ IF(UM.LT.U1) THEN ABTI0311
+ CALL QUANC8(FINT3,UM,U1,AERR,RERR,Y,ERROR(3),NOFU,FLAG(3)) ABTI0312
+ F = F-Y ABTI0313
+ ENDIF ABTI0314
+ ELSE ABTI0315
+ FLAG(3) = 0.0D0 ABTI0316
+ ENDIF ABTI0317
+ ISTOP = 0 ABTI0318
+ DO 5 IE=1,3 ABTI0319
+ IF(FLAG(IE).EQ.0.0D0) GOTO 5 ABTI0320
+ WRITE(6,100) IE,FLAG(IE),ERROR(IE) ABTI0321
+ IF(FLAG(IE)-AINT(FLAG(IE)).GT.0.5D-02) ISTOP = ISTOP+1 ABTI0322
+ 5 CONTINUE ABTI0323
+ IF(ISTOP.GT.0) STOP '*** EXECUTION ABORTED ***' ABTI0324
+ F = (2.0D0/3.141592653589793238D0)**2*F ABTI0325
+ RETURN ABTI0326
+ 100 FORMAT(' +++ FLAG(',I1,') =',F10.5,', ERROR =',D11.4,' +++') ABTI0327
+ END ABTI0328
+ SUBROUTINE QUANC8(FUN,A,B,ABSERR,RELERR,RESULT,ERREST,NOFUN,FLAG) ABTI0329
+C ABTI0330
+C ESTIMATE THE INTEGRAL OF FUN(X) FROM A TO B ABTI0331
+C TO A USER PROVIDED TOLERANCE. ABTI0332
+C AN AUTOMATIC ADAPTIVE ROUTINE BASED ON ABTI0333
+C THE 8-PANEL NEWTON-COTES RULE (G. FORSYTHE ET AL, 1977, P. 92) ABTI0334
+C ABTI0335
+C INPUT .. ABTI0336
+C ABTI0337
+C FUN THE NAME OF THE INTEGRAND FUNCTION SUBPROGRAM FUN(X). ABTI0338
+C A THE LOWER LIMIT OF INTEGRATION. ABTI0339
+C B THE UPPER LIMIT OF INTEGRATION.(B MAY BE LESS THAN A.) ABTI0340
+C RELERR A RELATIVE ERROR TOLERANCE. (SHOULD BE NON-NEGATIVE) ABTI0341
+C ABSERR AN ABSOLUTE ERROR TOLERANCE. (SHOULD BE NON-NEGATIVE) ABTI0342
+C ABTI0343
+C OUTPUT .. ABTI0344
+C ABTI0345
+C RESULT AN APPROXIMATION TO THE INTEGRAL HOPEFULLY SATISFYING THE ABTI0346
+C LEAST STRINGENT OF THE TWO ERROR TOLERANCES. ABTI0347
+C ERREST AN ESTIMATE OF THE MAGNITUDE OF THE ACTUAL ERROR. ABTI0348
+C NOFUN THE NUMBER OF FUNCTION VALUES USED IN CALCULATION OF RESULT.ABTI0349
+C FLAG A RELIABILITY INDICATOR. IF FLAG IS ZERO, THEN RESULT ABTI0350
+C PROBABLY SATISFIES THE ERROR TOLERANCE. IF FLAG IS ABTI0351
+C XXX.YYY , THEN XXX = THE NUMBER OF INTERVALS WHICH HAVE ABTI0352
+C NOT CONVERGED AND 0.YYY = THE FRACTION OF THE INTERVAL ABTI0353
+C LEFT TO DO WHEN THE LIMIT ON NOFUN WAS APPROACHED. ABTI0354
+C ABTI0355
+ DOUBLE PRECISION FUN, A, B, ABSERR, RELERR, RESULT, ERREST, FLAG ABTI0356
+ DOUBLE PRECISION W0,W1,W2,W3,W4,AREA,X0,F0,STONE,STEP,COR11,TEMP ABTI0357
+ DOUBLE PRECISION QPREV,QNOW,QDIFF,QLEFT,ESTERR,TOLERR ABTI0358
+ DOUBLE PRECISION QRIGHT(31),F(16),X(16),FSAVE(8,30),XSAVE(8,30) ABTI0359
+ DOUBLE PRECISION DABS,DMAX1 ABTI0360
+ INTEGER NOFUN ABTI0361
+ INTEGER LEVMIN,LEVMAX,LEVOUT,NOMAX,NOFIN,LEV,NIM,I,J ABTI0362
+C ABTI0363
+C *** STAGE 1 *** GENERAL INITIALIZATION ABTI0364
+C SET CONSTANTS. ABTI0365
+C ABTI0366
+ LEVMIN = 1 ABTI0367
+ LEVMAX = 30 ABTI0368
+ LEVOUT = 6 ABTI0369
+ NOMAX = 5000 ABTI0370
+ NOFIN = NOMAX - 8*(LEVMAX-LEVOUT+2**(LEVOUT+1)) ABTI0371
+C ABTI0372
+C TROUBLE WHEN NOFUN REACHES NOFIN ABTI0373
+C ABTI0374
+ W0 = 3956.0D0 / 14175.0D0 ABTI0375
+ W1 = 23552.0D0 / 14175.0D0 ABTI0376
+ W2 = -3712.0D0 / 14175.0D0 ABTI0377
+ W3 = 41984.0D0 / 14175.0D0 ABTI0378
+ W4 = -18160.0D0 / 14175.0D0 ABTI0379
+C ABTI0380
+C INITIALIZE RUNNING SUMS TO ZERO. ABTI0381
+C ABTI0382
+ FLAG = 0.0D0 ABTI0383
+ RESULT = 0.0D0 ABTI0384
+ COR11 = 0.0D0 ABTI0385
+ ERREST = 0.0D0 ABTI0386
+ AREA = 0.0D0 ABTI0387
+ NOFUN = 0 ABTI0388
+ IF (A .EQ. B) RETURN ABTI0389
+C ABTI0390
+C *** STAGE 2 *** INITIALIZATION FOR FIRST INTERVAL ABTI0391
+C ABTI0392
+ LEV = 0 ABTI0393
+ NIM = 1 ABTI0394
+ X0 = A ABTI0395
+ X(16) = B ABTI0396
+ QPREV = 0.0D0 ABTI0397
+ F0 = FUN(X0) ABTI0398
+ STONE = (B - A) / 16.0D0 ABTI0399
+ X(8) = (X0 + X(16)) / 2.0D0 ABTI0400
+ X(4) = (X0 + X(8)) / 2.0D0 ABTI0401
+ X(12) = (X(8) + X(16)) / 2.0D0 ABTI0402
+ X(2) = (X0 + X(4)) / 2.0D0 ABTI0403
+ X(6) = (X(4) + X(8)) / 2.0D0 ABTI0404
+ X(10) = (X(8) + X(12)) / 2.0D0 ABTI0405
+ X(14) = (X(12) + X(16)) / 2.0D0 ABTI0406
+ DO 25 J = 2, 16, 2 ABTI0407
+ F(J) = FUN(X(J)) ABTI0408
+ 25 CONTINUE ABTI0409
+ NOFUN = 9 ABTI0410
+C ABTI0411
+C *** STAGE 3 *** CENTRAL CALCULATION ABTI0412
+C REQUIRES QPREV,X0,X2,X4,...,X16,F0,F2,F4,...,F16. ABTI0413
+C CALCULATES X1,X3,...X15, F1,F3,...F15,QLEFT,QRIGHT,QNOW,QDIFF,AREA. ABTI0414
+C ABTI0415
+ 30 X(1) = (X0 + X(2)) / 2.0D0 ABTI0416
+ F(1) = FUN(X(1)) ABTI0417
+ DO 35 J = 3, 15, 2 ABTI0418
+ X(J) = (X(J-1) + X(J+1)) / 2.0D0 ABTI0419
+ F(J) = FUN(X(J)) ABTI0420
+ 35 CONTINUE ABTI0421
+ NOFUN = NOFUN + 8 ABTI0422
+ STEP = (X(16) - X0) / 16.0D0 ABTI0423
+ QLEFT = (W0*(F0 + F(8)) + W1*(F(1)+F(7)) + W2*(F(2)+F(6)) ABTI0424
+ 1 + W3*(F(3)+F(5)) + W4*F(4)) * STEP ABTI0425
+ QRIGHT(LEV+1)=(W0*(F(8)+F(16))+W1*(F(9)+F(15))+W2*(F(10)+F(14)) ABTI0426
+ 1 + W3*(F(11)+F(13)) + W4*F(12)) * STEP ABTI0427
+ QNOW = QLEFT + QRIGHT(LEV+1) ABTI0428
+ QDIFF = QNOW - QPREV ABTI0429
+ AREA = AREA + QDIFF ABTI0430
+C ABTI0431
+C *** STAGE 4 *** INTERVAL CONVERGENCE TEST ABTI0432
+C ABTI0433
+ ESTERR = DABS(QDIFF) / 1023.0D0 ABTI0434
+ TOLERR = DMAX1(ABSERR,RELERR*DABS(AREA)) * (STEP/STONE) ABTI0435
+ IF (LEV .LT. LEVMIN) GO TO 50 ABTI0436
+ IF (LEV .GE. LEVMAX) GO TO 62 ABTI0437
+ IF (NOFUN .GT. NOFIN) GO TO 60 ABTI0438
+ IF (ESTERR .LE. TOLERR) GO TO 70 ABTI0439
+C ABTI0440
+C *** STAGE 5 *** NO CONVERGENCE ABTI0441
+C LOCATE NEXT INTERVAL. ABTI0442
+C ABTI0443
+ 50 NIM = 2*NIM ABTI0444
+ LEV = LEV+1 ABTI0445
+C ABTI0446
+C STORE RIGHT HAND ELEMENTS FOR FUTURE USE. ABTI0447
+C ABTI0448
+ DO 52 I = 1, 8 ABTI0449
+ FSAVE(I,LEV) = F(I+8) ABTI0450
+ XSAVE(I,LEV) = X(I+8) ABTI0451
+ 52 CONTINUE ABTI0452
+C ABTI0453
+C ASSEMBLE LEFT HAND ELEMENTS FOR IMMEDIATE USE. ABTI0454
+C ABTI0455
+ QPREV = QLEFT ABTI0456
+ DO 55 I = 1, 8 ABTI0457
+ J = -I ABTI0458
+ F(2*J+18) = F(J+9) ABTI0459
+ X(2*J+18) = X(J+9) ABTI0460
+ 55 CONTINUE ABTI0461
+ GO TO 30 ABTI0462
+C ABTI0463
+C *** STAGE 6 *** TROUBLE SECTION ABTI0464
+C NUMBER OF FUNCTION VALUES IS ABOUT TO EXCEED LIMIT. ABTI0465
+C ABTI0466
+ 60 NOFIN = 2*NOFIN ABTI0467
+ LEVMAX = LEVOUT ABTI0468
+ FLAG = FLAG + (B - X0) / (B - A) ABTI0469
+ GO TO 70 ABTI0470
+C ABTI0471
+C CURRENT LEVEL IS LEVMAX. ABTI0472
+C ABTI0473
+ 62 FLAG = FLAG + 1.0D0 ABTI0474
+C ABTI0475
+C *** STAGE 7 *** INTERVAL CONVERGED ABTI0476
+C ADD CONTRIBUTIONS INTO RUNNING SUMS. ABTI0477
+C ABTI0478
+ 70 RESULT = RESULT + QNOW ABTI0479
+ ERREST = ERREST + ESTERR ABTI0480
+ COR11 = COR11 + QDIFF / 1023.0D0 ABTI0481
+C ABTI0482
+C LOCATE NEXT INTERVAL. ABTI0483
+C ABTI0484
+ 72 IF (NIM .EQ. 2*(NIM/2)) GO TO 75 ABTI0485
+ NIM = NIM/2 ABTI0486
+ LEV = LEV-1 ABTI0487
+ GO TO 72 ABTI0488
+ 75 NIM = NIM + 1 ABTI0489
+ IF (LEV .LE. 0) GO TO 80 ABTI0490
+C ABTI0491
+C ASSEMBLE ELEMENTS REQUIRED FOR THE NEXT INTERVAL. ABTI0492
+C ABTI0493
+ QPREV = QRIGHT(LEV) ABTI0494
+ X0 = X(16) ABTI0495
+ F0 = F(16) ABTI0496
+ DO 78 I = 1, 8 ABTI0497
+ F(2*I) = FSAVE(I,LEV) ABTI0498
+ X(2*I) = XSAVE(I,LEV) ABTI0499
+ 78 CONTINUE ABTI0500
+ GO TO 30 ABTI0501
+C ABTI0502
+C *** STAGE 8 *** FINALIZE AND RETURN ABTI0503
+C ABTI0504
+ 80 RESULT = RESULT + COR11 ABTI0505
+C ABTI0506
+C MAKE SURE ERREST NOT LESS THAN ROUNDOFF LEVEL. ABTI0507
+C ABTI0508
+ IF (ERREST .EQ. 0.0D0) RETURN ABTI0509
+ 82 TEMP = DABS(RESULT) + ERREST ABTI0510
+ IF (TEMP .NE. DABS(RESULT)) RETURN ABTI0511
+ ERREST = 2.0D0*ERREST ABTI0512
+ GO TO 82 ABTI0513
+ END ABTI0514
+ SUBROUTINE SETEPS(NEPS,NOS,OSC,JOS,EPSINF,WN,NAME) ABTI0515
+C ****************************************************************** ABTI0516
+C * * ABTI0517
+C * SET UP LONG-WAVELENGTH DIELECTRIC FUNCTIONS OF THE LAYERS FOR * ABTI0518
+C * THE PRESENT FREQUENCY WN (IN CM**-1) * ABTI0519
+C * * ABTI0520
+C ****************************************************************** ABTI0521
+ PARAMETER(LMAX=100) ABTI0522
+ CHARACTER NAME*10
+ DOUBLE PRECISION ARGMIN,ARGMAX,EPSINF,EPSMAC,OSC,THICK,WN,X ABTI0523
+ COMPLEX DENO,EPS,DENO1,DENO2 ABTI0524
+ DIMENSION NOS(NEPS),OSC(3,JOS),EPSINF(NEPS),NAME(NEPS) ABTI0525
+ COMMON/MULAYR/EPS(LMAX),THICK(LMAX),ARGMIN,ARGMAX,EPSMAC, ABTI0526
+ , LAYERS,NPER ABTI0527
+ J = 0 ABTI0528
+ DO 2 L=1,NEPS
+ M = NOS(L)/2
+ DENO1 = CMPLX(1.0D0,0.0D0)
+ DENO2 = CMPLX(1.0D0,0.0D0)
+ IF(NOS(L).GT.1) THEN
+ DO 1 K=1,M
+ J = J+1
+ DENO1 = DENO1*(OSC(1,J+M)**2-WN**2-CMPLX(0.0,WN*OSC(3,J+M)))
+ DENO2 = DENO2*(OSC(1,J )**2-WN**2-CMPLX(0.0,WN*OSC(3,J )))
+ 1 CONTINUE
+ EPS(L) = EPSINF(L)*DENO1/DENO2
+ ELSE IF(NAME(L).EQ.'metal') THEN
+ J = J+1
+ EPS(L) = -OSC(1,J)**2/( WN**2+CMPLX(0.0D0,WN*OSC(3,J)) )
+ ELSE
+ EPS(L) = EPSINF(L)
+ J = J+1
+ X = WN/OSC(1,J)
+ DENO = X*CMPLX(X,OSC(3,J))
+ IF(OSC(2,J).GE.0.0D0) DENO = 1.0-DENO
+ IF(ABS(DENO).EQ.0.0) DENO = EPSMAC
+ EPS(L) = EPS(L)+OSC(2,J)/DENO
+ ENDIF
+ 2 CONTINUE ABTI0538
+ IF(NEPS.EQ.LAYERS+1) THEN ABTI0540
+C THE SUBSTRATE IS A ANISOTROPIC UNIAXIAL MATERIAL ABTI0541
+ EPS(LAYERS) = SQRT(EPS(LAYERS)*EPS(LAYERS+1)) ABTI0542
+ IF(AIMAG(EPS(LAYERS)).LT.0.0) EPS(LAYERS) = -EPS(LAYERS) ABTI0543
+ ENDIF ABTI0544
+
+C **** log modification start
+ write (99, '(30g15.7)') wn, (eps(j), j = 1, neps)
+C **** log modification end
+
+ RETURN ABTI0545
+ END ABTI0546
+ DOUBLE PRECISION FUNCTION FUN(PHI) ABTI0547
+C ****************************************************************** ABTI0548
+C * * ABTI0549
+C * INTEGRAND OF THE EXPRESSION OF THE 1ST ORDER TERM IN THE * ABTI0550
+C * EXPANSION OF THE EELS INTEGRAL FOR A HOMOGENEOUS TARGET. * ABTI0551
+C * * ABTI0552
+C ****************************************************************** ABTI0553
+ LOGICAL USER,RATION ABTI0554
+ DOUBLE PRECISION ACOEF,BCOEF,CCOEF,COSPSI,DLIMF,ELLEPS,PHI,RU, ABTI0555
+ ,SINPHI,SINPSI,TANPSI,UM,WN ABTI0556
+ COMMON/PARAM/ACOEF,BCOEF,CCOEF,ELLEPS,COSPSI,SINPSI,TANPSI, ABTI0557
+ , RU,UM,DLIMF,WN,USER,RATION ABTI0558
+ SINPHI = SIN(PHI) ABTI0559
+ FUN = SQRT((1.0D0-ELLEPS+ELLEPS*SINPHI**2)*(1.0D0-SINPSI*SINPHI)* ABTI0560
+ * (1.0D0+SINPSI*SINPHI)) ABTI0561
+ RETURN ABTI0562
+ END ABTI0563
+ DOUBLE PRECISION FUNCTION FINT1(U) ABTI0564
+C ****************************************************************** ABTI0565
+C * * ABTI0566
+C * INTEGRATION OVER THE AZIMUTAL ANGLE FROM 0.0 TO PI * ABTI0567
+C * * ABTI0568
+C ****************************************************************** ABTI0569
+ LOGICAL RATION,USER ABTI0570
+ DOUBLE PRECISION ACOEF,BCOEF,CCOEF,COSPSI,DEN,DIF,DLIMF,E,ELLEPS, ABTI0571
+ ,PI,ROM,ROP,RU,SINPSI,SUM,SURLOS,U,UM,USURLO,T,TANPSI,WN ABTI0572
+ COMMON/PARAM/ACOEF,BCOEF,CCOEF,ELLEPS,COSPSI,SINPSI,TANPSI, ABTI0573
+ , RU,UM,DLIMF,WN,USER,RATION ABTI0574
+ DATA PI/3.141592653589793238D0/ ABTI0575
+ IF(U.EQ.0.0D0) THEN ABTI0576
+ FINT1 = 0.0D0 ABTI0577
+ RETURN ABTI0578
+ ENDIF ABTI0579
+ E = TANPSI*U ABTI0580
+ ROM = (1.0D0-E)**2+U**2 ABTI0581
+ ROP = (1.0D0+E)**2+U**2 ABTI0582
+ SUM = ROP+ROM ABTI0583
+ ROM = SQRT(ROM) ABTI0584
+ ROP = SQRT(ROP) ABTI0585
+ DIF = ROP-ROM ABTI0586
+ DEN = SQRT((2.0D0-DIF)*(2.0D0+DIF))*ROP*ROM ABTI0587
+ FINT1 = PI*U**2*(4.0D0*SUM-DIF**2*(SUM-ROP*ROM))/DEN**3 ABTI0588
+ IF(RATION) RETURN ABTI0589
+ IF(USER) THEN ABTI0590
+ FINT1 = FINT1*USURLO(RU*U,WN) ABTI0591
+ ELSE ABTI0592
+ FINT1 = FINT1*SURLOS(RU*U) ABTI0593
+ IF(DLIMF.GT.0.0D0) THEN ABTI0594
+ T = RU*U*DLIMF ABTI0595
+ FINT1 = FINT1*(1.D0+T*LOG(T/(T+0.26D0)))**2/(1.D0+1.40D0*T) ABTI0596
+ ENDIF ABTI0597
+ ENDIF ABTI0598
+ RETURN ABTI0599
+ END ABTI0600
+ DOUBLE PRECISION FUNCTION FINT2(U) ABTI0601
+C ****************************************************************** ABTI0602
+C * * ABTI0603
+C * INTEGRATION OVER THE AZIMUTAL ANGLE FROM 0.0 TO PHI < PI * ABTI0604
+C * * ABTI0605
+C ****************************************************************** ABTI0606
+ LOGICAL RATION,USER ABTI0607
+ DOUBLE PRECISION A,ARG,B,B2,C,CCOEF,COSPSI,DLIMF,ELLEPS,PHI, ABTI0608
+ ,PHINT,RU,SINPSI,SURLOS,U,UM,USURLO,T,TANPSI,WN,X ABTI0609
+ COMMON/PARAM/A,B,CCOEF,ELLEPS,COSPSI,SINPSI,TANPSI,RU,UM,DLIMF, ABTI0610
+ , WN,USER,RATION ABTI0611
+ IF(U.EQ.0.0D0) THEN ABTI0612
+ FINT2 = 0.0D0 ABTI0613
+ RETURN ABTI0614
+ ENDIF ABTI0615
+ B2 = B**2 ABTI0616
+ C = (1.0D0-ELLEPS)*(COSPSI*U)**2+(B-CCOEF)*(B+CCOEF) ABTI0617
+ IF(ABS(A*C).GT.1.0D-03*B2) THEN ABTI0618
+ X = (B-SQRT(B2-A*C))/A ABTI0619
+ ELSE ABTI0620
+ X = A*C/B2 ABTI0621
+ X = 0.5D0*C*(1.D0+0.25D0*X*(1.D0+0.5D0*X*(1.D0+0.625D0*X)))/B ABTI0622
+ ENDIF ABTI0623
+ ARG = X/U ABTI0624
+ IF(ABS(ARG).GT.1.0D0) ARG = SIGN(1.0D0,ARG) ABTI0625
+ PHI = ACOS(ARG) ABTI0626
+ FINT2 = PHINT(PHI,TANPSI,U) ABTI0627
+ IF(RATION) RETURN ABTI0628
+ IF(USER) THEN ABTI0629
+ FINT2 = FINT2*USURLO(RU*U,WN) ABTI0630
+ ELSE ABTI0631
+ FINT2 = FINT2*SURLOS(RU*U) ABTI0632
+ IF(DLIMF.GT.0.0D0) THEN ABTI0633
+ T = RU*U*DLIMF ABTI0634
+ FINT2 = FINT2*(1.D0+T*LOG(T/(T+0.26D0)))**2/(1.D0+1.40D0*T) ABTI0635
+ ENDIF ABTI0636
+ ENDIF ABTI0637
+ RETURN ABTI0638
+ END ABTI0639
+ DOUBLE PRECISION FUNCTION FINT3(U) ABTI0640
+C ****************************************************************** ABTI0641
+C * * ABTI0642
+C * INTEGRATION OVER THE AZIMUTAL ANGLE FROM PHI1 > 0 TO PHI2 < PI * ABTI0643
+C * * ABTI0644
+C ****************************************************************** ABTI0645
+ LOGICAL RATION,USER ABTI0646
+ DOUBLE PRECISION A,ARG,B,CCOEF,COSPSI,DLIMF,ELLEPS,PHI1,PHI2, ABTI0647
+ ,PHINT,SINPSI,SURLOS,RAC,RU,U,UM,USURLO,T,TANPSI,WN ABTI0648
+ COMMON/PARAM/A,B,CCOEF,ELLEPS,COSPSI,SINPSI,TANPSI,RU,UM,DLIMF, ABTI0649
+ , WN,USER,RATION ABTI0650
+ IF(U.EQ.0.0D0) THEN ABTI0651
+ FINT3 = 0.0D0 ABTI0652
+ RETURN ABTI0653
+ ENDIF ABTI0654
+ RAC = SIGN(1.0D0,A)*COSPSI*SQRT((1.0D0-ELLEPS)*A*(UM-U)*(UM+U)) ABTI0655
+ ARG = (B-RAC)/(U*A) ABTI0656
+ IF(ABS(ARG).GT.1.0D0) ARG = SIGN(1.0D0,ARG) ABTI0657
+ PHI2 = ACOS(ARG) ABTI0658
+ FINT3 = PHINT(PHI2,TANPSI,U) ABTI0659
+ ARG = (B+RAC)/(U*A) ABTI0660
+ IF(ABS(ARG).GT.1.0D0) ARG = SIGN(1.0D0,ARG) ABTI0661
+ PHI1 = ACOS(ARG) ABTI0662
+ FINT3 = FINT3-PHINT(PHI1,TANPSI,U) ABTI0663
+ IF(RATION) RETURN ABTI0664
+ IF(USER) THEN ABTI0665
+ FINT3 = FINT3*USURLO(RU*U,WN) ABTI0666
+ ELSE ABTI0667
+ FINT3 = FINT3*SURLOS(RU*U) ABTI0668
+ IF(DLIMF.GT.0.0D0) THEN ABTI0669
+ T = RU*U*DLIMF ABTI0670
+ FINT3 = FINT3*(1.D0+T*LOG(T/(T+0.26D0)))**2/(1.D0+1.40D0*T) ABTI0671
+ ENDIF ABTI0672
+ ENDIF ABTI0673
+ RETURN ABTI0674
+ END ABTI0675
+ DOUBLE PRECISION FUNCTION PHINT(PHI,A,U) ABTI0676
+C ****************************************************************** ABTI0677
+C * * ABTI0678
+C * EVALUATE THE INTEGRAL FROM ZERO TO PHI OF * ABTI0679
+C * * ABTI0680
+C * U 2 * ABTI0681
+C * ( ----------------------------- ) DPHI * ABTI0682
+C * 2 2 * ABTI0683
+C * (1 - A * U * COS(PHI)) + U * ABTI0684
+C * * ABTI0685
+C * FOR 0 <= PHI <= PI , U >= 0 AND A >= 0 * ABTI0686
+C * * ABTI0687
+C ****************************************************************** ABTI0688
+ DOUBLE PRECISION A,AI,AR,BI,BR,C,CPR,D,E,ESR,PI,PHI,QR,RI,RM,ROOT,ABTI0689
+ ,RP,RR,S,SPR,TM,TP,U,U2,X,ZETA,ZETAI,ZETAR,ZR ABTI0690
+ PI = 3.141592653589793238D0 ABTI0691
+ C = COS(PHI) ABTI0692
+ S = SIN(PHI) ABTI0693
+ U2 = U**2 ABTI0694
+ E = A*U ABTI0695
+ IF(U.LT.1.0D0 .AND. E.LT.1.0D-02*(1.0D0+U2)) GOTO 5 ABTI0696
+ RM = SQRT((1.0D0-E)**2+U2) ABTI0697
+ RP = SQRT((1.0D0+E)**2+U2) ABTI0698
+ ROOT = SQRT(RM*RP) ABTI0699
+ TM = 0.5D0*ATAN2(U,1.0D0-E) ABTI0700
+ TP = 0.5D0*ATAN2(U,1.0D0+E) ABTI0701
+ CPR = COS(TM+TP) ABTI0702
+ SPR = SIN(TM+TP) ABTI0703
+ IF(C.GE.0.0D0) GOTO 2 ABTI0704
+ IF(ABS(S).GT.1.0D-07) GOTO 1 ABTI0705
+ RR = -PI/ROOT*CPR ABTI0706
+ RI = PI/ROOT*SPR ABTI0707
+ AR = -RR+U*RI ABTI0708
+ AI = -RI-U*RR ABTI0709
+ GOTO 4 ABTI0710
+ 1 X = (1.0D0-C)/S ABTI0711
+ GOTO 3 ABTI0712
+ 2 X = S/(1.0D0+C) ABTI0713
+ 3 ZETA = SQRT(RM/RP) ABTI0714
+ ZETAR = -ZETA*SIN(TM-TP) ABTI0715
+ ZETAI = ZETA*COS(TM-TP) ABTI0716
+ BR = 0.5D0*LOG(((ZETAR+X)**2+ZETAI**2)/((ZETAR-X)**2+ZETAI**2)) ABTI0717
+ BI = ATAN2(ZETAI,ZETAR+X)-ATAN2(ZETAI,ZETAR-X) ABTI0718
+ RR = -1.0D0/ROOT*(BR*SPR-BI*CPR) ABTI0719
+ RI = -1.0D0/ROOT*(BI*SPR+BR*CPR) ABTI0720
+ D = E*S/((1.0D0-E*C)**2+U2) ABTI0721
+ AR = D*(1.0D0-E*C)-RR+U*RI ABTI0722
+ AI = -D*U-RI-U*RR ABTI0723
+ 4 QR = 1.0D0/(RM*RP)*(AR*(CPR-SPR)*(CPR+SPR)+2.0D0*AI*CPR*SPR) ABTI0724
+ PHINT = 0.5D0*(RI/U-QR) ABTI0725
+ RETURN ABTI0726
+ 5 ZR = 1.0D0+U2 ABTI0727
+ ESR = E/ZR ABTI0728
+ PHINT = U2/ZR**2*((( (4.D0/3.D0)*(2.0D0+C**2)*S*(5.0D0-3.0D0*U2)* ABTI0729
+ * ESR + (PHI+C*S)*(5.D0-U2))*ESR + 4.0D0*S)*ESR + PHI) ABTI0730
+ RETURN ABTI0731
+ END ABTI0732
+ DOUBLE PRECISION FUNCTION SURLOS(DK) ABTI0733
+C ****************************************************************** ABTI0734
+C * * ABTI0735
+C * EELS SURFACE LOSS FUNCTION FOR AN ARBITRARY MULTILAYERED TARGET* ABTI0736
+C * * ABTI0737
+C ****************************************************************** ABTI0738
+ PARAMETER(LMAX=100) ABTI0739
+ LOGICAL STATIC,ZERO ABTI0740
+ DOUBLE PRECISION ARG,ARGMIN,ARGMAX,CN,CNM1,D,DK,EPSMAC,SN,SNM1,TN ABTI0741
+ COMPLEX A,B,CSI,EPS,PNM2,PNM1,PN,PP,QNM2,QNM1,QN,QP,Z ABTI0742
+ DIMENSION ARG(LMAX) ABTI0743
+ COMMON/MULAYR/EPS(LMAX),D(LMAX),ARGMIN,ARGMAX,EPSMAC,LAYERS,NPER ABTI0744
+ ZERO(Z) = REAL(Z).EQ.0.0 .AND. AIMAG(Z).EQ.0.0 ABTI0745
+ LSTART = LAYERS-NPER+1 ABTI0746
+ STATIC = .TRUE. ABTI0747
+ N = 1 ABTI0748
+ 1 ARG(N) = DK*D(N) ABTI0749
+ IF(ARG(N).GT.ARGMAX .OR. ZERO(EPS(N))) GOTO 10 ABTI0750
+ IF(N.GE.LSTART .AND. ARG(N).GT.ARGMIN) STATIC = .FALSE. ABTI0751
+ N = N+1 ABTI0752
+ IF(N.LE.LAYERS) GOTO 1 ABTI0753
+C ABTI0754
+C *** PERIODIC CONTINUED FRACTION, PERIOD = NPER ABTI0755
+C ABTI0756
+ IF(NPER.GT.1) GOTO 2 ABTI0757
+ CSI = EPS(LAYERS) ABTI0758
+ GOTO 9 ABTI0759
+ 2 IF(STATIC) GOTO 5 ABTI0760
+ CN = COSH(ARG(LSTART)) ABTI0761
+ SN = SINH(ARG(LSTART)) ABTI0762
+ PNM1 = 1.0 ABTI0763
+ PN = CN ABTI0764
+ PP = EPS(LSTART)*SN ABTI0765
+ QNM1 = 0.0 ABTI0766
+ QN = SN/EPS(LSTART) ABTI0767
+ QP = PN ABTI0768
+ DO 3 N=LSTART+1,LAYERS ABTI0769
+ CNM1 = CN ABTI0770
+ SNM1 = SN ABTI0771
+ CN = COSH(ARG(N)) ABTI0772
+ SN = SINH(ARG(N)) ABTI0773
+ A = EPS(N)*SN ABTI0774
+ PP = CN*PP+A*PN ABTI0775
+ QP = CN*QP+A*QN ABTI0776
+ B = (EPS(N-1)/EPS(N))*(SN/SNM1) ABTI0777
+ A = CNM1*B+CN ABTI0778
+ PNM2 = PNM1 ABTI0779
+ PNM1 = PN ABTI0780
+ QNM2 = QNM1 ABTI0781
+ QNM1 = QN ABTI0782
+ PN = A*PNM1-B*PNM2 ABTI0783
+ QN = A*QNM1-B*QNM2 ABTI0784
+ 3 CONTINUE ABTI0785
+ IF(ZERO(QN)) GOTO 4 ABTI0786
+ A = 0.5*(PN-QP)/QN ABTI0787
+ B = SQRT(A**2+PP/QN) ABTI0788
+ PN = A-PN/QN ABTI0789
+ IF(ABS(PN+B).GT.ABS(PN-B)) B = -B ABTI0790
+ CSI = A+B ABTI0791
+ GOTO 9 ABTI0792
+ 4 A = QP-PN ABTI0793
+ IF(ZERO(A)) GOTO 12 ABTI0794
+ CSI = PP/A ABTI0795
+ GOTO 9 ABTI0796
+C ABTI0797
+C *** SMALL-DK LIMIT OF THE PERIODIC TAIL ABTI0798
+C ABTI0799
+ 5 PN = 0.0 ABTI0800
+ QN = 0.0 ABTI0801
+ DO 6 N=LSTART,LAYERS ABTI0802
+ PN = PN+D(N)*EPS(N) ABTI0803
+ QN = QN+D(N)/EPS(N) ABTI0804
+ 6 CONTINUE ABTI0805
+ IF(ZERO(QN)) GOTO 12 ABTI0806
+ CSI = SQRT(PN/QN) ABTI0807
+ IF(AIMAG(CSI).GT.0.0) GOTO 9 ABTI0808
+ IF(AIMAG(CSI).LT.0.0) THEN ABTI0809
+ CSI = -CSI ABTI0810
+ ELSE ABTI0811
+ IF(REAL(QN).LT.0.0) CSI = -CSI ABTI0812
+ ENDIF ABTI0813
+ 9 N = LSTART ABTI0814
+ GOTO 11 ABTI0815
+ 10 CSI = EPS(N) ABTI0816
+C ABTI0817
+C *** BACKWARD ALGORITHM ABTI0818
+C ABTI0819
+ 11 N = N-1 ABTI0820
+ IF(N.LE.0) GOTO 15 ABTI0821
+ IF(ARG(N).EQ.0.0D0) GOTO 11 ABTI0822
+ T = TANH(ARG(N)) ABTI0823
+ B = EPS(N)+CSI*T ABTI0824
+ IF(ZERO(B)) GOTO 13 ABTI0825
+ CSI = EPS(N)*(CSI+T*EPS(N))/B ABTI0826
+ GOTO 11 ABTI0827
+ 12 N = LSTART ABTI0828
+ 13 IF(N.LE.1) GOTO 14 ABTI0829
+ N = N-1 ABTI0830
+ CSI = EPS(N)/TANH(ARG(N)) ABTI0831
+ GOTO 11 ABTI0832
+ 14 SURLOS = 0.0D0 ABTI0833
+ GOTO 17 ABTI0834
+ 15 A = CSI+1.0 ABTI0835
+ IF(ZERO(A)) GOTO 16 ABTI0836
+ SURLOS = AIMAG(-2.0/A) ABTI0837
+ GOTO 17 ABTI0838
+ 16 SURLOS = 2.0D0/EPSMAC ABTI0839
+ 17 RETURN
+ END ABTI0841
+ DOUBLE PRECISION FUNCTION USURLO(DQ,WN) ABTI0842
+C ****************************************************************** ABTI0843
+C * * ABTI0844
+C * USER-SUPPLIED DIELECTRIC SURFACE LOSS FUNCTION AIMAG(G(DQ,WN)) * ABTI0845
+C * INPUT ARGUMENTS : * ABTI0846
+C * DQ : MODULUS OF THE TWO-DIMENSIONAL SURFACE WAVE VECTOR * ABTI0847
+C * (ANGSTROEM**-1) * ABTI0848
+C * WN : FREQUENCY (CM**-1) * ABTI0849
+C * * ABTI0850
+C ****************************************************************** ABTI0851
+ DOUBLE PRECISION DQ,WN ABTI0852
+ USURLO = 1.0D0 ABTI0853
+ RETURN ABTI0854
+ END ABTI0855
diff --git a/tests/SetEpsTestCasesFromScratch/f77/epsLog/seteps001.log b/tests/SetEpsTestCasesFromScratch/f77/epsLog/seteps001.log
new file mode 100644
index 0000000000000000000000000000000000000000..cdd1cfa4e76e9347f7b76b7104d26940d549e9d9
Binary files /dev/null and b/tests/SetEpsTestCasesFromScratch/f77/epsLog/seteps001.log differ
diff --git a/tests/SetEpsTestCasesFromScratch/f77/epsLog/seteps004.log b/tests/SetEpsTestCasesFromScratch/f77/epsLog/seteps004.log
new file mode 100644
index 0000000000000000000000000000000000000000..a851395f53c8fb965813dff69a69c879f8039f6a
Binary files /dev/null and b/tests/SetEpsTestCasesFromScratch/f77/epsLog/seteps004.log differ
diff --git a/tests/SetEpsTestCasesFromScratch/f77/epsLog/seteps006.log b/tests/SetEpsTestCasesFromScratch/f77/epsLog/seteps006.log
new file mode 100644
index 0000000000000000000000000000000000000000..58cff9d956f25913fa66096d764d8464719edfb9
Binary files /dev/null and b/tests/SetEpsTestCasesFromScratch/f77/epsLog/seteps006.log differ