subroutine fourn(data,nn,ndim,isign) c c Replaces DATA by its NDIM-dimensional discrete Fourier transform, c if ISIGN is input as 1. NN is an integer array of length NDIM, c containing the lengths of each dimension (number of complex values), c which must all be powers of 2. DATA is a real array of length twice c the product of these lengths, in which the data are stored as in a c multidimensional complex Fortran array. If ISIGN is input as -1, c DATA is replaced by its inverse transform times the product of the c lengths of all dimensions. From Press, W.H., Flannery, B.P., c Teukolsky, S.A., and Vetterling, W.T., 1986, Numerical Recipes, c Cambridge Univ. Press, p. 451-453. c real*8 wr,wi,wpr,wpi,wtemp,theta dimension nn(ndim),data(*) ntot=1 do 11 iidim=1,ndim ntot=ntot*nn(iidim) 11 continue nprev=1 do 18 iidim=1,ndim n=nn(iidim) nrem=ntot/(n*nprev) ip1=2*nprev ip2=ip1*n ip3=ip2*nrem i2rev=1 do 14 i2=1,ip2,ip1 if(i2.lt.i2rev)then do 13 i1=i2,i2+ip1-2,2 do 12 i3=i1,ip3,ip2 i3rev=i2rev+i3-i2 tempr=data(i3) tempi=data(i3+1) data(i3)=data(i3rev) data(i3+1)=data(i3rev+1) data(i3rev)=tempr data(i3rev+1)=tempi 12 continue 13 continue endif ibit=ip2/2 1 if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then i2rev=i2rev-ibit ibit=ibit/2 go to 1 endif i2rev=i2rev+ibit 14 continue ifp1=ip1 2 if(ifp1.lt.ip2)then ifp2=2*ifp1 theta=isign*6.28318530717959d0/(ifp2/ip1) wpr=-2.d0*dsin(0.5d0*theta)**2 wpi=dsin(theta) wr=1.d0 wi=0.d0 do 17 i3=1,ifp1,ip1 do 16 i1=i3,i3+ip1-2,2 do 15 i2=i1,ip3,ifp2 k1=i2 k2=k1+ifp1 tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1) tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2) data(k2)=data(k1)-tempr data(k2+1)=data(k1+1)-tempi data(k1)=data(k1)+tempr data(k1+1)=data(k1+1)+tempi 15 continue 16 continue wtemp=wr wr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi 17 continue ifp1=ifp2 go to 2 endif nprev=n*nprev 18 continue return end