R-Lum / RLumModel
1
//set_ODE_Rcpp_LM_OSL.cpp
2
//author: Johannes Friedrich, University of Bayreuth (Germany)
3
//version: 0.1.0 [2016-03-18]
4
//Function calculates the ODEs for LM-OSL for all quartz luminescence models iterativly
5
//
6

7

8
// [[Rcpp::depends(RcppArmadillo)]]
9
#include <RcppArmadillo.h>
10
using namespace Rcpp;
11

12
//  [[Rcpp::export(".set_ODE_Rcpp_LM_OSL")]]
13 4
List set_ODE_Rcpp_LM_OSL(double t, arma::vec n, Rcpp::List parameters) {
14

15
  //unpack parameters for ODEs
16

17 6
  arma::vec N = parameters["N"];
18 6
  arma::vec E = parameters["E"];
19 6
  arma::vec s = parameters["s"];
20 6
  arma::vec A = parameters["A"];
21 6
  arma::vec B = parameters["B"];
22 6
  arma::vec Th = parameters["Th"];
23 6
  arma::vec E_th = parameters["E_th"];
24

25 6
  double k_B = parameters["k_B"];
26

27 6
  double R = parameters["R"];
28 6
  double P = parameters["P"];
29 6
  double temp = parameters["temp"];
30 6
  double b = parameters["b"];
31 6
  double a = parameters["a"];
32 6
  arma::vec dn(N.size()+2);
33

34 6
  int j = 0;
35 6
  int jj = 0;
36 6
  for (std::size_t i = 0; i < N.size(); ++i){
37 6
    if (B[i] == 0)    {//use recombination propability of recombination centers to identify electron traps, because they had no recombination propability to recomibnation centers from conduction band
38 6
      j++;
39 6
      jj++;
40

41 6
      dn[i] = n[N.size()]*(N[i]-n[i])*A[i] - n[i]*P*a*t*Th[i]*exp(-E_th[i]/(k_B*(273+temp+b*t))) - n[i]*s[i]*exp(-E[i]/(k_B*(273+temp+b*t)));
42 2
    } else {//calculate recombination centers
43 6
      jj++;
44 6
      dn[i] = n[N.size()+1]*(N[i]-n[i])*A[i] - n[i]*s[i]*exp(-E[i]/(k_B*(273+temp+b*t))) - n[N.size()]*n[i]*B[i];
45
    }
46
  }
47

48
  //build sub-vectors for conduction/valenece band calculation
49 6
  arma::vec temp_dn1 = dn.subvec(0,j-1);
50 6
  arma::vec temp_dn2 = dn.subvec(j,jj-1);
51

52 6
  arma::vec temp_n = n[N.size()]*n.subvec(j,jj-1);
53 6
  arma::vec temp_B = B.subvec(j,jj-1);
54

55
  //conduction band
56 6
  dn[N.size()] = R - sum(temp_dn1) - sum(arma::trans(temp_n) * temp_B);
57

58
  //valence band
59 6
  dn[N.size()+1] = R - sum(temp_dn2) - sum(arma::trans(temp_n) * temp_B);
60

61 6
  return(Rcpp::List::create(dn));
62
}

Read our documentation on viewing source code .

Loading