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 5
List set_ODE_Rcpp_LM_OSL(double t, arma::vec n, Rcpp::List parameters) {
14

15
  //unpack parameters for ODEs
16

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

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

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

34 5
  int j = 0;
35 5
  int jj = 0;
36 5
  for (std::size_t i = 0; i < N.size(); ++i){
37 5
    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 5
      j++;
39 5
      jj++;
40

41 5
      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 5
      jj++;
44 5
      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 5
  arma::vec temp_dn1 = dn.subvec(0,j-1);
50 5
  arma::vec temp_dn2 = dn.subvec(j,jj-1);
51

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

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

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

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

Read our documentation on viewing source code .

Loading