esputr1993.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 esputr1993

!***********************************************************************************
!   DEFINITION OF PARAMETRS AND GLOBAL VARIABLES
!***********************************************************************************
    use esputr
    implicit none
    public esputr1993_init, esputr1993_if_initialized, esputr1993_check, esputr1993_deallocate
    public esputr1993_yn, esputr1993_yth, esputr1993_eth
    public esputr1993_getprojectileid, esputr1993_gettargetid
    public esputr1993_availablecombinations    
    integer, parameter, public :: esputr1993_first_col_width = 11
    private
    integer, save, private :: nheaderlines
    integer, parameter, private :: col_width = 9
    integer, parameter, private :: lines_per_target = 6
    integer, parameter, private :: nconstants = 3
    integer,save, private :: ntargets
    integer,save, private :: nprojectiles
    real(ESPUTR_DP), dimension(:,:,:), allocatable, save, private :: fitvalues
    character(ESPUTR1993_FIRST_COL_WIDTH), dimension(:), allocatable, save, private :: targets_names
    character(ESPUTR1993_FIRST_COL_WIDTH), dimension(:), allocatable, save, private :: projectiles_names
    
    contains

!***********************************************************************************
!   INITIALIZATION/DEINITIALIZATION ROUTINES  
!***********************************************************************************
    subroutine esputr1993_init(constantsFile, err)

        character(*), intent(in) :: constantsFile
        integer, intent(out) :: err
        integer, parameter :: UNIT = 3
        integer :: i, j, Nlines, res
        character(ESPUTR_MAX_LINE_LENGTH) :: line
        !lineCols(NPROJECTILES + 1) : data on all projectiles for a target
        character(ESPUTR1993_FIRST_COL_WIDTH), dimension(:), allocatable :: lineCols
        intrinsic trim        

        if(allocated(fitvalues)) then
            write(esputr_unit, *) "ERROR in esputr1993_init: esputr1993_init has already been called; &
                                    &call esputr1993_deallocate first"
            err = 100
            return
        end if
        
        err = 0
               
        nheaderlines = linesuntil(constantsfile, '===', err)
        if(err /= 0) then
            write(esputr_unit, *) "ERROR in esputr1993_init: error while determining header lines of " // trim(constantsfile)
            err = 300
            return
        end if

        ! minus the two lines for the projectiles
        nlines=numberoflines(constantsfile, err=err)-nheaderlines-2
        if(err /= 0) then
            write(esputr_unit, *) "ERROR in esputr1993_init: error while determining number of lines in "&
                // trim(constantsfile), err
            err = 300
            return
        end if
        !VK to produce error if lines are missing for the last target
        ntargets = nlines/lines_per_target 
        if (ntargets*lines_per_target<nlines)  ntargets=ntargets+1 !VK
        if(ntargets<1) then 
            write(esputr_unit, *) "ERROR in esputr1993_init: file ",trim(constantsfile)," contains no data"
            err = 300
            return
        end if

        open(unit, iostat=err, file=trim(constantsfile), status='old', action='read')
        
        ifioerr: if(err /= 0) then
            write(esputr_unit, *) "ERROR in esputr1993_init: cannot open file " // trim(constantsfile)
            err = 300
            close(unit)
            return
        else
            skipheader: do i = 1, nheaderlines
                read(unit, '(A)', iostat=err) line
                if(err /= 0) then                  
                    write(esputr_unit, *) "ERROR in esputr1993_init while reading the comment (first lines) in " &
                                        // trim(constantsfile)
                    err = 300
                    close(unit)
                    return
                end if
            end do skipheader
            
            ! read line with names of projectiles and allocate memory for arrays
            read(unit, '(A)', iostat=err) line
            if(err /= 0) goto 1100
            ! ceiling because of missing spaces at end of line
            nprojectiles = ceiling((len_trim(line)-esputr1993_first_col_width+1.)/col_width)
            if(nprojectiles<1) then
                write(esputr_unit, *) "ERROR in esputr1993_init: "//trim(constantsfile)//" contains no data"
                err = 300
                close(unit)
                return
            end if
            allocate(targets_names(ntargets), stat=err)
            if(err /= 0) goto 1200
            allocate(projectiles_names(nprojectiles), stat=err)
            if(err /= 0) goto 1200
            allocate(fitvalues(nconstants, ntargets, nprojectiles), stat=err)
            if(err /= 0) goto 1200
            allocate(linecols(nprojectiles + 1), stat=err)
            if(err /= 0) goto 1200
            backspace(unit, iostat=err)
            if(err /= 0) goto 1300
            read(unit, *, iostat=err) linecols
            if(err /= 0) goto 1100
            projectiles_names(:) = linecols(2:)
            
            read(unit, *, iostat=err) linecols(:)
            if(err /= 0) goto 1100
