esputr2001.f90

Go to the documentation of this file.

!    Copyright (c) 2016  Forschungszentrum Juelich GmbH
!                        Markus Brenneis, Vladislav Kotov 
!
!    This file is part of ESPUTR.
!
!    ESPUTR is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.
!
!    ESPUTR is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with ESPUTR.  If not, see <http://www.gnu.org/licenses/>.
!

module esputr2001

!***********************************************************************************
!   DEFINITION OF PARAMETRS AND GLOBAL VARIABLES
!***********************************************************************************
    use esputr
    implicit none
    public esputr2001_init, esputr2001_initn, esputr2001_initth, esputr2001_if_initialized, esputr2001_deallocate
    public esputr2001_check
    public esputr2001_getprojectiletargetidn, esputr2001_getprojectiletargetidsth
    public esputr2001_yn, esputr2001_yth
    public esputr2001_eth, esputr2001_availablecombinationsn, esputr2001_availablecombinationsth
    public esputr2001_getavailableenergyrange
    logical,public,save :: esputr2001_extrapolate_angular=.false.
    integer, parameter, public :: esputr2001_max_name_length = 16
    private
    integer, parameter, private :: nconstants_n = 5
    integer, parameter, private :: nconstants_th = 6
    integer,save, private :: ncombinations_n
    integer,save, private :: ncombinations_th
    real(ESPUTR_DP), dimension(:,:), allocatable, save, private :: fitvalues_n
    real(ESPUTR_DP), dimension(:,:), allocatable, save, private :: fitvalues_th
    character(ESPUTR2001_MAX_NAME_LENGTH), dimension(:,:), allocatable, save, private :: projectiles_targets_names_n
    character(ESPUTR2001_MAX_NAME_LENGTH), dimension(:,:), allocatable, save, private :: projectiles_targets_names_th
    integer, save, private :: nprojectilestargetsth = 0
    
    contains

