esputr_YTH.F

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/>.
!

#include "fintrf.h"

      subroutine mexfunction(nlhs, plhs, nrhs, prhs)
      use esputr
      use esputr1993
      use esputr2001
      implicit none
      integer :: nlhs 
      integer :: nrhs
      mwpointer, dimension(*) :: plhs
      mwpointer, dimension(*) :: prhs

      mwpointer :: mxgetpr, mxcreatedoublematrix, mxgetstring
      integer :: mxIsNumeric, mxIsChar
      mwsize :: mxgetm, mxgetn
      
      mwsize :: n_e0s, n_ths, n_projs, n_targs, n_combinations
      integer :: err, i, p, t, c, ind
      logical :: ifexist
      real(kind=8) :: Emin, Emax
      mwpointer :: e0s_pr, ths_pr, ys_pr
      real(kind=8), dimension(:), allocatable :: E0s, THs
      real(kind=8), dimension(:,:), allocatable :: Ys
      integer, parameter :: MAX_STR_LEN = 255
      character(MAX_STR_LEN) :: projsstr, targsstr
      character(MAX_STR_LEN), dimension(:), allocatable :: projs, targs
      mwsize :: projs_length, targs_length
      character(MAX_STR_LEN) :: version, fileTH
      
      if(nrhs < 5) then
          call mexerrmsgtxt('At least 5 arguments required.')
      elseif (mxisnumeric(prhs(1)) /= 1) then
          call mexerrmsgtxt('Input # 1 is not a numeric array.')
      elseif(mxisnumeric(prhs(2)) /= 1) then
          call mexerrmsgtxt('Input # 2 is not a numeric array.')
      elseif(mxischar(prhs(3)) /= 1) then
          call mexerrmsgtxt('Input # 3 must be a string.')
      elseif(mxischar(prhs(4)) /= 1) then
          call mexerrmsgtxt('Input # 4 must be a string.')
      elseif(mxischar(prhs(5)) /= 1) then
          call mexerrmsgtxt('Input # 5 must be a string.')
      elseif(nrhs > 5 .and. mxischar(prhs(6)) /= 1) then
          call mexerrmsgtxt('Input # 6 must be a string.')
      endif

      n_ths = mxgetm(prhs(1))
      if(mxgetn(prhs(1)) /= 1) then
          call mexerrmsgtxt('Input #1 must be a column vector.')
      endif

      n_e0s = mxgetn(prhs(2))
      if(mxgetm(prhs(2)) /= 1) then
          call mexerrmsgtxt('Input #2 must be a row vector.')
      endif

      n_projs = mxgetm(prhs(3))
      projs_length = mxgetn(prhs(3))
      n_targs = mxgetm(prhs(4))
      targs_length = mxgetn(prhs(4))
      
      n_combinations = max(n_projs, n_targs)
      
      allocate(projs(n_combinations))
      allocate(targs(n_combinations))
      
      ths_pr = mxgetpr(prhs(1))
      e0s_pr = mxgetpr(prhs(2))
      ! If the string array contains several rows, they are copied, one column at a time, into one long string array.
      err = mxgetstring(prhs(3), projsstr, max_str_len)
      projs = " "
      foreachprojectile: do p = 1, n_projs
          foreachchar: do c = 1, projs_length
              ind = (c-1)*n_projs+p
              projs(p)(c:c) = projsstr(ind:ind)
          end do foreachchar
      end do foreachprojectile
      projs(n_projs:n_combinations) = projs(n_projs)
      
      err = mxgetstring(prhs(4), targsstr, max_str_len)
      targs = " "
      foreachtarget: do t = 1, n_targs
          foreachchar2: do c = 1, targs_length
              ind = (c-1)*n_targs+t
              targs(t)(c:c) = targsstr(ind:ind)
          end do foreachchar2
      end do foreachtarget
      targs(n_targs:n_combinations) = targs(n_targs)
      
      err = mxgetstring(prhs(5), version, max_str_len)
      if(version.ne."2001".and.version.ne."1993") then
        call mexerrmsgtxt('unknown model: '//trim(version))
      end if

      if(version == "2001" .and. nrhs == 6) then
         err = mxgetstring(prhs(6), fileth, max_str_len)
      else if(version == "2001" .and. nrhs < 6) then
         call mexerrmsgtxt('6 arguments required')
      end if
      
      allocate(e0s(n_e0s))
      allocate(ths(n_ths))
      allocate(ys(n_ths,n_e0s))
      call mxcopyptrtoreal8(e0s_pr, e0s, n_e0s)
      call mxcopyptrtoreal8(ths_pr, ths, n_ths)
      
      if(version == "2001") then
          ! check if input file exists
          inquire(file=fileth, exist=ifexist, iostat=err) 
          if (.not.ifexist .or. err.ne.0)
     &               call mexerrmsgtxt('File not found: '//trim(fileth))
          call esputr2001_initth(fileth, err)
      end if
      if(err /= 0) then
          call esputr2001_deallocate(err)
          call mexerrmsgtxt('initialization error in ESPUTR')
      endif
      
      foreachcombination: do c=1, n_combinations
          call esputr_yth(ths, n_ths, e0s, n_e0s, projs(c),
     +                                       targs(c), version, ys, err)
          if(err /= 0) then
              if(err == 103 .or. err == 157) then
                  call mexerrmsgtxt('incorrect incident angle '//
     &                              '(it must be 0 <= theta <= pi/2)')
              else if(err == 152) then
                  call mexerrmsgtxt('projectile-target combination '//
     &                trim(projs(c))//'-'//trim(targs(c))//' not found')
              else if(err == 156) then
                  call mexerrmsgtxt('incident energy must be > 0')
              else if(err == 158) then
                  call esputr2001_initth(fileth, err)
                  call esputr2001_getavailableenergyrangebyname(
     &                     projs(c), targs(c), emin, emax, err)
                  write(projsstr,*) emin
                  write(targsstr,*) emax
                  call mexerrmsgtxt('no fiting parameters for all '//
     &     'requested energies for the projectile-target combination '//
     &         trim(projs(c))//'-'//trim(targs(c))//
     &         '; available energy range: ['//
     &         trim(projsstr)//','//trim(targsstr)//'] eV')
              else
                  call mexerrmsgtxt('error in ESPUTR')
              endif
              call esputr2001_deallocate(err)
          endif
          
          plhs(c) = mxcreatedoublematrix(n_ths,n_e0s,0)
          ys_pr = mxgetpr(plhs(c))
          call mxcopyreal8toptr(ys,ys_pr,n_e0s*n_ths)
      end do foreachcombination
      
      call esputr2001_deallocate(err)
      
      deallocate(e0s)
      deallocate(ys)
      
      end subroutine mexfunction