From 0fcb42d716a9d8e9d2976a80c2de02d3ebd8133e Mon Sep 17 00:00:00 2001
From: Karl-Michael Schindler <karl-michael.schindler@web.de>
Date: Thu, 10 Oct 2024 16:45:10 +0200
Subject: [PATCH] =?UTF-8?q?read=20yaml=20file=20for=20eels=20part=20should?=
 =?UTF-8?q?=20work=20=F0=9F=98=89?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

boson part to be done.
---
 source/f90/Makefile               |  10 +-
 source/f90/eels-boson.f90         | 146 ++++++++++++++++++++++++------
 source/f90/getOscillator.f90      |  56 ++++++++++++
 source/f90/oscillatorType.f90     |  21 ++++-
 source/f90/seteps.f90             |   2 +-
 tests/inputFiles/eelsbosin001.yml |   8 +-
 6 files changed, 201 insertions(+), 42 deletions(-)
 create mode 100644 source/f90/getOscillator.f90

diff --git a/source/f90/Makefile b/source/f90/Makefile
index bc645ae..1711e75 100644
--- a/source/f90/Makefile
+++ b/source/f90/Makefile
@@ -34,7 +34,7 @@ build/sicot_mod.mod: build/sintr_mod.mod
 
 bosonsubs = build/doboson.o build/respon.o build/sicot.o build/sintr.o build/rcffi.o build/o.o
 
-eelsmods = build/quanc8_mod.mod build/queels_mod.mod build/seteps_mod.mod build/param_mod.mod build/oscillatorType_mod.mod 
+eelsmods = build/quanc8_mod.mod build/queels_mod.mod build/seteps_mod.mod build/param_mod.mod build/oscillatorType_mod.mod  build/getOscillator_mod.mod
 build/doeels.o: doeels.f90 $(eelsmods)
 
 build/queels_mod.mod: build/quanc8_mod.mod build/fint1_mod.mod build/fint2_mod.mod build/fint3_mod.mod
@@ -75,19 +75,19 @@ build/sufr_dummy.mod:
 build/get_commandline_options.o: get_commandline_options.f90 $(getoptmods) $(getoptsubs)
 	$(FC) $(FFLAGS) -c -o $@ get_commandline_options.f90
 
-eelssubs = build/doeels.o build/usurlo.o build/quanc8.o build/fun.o build/queels.o build/fint1.o build/fint2.o build/fint3.o build/surlos.o build/seteps.o build/phint.o build/qrat.o build/param.o build/oscillatorType.o
+eelssubs = build/doeels.o build/usurlo.o build/quanc8.o build/fun.o build/queels.o build/fint1.o build/fint2.o build/fint3.o build/surlos.o build/seteps.o build/phint.o build/qrat.o build/param.o build/oscillatorType.o build/getOscillator.o
 
 eels-boson: eels-boson.f90 build/get_commandline_options.o $(getoptsubs) $(getoptmods) $(eelssubs) $(eelsmods) $(bosonsubs) $(bosonmods)
 	$(FC) $(FFLAGS) $(YAMLLIBS) -o build/$@ eels-boson.f90 build/get_commandline_options.o $(yamlsubs) $(getoptsubs) $(eelssubs) $(bosonsubs)
 
 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 param.f90 oscillatorType.f90 fun.f90
