From 3d054191ae17387f240fc2d3cc33f8a04bb4e96e Mon Sep 17 00:00:00 2001
From: kamischi <karl-michael.schindler@web.de>
Date: Tue, 1 Mar 2022 22:45:18 +0100
Subject: [PATCH] code improvements
- rename d to thick for uniform naming.
- remove unused thick and nper in subroutine seteps.
- use dcmplx for double precision
---
source/f90/doeels.f90 | 2 +-
source/f90/fint1.f90 | 9 ++++-----
source/f90/fint2.f90 | 9 ++++-----
source/f90/fint3.f90 | 9 ++++-----
source/f90/phint.f90 | 1 +
source/f90/quanc8.f90 | 14 +++++++-------
source/f90/queels.f90 | 14 +++++++-------
source/f90/seteps.f90 | 17 ++++++++---------
source/f90/surlos.f90 | 12 ++++++------
9 files changed, 42 insertions(+), 45 deletions(-)
diff --git a/source/f90/doeels.f90 b/source/f90/doeels.f90
index 273be7d..1d44ab6 100644
--- a/source/f90/doeels.f90
+++ b/source/f90/doeels.f90
@@ -177,7 +177,7 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
! 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, thick, layers, nper, mode)
+ if (.not. user) call seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
x = wn / (2.0d0 * ener * psia)
if (rational) then
diff --git a/source/f90/fint1.f90 b/source/f90/fint1.f90
index b1a1f7e..98702bb 100644
--- a/source/f90/fint1.f90
+++ b/source/f90/fint1.f90
@@ -1,6 +1,6 @@
module fint1_mod
contains
-double precision function fint1(u, eps, d, layers, nper, eps_size)
+double precision function fint1(u, eps, thick, layers, nper, eps_size)
! ******************************************************************
! * *
@@ -13,9 +13,8 @@ double precision function fint1(u, eps, d, layers, nper, eps_size)
implicit none
double precision, intent(in) :: u
-
- double precision, intent(in) :: d(eps_size)
double complex, intent(in) :: eps(eps_size)
+ double precision, intent(in) :: thick(eps_size)
integer, intent(in) :: layers, nper, eps_size
logical :: rational, user
@@ -29,7 +28,7 @@ double precision function fint1(u, eps, d, layers, nper, eps_size)
data pi / 3.141592653589793238d0 /
! write (*,*) 'fint1:'
-! write (*,*) 'd: ', size(d)
+! write (*,*) 'thick: ', size(thick)
! write (*,*) 'eps: ', size(eps)
if (u == 0.0d0) then
@@ -52,7 +51,7 @@ double precision function fint1(u, eps, d, layers, nper, eps_size)
if (user) then
fint1 = fint1 * usurlo(ru * u, wn)
else
- fint1 = fint1 * surlos(ru * u, eps, d, layers, nper)
+ fint1 = fint1 * surlos(ru * u, eps, thick, layers, nper)
if (dlimf > 0.0d0) then
t = ru * u * dlimf
fint1 = fint1 * (1.d0 + t * log(t / (t + 0.26d0)))**2 / (1.d0 + 1.40d0 * t)
diff --git a/source/f90/fint2.f90 b/source/f90/fint2.f90
index c6c2404..ab5562f 100644
--- a/source/f90/fint2.f90
+++ b/source/f90/fint2.f90
@@ -1,6 +1,6 @@
module fint2_mod
contains
-double precision function fint2(u, eps, d, layers, nper, eps_size)
+double precision function fint2(u, eps, thick, layers, nper, eps_size)
! ******************************************************************
! * *
@@ -13,9 +13,8 @@ double precision function fint2(u, eps, d, layers, nper, eps_size)
implicit none
double precision, intent(in) :: u
-
- double precision, intent(in) :: d(eps_size)
double complex, intent(in) :: eps(eps_size)
+ double precision, intent(in) :: thick(eps_size)
integer, intent(in) :: layers, nper, eps_size
logical :: rational, user
@@ -26,7 +25,7 @@ double precision function fint2(u, eps, d, layers, nper, eps_size)
ru, um, dlimf, wn, user, rational
! write (*,*) 'fint2:'
-! write (*,*) 'd: ', size(d)
+! write (*,*) 'thick: ', size(thick)
! write (*,*) 'eps: ', size(eps)
if (u == 0.0d0) then
@@ -53,7 +52,7 @@ double precision function fint2(u, eps, d, layers, nper, eps_size)
if (user) then
fint2 = fint2 * usurlo(ru * u, wn)
else
- fint2 = fint2 * surlos(ru * u, eps, d, layers, nper)
+ fint2 = fint2 * surlos(ru * u, eps, thick, layers, nper)
if (dlimf > 0.0d0) then
t = ru * u * dlimf
fint2 = fint2 * (1.d0 + t * log(t / (t + 0.26d0)))**2 / (1.d0 + 1.40d0 * t)
diff --git a/source/f90/fint3.f90 b/source/f90/fint3.f90
index 2c1a586..9ae66c9 100644
--- a/source/f90/fint3.f90
+++ b/source/f90/fint3.f90
@@ -1,6 +1,6 @@
module fint3_mod
contains
-double precision function fint3(u, eps, d, layers, nper, eps_size)
+double precision function fint3(u, eps, thick, layers, nper, eps_size)
! ******************************************************************
! * *
@@ -13,9 +13,8 @@ double precision function fint3(u, eps, d, layers, nper, eps_size)
implicit none
double precision, intent(in) :: u
-
- double precision, intent(in) :: d(eps_size)
double complex, intent(in) :: eps(eps_size)
+ double precision, intent(in) :: thick(eps_size)
integer, intent(in) :: layers, nper, eps_size
logical :: rational, user
@@ -26,7 +25,7 @@ double precision function fint3(u, eps, d, layers, nper, eps_size)
ru, um, dlimf, wn, user, rational
! write (*,*) 'fint3:'
-! write (*,*) 'd: ', size(d)
+! write (*,*) 'thick: ', size(thick)
! write (*,*) 'eps: ', size(eps)
if (u == 0.0d0) then
@@ -46,7 +45,7 @@ double precision function fint3(u, eps, d, layers, nper, eps_size)
if (user) then
fint3 = fint3 * usurlo(ru * u, wn)
else
- fint3 = fint3 * surlos(ru * u, eps, d, layers, nper)
+ fint3 = fint3 * surlos(ru * u, eps, thick, layers, nper)
if (dlimf > 0.0d0) then
t = ru * u * dlimf
fint3 = fint3 * (1.d0 + t * log(t / (t + 0.26d0)))**2 / (1.d0 + 1.40d0 * t)
diff --git a/source/f90/phint.f90 b/source/f90/phint.f90
index 271bc37..a089b8a 100644
--- a/source/f90/phint.f90
+++ b/source/f90/phint.f90
@@ -40,6 +40,7 @@ double precision function phint(phi, a, u)
root = sqrt(rm * rp)
cpr = cos(tm + tp)
spr = sin(tm + tp)
+ x = 0.0d0 ! ensure initialization
if (c >= 0.0d0) then
x = s / (1.0d0 + c)
elseif (abs(s) > 1.0d-07) then
diff --git a/source/f90/quanc8.f90 b/source/f90/quanc8.f90
index 6496e7c..e3c60b1 100644
--- a/source/f90/quanc8.f90
+++ b/source/f90/quanc8.f90
@@ -1,6 +1,6 @@
module quanc8_mod
contains
-subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, d, layers, nper)
+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.
@@ -41,7 +41,7 @@ subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, d
external fun
- double precision, intent(in) :: d(:)
+ double precision, intent(in) :: thick(:)
double complex, intent(in) :: eps(:)
integer, intent(in) :: layers, nper
@@ -56,7 +56,7 @@ subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, d
! set constants.
! write (*,*) 'quanc8:'
-! write (*,*) 'd: ', size(d)
+! write (*,*) 'thick: ', size(thick)
! write (*,*) 'eps: ', size(eps)
levmin = 1
@@ -90,7 +90,7 @@ subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, d
x0 = a
x(16) = b
qprev = 0.0d0
- f0 = fun(x0, eps, d, layers, nper, size(eps))
+ f0 = fun(x0, eps, thick, layers, nper, size(eps))
stone = (b - a) / 16.0d0
x(8) = (x0 + x(16)) / 2.0d0
x(4) = (x0 + x(8)) / 2.0d0
@@ -100,7 +100,7 @@ subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, d
x(10) = (x(8) + x(12)) / 2.0d0
x(14) = (x(12) + x(16)) / 2.0d0
do j = 2, 16, 2
- f(j) = fun(x(j), eps, d, layers, nper, size(eps))
+ f(j) = fun(x(j), eps, thick, layers, nper, size(eps))
enddo
nofun = 9
@@ -111,10 +111,10 @@ subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, d
! calculates x1, x3, ...x15, f1, f3, ...f15, qleft, qright, qnow, qdiff, area.
x(1) = (x0 + x(2)) / 2.0d0
- f(1) = fun(x(1), eps, d, layers, nper, size(eps))
+ 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.0d0
- f(j) = fun(x(j), eps, d, layers, nper, size(eps))
+ f(j) = fun(x(j), eps, thick, layers, nper, size(eps))
enddo
nofun = nofun + 8
step = (x(16) - x0) / 16.0d0
diff --git a/source/f90/queels.f90 b/source/f90/queels.f90
index aefd51f..2a7fe02 100644
--- a/source/f90/queels.f90
+++ b/source/f90/queels.f90
@@ -1,6 +1,6 @@
module queels_mod
contains
-subroutine queels(x, f, aerr, rerr, facru, eps, d, layers, nper)
+subroutine queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
! ******************************************************************
! * *
@@ -27,7 +27,7 @@ subroutine queels(x, f, aerr, rerr, facru, eps, d, layers, nper)
double precision, intent(in out) :: rerr
double precision, intent(in) :: facru
- double precision, intent(in) :: d(:)
+ double precision, intent(in) :: thick(:)
double complex, intent(in) :: eps(:)
integer, intent(in) :: layers, nper
@@ -43,7 +43,7 @@ subroutine queels(x, f, aerr, rerr, facru, eps, d, layers, nper)
! write (*,*) 'queels:'
! write (*,*) 'eps: ', size(eps)
-! write (*,*) 'd: ', size(d)
+! write (*,*) 'thick: ', size(thick)
f = 0.0d0
if (x <= 0.0d0) then
@@ -55,13 +55,13 @@ subroutine queels(x, f, aerr, rerr, facru, eps, d, layers, nper)
u1 = abs(ut)
u2 = ccoef + bcoef
if (ut > 0.0d0) then
- call quanc8(fint1, 0.0d0, u1, aerr, rerr, y, error(1), nofu, flag(1), eps, d, layers, nper)
+ 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, d, layers, nper)
+ 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
@@ -69,11 +69,11 @@ subroutine queels(x, f, aerr, rerr, facru, eps, d, layers, nper)
if (abs(acoef) > x * (1.0d0 - elleps) * bcoef) then
um = sqrt(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, d, layers, nper)
+ 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, d, layers, nper)
+ call quanc8(fint3, um, u1, aerr, rerr, y, error(3), nofu, flag(3), eps, thick, layers, nper)
f = f - y
endif
else
diff --git a/source/f90/seteps.f90 b/source/f90/seteps.f90
index 83a744a..b890f44 100644
--- a/source/f90/seteps.f90
+++ b/source/f90/seteps.f90
@@ -1,6 +1,6 @@
module seteps_mod
contains
-subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, thick, layers, nper, mode)
+subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
! ******************************************************************
! * *
@@ -19,8 +19,7 @@ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, thick, layers, nper, mo
character (len=10), intent(in) :: mode
double complex, intent(in out) :: eps(:)
- double precision, intent(in) :: thick(:)
- integer, intent(in) :: layers, nper
+ integer, intent(in) :: layers
double precision :: argmin, argmax, epsmac, x
double complex :: deno, nomi
@@ -39,8 +38,8 @@ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, thick, layers, nper, mo
j = 0
do l = 1, neps
m = nos(l)/2
- nomi = cmplx(1.0d0, 0.0d0)
- deno = cmplx(1.0d0, 0.0d0)
+ nomi = dcmplx(1.0d0, 0.0d0)
+ deno = dcmplx(1.0d0, 0.0d0)
if (mode == 'kurosawa') then
do k = 1, m
j = j + 1
@@ -49,15 +48,15 @@ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, thick, layers, nper, mo
! Furthermore the first term can be rewritten as
! (osc(1, j + m) - wn) * (osc(1, j + m) + wn)
! nomi = nomi * cmplx(osc(1, j + m)**2 - wn**2, - wn * osc(3, j + m))
- nomi = nomi * (osc(1, j + m)**2 - wn**2 - cmplx(0.0, wn * osc(3, j + m)))
- deno = deno * (osc(1, j )**2 - wn**2 - cmplx(0.0, wn * osc(3, j )))
+ nomi = nomi * (osc(1, j + m)**2 - wn**2 - dcmplx(0.0d0, wn * osc(3, j + m)))
+ deno = deno * (osc(1, j )**2 - wn**2 - dcmplx(0.0d0, wn * osc(3, j )))
enddo
eps(l) = epsinf(l) * nomi / deno
elseif (name(l) == 'metal') then
j = j + 1
! suggestion for replacement
! eps(l) = -osc(1, j)**2 / cmplx(wn**2, wn * osc(3, j))
- eps(l) = -osc(1, j)**2 / ( wn**2 + cmplx(0.0d0, wn * osc(3, j)) )
+ eps(l) = -osc(1, j)**2 / ( wn**2 + dcmplx(0.0d0, wn * osc(3, j)) )
else
eps(l) = epsinf(l)
! The original version had this additional loop. It seems, it has been removed
@@ -65,7 +64,7 @@ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, thick, layers, nper, mo
do k = 1, nos(l)
j = j + 1
x = wn / osc(1, j)
- deno = x * cmplx(x, osc(3, j))
+ deno = x * dcmplx(x, osc(3, j))
if (osc(2, j) >= 0.0d0) then
deno = 1.0 - deno
endif
diff --git a/source/f90/surlos.f90 b/source/f90/surlos.f90
index 4d819a6..3cef23c 100644
--- a/source/f90/surlos.f90
+++ b/source/f90/surlos.f90
@@ -1,6 +1,6 @@
module surlos_mod
contains
-double precision function surlos(dk, eps, d, layers, nper)
+double precision function surlos(dk, eps, thick, layers, nper)
! ******************************************************************
! * *
@@ -12,7 +12,7 @@ double precision function surlos(dk, eps, d, layers, nper)
double precision, intent(in) :: dk
double complex, intent(in) :: eps(:)
- double precision, intent(in) :: d(:)
+ double precision, intent(in) :: thick(:)
integer, intent(in) :: layers, nper
integer :: lmax
@@ -27,7 +27,7 @@ double precision function surlos(dk, eps, d, layers, nper)
zero(z) = (real(z) == 0.0) .and. (aimag(z) == 0.0)
! write (*,*) 'surlos:'
-! write (*,*) 'd: ', size(d)
+! write (*,*) 'thick: ', size(thick)
! write (*,*) 'eps: ', size(eps)
lmax = size(eps)
@@ -37,7 +37,7 @@ double precision function surlos(dk, eps, d, layers, nper)
skip = .false.
do n = 1, layers
- arg(n) = dk * d(n)
+ arg(n) = dk * thick(n)
if (arg(n) > argmax .or. zero(eps(n))) then
csi = eps(n)
skip = .true.
@@ -57,8 +57,8 @@ double precision function surlos(dk, eps, d, layers, nper)
pn = 0.0
qn = 0.0
do n = lstart, layers
- pn = pn + d(n) * eps(n)
- qn = qn + d(n) / eps(n)
+ pn = pn + thick(n) * eps(n)
+ qn = qn + thick(n) / eps(n)
enddo
if (zero(qn)) then
n = lstart
--
GitLab