!             line with AZ of projectiles is skiped
!             read(lineCols(2:), *, err=1300) SPUTER1993_PROJECTILES_MASSES
            
            foreachtarget: do i = 1, ntargets
                ! read name !and mass
                read(unit, *, iostat=err) linecols(1:2)
                if(err /= 0) goto 1300
                read(linecols(1:1), *, iostat=err) targets_names(i) !, SPUTER1993_TARGETS_MASSES(i)
                if(err /= 0) goto 1300

                ! skip Es
                read(unit, *, iostat=err) line
                if(err /= 0) goto 1300

                ! skip M2/M1
                read(unit, *, iostat=err) line
                if(err /= 0) goto 1300

                ! read E_tf
                read(unit, *, iostat=err) linecols
                if(err /= 0) goto 1300
                read(linecols(2:nprojectiles + 1), *, iostat=err) fitvalues(1,i,:)
                if(err /= 0) goto 1300

                ! read E_th
                read(unit, *, iostat=err) linecols
                if(err /= 0) goto 1300
                read(linecols(2:nprojectiles + 1), *, iostat=err) fitvalues(2,i,:)
                if(err /= 0) goto 1300

                ! read Q
                read(unit, *, iostat=err) linecols
                if(err /= 0) goto 1300
                read(linecols(2:nprojectiles + 1), *, iostat=err) fitvalues(3,i,:)
                if(err /= 0) goto 1300
                
                do j = 1, nprojectiles
                    if(err /= 0) then
                        write(esputr_unit, *) "ERROR in esputr1993_init: fitvalues are not sensible: ",&
                                                "nr ", i, ", value ", fitvalues(err,i,j)
                        close(unit)
                        return
                    end if
                end do
            end do foreachtarget
            
            close(unit)
            deallocate(linecols, stat=err)
            if(err /= 0) then
              write(esputr_unit, *) "WARNING from esputr1993_init: problem while deallocating memory"
              err = -200
            end if
           
            call esputr1993_check(res, err)
            if(res /=0 .or. err /=0) then
             err=400
             write(esputr_unit, *) "ERROR in esputr1993_init: the data are incorrect"
             return
            end if
          
            write(esputr_unit, *) "esputr1993_init: initialization from file "//trim(constantsfile)//" completed"
            write(esputr_unit, *) " Number of projectiles = ",nprojectiles, &
                                  ", Number of targets = ",ntargets
            return

1100        write(esputr_unit, *) "ERROR in esputr1993_init while reading the table header in " & 
                                // trim(constantsfile)
            err = 300
            close(unit)
            return             
1200        write(esputr_unit, *) "ERROR in esputr1993_init: could not allocate memory for ", ntargets,&
                " targets and ", nprojectiles, " projectiles with ", & 
                  nconstants, " parameters for each combination"
            err = 200
            close(unit)
            return
1300        write(esputr_unit, *) "ERROR in esputr1993_init while reading " // trim(constantsfile) &
                                 // " for target number ",i," (NAME = ",trim(targets_names(i)),")"
            err = 300
            close(unit)
            return
        end if ifioerr
        
    end subroutine esputr1993_init
    
    logical function esputr1993_if_initialized()
        esputr1993_if_initialized = allocated(fitvalues)
    end function esputr1993_if_initialized
    
    subroutine esputr1993_check(res, err)           

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

        integer :: itarg,iproj, n1,n2,m1,m2
        real(ESPUTR_DP) :: E_tf, E_th, Q
        intrinsic size

        err=0
        res=0
        if( (.not.allocated(fitvalues)) .or. (.not.allocated(targets_names)) &
             .or. (.not.allocated(projectiles_names)) ) then
          write(esputr_unit, *) "ERROR detected in esputr1993_check: some or all arrays are not allocated"
          res=1
          return
        end if 

        n1=size(fitvalues,2)
        n2=size(fitvalues,3)
        m1=size(targets_names,1)
        m2=size(projectiles_names,1)
        if(  ntargets /= n1 .or. nprojectiles /= n2 .or. &
             ntargets /= m1 .or. nprojectiles /= m2 ) then
          write(esputr_unit, *) "ERROR detected in esputr1993_check: inconsistent array size"
          write(esputr_unit, *) " NPROJECTILES , size(FITVALUES,1), size(TARGETS_NAMES,1) ",nprojectiles,n1,m1
          write(esputr_unit, *) " NTARGETS ,  size(FITVALUES,2), size(PROJECTILES_NAMES,1) ",ntargets,n2,m2
          res=2        
          return
        end if 
      
        do iproj=1,nprojectiles
          do itarg=1,ntargets
           e_tf=fitvalues(1,itarg,iproj)
           e_th=fitvalues(2,itarg,iproj) 
           q=fitvalues(3,itarg,iproj)
           if(e_tf <= 0) then
              res = 2
              write(esputr_unit, *) "ERROR detected in esputr1993_check: wrong E_tf"
              write(esputr_unit, *) " itarg, iproj ",itarg,iproj
              return  
           else if(e_th <= 0) then
              res = 3
              write(esputr_unit, *) "ERROR detected in esputr1993_check: wrong E_th"
              write(esputr_unit, *) " itarg, iproj ",itarg,iproj
              return  
           else if(q <= 0) then
              res = 4
              write(esputr_unit, *) "ERROR detected in esputr1993_check: wrong Q"
              write(esputr_unit, *) " itarg, iproj ",itarg,iproj
              return  
           end if 
          end do
        end do
    end subroutine esputr1993_check
    
    subroutine esputr1993_deallocate(err)

        integer, intent(out) :: err
        integer, dimension(3) :: errs
        err = 0
        errs = 0
        if(allocated(fitvalues)) deallocate(fitvalues, stat=errs(1))
        if(allocated(targets_names)) deallocate(targets_names, stat=errs(2))
        if(allocated(projectiles_names)) deallocate(projectiles_names, stat=errs(3))
        if(.not. all(errs==0)) then
            write(esputr_unit, *) "WARNING from esputr1993_deallocate: problem while deallocating memory"
            err = -200
        end if
        return
    end subroutine esputr1993_deallocate
    