+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 getOscillator.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
+	$(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 getOscillator.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
+	$(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 getOscillator.f90 fun.f90 -m doeels
 endif
 
 doboson-py: doboson.f90 sicot.f90 sintr.f90 rcffi.f90 respon.f90 o.f90
diff --git a/source/f90/eels-boson.f90 b/source/f90/eels-boson.f90
index bc3265d..d6a952b 100644
--- a/source/f90/eels-boson.f90
+++ b/source/f90/eels-boson.f90
@@ -8,9 +8,12 @@ program eels_boson
 ! *                                                                *
 ! ******************************************************************
 
+  use, intrinsic :: iso_fortran_env, only:  output_unit
+
   use get_commandline_options_mod
   use doeels_mod
-  use oscillator_mod
+  use oscillatorType_mod
+  use getOscillator_mod
   use fortran_yaml_c
 
   implicit none
@@ -21,9 +24,12 @@ program eels_boson
   character(:), allocatable :: err
     
   class(type_node), pointer :: root
+  type(type_error), allocatable :: io_err
+  class(type_list), pointer :: layerList, termList
+  class(type_list_item), pointer :: layerListItem, termListItem
 
   double precision :: e0, theta, phia, phib, wmin, wmax, dw
-  integer :: i, j, jos, k, l, layers, neps, nper, nw
+  integer :: i, j, jos, l, layers, neps, nper, nw
   logical :: user
   character (len = 72) :: comment(2)
   character (len = 10) :: contrl
@@ -40,7 +46,6 @@ program eels_boson
 
   double precision :: asym, emax, emin, gauss, t, width
   logical :: debug
-  integer :: ioStatus
   
   character (len = :), allocatable :: eelsin_name, bosonout_name
 
@@ -56,23 +61,73 @@ program eels_boson
     stop 1
   endif
 
+  call file%dump(unit = output_unit, indent = 0)
+
   root => file%root
   
-! *** read spectrometer parameters
-  open(unit = 11, file = eelsin_name)
+  select type (root)
+    class is (type_dictionary)
+      
+! *** assign spectrometer parameters
 ! impact energy (ev)
-  read(11, *) e0
+  e0 = root%get_real('E0', error = io_err)
+  if (allocated(io_err)) then
+    print *, io_err%message
+    stop 1
+  endif
 ! incidence angle (%)
-  read(11, *) theta
+  theta = root%get_real('Theta', error = io_err)
+  if (allocated(io_err)) then
+    print *, io_err%message
+    stop 1
+  endif
 ! angular apertures of the elliptic detector (%)
-  read(11, *) phia
-  read(11, *) phib
+  phia = root%get_real('PhiA', error = io_err)
+  if (allocated(io_err)) then
+    print *, io_err%message
+    stop 1
+  endif
+  phib = root%get_real('PhiB', error = io_err)
+  if (allocated(io_err)) then
+    print *, io_err%message
+    stop 1
+  endif
 ! energy-loss interval and step size (cm**-1)
-  read(11, *) wmin
-  read(11, *) wmax
-  read(11, *) dw
+  wmin = root%get_real('wMin', error = io_err)
+  if (allocated(io_err)) then
+    print *, io_err%message
+    stop 1
+  endif
+  wmax = root%get_real('wMax', error = io_err)
+  if (allocated(io_err)) then
+    print *, io_err%message
+    stop 1
+  endif
+  dw = root%get_real('dw', error = io_err)
+  if (allocated(io_err)) then
+    print *, io_err%message
+    stop 1
+  endif
 ! comment lines
-  read(11, '(a72)') (comment(k), k = 1, 2)
+  comment(1) = root%get_string("Comment1", error = io_err)
+  if (allocated(io_err)) then
+    print *, io_err%message
+    stop 1
+  endif
+  comment(2) = root%get_string("Comment2", error = io_err)
+  if (allocated(io_err)) then
+    print *, io_err%message
+    stop 1
+  endif
+
+    layerList => root%get_list('Layers', required = .true., error = io_err)
+    if (allocated(io_err)) then
+      print *, io_err%message
+      stop 1
+    endif
+
+    print'(1x,a,i3)',"Number of layers: ", layerList%size()
+  end select
 
   write(*,*) 'program eels (September 2020), version: ', version
   write(*,'(a, f6.2, a, f5.1, a, f5.2, a, f5.2, a)') ' e0 = ', e0,       &
@@ -91,8 +146,9 @@ program eels_boson
 
 ! *** read target specifications
 
-  read(11, *) layers, nper
+  layers = layerList%size()
   user = layers == 0
+  nper = 1
   if (.not. user) then
     neps = layers
     if (nper == -1) then
@@ -106,21 +162,53 @@ program eels_boson
     write(6, 101) layers, nper
     jos = 0
     allocate (name(neps), thick(neps), epsinf(neps), nos(neps))
+    layerListItem => layerList%first
+    call layerListItem%node%dump(unit = output_unit, indent = 0)
+
     do l = 1, neps
       if (l <= layers) then
-        read(11, 102) name(l), thick(l)
+        write (*,*)
+        write (*,*) '*** Layer:', l, 'from', layers, 'layers. ***'
+        write (*,*)
+        select type (layer => layerListItem%node)
+        class is (type_dictionary)
+          name(l) = layer%get_string('Name', error = io_err)
+          if (allocated(io_err)) then
+            print *, io_err%message
+            stop 1
+          endif
+          thick(l) = layer%get_real('thickness', error = io_err)
+          if (allocated(io_err)) then
+            print *, io_err%message
+            stop 1
+          endif
+          epsinf(l) = layer%get_real('epsInfinity', error = io_err)
+          if (allocated(io_err)) then
+            print *, io_err%message
+            stop 1
+          endif
+          termList => layer%get_list('Terms', required = .true., error = io_err)
+          if (allocated(io_err)) then
+            print *, io_err%message
+            stop 1
+          endif
+          nos(l) = termList%size()
+        end select
+        write (*,*) 'Name: ', name(l)
+        write (*,*) 'Thickness: ', thick(l)
+        write (*,*) 'epsInfinity: ', epsinf(l)
+        write (*,*) 'Number of Oscillators: ', nos(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
+        termListItem => termList%first
         do j = 1, nos(l)
           if (.not. allocated(osc)) then
             allocate(osc(4, 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)))
@@ -138,13 +226,14 @@ program eels_boson
             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
+          call getOscillator(termListItem,  oscType, osc, jos)
+!          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)
@@ -158,17 +247,15 @@ program eels_boson
           else
             write(6, 108) (osc(i, jos), i = 1, 4)
           endif
+          termListItem => termListItem%next
         enddo
       endif
+      layerListItem => layerListItem%next
     enddo
     write(*,*)
-    read(11, 102, IOSTAT = ioStatus) contrl
-    if (ioStatus /= 0) then
-      contrl = ''
-    endif
   endif
 
-  close(unit = 11)
+!  close(unit = 11)
 
   nw = 1 + int((wmax - wmin) / dw)
   allocate (wn_array(nw), f(nw))
@@ -227,7 +314,6 @@ program eels_boson
 
 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, f9.4)
 105 format(31x, f10.4, f12.4, f10.4, f9.4, f9.4)
