Skip to content
Snippets Groups Projects
Commit f32feb91 authored by Karl-Michael Schindler's avatar Karl-Michael Schindler
Browse files

import new seteps routine from branch newSeteps.

parent 01cf5a48
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
source/f90/CHANGELOG
+ 2
− 0
View file @ f32feb91
History
2025 07 09 use new seteps function in doeels, eels, eels-boson.
2025 04 17 Makefile also creates eels-boson
2014 07 22 separate file IO and create subroutine calls.
......
source/f90/Makefile
+ 25
− 5
View file @ f32feb91
......@@ -16,6 +16,9 @@ endif
# python3 interpreter
PYTHON3 = python3
# Platform, i.e. operating system
OS := $(shell echo `uname -s`)
# implicit rules for compiling fortran files
%.o: %.f90 ; $(FC) $(FFLAGS) -c $<
%_mod.mod: %.f90 ; $(FC) $(FFLAGS) -c $<
......@@ -31,7 +34,7 @@ bosonsubs = doboson.o respon.o sicot.o sintr.o rcffi.o o.o
boson: boson.f90 change_working_dir.o $(bosonsubs) $(bosonmods)
$(FC) $(FFLAGS) -o bosonf90 boson.f90 change_working_dir.o $(bosonsubs)
eelsmods = quanc8_mod.mod queels_mod.mod seteps_mod.mod param_mod.mod
eelsmods = quanc8_mod.mod queels_mod.mod seteps_mod.mod param_mod.mod oscillatorType_mod.mod
doeels.o: doeels.f90 $(eelsmods)
queels_mod.mod: quanc8_mod.mod fint1_mod.mod fint2_mod.mod fint3_mod.mod
......@@ -40,6 +43,7 @@ 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: oscillatorType_mod.mod
getoptsubs = getopt.o text.o system.o constants.o date_and_time.o kinds.o dummy_variables.o
......@@ -66,7 +70,7 @@ sufr_kinds.mod:
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
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 oscillatorType.o
eels: eels.f90 change_working_dir.o $(eelssubs) $(eelsmods)
$(FC) $(FFLAGS) -o eelsf90 eels.f90 change_working_dir.o $(eelssubs)
......@@ -78,11 +82,27 @@ eels-boson: eels-boson.f90 get_commandline_options.o $(getoptsubs) $(getoptmods)
pylibs: doeels-py doboson-py
doeels-py: doeels.f90 quanc8.f90 surlos.f90 fint1.f90 fint2.f90 fint3.f90 queels.f90 phint.f90 qrat.f90 usurlo.f90
$(PYTHON3) -m numpy.f2py -c seteps.f90 doeels.f90 quanc8.f90 surlos.f90 fint1.f90 fint2.f90 fint3.f90 queels.f90 phint.f90 qrat.f90 usurlo.f90 -m doeels
doeels-py: doeels.f90 quanc8.f90 surlos.f90 fint1.f90 fint2.f90 fint3.f90 queels.f90 phint.f90 qrat.f90 usurlo.f90 param.f90 oscillatorType.f90 fun.f90
ifeq ($(OS),Darwin)
$(PYTHON3) -m numpy.f2py -c seteps.f90 doeels.f90 quanc8.f90 surlos.f90 fint1.f90 fint2.f90 fint3.f90 queels.f90 phint.f90 qrat.f90 usurlo.f90 param.f90 oscillatorType.f90 fun.f90 -m doeels
install_name_tool -change @rpath/libgfortran.5.dylib /opt/local/lib/libgcc/libgfortran.5.dylib doeels.cpython-312-darwin.so
else
$(PYTHON3) -m numpy.f2py -c seteps.f90 doeels.f90 quanc8.f90 surlos.f90 fint1.f90 fint2.f90 fint3.f90 queels.f90 phint.f90 qrat.f90 usurlo.f90 param.f90 oscillatorType.f90 fun.f90 -m doeels
endif
doboson-py: doboson.f90 sicot.f90 sintr.f90 rcffi.f90 respon.f90 o.f90
$(PYTHON3) -m numpy.f2py -c o.90 sintr.f90 rcffi.f90 respon.f90 sicot.f90 doboson.f90 -m doboson
ifeq ($(OS),Darwin)
$(PYTHON3) -m numpy.f2py -c o.f90 sintr.f90 rcffi.f90 respon.f90 sicot.f90 doboson.f90 -m doboson
install_name_tool -change @rpath/libgfortran.5.dylib /opt/local/lib/libgcc/libgfortran.5.dylib doboson.cpython-312-darwin.so
else
$(PYTHON3) -m numpy.f2py -c o.f90 sintr.f90 rcffi.f90 respon.f90 sicot.f90 doboson.f90 -m doboson
endif
mytest:
echo $(OS)
ifeq ( $(OS), Darwin )
echo $(OS)
endif
clean:
rm -f *.o
......
source/f90/TODO
+ 0
− 2
View file @ f32feb91
- usage output when called without parameter or with options -h --help
- options for specifying input and output files.
- Makefile: make tests
\ No newline at end of file
source/f90/doeels.f90
+ 13
− 8
View file @ f32feb91
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)
layers, neps, nper, name, name_size, thick, epsinf, nos, oscType, osc, osc_size, &
contrl, wn_array, debug, f_array, wn_array_size)
! ******************************************************************
! * *
......@@ -24,10 +24,10 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
! * thick: layer thickness *
! * epsinf: epsilon at infinite frequency *
! * nos: number of oscillators *
! * oscType: oscillator type: Lorentz, Kurosawa, drude *
! * 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 *
......@@ -47,10 +47,10 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
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) :: thick(name_size), epsinf(name_size), osc(4, osc_size)
character (len = 10) :: contrl
integer, intent(in) :: comment_size, name_size, osc_size, wn_array_size, oscType(:)
integer, intent(in) :: 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
......@@ -113,10 +113,15 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
nout = 1 + nw / 20
ener = 8065 * e0
psia = phia / 180 * pi
! *** kms
psii = theta / 180 * pi
cospsi = dcos(psii)
sinpsi = dsin(psii)
tanpsi = dtan(psii)
! *** f2023 replacement suggestion
! cospsi = dcosd(phia)
! sinpsi = dsind(phia)
! tanpsi = dtand(phia)
prefac = dsqrt(255500 / e0)/(137 * cospsi)
facru = psia / cospsi * dsqrt(0.2624664d0 * e0)
elleps = (1.0d0 - phia / phib) * (1.0d0 + phia / phib)
......@@ -192,7 +197,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, layers, mode)
if (.not. user) call seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
x = wn / (2 * ener * psia)
if (rational) then
......
source/f90/eels-boson.f90
+ 30
− 14
View file @ f32feb91
......@@ -10,6 +10,7 @@ program eels_boson
use get_commandline_options_mod
use doeels_mod
use oscillator_mod
implicit none
......@@ -19,16 +20,16 @@ program eels_boson
integer :: i, j, jos, k, l, layers, neps, nper, nw
logical :: user
character (len = 72) :: comment(2)
character (len = 10) :: contrl, mode
character (len = 10) :: contrl
double precision, allocatable :: xout(:), yout(:)
integer :: nout
double precision, allocatable :: thick(:), epsinf(:), osc(:, :)
integer, allocatable :: nos(:)
integer, allocatable :: nos(:), oscType(:), tmp_oscType(:)
character (len = 10), allocatable :: name(:)
double precision, allocatable :: wn_array(:), f(:)
integer :: old_size_1, old_size_2
integer :: old_size, old_size_1, old_size_2
double precision, allocatable :: tmp_osc(:, :)
double precision :: asym, emax, emin, gauss, t, width
......@@ -79,14 +80,14 @@ program eels_boson
! *** read target specifications
read(11, *) layers, nper, mode
read(11, *) layers, nper
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'
write(*,*) 'The substrate is an anisotropic uniaxial material'
endif
if (layers < 0 .or. nper < 1 .or. nper > layers) then
stop '*** invalid target specifications ***'
......@@ -106,7 +107,7 @@ program eels_boson
else
do j = 1, nos(l)
if (.not. allocated(osc)) then
allocate(osc(3, nos(l)))
allocate(osc(4, nos(l)))
else if (j == 1) then
! call extend3(osc, nos(j))
old_size_1 = size(osc, 1)
......@@ -116,20 +117,35 @@ program eels_boson
deallocate(osc)
call move_alloc(tmp_osc, osc)
endif
if (.not. allocated(oscType)) then
allocate(oscType(nos(l)))
else if (j == 1) then
old_size = size(oscType)
allocate(tmp_oscType(old_size + nos(l)))
tmp_oscType(1:old_size) = oscType
deallocate(oscType)
call move_alloc(tmp_oscType, oscType)
endif
jos = jos + 1
read(11, *) (osc(k, jos), k = 1, 3)
oscType(jos) = Lorentz
if (osc(2, jos) == -1.0) then ! Drude term
osc(2, jos) = osc(3, jos)
osc(3, jos) = 0.0
oscType(jos) = Drude
endif
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)
write(6, 104) l, name(l), thick(l), epsinf(l), (osc(i, jos), i = 1, 4)
else
write(6, 105) epsinf(l), (osc(i, jos), i = 1, 3)
write(6, 105) epsinf(l), (osc(i, jos), i = 1, 4)
endif
else
write(6, 108) (osc(i, jos), i = 1, 3)
write(6, 108) (osc(i, jos), i = 1, 4)
endif
enddo
endif
......@@ -148,8 +164,8 @@ program eels_boson
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))
layers, neps, nper, name, size(name), thick, epsinf, nos, oscType, osc, size(osc, 2),&
contrl, wn_array, debug, f, size(wn_array))
open(unit = 13, file = bosin_name)
! target temperature (Kelvin)
......@@ -202,9 +218,9 @@ program eels_boson
'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)
104 format(1x, i3, 2x, a10, g15.3, f10.4, f12.4, f10.4, f9.4, f9.4)
105 format(31x, f10.4, f12.4, f10.4, f9.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)
108 format(41x, f12.4, f10.4, f9.4, f9.4)
end program eels_boson
source/f90/eels.f90
+ 30
− 14
View file @ f32feb91
......@@ -9,6 +9,7 @@ program eels
! ******************************************************************
use doeels_mod
use oscillator_mod
implicit none
......@@ -16,14 +17,14 @@ program eels
integer :: i, j, jos, k, l, layers, neps, nper, nw
logical :: user
character (len = 72) :: comment(2)
character (len = 10) :: contrl, mode
character (len = 10) :: contrl
double precision, allocatable :: thick(:), epsinf(:), osc(:, :)
integer, allocatable :: nos(:)
integer, allocatable :: nos(:), oscType(:), tmp_oscType(:)
character (len = 10), allocatable :: name(:)
double precision, allocatable :: wn_array(:), f(:)
integer :: old_size_1, old_size_2
integer :: old_size, old_size_1, old_size_2
double precision, allocatable :: tmp_osc(:, :)
logical :: debug
integer :: ioStatus
......@@ -65,14 +66,14 @@ program eels
! *** read target specifications
read(11, *) layers, nper, mode
read(11, *) layers, nper
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'
write(*,*) 'The substrate is an anisotropic uniaxial material'
endif
if (layers < 0 .or. nper < 1 .or. nper > layers) then
stop '*** invalid target specifications ***'
......@@ -92,7 +93,7 @@ program eels
else
do j = 1, nos(l)
if (.not. allocated(osc)) then
allocate(osc(3, nos(l)))
allocate(osc(4, nos(l)))
else if (j == 1) then
! call extend3(osc, nos(j))
old_size_1 = size(osc, 1)
......@@ -102,20 +103,35 @@ program eels
deallocate(osc)
call move_alloc(tmp_osc, osc)
endif
if (.not. allocated(oscType)) then
allocate(oscType(nos(l)))
else if (j == 1) then
old_size = size(oscType)
allocate(tmp_oscType(old_size + nos(l)))
tmp_oscType(1:old_size) = oscType
deallocate(oscType)
call move_alloc(tmp_oscType, oscType)
endif
jos = jos + 1
read(11, *) (osc(k, jos), k = 1, 3)
oscType(jos) = Lorentz
if (osc(2, jos) == -1.0) then ! Drude term
osc(2, jos) = osc(3, jos)
osc(3, jos) = 0.0
oscType(jos) = Drude
endif
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)
write(6, 104) l, name(l), thick(l), epsinf(l), (osc(i, jos), i = 1, 4)
else
write(6, 105) epsinf(l), (osc(i, jos), i = 1, 3)
write(6, 105) epsinf(l), (osc(i, jos), i = 1, 4)
endif
else
write(6, 108) (osc(i, jos), i = 1, 3)
write(6, 108) (osc(i, jos), i = 1, 4)
endif
enddo
endif
......@@ -134,8 +150,8 @@ program eels
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))
layers, neps, nper, name, size(name), thick, epsinf, nos, oscType, osc, size(osc, 2),&
contrl, wn_array, debug, f, size(wn_array))
open (unit = 12, file = 'eelsou')
write (12, 207) e0, theta, phia, phib
......@@ -151,10 +167,10 @@ program eels
'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)
104 format(1x, i3, 2x, a10, g15.3, f10.4, f12.4, f10.4, f9.4, f9.4)
105 format(31x, f10.4, f12.4, f10.4, f9.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)
108 format(41x, f12.4, f10.4, f9.4, f9.4)
207 format('e0 =', f6.2, ' theta =', f5.1, ' phia =', f5.2, ' phib =', f5.2)
end program eels
source/f90/oscillatorType.f90 0 → 100644
+ 4
− 0
View file @ f32feb91
module oscillator_mod
implicit none
integer, parameter :: Lorentz = 1, Kurosawa = 2, Drude = 3
end module oscillator_mod
source/f90/seteps.f90
+ 34
− 65
View file @ f32feb91
module seteps_mod
contains
subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
subroutine seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
! ******************************************************************
! * *
......@@ -9,83 +9,52 @@ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
! * *
! * neps: number of epsilon values *
! * nos: number of oscillators *
! * oscType: oscillator type: Lorentz, Kurosawa, drude *
! * osc: oscillator parameters, wTO, Q, lambda *
! * epsinf: epsilon at infinite frequency *
! * wn: frequency *
! * name: layer names *
! * wn: frequency *
! * eps: epsilon *
! * layers: number of layers *
! * mode: 'kurosawa' for kurosawa model *
! * *
! ******************************************************************
implicit none
integer, intent(in) :: neps
integer, intent(in) :: nos(:)
use oscillator_mod
integer, intent(in) :: neps, nos(:), oscType(:)
double precision, intent(in) :: osc(:, :), epsinf(:), wn
character (len=10), intent(in) :: name(:)
double complex, intent(in out) :: eps(:)
integer, intent(in) :: layers
character (len=10), intent(in) :: mode
double precision :: x
double complex :: deno, nomi
integer :: j, k, l, m
! write (*,*) 'seteps:'
! write (*,*) 'nos: ', size(nos)
! write (*,*) 'osc: ', size(osc)
! write (*,*) 'epsinf: ', size(epsinf)
! write (*,*) 'name: ', size(name)
! write (*,*) 'eps: ', size(eps)
! write (*,*) 'thick: ', size(thick)
integer :: i, j, k
double complex :: one, nomi, deno, addLorentz, addDrude
double precision :: wn2, b
j = 0
do l = 1, neps
m = nos(l)/2
nomi = dcmplx(1.0d0, 0.0d0)
deno = dcmplx(1.0d0, 0.0d0)
if (mode == 'kurosawa') then
do k = 1, m
j = j + 1
! since osc(1,*) and wn are real, the following should be equivalent
! Check required.
! Furthermore the first term can be rewritten as
! (osc(1, j + m) - wn) * (osc(1, j + m) + wn)
! nomi = nomi * cmplx(osc(1, j + m)**2 - wn**2, - wn * osc(3, j + m))
nomi = nomi * (osc(1, j + m)**2 - wn**2 - dcmplx(0.0d0, wn * osc(3, j + m)))
deno = deno * (osc(1, j )**2 - wn**2 - dcmplx(0.0d0, wn * osc(3, j )))
enddo
eps(l) = epsinf(l) * nomi / deno
elseif (name(l) == 'metal') then
j = j + 1
! suggestion for replacement
! eps(l) = -osc(1, j)**2 / cmplx(wn**2, wn * osc(3, j))
eps(l) = -osc(1, j)**2 / ( wn**2 + dcmplx(0.0d0, wn * osc(3, j)) )
else
eps(l) = epsinf(l)
! The original version had this additional loop. It seems, it has been removed
! because all cases of nos(l) > 1 are treated in case 1 above
do k = 1, nos(l)
j = j + 1
if (osc(1, j) > 1d-3) then ! prevent division by zero
x = wn / osc(1, j)
else
x = wn * 1d3
endif
deno = x * dcmplx(x, osc(3, j))
if (osc(2, j) >= 0.0d0) then
deno = 1.0d0 - deno
endif
if (cdabs(deno) == 0.0d0) then ! replace 0 by machine epsilon
! if deno is always > 0 then this would do it:
! deno = cdmax(deno, epsilon(1.0d0) / 2)
deno = epsilon(1.0d0) / 2
endif
eps(l) = eps(l) + osc(2, j) / deno
enddo
endif
one = dcmplx(1.0, 0.0)
wn2 = wn**2
k = 1
do i = 1, neps
nomi = one
deno = one
addLorentz = dcmplx(0.0, 0.0)
addDrude = dcmplx(0.0, 0.0)
do j = 1, nos(i)
select case (oscType(i))
case (Lorentz)
addLorentz = addLorentz + osc(1,k)**2 * osc(2,k) / dcmplx(osc(1,k)**2 - wn2, -wn*osc(1,k)*osc(3,k))
case (Kurosawa)
b = wn / osc(1,k)
nomi = nomi * osc(3,k)**2 * (1.0 - b * dcmplx( b, osc(4,k)/osc(3,k)) )
deno = deno * (osc(1,k)**2 - wn * dcmplx(wn, osc(2,k) ) )
case (Drude)
addDrude = addDrude - osc(1,k)**2 / dcmplx(wn2, wn*osc(1,j)*osc(2,j))
case default
write (*,'(A)') '*** Error in oscillator type! ***'
end select
k = k + 1
eps(i) = epsinf(i) * (nomi / deno + addDrude) + addLorentz
enddo
enddo
if (neps == layers + 1) then
! the substrate is a anisotropic uniaxial material
eps(layers) = cdsqrt(eps(layers) * eps(layers + 1))
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment