eiram_octave.cpp

Go to the documentation of this file.

//    Copyright (c) 2016  Forschungszentrum Juelich GmbH
//                        Markus Brenneis, Vladislav Kotov 
//
//    This file is part of EIRAM.
//
//    EIRAM 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.
//
//    EIRAM 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 EIRAM.  If not, see <http://www.gnu.org/licenses/>.

#include <octave/oct.h>
#include <cstring>
#include <string>
#include "eiram_octave.hpp"

DEFUN_DLD (eiram, args, nargout,
"\n\
 Calculate reaction rates and cross-sections described by polynomial fit formulas \n\
 (interface to fortran module EIRAM) \n\
\n\
  [Y1 Y2 .. ] = eiram(X1, [X2,] name, type, index, [path]) \n\
\n\
   Y1, Y2 ... : vectors or matrices with the results, \n\
                each Y1, Y2 ..., corresponds to one reaction defined in INDEX, \n\
                each column in Y1, Y2 ... corresponds to one value of X1 \n\
   X1    : row vector, first argument of the fit \n\
   X2    : optional,  column vector, second argument of the fit \n\
   name  : string, name of the file with the data-set ('HYDHEL', 'AMJUEL' etc.) \n\
   type  : string, type of the reaction ('H.1','H.3','H.4' etc.) \n\
   index : array of strings, index of the reaction ('3.1.8', '2.1.5' etc.) \n\
   path  : optional, string, complete path to the file defined by NAME, \n\
           if PATH is omitted than it is assumed that the data-set is in the current folder \n\
\n\
   Units : energy and temperature in [eV], density in [cm^-3], \n\
           cross sections in cm^2, collision rates in cm^3/s,  \n\
           energy loss rates in cm^3/s*eV \n\
\n\
  Examples: \n\
  % Charge-exchange H + p, Ebeam=10 eV  \n\
    svcx = eiram(10,[1:100],'HYDHEL','H.3','3.1.8'); \n\
\n\
  % Momentum transfer cross-sections for elastic collisions p+H, p+He, p+H2 \n\
   [sel_D sel_He sel_D2] = eiram([0.1:10],'AMJUEL','H.1',['0.1D'; '0.2D'; '0.3D'],'.../eiram/data/'); \n\
\n\
  % Ionization rate of H atoms, electron density 1e12 cm-3 and 1e14 cm-3  \n\
    svi = eiram([1e12 1e14],[0.1:0.1:10 10+1:1:100],'AMJUEL','H.4','2.1.5');"){
    octave_value_list retval;
    int nargin = args.length();
    
    if(nargin < 4) { error("At least 4 input arguments required"); return retval;}

    if(!args(0).is_numeric_type()) { error("1st argument must be of numeric type"); return retval;}

    int numberOfArgs;

    if(args(1).is_numeric_type()) {  numberOfArgs = 2; }
    else {  numberOfArgs = 1; }

    if ( (nargin-numberOfArgs) < 3 ) { error("Too few input arguments"); }

    NDArray Arg1s = args(0).array_value();
    double* Arg1s_v = Arg1s.fortran_vec();      

    if(error_state) {return retval;}
    octave_idx_type n_Arg1s = Arg1s.numel();

    std::string fileName = args(numberOfArgs).char_matrix_value().row_as_string(0);
    std::string reactionType = args(numberOfArgs+1).char_matrix_value().row_as_string(0);
    charMatrix reactionIndex = args(numberOfArgs+2).char_matrix_value();
    if(error_state) {return retval;}

    std::string filePath;
    if( (nargin-numberOfArgs) > 3 ) {filePath = args(numberOfArgs+3).char_matrix_value().row_as_string(0);}

    int err = 0;
            
    if(error_state) {return retval;} 

     F77_XFCN(eiram_load_wrapper, FORTSUB, (filePath.c_str(), fileName.c_str(), err, filePath.length(), fileName.length()));
            
    if(err != 0) {
       char str[1024];           
      switch(err){
        case 310 : 
          sprintf(str,"Input file %s%s is not found",filePath.c_str(),fileName.c_str());
          error(str);
          break;
        case 300 :
      sprintf(str,"Cannot read input file %s%s",filePath.c_str(),fileName.c_str());
          error(str);
          break;
        default :
          error("error during initialization of EIRAM");
      }
      F77_XFCN(eiram_deallocate_wrapper, FORTSUB, (err));
      return retval;
    }

    octave_idx_type n_reactions = reactionIndex.rows();

    for(int i=0;i<n_reactions;i++) {    
        std::string index=reactionIndex.row_as_string(i);

        if( numberOfArgs == 1) {
         Matrix Ys(n_Arg1s, 1);

         F77_XFCN(eiram_matlab_calc1, FORTSUB,
                 (Ys.fortran_vec(),
                  Arg1s_v, n_Arg1s,
                  fileName.c_str(),reactionType.c_str(),
                  index.c_str(),err
                  F77_CHAR_ARG_LEN(fileName.length())
                  F77_CHAR_ARG_LEN(reactionType.length())
                  F77_CHAR_ARG_LEN(index.length())
          )); 
         if(err == 0) {retval(i) = Ys;}

        }
    else if( numberOfArgs == 2) {  
         NDArray Arg2s=args(1).array_value();
         double* Arg2s_v = Arg2s.fortran_vec();  
         octave_idx_type n_Arg2s = Arg2s.numel();
         Matrix Ys(n_Arg2s, n_Arg1s);     

         F77_XFCN(eiram_matlab_calc2, FORTSUB,
                 (Ys.fortran_vec(),
                  Arg1s_v, Arg2s_v,n_Arg1s,n_Arg2s,
                  fileName.c_str(),reactionType.c_str(),
                  index.c_str(),err
                  F77_CHAR_ARG_LEN(fileName.length())
                  F77_CHAR_ARG_LEN(reactionType.length())
                 F77_CHAR_ARG_LEN(index.length())
               ));
         if(err == 0) {retval(i) = Ys;}
        }

       switch(err){
         case 11 : 
              error("error : negative argument");
              break;
         case 12 :
              error("error : negative 1st argument");
              break;
         case 13 :
              error("error : negative 2nd argument");
              break;
         case 110 :
       char str[128];           
       sprintf(str,"error: requested reaction %s %s %s is not found",fileName.c_str(),reactionType.c_str(),index.c_str());
           error(str);
           break; 
         case 120 :
           error("error : double fit is requested for reaction described by a single fit");
           break; 
         case 130 :
           error("error : single fit is requested for reaction described by a double fit");
           break; 
         default :
          if(err != 0) {error("error during execution of EIRAM");}
       }      
       if(err != 0) {
          F77_XFCN(eiram_deallocate_wrapper, FORTSUB, (err));
          return retval; }
    } //for(int i=0;i<n_reactions;i++)
   
    F77_XFCN(eiram_deallocate_wrapper, FORTSUB, (err));
    return retval;
}