!***********************************************************************************
!   INITIALIZATION/DEINITIALIZATION ROUTINES  
!***********************************************************************************
    
    subroutine esputr2001_init(fileNName, fileThName, err)

        character(*), intent(in) :: fileNName
        character(*), intent(in) :: fileThName
        integer, intent(out) :: err
 
        integer :: res

        err = 0
     ! read data for normal incidence
        call esputr2001_initn(filenname, err)
     ! read data for angular dependence
        if(err /= 0) return
        call esputr2001_initth(filethname, err)
        if(err /= 0) return
        call esputr2001_check(res, err)
        if( res /=0 ) then
          write(esputr_unit, *) "ERROR in esputr2001_init: created module is invalid"
          err=400
        end if
    end subroutine esputr2001_init
    
    subroutine esputr2001_initn(fileNName, err)

        character(*), intent(in) :: fileNName
        integer, intent(out) :: err
        err = 0
        if(allocated(fitvalues_n)) then
            write(esputr_unit, *) "ERROR in esputr2001_initN: esputr2001_initN has already been called; &
                                    &call esputr2001_deallocate first."
            err = 100
            return
        end if
        call init2001_(filenname, fitvalues_n, nconstants_n, & 
                        projectiles_targets_names_n, ncombinations_n, err)
        if(err == 0) then
           write(esputr_unit, *) " Number of projectile-target combinations = ", size(fitvalues_n, 2)
        end if
    end subroutine esputr2001_initn
    
    subroutine esputr2001_initth(fileTHName, err)

        character(*), intent(in) :: fileTHName
        integer, intent(out) :: err
        err = 0
        if(allocated(fitvalues_th)) then
            write(esputr_unit, *) "ERROR in esputr2001_initTH: esputr2001_initTH has already been called; &
                                    &call esputr2001_deallocate first."
            err = 100
            return
        end if
        call init2001_(filethname, fitvalues_th, nconstants_th, &
                        projectiles_targets_names_th, ncombinations_th, err)
        if(err == 0) then
           nprojectilestargetsth = numberofprojtargcombinations()
           write(esputr_unit, *) " Number of projectile-target combinations = ", nprojectilestargetsth
           write(esputr_unit, *) "  Number of projectile-target-energy combinations = ", size(projectiles_targets_names_th, 1)
        end if
        
        contains
        
        integer function numberofprojtargcombinations()
            character(ESPUTR_MAX_LINE_LENGTH) :: currentCombination
            character(ESPUTR_MAX_LINE_LENGTH) :: lastCombination
            integer :: i
            lastcombination = ""
            numberofprojtargcombinations = 0
            do i = 1, size(projectiles_targets_names_th, 1)
                currentcombination = projectiles_targets_names_th(i, 1) // projectiles_targets_names_th(i, 2)
                if(currentcombination /= lastcombination) then
                    numberofprojtargcombinations = numberofprojtargcombinations + 1
                    lastcombination = currentcombination
                end if
            end do
        end function numberofprojtargcombinations
    end subroutine esputr2001_initth
    
    subroutine init2001_(fileName, fitvalues, nConstants, projTargNames, nCombinations, err)

        character(*), intent(in) :: fileName
        real(ESPUTR_DP), dimension(:,:), allocatable, intent(inout) :: fitvalues
        integer, intent(in) :: nConstants
        character(ESPUTR2001_MAX_NAME_LENGTH), dimension(:,:), allocatable, intent(inout) :: projTargNames
        integer,intent(out) :: nCombinations
        integer, intent(out) :: err
    
        integer, parameter :: UNIT = 3
        integer :: nLines
        integer :: i, j, err_
        character(ESPUTR_MAX_LINE_LENGTH) :: line
        intrinsic trim, size

        i = 0
        j = 1
        
        ! -1 for the header line
        ncombinations = numberoflines(filename, .true., err=err)
        if (err /= 0) then
            write(esputr_unit, *) "ERROR in esputr2001_init: error while determining number of lines in "//trim(filename)
            err=300
            return
        end if
        if (ncombinations<1) then
            write(esputr_unit, *) "ERROR in esputr2001_init: "//trim(filename)//" contains no data"
            err=300
            return
        end if
        nlines = numberoflines(filename, .false., err)
        if (err /= 0) then
            write(esputr_unit, *) "ERROR in esputr2001_init: error while determining number of lines in "//trim(filename)
            err=300
            return
        end if
        
        allocate(fitvalues(nconstants,ncombinations), stat=err_); if(err_/=0) goto 3200
        allocate(projtargnames(ncombinations,2), stat=err_); if(err_/=0) goto 3200
        
        open(unit, iostat=err, file=filename, status='old', action='read')
        if (err /= 0 ) goto 3300        

        do i = 1, nlines
            read(unit, '(A)', iostat=err) line
            if (err /= 0 ) goto 3350    
            if(line(1:1) /= esputr_comment_start) then
                read(line, *, iostat=err) projtargnames(j,:), fitvalues(1:nconstants,j)
                if (err /= 0 ) goto 3350    
                j = j+1
            end if
        end do
        
        close(unit)
        
        write(esputr_unit, *) "esputr2001_init: initialization from file "//trim(filename)//" completed"
        return
        
3200    write(esputr_unit, *) "ERROR in esputr2001_init: could not allocate memory for data from file "//trim(filename)
        write(esputr_unit, *) " Number of entries = ",ncombinations,&
                              " number of parameter for each entry = ",nconstants
        err = 200
        return
3300    write(esputr_unit, *) "ERROR in esputr2001_init: cannot open "//trim(filename)
        err = 300
        return