diff --git a/source/f90/getOscillator.f90 b/source/f90/getOscillator.f90
new file mode 100644
index 0000000..73f403c
--- /dev/null
+++ b/source/f90/getOscillator.f90
@@ -0,0 +1,56 @@
+module getOscillator_mod
+  contains
+  subroutine getOscillator(termListItem,  oscType, osc, jos)
+
+    use, intrinsic :: iso_fortran_env, only:  output_unit
+    use oscillatorType_mod
+    use fortran_yaml_c
+
+    class(type_list_item), pointer, intent(in) :: termListItem
+    integer, allocatable, intent(in out) :: oscType(:)
+    double precision, allocatable, intent(in out) :: osc(:, :)
+    integer, intent(in) :: jos
+
+    type(type_error), allocatable :: io_err
+    character(:), allocatable :: oscTypeString
+
+    write (*,*) 'Oscillator term', jos, ':'
+    call termListItem%node%dump(unit = output_unit, indent = 0)
+
+    select type (oscillatorTerm => termListItem%node)
+    class is (type_dictionary)
+      oscTypeString = oscillatorTerm%get_string('Type', error = io_err)
+      if (allocated(io_err)) then
+        print *, io_err%message
+        stop 1
+      endif
+      oscType(jos) = convert(oscTypeString)
+      select case (oscType(jos))
+      case (Lorentz)
+        osc(1,jos) = oscillatorTerm%get_real('wTO', error = io_err)
+        osc(2,jos) = oscillatorTerm%get_real('Q', error = io_err)
+        osc(3,jos) = oscillatorTerm%get_real('Gamma', error = io_err)
+        osc(4,jos) = 0.0
+      case (Kurosawa)
+        osc(1,jos) = oscillatorTerm%get_real('wTO', error = io_err)
+        osc(2,jos) = oscillatorTerm%get_real('GammaTO', error = io_err)
+        osc(3,jos) = oscillatorTerm%get_real('wLO', error = io_err)
+        osc(4,jos) = oscillatorTerm%get_real('GammaTO', error = io_err)
+      case (Drude)
+        osc(1,jos) = oscillatorTerm%get_real('wTO', error = io_err)
+        osc(2,jos) = oscillatorTerm%get_real('lambda', error = io_err)
+        osc(3,jos) = 0.0
+        osc(4,jos) = 0.0
+      case (extendedDrude)
+        osc(1,jos) = oscillatorTerm%get_real('wTO', error = io_err)
+        osc(2,jos) = oscillatorTerm%get_real('lambda', error = io_err)
+        osc(3,jos) = 0.0
+        osc(4,jos) = 0.0
+      end select
+    end select
+
+    write (*,*) 'Oscillator: ', oscTypeString, osc(1,jos), osc(2,jos), osc(3,jos), osc(4,jos)
+    return
+
+  end subroutine getOscillator
+end module getOscillator_mod
diff --git a/source/f90/oscillatorType.f90 b/source/f90/oscillatorType.f90
index aa017ed..6770548 100644
--- a/source/f90/oscillatorType.f90
+++ b/source/f90/oscillatorType.f90
@@ -1,4 +1,21 @@
-module oscillator_mod
+module oscillatorType_mod
   implicit none
   integer, parameter :: Lorentz = 1, Kurosawa = 2, Drude = 3, extendedDrude = 4