!***********************************************************************************
!   GET INDEX
!***********************************************************************************

    integer function esputr1993_getprojectileid(proj, err)

        character(*), intent(in) :: proj
        integer, intent(out) :: err
        intrinsic trim
        err = 0
        
        if(.not. allocated(projectiles_names)) then
            write(esputr_unit, *) "ERROR in esputr1993_getProjectileId: arrays are not initialized" 
            write(esputr_unit, *) " Call esputr1993_init first!"
            esputr1993_getprojectileid = 0
            err =  50
            return
        end if
        
        esputr1993_getprojectileid = getindexof(proj, projectiles_names)
        
        if(esputr1993_getprojectileid <= 0) then
            write(esputr_unit, *) "ERROR in esputr1993_getProjectileId: projectile "//trim(proj)//" not found"
            err = 101
        end if
    end function esputr1993_getprojectileid
    

    integer function esputr1993_gettargetid(targ, err)

        character(*), intent(in) :: targ
        integer, intent(out) :: err
        intrinsic trim
        err = 0
        
        if(.not. allocated(targets_names)) then
            write(esputr_unit, *) "ERROR in esputr1993_getTargetId:  arrays are not initialized" 
            write(esputr_unit, *) " Call esputr1993_init first!"
            esputr1993_gettargetid = 0
            err =   50
            return
        end if
        
        esputr1993_gettargetid = getindexof(targ, targets_names)
        
        if(esputr1993_gettargetid <=0) then
            write(esputr_unit, *) "ERROR in esputr1993_getTargetId: target "//trim(targ)//" not found"
            err = 102
        end if
    end function esputr1993_gettargetid

!***********************************************************************************
!   CALCULATION OF SPUTTERING YIELDS
!***********************************************************************************    
    real(ESPUTR_DP) function esputr1993_yn(E0, proj_id, targ_id, err)

        real(ESPUTR_DP), intent(in)     :: E0
        integer, intent(in)  :: proj_id 
        integer, intent(in) ::  targ_id
        integer, intent(out) :: err
        
        real(ESPUTR_DP) :: Q, E_th, E_tf, EthOverE0, s_n_KrC, eps

        err = 0
        
        if(.not. allocated(fitvalues)) then
            write(esputr_unit,*) "ERROR in esputr1993_yn: arrays are not initialized"
            write(esputr_unit,*) " Call esputr1993_init first"
            esputr1993_yn = -1
            err =  50
            return
        end if
        if(targ_id < 1 .or. targ_id > ntargets .or.&
           proj_id < 1 .or. proj_id > nprojectiles) then
            write(esputr_unit, *) "ERROR in esputr1993_yn: targ_id or proj_id out of bounds", targ_id, proj_id
            esputr1993_yn = -1.
            err = 100
            return
        end if
        
        e_th = fitvalues(2, targ_id, proj_id)
                
        if(e0 > e_th) then
         e_tf = fitvalues(1, targ_id, proj_id)
         q    = fitvalues(3, targ_id, proj_id)
         eps = e0/e_tf
         ethovere0 = e_th / e0        
         s_n_krc = .5 * log(1 + 1.2288*eps) / (eps + .1728*sqrt(eps) + .008*eps**0.1504)        
         esputr1993_yn = q * s_n_krc * (1 - ethovere0**(2./3.)) * (1 - ethovere0)**2
        else 
         esputr1993_yn = 0.
        end if        

    end function esputr1993_yn
    
    real(ESPUTR_DP) function esputr1993_yth(theta, err)

        real(ESPUTR_DP), intent(in)  :: theta
        integer, intent(out) :: err
        
        real(ESPUTR_DP) :: cosTheta
        real(ESPUTR_DP), parameter :: f = 2., cos_th_opt= cos(75.*esputr_pi/180.) ! th_opt = 75.
                                                                                  ! Y_th(th_opt) = 3.3902

        err = 0
        
        if(theta < 0 .or. theta > esputr_pi2) then
            write(esputr_unit, *) "ERROR detected in esputr1993_yth: incorrect incident angle theta (radian) ", theta
            write(esputr_unit, *) " Incident angle cannot be < 0 or > pi/2"
            err = 103
            esputr1993_yth = -1
            return
        end if
        
        costheta = cos(theta)
        
        if(costheta > 0) then
            esputr1993_yth = costheta**(-f) * exp(f * (1-1./costheta) * cos_th_opt)
        else
            esputr1993_yth = 0
        end if
    end function esputr1993_yth