3350    write(esputr_unit, *) "ERROR in esputr2001_init while reading line", i, "of "//trim(filename)
        call esputr2001_deallocate(err)
        err = 300
        close(unit)
        return
        
    end subroutine init2001_
    
    subroutine esputr2001_check(res, err)

        integer, intent(out) :: res
        integer, intent(out) :: err

        integer :: i,n,m,nn,mm
        character(ESPUTR_MAX_LINE_LENGTH) :: currentCombination, lastCombination
        real(ESPUTR_DP) :: lastEnergy,lambda,q,Eth,Etf,E0,Esp
        intrinsic size, trim
        
        err = 0
        res = 0

        if ( (.not.allocated(fitvalues_n)).or.(.not.allocated(fitvalues_th)) .or. & 
             (.not.allocated(projectiles_targets_names_n)) .or. &
             (.not.allocated(projectiles_targets_names_th)) ) then
           write(esputr_unit, *) "ERROR detected in esputr2001_check: some or all arrays are not allocated"
           res = 1
           return
        end if

        n=size(fitvalues_n,2)
        m=size(fitvalues_th,2)
        nn=size(projectiles_targets_names_n,1)
        mm=size(projectiles_targets_names_th,1)
        if ( ncombinations_n /= n .or. ncombinations_th /= m .or. & 
             ncombinations_n /= nn .or. ncombinations_th /= mm ) then
          write(esputr_unit, *) "ERROR detected in esputr2001_check: inconsistent array size"
          write(esputr_unit, *) " NCOMBINATIONS_N, size(FITVALUES_N), size(PROJECTILES_TARGETS_NAMES_N) ",&
                                  ncombinations_n, n, nn
          write(esputr_unit, *) " NCOMBINATIONS_N, size(FITVALUES_TH), size(PROJECTILES_TARGETS_NAMES_TH) ",&
                                  ncombinations_th, m, mm
          res=2
          return
        end if  

        !cycle over the normal incidence data
        do i = 1, ncombinations_n
         lambda=fitvalues_n(1,i) 
         q=fitvalues_n(2,i) 
         eth=fitvalues_n(4,i)     
         etf=fitvalues_n(5,i)     
         if(lambda <= 0.) then
           write(esputr_unit, *) "ERROR detected in esputr2001_check: wrong lambda = ",lambda
           write(esputr_unit, *)  " pos, proj, targ ",i, &
                                    trim(projectiles_targets_names_n(i,1)),' ',trim(projectiles_targets_names_n(i,2))
           res = 3
           return
         else if(q <= 0. ) then
                 write(esputr_unit, *) "ERROR detected in esputr2001_check: wrong q = ", q
                 write(esputr_unit, *)  " pos, proj, targ ",i, &
                                          trim(projectiles_targets_names_n(i,1)),' ',trim(projectiles_targets_names_n(i,2))
                 res = 4
                 return
         else if(eth <= 0.) then
                 write(esputr_unit, *) "ERROR detected in esputr2001_check: wrong E_th = ", eth
                 write(esputr_unit, *)  " pos, proj, targ ",i, &
                                          trim(projectiles_targets_names_n(i,1)),' ',trim(projectiles_targets_names_n(i,2))
                 res = 5
                 return
         else if(etf < 0.) then
                 write(esputr_unit, *) "ERROR detected in esputr2001_check: wrong eps (E_tf) = ", etf
                 write(esputr_unit, *)  " pos, proj, targ ",i, &
                                          trim(projectiles_targets_names_n(i,1)),' ',trim(projectiles_targets_names_n(i,2))
                 res = 6
                 return
         end if
        end do 

        !cycle over the angular dependence data
        lastcombination = ""
        lastenergy = -1.
        do i = 1, ncombinations_th
         currentcombination = projectiles_targets_names_th(i,1) // projectiles_targets_names_th(i,2)
         e0=fitvalues_th(1,i)
         esp=fitvalues_th(6,i)
         if(e0 < 0.) then
            write(esputr_unit, *) "ERROR detected in esputr2001_check: wrong E0 = ", e0
            write(esputr_unit, *)  " pos, proj, targ ",i, &
                                    trim(projectiles_targets_names_th(i,1)),' ',trim(projectiles_targets_names_th(i,2))
            res = 7
            return
          else if(esp < 0.) then
                  write(esputr_unit, *) "ERROR detected in esputr2001_check: wrong E_sp = ", esp
                  write(esputr_unit, *)  " pos, proj, targ ",i, &
                                           trim(projectiles_targets_names_th(i,1)),' ',trim(projectiles_targets_names_th(i,2)) 
                 res = 8
                 return
         else if( currentcombination == lastcombination .and. e0 <= lastenergy ) then
                  write(esputr_unit, *) "ERROR detected in esputr2001_check: ", &
                                        "incident energy for angular dependence is not monotonically increasing "
                  write(esputr_unit, *)  " pos, proj, targ ",i, &
                                           trim(projectiles_targets_names_th(i,1)),' ',trim(projectiles_targets_names_th(i,2)) 
                  res = 9
                  return
         end if
         lastenergy = fitvalues_th(1,i)
         lastcombination = currentcombination
        end do

    end subroutine esputr2001_check
    
    logical function esputr2001_if_initialized()
        esputr2001_if_initialized = allocated(fitvalues_n) .and. allocated(fitvalues_th)
    end function esputr2001_if_initialized
    
    subroutine esputr2001_deallocate(err)

        integer, intent(out) :: err
        integer, dimension(4) :: errs
        err = 0
        errs = 0
        if(allocated(fitvalues_n)) deallocate(fitvalues_n, stat=errs(1))
        if(allocated(fitvalues_th)) deallocate(fitvalues_th, stat=errs(2))
        if(allocated(projectiles_targets_names_n)) deallocate(projectiles_targets_names_n, stat=errs(3))
        if(allocated(projectiles_targets_names_th)) deallocate(projectiles_targets_names_th, stat=errs(4))
        if(.not. all(errs==0)) then
            write(esputr_unit, *) "WARNING from esputr2001_deallocate: problem while deallocating memory"
            err = -200
        end if
        return
    end subroutine esputr2001_deallocate