-end module oscillator_mod
+
+  contains  
+  integer function convert(oscillatorTypeString)
+    character(:), allocatable :: oscillatorTypeString
+
+    if      (oscillatorTypeString == 'Lorentz') then
+      convert = Lorentz
+    else if (oscillatorTypeString == 'Kurosawa') then
+      convert = Kurosawa
+    else if (oscillatorTypeString == 'Drude') then
+      convert = Drude
+    else if (oscillatorTypeString == 'extendedDrude') then
+      convert = extendedDrude
+    endif
+    
+    return
+  end function convert
+end module oscillatorType_mod
diff --git a/source/f90/seteps.f90 b/source/f90/seteps.f90
index 94874f7..8b6265a 100644
--- a/source/f90/seteps.f90
+++ b/source/f90/seteps.f90
@@ -18,7 +18,7 @@ subroutine seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
 ! *                                                                *
 ! ******************************************************************
 
-  use oscillator_mod
+  use oscillatorType_mod
 
   integer, intent(in) :: neps, nos(:), oscType(:)
   double precision, intent(in) :: osc(:, :), epsinf(:), wn
diff --git a/tests/inputFiles/eelsbosin001.yml b/tests/inputFiles/eelsbosin001.yml
index 483fb99..8fdfcd9 100644
--- a/tests/inputFiles/eelsbosin001.yml
+++ b/tests/inputFiles/eelsbosin001.yml
@@ -21,9 +21,9 @@ Layers:
          Gamma: 0.05}
       ,
         {Type: Lorentz,
-         wTO: 269.0,
-         Q:    16.000,
-         Gamma: 0.05}
+         wTO: 279.0,
+         Q:    17.000,
+         Gamma: 0.06}
       ]
   }
 - { Name: "Platinum",
@@ -31,7 +31,7 @@ Layers:
     epsInfinity: 8.90,
     Terms:
       [
-        {TypeOfTerm: Drude,
+        {Type: Drude,
          wTO: 160000.0,
          lambda: 0.012}
       ]
-- 
GitLab