esputr_YN.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_projs, n_targs, n_combinations
      integer :: err, i, p, t, c, ind
      mwpointer :: e0s_pr, ys_pr
      real(kind=8), dimension(:), allocatable :: E0s
      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, fileN
      logical :: ifexist
      integer :: io
      intrinsic trim

      if (nrhs /= 5) then
          call mexerrmsgtxt('5 arguments required.')
      elseif (mxisnumeric(prhs(1)) /= 1) then
          call mexerrmsgtxt('Input # 1 is not a numeric array.')
      elseif(mxischar(prhs(2)) /= 1) then
          call mexerrmsgtxt('Input # 2 must be a string.')
      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.')
      endif

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

      n_projs = mxgetm(prhs(2))
      projs_length = mxgetn(prhs(2))
      n_targs = mxgetm(prhs(3))
      targs_length = mxgetn(prhs(3))
      
      n_combinations = max(n_projs, n_targs)
           
      allocate(projs(n_combinations))
      allocate(targs(n_combinations))
      
      e0s_pr = mxgetpr(prhs(1))
      ! If the string array contains several rows, they are copied, one
      ! column at a time, into one long string array. So we must
      ! "transpose" the strings here.
      err = mxgetstring(prhs(2), 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(3), 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(4), version, max_str_len)
      err = mxgetstring(prhs(5), filen, max_str_len)
      
      allocate(e0s(n_e0s))
      allocate(ys(n_e0s))
      call mxcopyptrtoreal8(e0s_pr, e0s, n_e0s)
            
      esputr_unit=0 !VK in order to print ESPUTR messages to the terminal

      !VK check if input file exist
      inquire(file=filen,exist=ifexist,iostat=io) 
      if (.not.ifexist.or.io.ne.0) 
     f call mexerrmsgtxt('File not found: '//trim(filen))

      if(version == "1993") then
          call esputr1993_init(filen, err)
      else if(version == "2001") then
          call esputr2001_initn(filen, err)
      else
          call mexerrmsgtxt('unknown model: '//trim(version))
      end if
      if(err /= 0) then
          call esputr1993_deallocate(err)
          call esputr2001_deallocate(err)
          call mexerrmsgtxt('initialization error in ESPUTR')
      endif
      
      foreachcombination: do c=1, n_combinations
          call esputr_yn(e0s, n_e0s, projs(c), targs(c), version,
     &                                                          ys, err)
          if(err /= 0) then
              if(err == 101) then
                  call mexerrmsgtxt('projectile '//trim(projs(c))//
     &                                                     ' not found')
              else if(err == 102) then
                  call mexerrmsgtxt('target '//trim(targs(c))//
     &                                                     ' not found')
              else if(err == 151) then
                  call mexerrmsgtxt('projectile-target combination '//
     &                trim(projs(c))//'-'//trim(targs(c))//' not found')
              else
                  call mexerrmsgtxt('error in ESPUTR')
              endif
              call esputr1993_deallocate(err)
              call esputr2001_deallocate(err)
          endif          
          plhs(c) = mxcreatedoublematrix(1,n_e0s,0)
          ys_pr = mxgetpr(plhs(c))
          call mxcopyreal8toptr(ys,ys_pr,n_e0s)
      end do foreachcombination
      
      call esputr1993_deallocate(err)
      call esputr2001_deallocate(err)
      
      deallocate(e0s)
      deallocate(ys)
      
      end subroutine mexfunction