!***********************************************************************************
!   AUXILIARY FUNCTIONS
!***********************************************************************************
    
    real(ESPUTR_DP) function esputr1993_eth(proj_id, targ_id, err)

        integer, intent(in)  :: proj_id 
        integer, intent(in) ::  targ_id
        integer, intent(out) :: err

        if(.not. allocated(fitvalues)) then
            esputr1993_eth = -1
            write(esputr_unit, *) "ERROR in esputr1993_Eth: arrays are not initialized" 
            write(esputr_unit, *) "  Call esputr1993 first!"
            err = 50
        end if

        if(targ_id < 1 .or. targ_id > ntargets .or.&
           proj_id < 1 .or. proj_id > nprojectiles) then
            write(esputr_unit, *) "ERROR in esputr1993_Eth: targ_id or proj_id out of bounds", targ_id, proj_id
            esputr1993_eth = -1.
            err = 100
            return
        end if

        esputr1993_eth = fitvalues(2, targ_id, proj_id)
    end function esputr1993_eth
    
    subroutine esputr1993_availablecombinations(combinations, err)

        character(ESPUTR1993_FIRST_COL_WIDTH), dimension(:,:), allocatable, intent(inout) :: combinations
        integer, intent(out) :: err
        integer ::  t, p, i

        err = 0

        if(.not.allocated(projectiles_names) .or. .not.allocated(targets_names)) then
            err = 50
            write(esputr_unit, *) "ERROR in esputr1993_availableCombinations: arrays are not initialized" 
            write(esputr_unit, *) "  Call esputr1993_init first!"
            return
        end if

        allocate(combinations(ntargets * nprojectiles, 2), stat = err)
        if(err /= 0) then
            write(esputr_unit, *) "ERROR in esputr1993_availableCombinations: error while allocating memory", err
            err = 200
            return
        end if
        do t = 1, ntargets
            do p = 1, nprojectiles
                i = (t-1)*nprojectiles + p
                combinations(i, 1) = projectiles_names(p)
                combinations(i, 2) = targets_names(t)
            end do
        end do
    end subroutine esputr1993_availablecombinations
    
    integer function linesuntil(filePath, string, err)

        character(*), intent(in) :: filePath
        character(*), intent(in) :: string
        integer, intent(out) :: err
        character(len(string)) :: line
        integer, parameter :: UNIT = 50
        intrinsic index,trim
        err = 0
        linesuntil = 0
        
        open(unit, iostat=err, file=trim(filepath), status='old', action='read')
        if(err /= 0) then
            write(esputr_unit, *), "ERROR in linesUntil: while opening file", trim(filepath)
            close(unit)
            err = 300
            return
        end if
        
        do
            read(unit, *, iostat=err, end=9100) line
            if(err /= 0) then
                write(esputr_unit, *), "ERROR in linesUntil: while reading file", trim(filepath)
                close(unit)
                err = 300
                return
            end if
            linesuntil = linesuntil + 1
            if(len_trim(line)>0) then
             if(index(line,string)>0) then
                 close(unit)
                 return
             end if
            end if
        end do
9100    linesuntil = 0
        close(unit)
    end function linesuntil
    
    pure integer function getindexof(subj, list)

        character(*), intent(in) :: subj 
        character(*), dimension(:), intent(in) :: list
        integer :: i
        intrinsic size
        getindexof = 0
        do i = 1, size(list, 1)
            if(list(i) == subj) then
                getindexof = i
                exit
            end if
        end do
    end function getindexof
    
end module esputr1993