final commit for ape 3.0

[ape.git] / src / rTrait.c

diff --git a/src/rTrait.c b/src/rTrait.c
index 6441e161b9c07413bb97e2cb4cdba48edce20788..43dc73ad85f3163b65e49333889a2088d30a06a5 100644 (file)
--- a/src/rTrait.c
+++ b/src/rTrait.c
@@ -1,6 +1,6 @@
-/* rTrait.c       2010-01-11 */
+/* rTrait.c       2011-06-25 */
 
-/* Copyright 2010 Emmanuel Paradis */
+/* Copyright 2010-2011 Emmanuel Paradis */
 
 /* This file is part of the R-package `ape'. */
 /* See the file ../COPYING for licensing issues. */
@@ -12,17 +12,26 @@ void rTraitCont(int *model, int *Nedge, int *edge1, int *edge2, double *el,
 {
 /* The tree must be in pruningwise order */
        int i;
+       double alphaT, M, S;
 
        switch(*model) {
        case 1 : for (i = *Nedge - 1; i >= 0; i--) {
                        GetRNGstate();
-                       x[edge2[i]] = x[edge1[i]] + el[i] * sigma[i] * norm_rand();
+                       x[edge2[i]] = x[edge1[i]] + sqrt(el[i]) * sigma[i] * norm_rand();
                        PutRNGstate();
                }
                break;
        case 2 : for (i = *Nedge - 1; i >= 0; i--) {
+                       if (alpha[i]) {
+                               alphaT = alpha[i] * el[i];
+                               M = exp(-alphaT);
+                               S = sigma[i] * sqrt((1 - exp(-2 * alphaT))/(2 * alpha[i]));
+                       } else { /* same than if (alpha[i] ==0) */
+                               M = 1;
+                               S = sqrt(el[i]) * sigma[i];
+                       }
                        GetRNGstate();
-                       x[edge2[i]] = x[edge1[i]] + (sigma[i]*norm_rand() - alpha[i]*(x[edge1[i]] - theta[i])) * el[i];
+                       x[edge2[i]] = x[edge1[i]] * M +  theta[i] * (1 - M) + S * norm_rand();
                        PutRNGstate();
                }
                break;