I want to integrate a 3d array over space in r, theta and phi (spherical polar). For 1d I use Simpson's 1/3rd rule but I am confused about that for 3d. Also, would you like to suggest any other method for integration or subroutine? I am using Fortran 95.
I have written the Fortran code for integration in 3d, I thought I should share with you. The code for calculating integration of a function is 3 dimension is:
!This program uses Simpson's 1/3 method to calulate volume
integral in r,theta & phi.
program SimpsonInteg3d
implicit none
integer::i,j,k
integer, parameter :: N=10,M=360,L=180
integer, parameter:: rmin=0,rmax=N,phimin=0,phimax=M,&
thetamin=0,thetamax=L
double precision,&
dimension(rmin:rmax,thetamin:thetamax,phimin:phimax)::U
real*8, parameter :: pi = 4*atan(1.0),dr=1./N,&
dtheta=pi/(L),dphi=2*pi/M
real*8 :: r(rmin:rmax)=(/(i*dr,i=rmin,rmax)/),&
theta(thetamin:thetamax)=(/(j*dtheta,j=thetamin,thetamax)/),&
p(phimin:phimax)=(/(k*dphi,k=phimin,phimax)/)
real*8::intg
do i=rmin,rmax
do j=thetamin, thetamax
do k=phimin,phimax
!The function which has to be integrated.
U(i,j,k)=r(i)* (sin((p(k)))**2) *sin(theta(j))
enddo
enddo
enddo
call Integration(Intg,U,r,theta,p)
print*,"Integration of function U using simpson's 1/3=", Intg
end program
!===============================================================!
!Subroutine for calculating integral of a function in 3d.
subroutine Integration(Intg,U,r,theta,p)
implicit none
integer::i,j,k
integer, parameter :: N=10,M=360,L=180
integer, parameter ::rmin=0,rmax=N,&
phimin=0,phimax=M,thetamin=0,thetamax=L
double precision,&
dimension(rmin:rmax,thetamin:thetamax,phimin:phimax):: U
real*8::
r(rmin:rmax),theta(thetamin:thetamax),p(phimin:phimax),Intg,Ia
double precision,dimension(rmin:rmax)::Itheta
real*8, parameter :: pi = 4*atan(1.0),dr=1./N,&
dtheta=pi/(L),dphi=2*pi/M
Intg=0
Ia=0
do i=rmin+1,rmax-1
call Integtheta(Itheta,i,U,r,theta,p)
if(mod(i,2).eq.0) then
Ia = Ia + 2*Itheta(i)*r(i)**2
else
Ia = Ia + 4*Itheta(i)*r(i)**2
endif
end do
call Integtheta(Itheta,rmin,U,r,theta,p)
call Integtheta(Itheta,rmax,U,r,theta,p)
Intg=(dr/3)*(Itheta(rmin)+Itheta(rmax)+ Ia)
end subroutine Integration
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Subroutine for calculating integral of U along theta and phi
subroutine Integtheta(Itheta,i,U,r,theta,p)
implicit none
integer::i,j,k
integer, parameter :: N=10,M=360,L=180
integer, parameter ::rmin=0,rmax=N,&
phimin=0,phimax=M,thetamin=0,thetamax=L
double precision,&
dimension(rmin:rmax,thetamin:thetamax,phimin:phimax)::U
real*8:: r(rmin:rmax),theta(thetamin:thetamax),p(phimin:phimax)
double precision,dimension(rmin:rmax)::Itheta,Itha
double precision,dimension(rmin:rmax,thetamin:thetamax)::Ip
real*8, parameter :: pi = 4*atan(1.0),dr=1./N,&
dtheta=pi/(L),dphi=2*pi/M
Itheta(i)=0
Itha(i)=0
do j=thetamin+1,thetamax-1
call Integphi(Ip,i,j,U,r,theta,p)
if(mod(j,2).eq.0) then
Itha(i) = Itha(i) + 2*Ip(i,j)*sin(theta(j))
else
Itha(i) = Itha(i) + 4*Ip(i,j)*sin(theta(j))
endif
end do
call Integphi(Ip,i,thetamin,U,r,theta,p)
call Integphi(Ip,i,thetamax,U,r,theta,p)
Itheta(i)=(dtheta/3)*(Ip(i,thetamin)+Ip(i,thetamax)+ Itha(i))
end subroutine Integtheta
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Subroutine for calculating integral of U along phi
subroutine Integphi(Ip,i,j,U,r,theta,p)
implicit none
integer::i,j,k
integer, parameter :: N=10,M=360,L=180
integer, parameter ::rmin=0,rmax=N,&
phimin=0,phimax=M,thetamin=0,thetamax=L
double precision,&
dimension(rmin:rmax,thetamin:thetamax,phimin:phimax)::U
real*8:: r(rmin:rmax),theta(thetamin:thetamax),p(phimin:phimax)
double precision,dimension(rmin:rmax,thetamin:thetamax)::Ip,Ipa
real*8, parameter :: pi = 4*atan(1.0),dr=1./N,&
dtheta=pi/(L),dphi=2*pi/M
Ipa(i,j)=0
do k=phimin+1,phimax-1
if(mod(k,2).eq.0) then
Ipa(i,j) = Ipa(i,j) + 2*U(i,j,k)
else
Ipa(i,j)= Ipa(i,j) + 4*U(i,j,k)
endif
end do
Ip(i,j)=(dphi/3)*(U(i,j,phimin)+U(i,j,phimax)+ Ipa(i,j))
end subroutine Integphi
It calculates the integration of the function U along phi first and then uses the function Ip to calculate integral along theta. Then finally the function Itheta is used to calculate integration along r.