!***********************************************************************************
!   GET INDEX
!***********************************************************************************
    
    integer function esputr2001_getprojectiletargetidn(proj, targ, err)

        character(*), intent(in) :: proj
        character(*), intent(in) :: targ
        integer, intent(out) :: err
        intrinsic trim
        err = 0
        
        if(.not. allocated(projectiles_targets_names_n)) then
            write(esputr_unit, *) "ERROR in esputr2001_getProjectileTargetIdN: arrays are not initialized" 
            write(esputr_unit, *) "  Call esputr2001_init first!"
            esputr2001_getprojectiletargetidn = 0
            err =   50
            return
        end if
        
        esputr2001_getprojectiletargetidn = getindexof2(proj, targ, projectiles_targets_names_n)
        
        if(esputr2001_getprojectiletargetidn == 0) then
            write(esputr_unit, *) "ERROR in esputr2001_getProjectileTargetIdN: projectile-target combination "&
                                 //trim(proj)//"-"//trim(targ)// " is not found (normal incidence)"
            err = 151
        end if
    end function esputr2001_getprojectiletargetidn
    
    subroutine esputr2001_getprojectiletargetidsth(proj, targ, thLower, thUpper, err)

        character(*), intent(in) :: proj
        character(*), intent(in) :: targ
        integer, intent(out) :: thLower 
        integer, intent(out) :: thUpper 
        integer, intent(out) :: err
        intrinsic trim
        err = 0
        
        if(.not. allocated(projectiles_targets_names_th)) then
            write(esputr_unit, *) "ERROR in esputr2001_getProjectileTargetIdsTh: arrays are not initialized"
            write(esputr_unit, *) "  Call esputr2001_init first!"
            thlower = 0
            thupper = 0
            err =  50
            return
        end if
        
        thlower = getindexof2(proj, targ, projectiles_targets_names_th, thupper)
        
        if(thlower == 0 .or. thupper == 0) then
            write(esputr_unit, *) "ERROR in esputr2001_getProjectileTargetIdTh: projectile-target combination "&
                                 //trim(proj)//"-"//trim(targ)// " is not found (angular dependence)"
            err = 152
        end if
    end subroutine esputr2001_getprojectiletargetidsth

