3D integration

#pragma rtGlobals=1     // Use modern global access method.

Function do3dIntegration(xmin,xmax,ymin,ymax,zmin,zmax)
    Variable xmin,xmax,ymin,ymax,zmin,zmax 
    Variable/G globalXmin=xmin
    Variable/G globalXmax=xmax
    Variable/G globalYmin=ymin
    Variable/G globalYmax=ymax
    Variable/G globalY, globalZ
    return  Integrate1D(userFunction3,zmin,zmax,1)      // Romberg integration
End

Function userFunction3(inZ)
    Variable inZ   
    NVAR globalZ=globalZ
    globalZ=inZ
    NVAR globalYmin=globalYmin
    NVAR globalYmax=globalYmax 
    return integrate1D(userFunction2,globalYmin,globalYmax,1)   // Romberg integration
End

Function userFunction2(inY)
    Variable inY   
    NVAR globalY=globalY
    globalY=inY
    NVAR globalXmin=globalXmin
    NVAR globalXmax=globalXmax 
    return integrate1D(userFunction1,globalXmin,globalXmax,1)   // Romberg integration
End

Function userFunction1(inX)
    Variable inX   
    NVAR globalY=globalY
    NVAR globalZ=globalZ   
    return ((inX^2+globalY^2+globalZ^2)^(3/2)) * exp(-(inX^2+globalY^2+globalZ^2))
//  r^3 exp(-r^2)  --> 2*pi * u^2 *exp(-u) --> 4 pi integral
//  r exp(-r^2)     --> 2*pi * u    *exp(-u) --> 2 pi integral
End

•printf "%8.6f\r", do3dIntegration(-5,5,-5,5,-5,5)
  12.566372