function plotqgaussian()

x = -4.1:.01:4.1;

figure; hold on;
plot(x, qgaussian(x,0,1), '-r','LineWidth',2)    
plot(x, qgaussian(x,1,1), '-k','LineWidth',2) 
plot(x, qgaussian(x,2,1), '-b','LineWidth',2) 
plot(x, qgaussian(x,2,3), '-g','LineWidth',2) 

function y = qgaussian(x,q,beta)

if q == 1;
    Cq = sqrt(pi);
elseif q < 1
    Cq = 2*sqrt(pi)/((3-q)*sqrt(1-q)) * exp(gammaln(1/(1-q)) - gammaln((3-q)/(2*(1-q))));
else 
    Cq = sqrt(pi)/sqrt(q-1) * exp(gammaln((3-q)/(2*(q-1))) - gammaln(1/(q-1)));    
end

y = sqrt(beta)/Cq * qexp(-beta*x.^2,q);

if q < 1
    y( abs(x) > (1/sqrt(beta*(1-q))) ) = 0;
end

function y = qexp(x,q)

if q == 1
    y = exp(x);
else
    y = (1+(1-q)*x).^(1/(1-q));
end