!***********************************************************************************
!   CALCULATION OF SPUTTERING YIELDS
!***********************************************************************************
       
    real(ESPUTR_DP) function esputr2001_yn(E0, projectileTarget_id, err)

        real(ESPUTR_DP),    intent(in)  :: E0
        integer, intent(in)  :: projectileTarget_id
        integer, intent(out) :: err
        
        real(ESPUTR_DP) :: lambda, q, mu, E_th, E_tf, &
                           powerMuTerm, denumerator, numerator, eps
        
        err = 0
        
        if(.not. allocated(fitvalues_n)) then
            write(esputr_unit, *) "ERROR in esputr2001_yn: arrays are not initialized"
            write(esputr_unit, *) "  Call esputr2001_init first!"
            esputr2001_yn = -1
            err =  50
            return
        end if
        if(projectiletarget_id < 1 .or. projectiletarget_id > ncombinations_n ) then
            write(esputr_unit, *) "ERROR in esputr2001_yn: projectileTarget_id is out of bounds", projectiletarget_id
            esputr2001_yn = -1
            err = 100
            return
        end if

        e_th   = fitvalues_n(4,projectiletarget_id)
               
        if (e0>e_th) then
         lambda = fitvalues_n(1,projectiletarget_id)
         q      = fitvalues_n(2,projectiletarget_id)
         mu     = fitvalues_n(3,projectiletarget_id)
         e_tf   = fitvalues_n(5,projectiletarget_id)
         eps = e0/e_tf
         powermuterm = (e0/e_th - 1)**mu        
         numerator   = powermuterm * log(1 + 1.2288*eps)
         denumerator = lambda + powermuterm * (eps + .1728*sqrt(eps) + .008*eps**.1504)       
         esputr2001_yn = 0.5 * q * numerator/denumerator
        else
         esputr2001_yn = 0
        end if

    end function esputr2001_yn
    
    real(ESPUTR_DP) function esputr2001_yth(E0, theta, projTargStartId, projTargEndId, err)
        !Incident energy, eV 
        real(ESPUTR_DP), intent(in) :: E0
        real(ESPUTR_DP), intent(in) :: theta
        integer, intent(in) :: projTargStartId
        integer, intent(in) :: projTargEndId
        integer, intent(out) :: err
        real(ESPUTR_DP) :: E_sp, f, c, b, cosTerm, theta0, S, dE, E00
        integer :: i1, i2
        
        err = 0
        
        if(.not. allocated(fitvalues_th)) then
            write(esputr_unit, *) "ERROR in esputr2001_yth: arrays are not initialized"
            write(esputr_unit, *) "Call esputr2001_init first"
            esputr2001_yth = -1
            err =  50
            return
        end if

        if(projtargstartid < 1) then
            write(esputr_unit, *) "ERROR in esputr2001_yth: projTargStartId is out of bounds ",projtargstartid
            esputr2001_yth = -1
            err = 153
            return
        end if

        if(projtargendid > ncombinations_th) then
            write(esputr_unit, *) "ERROR in esputr2001_yth: projTargEndId is out of bounds ",projtargendid
            esputr2001_yth = -1
            err = 154
            return
        end if

        if(projtargstartid > projtargendid ) then
            write(esputr_unit, *) "ERROR in esputr2001_yth: & 
                                  &projTargStartId and projTargEndId are inconsistent ",&
                                   projtargstartid,projtargendid
            esputr2001_yth = -1
            err = 155
            return
        end if

        if(e0 < 0) then
          write(esputr_unit, *) "ERROR in esputr2001_yth_: incident energy cannot be < 0"
          err = 156
          esputr2001_yth = -1
          return
        end if

        if(theta < 0 .or. theta > esputr_pi2) then
          write(esputr_unit, *) "ERROR in esputr2001_yth: incorrect incident angle theta (radian) ", theta
          write(esputr_unit, *) " Incident angle cannot be < 0 or > pi/2"
          err = 157
          esputr2001_yth = -1
          return
        end if
        
        !Find element of the array 
        i1=findindex(projtargstartid, projtargendid, e0, fitvalues_th(1,:), err) 
        e00=e0
        if(esputr2001_extrapolate_angular) then
         if(err == 21) then
          err=0 
          e00=fitvalues_th(1,i1)
         end if
        else
         if(err == 21) goto 4102
        end if
        if(err /= 0) goto 4100

        !VK linear interpolation        
        !Here it is assumed that E0 in the table is always monotonically increasing 
        i2=i1+1
        s=fitvalues_th(1,i2)-fitvalues_th(1,i1)
        s=1./s
        de=e00-fitvalues_th(1,i1)    

        f=fitvalues_th(2,i1)+(fitvalues_th(2,i2)-fitvalues_th(2,i1))*s*de
        b=fitvalues_th(3,i1)+(fitvalues_th(3,i2)-fitvalues_th(3,i1))*s*de
        c=fitvalues_th(4,i1)+(fitvalues_th(4,i2)-fitvalues_th(4,i1))*s*de
        e_sp=fitvalues_th(6,i1)+(fitvalues_th(6,i2)-fitvalues_th(6,i1))*s*de
                      
        theta0 = esputr_pi - acos(sqrt(e_sp/(e00+e_sp))) !VK always <= Pi/2 if both E_sp>=0 and E0>0
        
        costerm = cos((theta*esputr_pi2/theta0)**c)
        
        if(costerm > 0) then
            esputr2001_yth = costerm**(-f) * exp(b * (1-1/costerm))
        else
            esputr2001_yth = 0
        end if

        return
        
