X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=src%2FrTrait.c;h=43dc73ad85f3163b65e49333889a2088d30a06a5;hb=dc2dda7eb763b4e140f1e5adb7fa8bfaa2661f6d;hp=11cc63e2a0bdb2a66105873e690ab1e03d115074;hpb=f295ab19440298e543db5a270e54f10a84382197;p=ape.git

diff --git a/src/rTrait.c b/src/rTrait.c
index 11cc63e..43dc73a 100644
--- a/src/rTrait.c
+++ b/src/rTrait.c
@@ -1,6 +1,6 @@
-/* rTrait.c       2010-07-26 */
+/* 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,6 +12,7 @@ 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--) {
@@ -21,14 +22,16 @@ void rTraitCont(int *model, int *Nedge, int *edge1, int *edge2, double *el,
 		}
 		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]] - alpha[i]*el[i]*(x[edge1[i]] - theta[i]) + sqrt(el[i])*sigma[i]*norm_rand();
-			PutRNGstate();
-		}
-		break;
-	case 3 : for (i = *Nedge - 1; i >= 0; i--) {
-			GetRNGstate();
-			x[edge2[i]] = x[edge1[i]] - (1 - exp(alpha[i]*el[i]))*(x[edge1[i]] - theta[i]) + sqrt(el[i])*sigma[i]*norm_rand();
+			x[edge2[i]] = x[edge1[i]] * M +  theta[i] * (1 - M) + S * norm_rand();
 			PutRNGstate();
 		}
 		break;