Error while converting a logic in C++ to Excel VBA

[Juggler_IN]

I have the following c++ code for generating q Gaussian distribution:

<rcpp.h><math.h><iostream.h>
double qnormal_x,qnormal_y,qnormal_z;
double expq(double q, double w){
if(q==1.0){
return(exp(w));
}
else{
return (expl(log(1.0+(1.0-q)*w)/(1.0-q)));
}
}
double lnq(double q, double w){
if(q==1.0){
return(log(w));
}
else{
return ((exp(log(w)*(1.0-q))-1.0)/(1.0-q));
}
}
void setseed_qnormal(double v0, double z0){
qnormal_x = sqrt(1-v0*v0);
qnormal_y = v0;
qnormal_z = z0;
}
double Q8(double w, double v){
return(8*w*v*(((16.0*w*w-24.0)*w*w+10.0)*w*w-1.0));
}
double P8(double w){
return((((128.0*w*w-256.0)*w*w+160.0)*w*w-32.0)*w*w+1.0);
}
double f(double z){
return(1.0-fabs(1.0-1.99999*z));
}
void qnormal(double q){
double qq;
qnormal_y = Q8(qnormal_x,qnormal_y);
qnormal_x = P8(qnormal_x);
qq = (q+1.0)/(3.0-q);
qnormal_z = f(expq(qq,-qnormal_z*qnormal_z*0.5));
qnormal_z = sqrt(-2.0*lnq(qq,qnormal_z));
}  


// [[Rcpp::export]]
NumericVector Chaotic(int n,double q,double v0, double z0){
  double eta;
  double xi;
  NumericVector ux(n);
setseed_qnormal(v0,z0);
for(int i=0;i<n;i++){
qnormal(q);
xi = qnormal_x*qnormal_z;
eta = qnormal_y*qnormal_z;
ux[i]= xi;
}
return ux;
}

My VBA conversion of the above code is:

Public Function Chaotic(Optional ByVal q As Double = 1#) As Double

    If q > 3 Or q = 3 Then Exit Function


    Dim u1, u2, qq, z As Double
    Dim qnormal_x As Double, qnormal_y As Double, qnormal_z As Double


    u1 = Application.WorksheetFunction.NormInv(RND(), 0, 1) 'v0
    u2 = Application.WorksheetFunction.NormInv(RND(), 0, 1) 'z0


    qnormal_x = Sgn((1 - u1 * u1)) * Sqr(Abs(1 - u1 * u1))
    qnormal_y = u1
    qnormal_z = u2
    qnormal_y = (8 * qnormal_x * qnormal_y * (((16 * qnormal_x * qnormal_x - 24) * qnormal_x * qnormal_x + 10) * qnormal_x * qnormal_x - 1))
    qnormal_x = ((((128 * qnormal_x * qnormal_x - 256) * qnormal_x * qnormal_x + 160) * qnormal_x * qnormal_x - 32) * qnormal_x * qnormal_x + 1)
    qq = (1 + q) / (3 - q)
    qnormal_z = (1 - Abs(1 - 1.99999 * (qEXP(qq, -qnormal_z * qnormal_z * 0.5))))
    If qq = 1 Then z = Log(qnormal_z) Else z = ((Exp(Sgn(qnormal_z) * Log(Abs(qnormal_z)) * (1 - qq)) - 1) / (1 - qq))
    Chaotic = Sgn(-2 * z) * Sqr(Abs(-2 * z))


End Function

But the conversion is not giving the desired output in terms of mean and std. deviation as per wikipedia page

https://en.wikipedia.org/wiki/Q-Gaussian_distribution

I am not able to figure out the gap in my conversion.</n;i++){
</iostream.h></math.h></rcpp.h>