4100    write(esputr_unit, *) "ERROR in esputr2001_yth: failed to find fitting parameters for this energy"
        write(esputr_unit, *) " E0, E_min, E_max ",e0, &
                                fitvalues_th(1,projtargstartid),fitvalues_th(1,projtargendid)
        esputr2001_yth = -1
        err = 400
        return
4102    write(esputr_unit, *) "ERROR in esputr2001_yth: no data available for this energy"
        write(esputr_unit, *) " E0, E_min, E_max ",e0, &
                                fitvalues_th(1,projtargstartid),fitvalues_th(1,projtargendid)
        esputr2001_yth = -1
        err = 158
        return
    end function esputr2001_yth

!***********************************************************************************
!   AUXILIARY FUNCTIONS
!***********************************************************************************

    real(ESPUTR_DP) function esputr2001_eth(projectileTarget_id, err)

        integer, intent(in) :: projectileTarget_id
        integer, intent(out) :: err
        if(.not. allocated(fitvalues_n)) then
            esputr2001_eth = -1
            write(esputr_unit, *) "ERROR in esputr2001_Eth: arrays are not initialized" 
            write(esputr_unit, *) "  Call esputr2001_init first!"
            err = 50
            return
        end if
        esputr2001_eth = fitvalues_n(4, projectiletarget_id)
    end function esputr2001_eth
    
    subroutine esputr2001_availablecombinationsn(combinations, err)

        character(ESPUTR2001_MAX_NAME_LENGTH), dimension(:,:), allocatable, intent(inout) :: combinations
        integer, intent(out) :: err
        err = 0
        if(.not.allocated(projectiles_targets_names_n)) then
            err = 50
            write(esputr_unit, *) "ERROR in esputr2001_availableCombinationsN: arrays are not initialized" 
            write(esputr_unit, *) "  Call esputr2001_init first!"
            return
        end if
        allocate(combinations(size(projectiles_targets_names_n,1), size(projectiles_targets_names_n,2)), stat = err)
        if(err /= 0) then
            write(esputr_unit, *) "ERROR in esputr2001_availableCombinationsN: error while allocating memory", err
            err = 200
            return
        end if
        combinations = projectiles_targets_names_n
    end subroutine esputr2001_availablecombinationsn
    
    subroutine esputr2001_availablecombinationsth(combinations, err)

        character(ESPUTR2001_MAX_NAME_LENGTH*2), dimension(:,:), allocatable, intent(inout) :: combinations
        integer, intent(out) :: err
        character(ESPUTR_MAX_LINE_LENGTH) :: currentCombination, lastCombination
        integer :: i, numberOfCombination
        
        if(.not.allocated(projectiles_targets_names_th)) then
            err = 50
            write(esputr_unit, *) "ERROR in esputr2001_availableCombinationsTH: arrays are not initialized" 
            write(esputr_unit, *) "  Call esputr2001_init first!"
            return 
        end if
        err = 0
        
        allocate(combinations(nprojectilestargetsth, 2), stat = err)
        if(err /= 0) then
            write(esputr_unit, *) "ERROR in esputr2001_availableCombinationsTH: error while allocating memory", err
            err = 200
            return
        end if
        
        numberofcombination = 0
        lastcombination = ""
        do i = 1, size(projectiles_targets_names_th, 1)
            currentcombination = projectiles_targets_names_th(i, 1) // projectiles_targets_names_th(i, 2)
            if(currentcombination /= lastcombination) then
                numberofcombination = numberofcombination + 1
                combinations(numberofcombination,:) = (/ projectiles_targets_names_th(i, 1), projectiles_targets_names_th(i, 2) /)
                lastcombination = currentcombination
            end if
        end do
    end subroutine esputr2001_availablecombinationsth
    
    subroutine esputr2001_getavailableenergyrange(projTargStartId, projTargEndId, Emin, Emax, err)

        integer, intent(in) :: projTargStartId
        integer, intent(in) :: projTargEndId
        real(kind=ESPUTR_DP), intent(out) :: Emin
        real(kind=ESPUTR_DP), intent(out) :: Emax
        integer, intent(out) :: err
        err = 0
        if(.not.allocated(fitvalues_th)) then
            err = 50
            write(esputr_unit, *) "ERROR in esputr2001_getAvailableEnergyRange: arrays are not initialized" 
            write(esputr_unit, *) "  Call esputr2001_init first!"
            return
        end if
        if(size(fitvalues_th,2) < projtargendid .or. &
          0 >= projtargstartid .or. &
          projtargstartid > projtargendid) then
            write(esputr_unit, *) "ERROR in esputr2001_getAvailableEnergyRange: lower/upper bound is invalid",&
                projtargstartid, projtargstartid
            err = 100
            return
        end if
        emin = fitvalues_th(1, projtargstartid)
        emax = fitvalues_th(1, projtargendid)
    end subroutine esputr2001_getavailableenergyrange
    
    integer function getindexof2(subj1, subj2, list, lastIndex)

        character(*), intent(in) :: subj1, subj2
        character(*), dimension(:,:), intent(in) :: list
        integer, optional, intent(inout) :: lastIndex
        integer :: i, j, n
        intrinsic size,present 
        getindexof2 = 0
        if(present(lastindex)) lastindex = 0
        n=size(list, 1)
        foreachitem: do i = 1, n
            match: if(list(i,1) == subj1 .and. list(i,2) == subj2) then
                getindexof2 = i
                mustfindlast: if(present(lastindex)) then
                    lastindex = i
                    findlast: do j = i+1, n
                        if(list(j,1) == subj1 .and. list(j,2) == subj2) then
                            lastindex = j
                        else
                            exit findlast
                        end if
                    end do findlast
                end if mustfindlast
                exit foreachitem
            end if match
        end do foreachitem
    end function getindexof2
    
    integer function findindex(loBound, upBound, x, arr, err)

        integer, intent(in) :: loBound
        integer, intent(in) :: upBound
        real(ESPUTR_DP), intent(in) :: x
        real(ESPUTR_DP), dimension(:), intent(in) :: arr
        integer :: err

        integer :: i1, i2, i
        
        err = 0
        
        i1 = lobound
        i2 = upbound
        
        if(arr(lobound) > x) then
            err=21   
            findindex=lobound
            return
        end if
        if(arr(upbound) < x) then
            err=21
            findindex=upbound
            return 
        end if
        
        do
            i = (i1 + i2) / 2
            if(i1+1 == i2) then
                exit
            else if(i1 == i2) then
                if(i1 == upbound) then
                    i1 = i1 - 1
                end if
                exit
            else if(arr(i) < x) then
                i1 = i
            else
                i2 = i
            end if
        end do
        findindex=i1
    end function findindex
    
end module esputr2001