## 3D integration

Average rating
(1 vote)

Igor provides a function "Integrate1D()" for integrating functions of one-dimensional arguments. The associated help file gives a method for using it to integrate functions of two-dimensional arguments, with the statement "This method can be extended to higher dimensions." I found the effort to be not entirely trivial, and thought to save others the work with the following code. I have shown an interestingly complex user function of (x,y,z) that has a simple analytic result for large integration limits after converting to a radial coordinate system

```#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```

The exact result should be 4 *pi.