c................................................................... c subroutine set_gravity3d c c................................................................... include 'von_3defs.h' save c......Prepare fixed terms for calls to FFT 3--D potential solver c......This routine only called once during run. k=n2 dx=dx3b(1) C..... compute fixed constants for the Poisson solver do i=-k,k-1 do j=-k,k-1 do l=-k,k-1 Glm(i,j,l)=-1./sqrt + (eps**2+(dx*real(i))**2+(dx*real(j))**2+(dx*real(l))**2) end do end do end do return end c................................................................... c subroutine gravity3d c c................................................................... c......Use 3D Fourier Convolution theorem to solve for gravitational c......potential phi(i,j,l). include 'von_3defs.h' c......define local variables real*8 Mlm(-n2:n2-1,-n2:n2-1,-n2:n2-1) real*8 fac,data1(n2*2,n2*2,n2*2),data2(n2*2,n2*2,n2*2) complex*16 spec1(n2,n2*2,n2*2),speq1(n2*2,n2*2), + spec2(n2,n2*2,n2*2),speq2(n2*2,n2*2), + zpec1(n2*2*n2*2*n2),zpeq1(n2*2*n2*2), + zpec2(n2*2*n2*2*n2),zpeq2(n2*2*n2*2) equivalence (data1,spec1,zpec1), (data2,spec2,zpec2), + (speq1,zpeq1), (speq2,zpeq2) save k=n2 n=n2*2 do i=-k,k-1 do j=-k,k-1 do l=-k,k-1 if((i.lt.0).or.(j.lt.0).or.(l.lt.0))then Mlm(i,j,l)=0. else Mlm(i,j,l)=d(i+1,j+1,l+1)*dx**3 endif end do end do end do do i=1,n do j=1,n do l=1,n if (i.le.n/2.and.j.le.n/2.and.l.le.n/2) then data1(i,j,l)=Glm(i-1,j-1,l-1) data2(i,j,l)=Mlm(i-1,j-1,l-1) elseif (i.gt.n/2.and.j.le.n/2.and.l.le.n/2) then data1(i,j,l)=Glm(i-1-n,j-1,l-1) data2(i,j,l)=Mlm(i-1-n,j-1,l-1) elseif (i.le.n/2.and.j.le.n/2.and.l.gt.n/2) then data1(i,j,l)=Glm(i-1,j-1,l-1-n) data2(i,j,l)=Mlm(i-1,j-1,l-1-n) elseif (i.gt.n/2.and.j.le.n/2.and.l.gt.n/2) then data1(i,j,l)=Glm(i-1-n,j-1,l-1-n) data2(i,j,l)=Mlm(i-1-n,j-1,l-1-n) elseif (i.le.n/2.and.j.gt.n/2.and.l.gt.n/2) then data1(i,j,l)=Glm(i-1,j-1-n,l-1-n) data2(i,j,l)=Mlm(i-1,j-1-n,l-1-n) elseif (i.gt.n/2.and.j.gt.n/2.and.l.gt.n/2) then data1(i,j,l)=Glm(i-1-n,j-1-n,l-1-n) data2(i,j,l)=Mlm(i-1-n,j-1-n,l-1-n) elseif (i.le.n/2.and.j.gt.n/2.and.l.le.n/2) then data1(i,j,l)=Glm(i-1,j-1-n,l-1) data2(i,j,l)=Mlm(i-1,j-1-n,l-1) elseif (i.gt.n/2.and.j.gt.n/2.and.l.le.n/2) then data1(i,j,l)=Glm(i-1-n,j-1-n,l-1) data2(i,j,l)=Mlm(i-1-n,j-1-n,l-1) else write(*,*) 'Error in potential grid' stop end if end do end do end do c.. call the FFT routines and perform convolution call rlft3(data1,speq1,n,n,n,1) call rlft3(data2,speq2,n,n,n,1) fac=2./(N*N*N) do j=1,n*n*n/2 zpec1(j)=fac*zpec1(j)*zpec2(j) end do do j=1,n*n zpeq1(j)=fac*zpeq1(j)*zpeq2(j) end do call rlft3(data1,speq1,n,n,n,-1) do i=1,k do j=1,k do l=1,k phi(i,j,l)=data1(i,j,l) end do end do end do return end c............................................................................ subroutine rlft3(data,speq,nn1,nn2,nn3,isign) c............................................................................ implicit real*8(a-h,o-z), integer*4(i-n) integer*4 isign,nn1,nn2,nn3 complex*16 data(nn1/2,nn2,nn3),speq(nn2,nn3) c.. uses fourn c.. From Numerical Recipes, 2nd edition c---------------------------------------------------------------------------- c.. Given a two- or three- dimensional real array data whose dimensions c.. are nn1, nn2, nn3 (where nn3 is 1 for the case of a two-dimensional c.. array), this routine returns (for isign=1) the complex fast fourier c.. transform as two complex arrays: On output, data contains the zero c.. and positive frequency values for the first frequency component, c.. while speq contains the Nyquist critical frequency values of the c.. first frequency component. Second (and third) frequency components c.. are store for zero, positive, and negative frequencies, in standard c.. wrap-around order. For isign=-1, the inverse transform (times nn1*nn2 c.. *nn3/2 as a constant multiplicative factor) is performed, with c.. output data (viewed as a real array) deriving from input data (viewed c.. as complex) and speq. The dimensions nn1, nn2, nn3 must always be c.. integer powers of 2. c----------------------------------------------------------------------------- integer*4 i1,i2,i3,j1,j2,j3,nn(3) double precision theta,wi,wpi,wpr,wr,wtemp complex*16 c1,c2,h1,h2,w c1=cmplx(0.5,0.0) c2=cmplx(0.0,-0.5*isign) theta=6.28318530717959d0/dble(isign*nn1) wpr=-2.0d0*sin(0.5d0*theta)**2 wpi=sin(theta) nn(1)=nn1/2 nn(2)=nn2 nn(3)=nn3 if(isign.eq.1)then call fourn(data,nn,3,isign) do i3=1,nn3 do i2=1,nn2 speq(i2,i3)=data(1,i2,i3) end do end do end if do i3=1,nn3 j3=1 if (i3.ne.1) j3=nn3-i3+2 wr=1.0d0 wi=0.0d0 do i1=1,nn1/4+1 j1=nn1/2-i1+2 do i2=1,nn2 j2=1 if (i2.ne.1) j2=nn2-i2+2 if(i1.eq.1) then h1=c1*(data(1,i2,i3)+conjg(speq(j2,j3))) h2=c2*(data(1,i2,i3)-conjg(speq(j2,j3))) data(1,i2,i3)=h1+h2 speq(j2,j3)=conjg(h1-h2) else h1=c1*(data(i1,i2,i3)+conjg(data(j1,j2,j3))) h2=c2*(data(i1,i2,i3)-conjg(data(j1,j2,j3))) data(i1,i2,i3)=h1+w*h2 data(j1,j2,j3)=conjg(h1-w*h2) end if end do wtemp=wr wr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi w=cmplx(sngl(wr),sngl(wi)) end do end do if(isign.eq.-1)then call fourn(data,nn,3,isign) endif return end c.................................................................. c subroutine fourn(data,nn,ndim,isign) c c................................................................. implicit real*8(a-h,o-z), integer*4(i-n) integer*4 isign,ndim,nn(ndim) real*8 data(*) integer*4 i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2, + ip3,k1,k2,n,nprev,nrem,ntot real*8 tempi,tempr double precision theta,wi,wpi,wpr,wr,wtemp ntot=1 do idim=1,ndim ntot=ntot*nn(idim) end do nprev=1 do idim=1,ndim n=nn(idim) nrem=ntot/(n*nprev) ip1=2*nprev ip2=ip1*n ip3=ip2*nrem i2rev=1 do i2=1,ip2,ip1 if(i2.lt.i2rev) then do i1=i2,i2+ip1-2,2 do 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 end do end do end if ibit=ip2/2 1 if((ibit.ge.ip1).and.(i2rev.gt.ibit))then i2rev=i2rev-ibit ibit=ibit/2 goto 1 endif i2rev=i2rev+ibit end do ifp1=ip1 2 if(ifp1.lt.ip2)then ifp2=2*ifp1 theta=isign*6.28318530717959d0/(ifp2/ip1) wpr=-2.d0*sin(0.5d0*theta)**2 wpi=sin(theta) wr=1.d0 wi=0.d0 do i3=1,ifp1,ip1 do i1=i3,i3+ip1-2,2 do 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 end do end do wtemp=wr wr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi end do ifp1=ifp2 goto 2 endif nprev=n*nprev end do return end c..................End Routine........................................