c 270107 c try to speed the search of nearby pairs c with 2D partitioning c-------------------------------------------------------------- c 220605 c if(encmin.eq.0) then nerby gravity +impacts ignored c better to do something with impacts then! c fft-array definations moved to ssbox.inc 31/05/2003 c contains individual eromax 6/4/2003 c gravinew567: c gmet=5 old c gmet=6 da/dt c gmet=7 gmet=6 + impact list separately c gmet=8 only far-gravity c uses also glim: delta > glim*ssum c gravinew2 -> evaluate the mutual gravity of c potentially colliding particles at dimpacts c -> excluded here c gravinew -> evaluate also dA/dt for nearby gravity c remove all except gmet=5 c include non-colliding tracer-particles c in this case use encmin=eromax + non-zero softening c ~hsalo/POTSDAM/SSBOX/OCT2000/ssgravi_d2.f c correct bug in defination of dyp with gmet=5 c add individual version-number to gidl and ggidl-files c remove velocity,xgrid from gidl output c add ss-output to gidl output c update ggidl-output (gmet=5) c _d2 -> include gravity from replicas (to gmet=5) c determined by mxgravi,mygravi c fz on particles outside grid: linear interpolation c 17.10.2000 c _d3 -> include possible gravity softening c esoft=max(encsoft,ssum) c*********************************************************** SUBROUTINE GRAVIGRID_mc(ival) c*********************************************************** c calculate separately grid-gravity and c that due to nearby particles: HS 11/08/98 c every step: c 1) calculate direct forces c from particles with DIST < ENCMIN c 2) add contribution of gravity c from particles with ENCMIN < DIST < ENCLIM c using FFT on xyz-grid: nxg*nyg*nyz cells c inclined grid c 3) add vertical gravity due to infinite sheet c ival determines test-output c ival=0 - initialization c ival=1 - force profiles c ival=2 - force profiles + individual forces c IVAL=2 -> output forces to DGIDL-file, UNIT 30 c R,AGRAV_near c V,AGRAV_tot c xgrid,ygrid,zgrid, fxgrid,fygrid,fzgrid c -> calculation of gravity can be checked c IVAL=1,2 -> output grid-gravity cuts to DGGIDL, UNIT 31 c for 20 locations along y,z=0 c for 10 location along x,z=0 c for 60 locations along x,y=0 c cpugrav1 = CPU spent in sorting c cpugrav2 = CPU spent in tabulating distant gravity c cpugrav3 = CPU spent in calculation of nearby gravity c cpugrav4 = CPU spent in interpolation of distant gravity C*********************************************************** INCLUDE 'ssbox.inc' c for near gravity real xvar(nstar2),yvar(nstar2) dimension indx(nstar2) c for tests real xtest(160),ytest(160),ztest(160) real fxapu31(160),fyapu31(160),fzapu31(160) dimension agrav_sheet(nstar) c for FFT c 3D: maximum true grid-size 512*256*8 c -> 176 mega c parameter (nxmax=512,nymax=256,nzmax=8) c parameter (nxmax=128,nymax=128,nzmax=32) c parameter (nzmax2=2*nzmax,ndatamax=4*nxmax*nymax*nzmax) c parameter declaration moved to include-file 310503 real rhogrid3(nxmax,nymax,nzmax) real fxgrid3(nxmax,nymax,nzmax) real fygrid3(nxmax,nymax,nzmax) real fzgrid3(nxmax,nymax,nzmax) real frrho3(nxmax,nymax,nzmax2) real firho3(nxmax,nymax,nzmax2) real distx3(nxmax,nymax,nzmax2) real disty3(nxmax,nymax,nzmax2) real distz3(nxmax,nymax,nzmax2) real frdist3(nxmax,nymax,nzmax2) real fidist3(nxmax,nymax,nzmax2) real data(ndatamax) dimension nnn(3) c 2d-partitioning in ssbox dimension pottab_mc(nstar2) dimension n_pottab_mc(nstar2) dimension m_pottab_mc(nstar2) VVAK=(KAPPA2-4.D0*OMEGA2)/OMEGA*XL CALL CPU(TIME1) encmin2=encmin**2 tgrid=0.d0 svak=bvak/2.d0/xl c interval for f_checks (in steps) i200=20 c*********************************************** c*********************************************** c TABULATE DISTANT GRAVITY c forces from xp,yp,zp density with FFT c*********************************************** c*********************************************** c encmin = limit for nearby gravity c enclim = limit for near+far gravity c encmin < enclim -> use fft for the zone from encmin to enclim c if(encmin.lt.enclim) then c CALL CPU(TTIME1) svak=bvak/2.d0/xl n128x=nxg n128y=nyg n128z=2*nzg n64x=n128x/2 n64y=n128y/2 n64z=n128z/2 dxg=2.d0*xl/n128x dyg=2.d0*yl/n128y c----------------------------------------------- c vertical extent defined by fzrange c fzrange > 0 use zmax=fzrange*stdev(zz) c fzrange < 0 use zmax=abs(fzrange) apu=0.d0 apu2=0.d0 do 400 i=1,lkm apu2=apu2+r(i,3)**2 apu=apu+r(i,3) 400 continue apu=apu/lkm apu2=sqrt(apu2/lkm) apu_zdisp=apu2 apu_zmean=apu if(fzrange.gt.0) then zmax=fzrange*apu2 endif if(fzrange.lt.0) zmax=-fzrange zmax_stamp=zmax zmin=-zmax dzg=(zmax-zmin)/nzg c use this variable for storing zmax (output to dgidl and dggidl) tgrid=zmax c-------------------------------------------------------------- c moved here -> above test-quantities are always calculated if(encmin.lt.enclim) then CALL CPU(TTIME1) c----------------------------------------------- c clear density arrays if(time.lt.period/20.) then CALL CPU(TTTIME1) endif do 401 i=1,n128x do 402 j=1,n128y do 403 k=1,n128z rhogrid3(i,j,k)=0.d0 distx3(i,j,k)=0. disty3(i,j,k)=0. distz3(i,j,k)=0. 403 continue 402 continue 401 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-GRID',time/period, + tttime2-tttime1 write(6,*) rhogrid3(1,1,nzg/2),rhogrid3(1,1,nzg/2+1) CALL CPU(TTTIME1) endif c------------------------------------------------------- c calculate density rhogrid3 c cic=0 ngp c cic=2 cic c ngp-assignment if(cic.lt.1) then do 410 i=1,lkm xp=r(i,1) yp=r(i,2)-svak*xp zp=r(i,3) if(yp.gt.yl) yp=yp-2.d0*yl if(yp.lt.-yl) yp=yp+2.d0*yl if(yp.gt.yl) yp=yp-2.d0*yl if(yp.lt.-yl) yp=yp+2.d0*yl if(xp.gt.xl) xp=xp-2.d0*xl if(xp.lt.-xl) xp=xp+2.d0*xl ixp=(xp+xl)/dxg+1 iyp=(yp+yl)/dyg+1 izp=(zp-zmin)/dzg+1 if(izp.ge.1.and.izp.le.n64z) + rhogrid3(ixp,iyp,izp)= + rhogrid3(ixp,iyp,izp)+gg*xmass(i) 410 continue endif c cic-assignment if(cic.gt.1) then do 411 i=1,lkm xp=r(i,1) yp=r(i,2)-svak*xp zp=r(i,3) if(yp.gt.yl) yp=yp-2.d0*yl if(yp.lt.-yl) yp=yp+2.d0*yl if(yp.gt.yl) yp=yp-2.d0*yl if(yp.lt.-yl) yp=yp+2.d0*yl if(xp.gt.xl) xp=xp-2.d0*xl if(xp.lt.-xl) xp=xp+2.d0*xl px=(xp+xl)/2./xl*nxg+n128x-0.5 py=(yp+yl)/2./yl*nyg+n128y-0.5 pz=(zp-zmin)/(zmax-zmin)*nzg+n128z-0.5 iix=int(px)-n128x iiy=int(py)-n128y iiz=int(pz)-n128z px=px-n128x py=py-n128y pz=pz-n128z dx=px-iix dy=py-iiy dz=pz-iiz tx=1.-dx ty=1.-dy tz=1.-dz ix1=iix+1 ix2=iix+2 iy1=iiy+1 iy2=iiy+2 iz1=iiz+1 iz2=iiz+2 if(ix1.gt.n128x) ix1=ix1-n128x if(ix2.gt.n128x) ix2=ix2-n128x if(iy1.gt.n128y) iy1=iy1-n128y if(iy2.gt.n128y) iy2=iy2-n128y if(ix1.lt.1) ix1=ix1+n128x if(ix2.lt.1) ix2=ix2+n128x if(iy1.lt.1) iy1=iy1+n128y if(iy2.lt.1) iy2=iy2+n128y if(ix2.gt.n128x.or.ix1.lt.1) then write(6,*) 'x-outgrid1',time/period,i, + r(i,1),ix1,ix2,n128x goto 199 endif if(iy2.gt.n128y.or.iy1.lt.1) then write(6,*) 'y-outgrid1',time/period,i, + r(i,2),iy1,iy2,n128y goto 199 endif c mass falling outside the grid if(iz2.gt.nzg) tz=0. if(iz2.lt.1) tz=0. if(iz1.gt.nzg) dz=0. if(iz1.lt.1) dz=0. c these just for convenience if(iz1.lt.1) iz1=1 if(iz2.gt.nzg) iz2=nzg if(iz2.lt.1) iz2=1 if(iz1.gt.nzg) iz1=nzg rhogrid3(ix2,iy2,iz2)= + rhogrid3(ix2,iy2,iz2)+gg*xmass(i)*dx*dy*dz rhogrid3(ix2,iy2,iz1)= + rhogrid3(ix2,iy2,iz1)+gg*xmass(i)*dx*dy*tz rhogrid3(ix2,iy1,iz2)= + rhogrid3(ix2,iy1,iz2)+gg*xmass(i)*dx*ty*dz rhogrid3(ix2,iy1,iz1)= + rhogrid3(ix2,iy1,iz1)+gg*xmass(i)*dx*ty*tz rhogrid3(ix1,iy2,iz2)= + rhogrid3(ix1,iy2,iz2)+gg*xmass(i)*tx*dy*dz rhogrid3(ix1,iy2,iz1)= + rhogrid3(ix1,iy2,iz1)+gg*xmass(i)*tx*dy*tz rhogrid3(ix1,iy1,iz2)= + rhogrid3(ix1,iy1,iz2)+gg*xmass(i)*tx*ty*dz rhogrid3(ix1,iy1,iz1)= + rhogrid3(ix1,iy1,iz1)+gg*xmass(i)*tx*ty*tz 199 continue 411 continue endif if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-RHO',tttime2-tttime1,cic write(6,*) rhogrid3(1,1,nzg/2),rhogrid3(1,1,nzg/2+1) CALL CPU(TTTIME1) endif c------------------------------------------------------- c construct distx3,disty3,distz3 c this makes sure that forces are symmetric if only one box is used i_yla=n64x-1 j_yla=n64y-1 if(mxgravi.gt.0) i_yla=n64x if(mygravi.gt.0) j_yla=n64y do 420 i=-n64x+1,i_yla do 421 j=-n64y+1,j_yla dx=i*dxg dy=j*dyg+svak*dx apu2=dx*dx+dy*dy do 422 k=-n64z+1,n64z-1 dz=k*dzg if(apu2.gt.encmin2.and.apu2.le.enclim2) then apu=dx*dx+dy*dy+dz*dz ii=i+1 if(ii.lt.1) ii=ii+n128x jj=j+1 if(jj.lt.1) jj=jj+n128y kk=k+1 if(kk.lt.1) kk=kk+n128z apu=apu+encsoft2 apu3=1./sqrt(apu)/apu apux=dx*apu3 apuy=dy*apu3 apuz=dz*apu3 distx3(ii,jj,kk)=apux disty3(ii,jj,kk)=apuy distz3(ii,jj,kk)=apuz endif c _d2: if(mxgravi.gt.0.or.mygravi.gt.0) then do 1420 imx=-mxgravi,mxgravi do 1421 imy=-mygravi,mygravi if(imx.ne.0.or.imy.ne.0) then dx=i*dxg+2.*xl*imx dy=j*dyg+svak*dx+2.*yl*imy apu=dx*dx+dy*dy+dz*dz apu2=dx*dx+dy*dy if(apu2.gt.encmin2. + and.apu.le.enclim2) then apu=apu+encsoft2 apu3=1./sqrt(apu)/apu apux=dx*apu3 apuy=dy*apu3 apuz=dz*apu3 distx3(ii,jj,kk)= + distx3(ii,jj,kk)+apux disty3(ii,jj,kk)= + disty3(ii,jj,kk)+apuy distz3(ii,jj,kk)= + distz3(ii,jj,kk)+apuz endif endif 1421 continue 1420 continue endif 422 continue 421 continue 420 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-DIST',tttime2-tttime1 CALL CPU(TTTIME1) endif c------------------------------------------------------- c FFT of rhogrid3 ndim=3 nnn(1)=n128x nnn(2)=n128y nnn(3)=n128z lask=0 do 430 k=1,n128z do 431 j=1,n128y do 432 i=1,n128x lask=lask+1 data(lask)=rhogrid3(i,j,k) lask=lask+1 data(lask)=0. 432 continue 431 continue 430 continue ISIGN=+1 CALL FOURN2007(DATA,NNN,NDIM,ISIGN) lask=0 do 433 k=1,n128z do 434 j=1,n128y do 435 i=1,n128x lask=lask+1 frrho3(i,j,k)=data(lask) lask=lask+1 firho3(i,j,k)=data(lask) 435 continue 434 continue 433 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-FFT(RHO)',tttime2-tttime1 CALL CPU(TTTIME1) endif c---------------------------------- c only vertical forces? c 09/10/2000 <- 24/1/2000: not done if encz=1 if(encz.eq.1) then do 793 k=1,n64z do 794 j=1,n128y do 795 i=1,n128x fxgrid3(i,j,k)=0. fygrid3(i,j,k)=0. 795 continue 794 continue 793 continue endif if(encz.ne.1) then c------------------------------------------------------- c CALCULATION OF FXGRID3 c------------------------------------------------------- c FFT of distx3 lask=0 do 440 k=1,n128z do 441 j=1,n128y do 442 i=1,n128x lask=lask+1 data(lask)=distx3(i,j,k) lask=lask+1 data(lask)=0. 442 continue 441 continue 440 continue ISIGN=+1 CALL FOURN2007(DATA,NNN,NDIM,ISIGN) lask=0 do 443 k=1,n128z do 444 j=1,n128y do 445 i=1,n128x lask=lask+1 frdist3(i,j,k)=data(lask) lask=lask+1 fidist3(i,j,k)=data(lask) 445 continue 444 continue 443 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-FFT(distx3)',tttime2-tttime1 CALL CPU(TTTIME1) endif c------------------------------------------------------- c multiply fft(rhogrid3) by fft(dist3) c and insert into data lask=0 do 450 k=1,n128z do 451 j=1,n128y do 452 i=1,n128x lask=lask+1 data(lask)=frdist3(i,j,k)*frrho3(i,j,k)- + fidist3(i,j,k)*firho3(i,j,k) lask=lask+1 data(lask)=frdist3(i,j,k)*firho3(i,j,k)+ + fidist3(i,j,k)*frrho3(i,j,k) 452 continue 451 continue 450 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-MULTIPLY',tttime2-tttime1 CALL CPU(TTTIME1) endif c------------------------------------------------------- c obtain FX by inverse transform c just the real part isign=-1 CALL FOURN2007(DATA,NNN,NDIM,ISIGN) lask=0 do 453 k=1,n64z do 454 j=1,n128y do 455 i=1,n128x lask=lask+1 fxgrid3(i,j,k)=-data(lask)/n128x/n128y/n128z lask=lask+1 455 continue 454 continue 453 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-FXGRID',tttime2-tttime1 CALL CPU(TTTIME1) endif c------------------------------------------------------- c CALCULATION OF FYGRID3 c------------------------------------------------------- c FFT of disty3 lask=0 do 460 k=1,n128z do 461 j=1,n128y do 462 i=1,n128x lask=lask+1 data(lask)=disty3(i,j,k) lask=lask+1 data(lask)=0. 462 continue 461 continue 460 continue ISIGN=+1 CALL FOURN2007(DATA,NNN,NDIM,ISIGN) lask=0 do 463 k=1,n128z do 464 j=1,n128y do 465 i=1,n128x lask=lask+1 frdist3(i,j,k)=data(lask) lask=lask+1 fidist3(i,j,k)=data(lask) 465 continue 464 continue 463 continue c------------------------------------------------------- c multiply fft(rhogrid3) by fft(dist3) c and insert into data lask=0 do 470 k=1,n128z do 471 j=1,n128y do 472 i=1,n128x lask=lask+1 data(lask)=frdist3(i,j,k)*frrho3(i,j,k)- + fidist3(i,j,k)*firho3(i,j,k) lask=lask+1 data(lask)=frdist3(i,j,k)*firho3(i,j,k)+ + fidist3(i,j,k)*frrho3(i,j,k) 472 continue 471 continue 470 continue c------------------------------------------------------- c obtain FY by inverse transform c just the real part isign=-1 CALL FOURN2007(DATA,NNN,NDIM,ISIGN) lask=0 do 473 k=1,n64z do 474 j=1,n128y do 475 i=1,n128x lask=lask+1 fygrid3(i,j,k)=-data(lask)/n128x/n128y/n128z lask=lask+1 475 continue 474 continue 473 continue endif c encz ne 1? c------------------------------------------------------- c CALCULATION OF FZGRID3 c------------------------------------------------------- c FFT of distz3 lask=0 do 480 k=1,n128z do 481 j=1,n128y do 482 i=1,n128x lask=lask+1 data(lask)=distz3(i,j,k) lask=lask+1 data(lask)=0. 482 continue 481 continue 480 continue ISIGN=+1 CALL FOURN2007(DATA,NNN,NDIM,ISIGN) lask=0 do 483 k=1,n128z do 484 j=1,n128y do 485 i=1,n128x lask=lask+1 frdist3(i,j,k)=data(lask) lask=lask+1 fidist3(i,j,k)=data(lask) 485 continue 484 continue 483 continue c------------------------------------------------------- c multiply fft(rhogrid3) by fft(dist3) c and insert into data lask=0 do 490 k=1,n128z do 491 j=1,n128y do 492 i=1,n128x lask=lask+1 data(lask)=frdist3(i,j,k)*frrho3(i,j,k)- + fidist3(i,j,k)*firho3(i,j,k) lask=lask+1 data(lask)=frdist3(i,j,k)*firho3(i,j,k)+ + fidist3(i,j,k)*frrho3(i,j,k) 492 continue 491 continue 490 continue c------------------------------------------------------- c obtain FZ by inverse transform c just the real part isign=-1 CALL FOURN2007(DATA,NNN,NDIM,ISIGN) lask=0 do 493 k=1,n64z do 494 j=1,n128y do 495 i=1,n128x lask=lask+1 fzgrid3(i,j,k)=-data(lask)/n128x/n128y/n128z lask=lask+1 495 continue 494 continue 493 continue CALL CPU(TTIME2) cpugrav2=cpugrav2+ttime2-ttime1 c encmin = limit for nearby gravity c enclim = limit for near+far gravity c encmin < enclim -> use fft for the zone from encmin to enclim c encmin exclude gravity of particles in impact list c**************************************************** c**************************************************** raja1_sum=0. raja1_lkm=0. raja1_max=0. ix=0 if(encmin.gt.1d-6) then c if(simustep.eq.1) then c write(97,*) 'ii,jj,nntemp,mmtemp,ssum,ero,dero,raja1_use' c write(97,*) 'dx,dy,dz,dvx,dvy,dvz' c endif c--------------------------------------------------- CALL CPU(TTIME1) call tab_init CALL CPU(TTIME2) cpugrav1=cpugrav1+ttime2-ttime1 c calculate gravity + find nearby pairs CALL CPU(TTIME1) ix=0 rajatmp=encmin c target particles c npotsum1 = checked for distance c npotsum2 = accepted to contribute npotsum1=0 npotsum2=0 npotsum3=0 npotsum4=0 do 20 ii=1,lkm-1 c choose particles jj potentially contributing c to gravity at particle ii npot=0 xt=r(ii,1) yt=r(ii,2) n_mc1=0 n_mc2=0 m_mc1=0 m_mc2=0 if(xt.gt.(xl-encmin)) n_mc1=-1 if(xt.lt.(-xl+encmin)) n_mc2=1 if(yt.gt.(yl-encmin)) m_mc1=-1 if(yt.lt.(-yl+encmin)) m_mc2=1 do 9008 n_mc=n_mc1,n_mc2 do 9009 m_mc=m_mc1,m_mc2 do 9010 imc=1,mx_mc*my_mc if(numtab_mc(imc).le.0) goto 9010 xt=r(ii,1)+n_mc*avak yt=r(ii,2)+n_mc*bvak+m_mc*cvak dx=cxtab_mc(imc)-xt dy=cytab_mc(imc)-yt d2=dx**2+dy**2 if(d2.gt.dtab2_mc(imc)) goto 9010 do 9020 jmc=1,numtab_mc(imc) npot=npot+1 k=mintab_mc(imc)+(jmc-1) pottab_mc(npot)=indtab_mc(k) n_pottab_mc(npot)=n_mc m_pottab_mc(npot)=m_mc 9020 continue 9010 continue 9009 continue 9008 continue npotsum1=npotsum1+npot do 30 j=1,npot jj=pottab_mc(j) if(jj.le.ii) goto 30 nntemp=n_pottab_mc(j) mmtemp=m_pottab_mc(j) dx=r(ii,1)-r(jj,1)+nntemp*avak dy=r(ii,2)-r(jj,2)+nntemp*bvak+mmtemp*cvak dz=r(ii,3)-r(jj,3) ddist2=dx*dx+dy*dy ddist=ddist2+dz*dz c these have to be here ssum=s(ii)+s(jj) dist=sqrt(ddist) ero=dist-ssum c******************************************* c nearby gravitational forces if(ddist2.le.encmin2) then npotsum2=npotsum2+1 dist1=dist if(dist1.le.ssum*glim) dist1=ssum*glim if(encsoft.gt.0) dist1=sqrt(dist1**2+encsoft2) dvx=rp(ii,1)-rp(jj,1) dvy=rp(ii,2)-rp(jj,2)+nntemp*vvak dvz=rp(ii,3)-rp(jj,3) c------------------------------------------------------------------ c gravinew567 c gmet=8: no nearby gravity c this is calculated in dimpacts if(gmet.eq.8) goto 9876 c-------------------------------- c individual eromax c used if particles have experienced enough impacts c (together more than coll_eropart impacts since last update) c (coll_eropart negative -> do not use) c OTHERWISE MAKE THE SAME WITH EROPART_OLD,EROCOLL_OLD C LAST REMEDY IS TO USE RAJA1 raja1_use=raja1 raja1_use_ori=0. if(coll_eropart.ge.0) then C normal case c----------------------------------------------------- if(collpart(ii)+collpart(jj).ge. + coll_eropart) then raja1_use=eropart(ii)*SAFE_EROPART raja1_use2=eropart(JJ)*SAFE_EROPART if(RAJA1_USE2.gt.raja1_use) + raja1_use=RAJA1_USE2 raja1_use_ori=raja1_use if(raja1_use.lt.raja1min) + raja1_use=raja1min goto 1234 endif C from previous step c this helps to avoid spurious spikes in ix? c----------------------------------------------------- if(collpart_old(ii)+collpart_old(jj).ge. + coll_eropart) then raja1_use=eropart_old(ii)*SAFE_EROPART raja1_use2=eropart_old(JJ)*SAFE_EROPART if(RAJA1_USE2.gt.raja1_use) + raja1_use=RAJA1_USE2 raja1_use_ori=raja1_use if(raja1_use.lt.raja1min) + raja1_use=raja1min endif endif 1234 continue raja1_sum=raja1_sum+raja1_use raja1_lkm=raja1_lkm+1. if(raja1_use_ori.gt.raja1_max) + raja1_max=raja1_use_ori c gmet=7: exclude pairs which are inserted to collision list if(gmet.eq.7) then if(ero.lt.raja1_use) then dero=EROVAK1*sqrt(dvx**2+dvy**2+dvz**2)*dt+ + EROVAK2*raja1_use if(ero.lt.dero) then c write(6,*) '*gmet7:excl',time,ii,jj,dist goto 9876 endif endif endif c------------------------------------------------------------------ FAPU=GG/dist1/dist1/dist1 FAPU1=-FAPU*XMASS(JJ) FAPU2= FAPU*XMASS(II) agrav(ii,1)=agrav(ii,1)+fapu1*dx agrav(ii,2)=agrav(ii,2)+fapu1*dy agrav(ii,3)=agrav(ii,3)+fapu1*dz agrav(jj,1)=agrav(jj,1)+fapu2*dx agrav(jj,2)=agrav(jj,2)+fapu2*dy agrav(jj,3)=agrav(jj,3)+fapu2*dz c energy due to over-border gravity c s=-shear c but nntemp=-1 -> negative sign if(nntemp.ne.0) then enes_grav=enes_grav-xmass(ii)*2.*xl*shear* + fapu1*dy*dt e_obg_size(is(ii))=e_obg_size(is(ii))- + xmass(ii)*2.*xl*shear*fapu1*dy*dt endif c gmet=5: do not calculate dA/dt if(gmet.eq.5) goto 9876 c------------------------------------------------------------------ c dA/dt only for gmet=6 and gmet=7 cv=dvx*dx+dvy*dy+dvz*dz apux=dvx-3.d0*cv*dx/dist1/dist1 apuy=dvy-3.d0*cv*dy/dist1/dist1 apuz=dvz-3.d0*cv*dz/dist1/dist1 aagrav(ii,1)=aagrav(ii,1)+fapu1*apux aagrav(ii,2)=aagrav(ii,2)+fapu1*apuy aagrav(ii,3)=aagrav(ii,3)+fapu1*apuz aagrav(jj,1)=aagrav(jj,1)+fapu2*apux aagrav(jj,2)=aagrav(jj,2)+fapu2*apuy aagrav(jj,3)=aagrav(jj,3)+fapu2*apuz c no gravity or dA/dt for this pair 9876 continue c--------------------------------------------------- c dist use fft for the zone from encmin to enclim if(encmin.lt.enclim) then CALL CPU(TTIME1) c**************************************************** c first normal particles c**************************************************** do 60 i=1,lkm fxapu=0. fyapu=0. fzapu=0. x0p=r(i,1) y0p=r(i,2)-svak*x0p z0p=r(i,3) if(y0p.gt.yl) y0p=y0p-2.*yl if(y0p.lt.-yl) y0p=y0p+2.*yl px=(x0p+xl)/2./xl*nxg+n128x-0.5 py=(y0p+yl)/2./yl*nyg+n128y-0.5 pz=(z0p-zmin)/(zmax-zmin)*nzg+n128z-0.5 iix=int(px)-n128x iiy=int(py)-n128y iiz=int(pz)-n128z px=px-n128x py=py-n128y pz=pz-n128z dx=px-iix dy=py-iiy dz=pz-iiz tx=1.-dx ty=1.-dy tz=1.-dz ix1=iix+1 ix2=iix+2 iy1=iiy+1 iy2=iiy+2 iz1=iiz+1 iz2=iiz+2 if(ix1.gt.n128x) ix1=ix1-n128x if(ix2.gt.n128x) ix2=ix2-n128x if(iy1.gt.n128y) iy1=iy1-n128y if(iy2.gt.n128y) iy2=iy2-n128y if(ix1.lt.1) ix1=ix1+n128x if(ix2.lt.1) ix2=ix2+n128x if(iy1.lt.1) iy1=iy1+n128y if(iy2.lt.1) iy2=iy2+n128y if(ix2.gt.n128x.or.ix1.lt.1) then write(6,*) 'x-outgrid',time/period,i, + r(i,1),ix1,ix2,n128x goto 99 endif if(iy2.gt.n128y.or.iy1.lt.1) then write(6,*) 'y-outgrid',time/period,i, + r(i,2),iy1,iy2,n128y goto 99 endif c mass falling outside the grid c 9.5.2002 if(iz2.gt.nzg) tz=0. if(iz2.lt.1) tz=0. if(iz1.gt.nzg) dz=0. if(iz1.lt.1) dz=0. if(iz1.lt.1) iz1=1 if(iz2.lt.1) iz2=1 if(iz1.gt.nzg) iz1=nzg if(iz2.gt.nzg) iz2=nzg fxapu= + dx*dy*dz*fxgrid3(ix2,iy2,iz2)+ + dx*dy*tz*fxgrid3(ix2,iy2,iz1)+ + dx*ty*dz*fxgrid3(ix2,iy1,iz2)+ + dx*ty*tz*fxgrid3(ix2,iy1,iz1)+ + tx*dy*dz*fxgrid3(ix1,iy2,iz2)+ + tx*dy*tz*fxgrid3(ix1,iy2,iz1)+ + tx*ty*dz*fxgrid3(ix1,iy1,iz2)+ + tx*ty*tz*fxgrid3(ix1,iy1,iz1) fyapu= + dx*dy*dz*fygrid3(ix2,iy2,iz2)+ + dx*dy*tz*fygrid3(ix2,iy2,iz1)+ + dx*ty*dz*fygrid3(ix2,iy1,iz2)+ + dx*ty*tz*fygrid3(ix2,iy1,iz1)+ + tx*dy*dz*fygrid3(ix1,iy2,iz2)+ + tx*dy*tz*fygrid3(ix1,iy2,iz1)+ + tx*ty*dz*fygrid3(ix1,iy1,iz2)+ + tx*ty*tz*fygrid3(ix1,iy1,iz1) fzapu= + dx*dy*dz*fzgrid3(ix2,iy2,iz2)+ + dx*dy*tz*fzgrid3(ix2,iy2,iz1)+ + dx*ty*dz*fzgrid3(ix2,iy1,iz2)+ + dx*ty*tz*fzgrid3(ix2,iy1,iz1)+ + tx*dy*dz*fzgrid3(ix1,iy2,iz2)+ + tx*dy*tz*fzgrid3(ix1,iy2,iz1)+ + tx*ty*dz*fzgrid3(ix1,iy1,iz2)+ + tx*ty*tz*fzgrid3(ix1,iy1,iz1) c set grid-force to zero if outside grid c add later that due to a constant-density sheet if(abs(z0p).gt.zmax) fzapu=0. 99 continue agrav(i,1)=agrav(i,1)+fxapu agrav(i,2)=agrav(i,2)+fyapu agrav(i,3)=agrav(i,3)+fzapu if(encz.eq.1) agrav(i,1)=0. if(encz.eq.1) agrav(i,2)=0. 60 continue c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c add checks: sum of vertical gravity should be zero c contribution from near+far particles: fz_check2 if(mod(simustep,i200).eq.0) then fx_check2=0. fy_check2=0. fz_check2=0. axmin=1d8 aymin=1d8 azmin=1d8 axmax=-1d8 aymax=-1d8 azmax=-1d8 n_zover=0 do 102 i=1,lkm fx_check2=fx_check2+xmass(i)*agrav(i,1) fy_check2=fy_check2+xmass(i)*agrav(i,2) fz_check2=fz_check2+xmass(i)*agrav(i,3) if(agrav(i,1).lt.axmin) axmin=agrav(i,1) if(agrav(i,2).lt.aymin) aymin=agrav(i,2) if(agrav(i,3).lt.azmin) azmin=agrav(i,3) if(agrav(i,1).gt.axmax) axmax=agrav(i,1) if(agrav(i,2).gt.aymax) aymax=agrav(i,2) if(agrav(i,3).gt.azmax) azmax=agrav(i,3) if(abs(r(i,3)).gt.zmax) n_zover=n_zover 102 continue fx_check2=fx_check2/tmass/omega2 fy_check2=fy_check2/tmass/omega2 fz_check2=fz_check2/tmass/omega2 axmin2=axmin/omega2 aymin2=aymin/omega2 azmin2=azmin/omega2 axmax2=axmax/omega2 aymax2=aymax/omega2 azmax2=azmax/omega2 write(6,1104) 'f_check2:',orb,fx_check2,fy_check2,fz_check2 write(6,1105) 'f_check2: zmean, zdisp, zmax, nov', + apu_zmean,apu_zdisp,zmax,n_zover write(6,1104) 'f_amin2: ',orb,axmin2,aymin2,azmin2 write(6,1104) 'f_amax2: ',orb,axmax2,aymax2,azmax2 endif c**************************************************** c tracer-particles c**************************************************** if(n_tracer.gt.0) then do 160 i=1,n_tracer x0p=r_tracer(i,1) y0p=r_tracer(i,2)-svak*x0p z0p=r_tracer(i,3) if(y0p.gt.yl) y0p=y0p-2.*yl if(y0p.lt.-yl) y0p=y0p+2.*yl px=(x0p+xl)/2./xl*nxg+n128x-0.5 py=(y0p+yl)/2./yl*nyg+n128y-0.5 pz=(z0p-zmin)/(zmax-zmin)*nzg+n128z-0.5 iix=int(px)-n128x iiy=int(py)-n128y iiz=int(pz)-n128z px=px-n128x py=py-n128y pz=pz-n128z dx=px-iix dy=py-iiy dz=pz-iiz tx=1.-dx ty=1.-dy tz=1.-dz ix1=iix+1 ix2=iix+2 iy1=iiy+1 iy2=iiy+2 iz1=iiz+1 iz2=iiz+2 if(ix1.gt.n128x) ix1=ix1-n128x if(ix2.gt.n128x) ix2=ix2-n128x if(iy1.gt.n128y) iy1=iy1-n128y if(iy2.gt.n128y) iy2=iy2-n128y if(ix1.lt.1) ix1=ix1+n128x if(ix2.lt.1) ix2=ix2+n128x if(iy1.lt.1) iy1=iy1+n128y if(iy2.lt.1) iy2=iy2+n128y c mass falling outside the grid if(iz2.gt.nzg) tz=0. if(iz2.lt.1) tz=0. if(iz1.gt.nzg) dz=0. if(iz1.lt.1) dz=0. if(iz1.lt.1) iz1=1 if(iz2.lt.1) iz2=1 if(iz1.gt.nzg) iz1=nzg if(iz2.gt.nzg) iz2=nzg fxapu= + dx*dy*dz*fxgrid3(ix2,iy2,iz2)+ + dx*dy*tz*fxgrid3(ix2,iy2,iz1)+ + dx*ty*dz*fxgrid3(ix2,iy1,iz2)+ + dx*ty*tz*fxgrid3(ix2,iy1,iz1)+ + tx*dy*dz*fxgrid3(ix1,iy2,iz2)+ + tx*dy*tz*fxgrid3(ix1,iy2,iz1)+ + tx*ty*dz*fxgrid3(ix1,iy1,iz2)+ + tx*ty*tz*fxgrid3(ix1,iy1,iz1) fyapu= + dx*dy*dz*fygrid3(ix2,iy2,iz2)+ + dx*dy*tz*fygrid3(ix2,iy2,iz1)+ + dx*ty*dz*fygrid3(ix2,iy1,iz2)+ + dx*ty*tz*fygrid3(ix2,iy1,iz1)+ + tx*dy*dz*fygrid3(ix1,iy2,iz2)+ + tx*dy*tz*fygrid3(ix1,iy2,iz1)+ + tx*ty*dz*fygrid3(ix1,iy1,iz2)+ + tx*ty*tz*fygrid3(ix1,iy1,iz1) fzapu= + dx*dy*dz*fzgrid3(ix2,iy2,iz2)+ + dx*dy*tz*fzgrid3(ix2,iy2,iz1)+ + dx*ty*dz*fzgrid3(ix2,iy1,iz2)+ + dx*ty*tz*fzgrid3(ix2,iy1,iz1)+ + tx*dy*dz*fzgrid3(ix1,iy2,iz2)+ + tx*dy*tz*fzgrid3(ix1,iy2,iz1)+ + tx*ty*dz*fzgrid3(ix1,iy1,iz2)+ + tx*ty*tz*fzgrid3(ix1,iy1,iz1) c set grid-force to zero if outside grid c add later that due to constant density sheet if(abs(z0p).gt.zmax) fzapu=0. agrav_tracer(i,1)=fxapu agrav_tracer(i,2)=fyapu agrav_tracer(i,3)=fzapu if(encz.eq.1) agrav_tracer(i,1)=0. if(encz.eq.1) agrav_tracer(i,2)=0. 160 continue c n_tracer gt 0 endif c**************************************************** c**************************************************** c output slices of force to unit 31 c ival=1 or ival=2 c**************************************************** c**************************************************** if(ival.eq.1.or.ival.eq.2) then if(mod(simustep,5).eq.0) then n30=90 write(31) time,omega,tgrid,igrav write(31) nxg,nyg,n30 write(31) xl,yl,enclim,encmin do 61 i=1,20 xtest(i)=-xl+0.05*(i-0.5)*2.*xl ytest(i)=0. ztest(i)=0. 61 continue do 62 i=1,10 xtest(20+i)=0. ytest(20+i)=-yl+0.10*(i-0.5)*2.*yl ztest(20+i)=0. 62 continue do 621 i=1,60 xtest(30+i)=0. ytest(30+i)=0. ztest(30+i)=zmin*2+(i-0.5)*(zmax-zmin)*2/30. 621 continue do 63 i=1,n30 x0p=xtest(i) y0p=ytest(i)-svak*x0p z0p=ztest(i) if(y0p.gt.yl) y0p=y0p-2.*yl if(y0p.lt.-yl) y0p=y0p+2.*yl px=(x0p+xl)/2./xl*nxg+n128x-0.5 py=(y0p+yl)/2./yl*nyg+n128y-0.5 pz=(z0p-zmin)/(zmax-zmin)*nzg+n128z-0.5 iix=int(px)-n128x iiy=int(py)-n128y iiz=int(pz)-n128z px=px-n128x py=py-n128y pz=pz-n128z dx=px-iix dy=py-iiy dz=pz-iiz tx=1.-dx ty=1.-dy tz=1.-dz ix1=iix+1 ix2=iix+2 iy1=iiy+1 iy2=iiy+2 iz1=iiz+1 iz2=iiz+2 if(ix1.gt.n128x) ix1=ix1-n128x if(ix2.gt.n128x) ix2=ix2-n128x if(iy1.gt.n128y) iy1=iy1-n128y if(iy2.gt.n128y) iy2=iy2-n128y if(ix1.lt.1) ix1=ix1+n128x if(ix2.lt.1) ix2=ix2+n128x if(iy1.lt.1) iy1=iy1+n128y if(iy2.lt.1) iy2=iy2+n128y c mass falling outside the grid if(iz2.gt.nzg) tz=0. if(iz2.lt.1) tz=0. if(iz1.gt.nzg) dz=0. if(iz1.lt.1) dz=0. if(iz1.lt.1) iz1=1 if(iz2.lt.1) iz2=1 if(iz1.gt.nzg) iz1=nzg if(iz2.gt.nzg) iz2=nzg fxapu= + dx*dy*dz*fxgrid3(ix2,iy2,iz2)+ + dx*dy*tz*fxgrid3(ix2,iy2,iz1)+ + dx*ty*dz*fxgrid3(ix2,iy1,iz2)+ + dx*ty*tz*fxgrid3(ix2,iy1,iz1)+ + tx*dy*dz*fxgrid3(ix1,iy2,iz2)+ + tx*dy*tz*fxgrid3(ix1,iy2,iz1)+ + tx*ty*dz*fxgrid3(ix1,iy1,iz2)+ + tx*ty*tz*fxgrid3(ix1,iy1,iz1) fyapu= + dx*dy*dz*fygrid3(ix2,iy2,iz2)+ + dx*dy*tz*fygrid3(ix2,iy2,iz1)+ + dx*ty*dz*fygrid3(ix2,iy1,iz2)+ + dx*ty*tz*fygrid3(ix2,iy1,iz1)+ + tx*dy*dz*fygrid3(ix1,iy2,iz2)+ + tx*dy*tz*fygrid3(ix1,iy2,iz1)+ + tx*ty*dz*fygrid3(ix1,iy1,iz2)+ + tx*ty*tz*fygrid3(ix1,iy1,iz1) fzapu= + dx*dy*dz*fzgrid3(ix2,iy2,iz2)+ + dx*dy*tz*fzgrid3(ix2,iy2,iz1)+ + dx*ty*dz*fzgrid3(ix2,iy1,iz2)+ + dx*ty*tz*fzgrid3(ix2,iy1,iz1)+ + tx*dy*dz*fzgrid3(ix1,iy2,iz2)+ + tx*dy*tz*fzgrid3(ix1,iy2,iz1)+ + tx*ty*dz*fzgrid3(ix1,iy1,iz2)+ + tx*ty*tz*fzgrid3(ix1,iy1,iz1) c _d2 c if(abs(z0p).gt.zmax) fzapu=fzapu*abs(z0p)/zmax c _tracer c extrapolation: include forces beyond encmin c encmin < enclim if(abs(z0p).gt.zmax) then finf=2.*pii*gg*tmass/4./xl/yl fac=(1.+(encmin/z0p)**2)**(-.5) fzapu=-finf*fac*z0p/abs(z0p) endif if(abs(z0p).lt.zmax) then finf=2.*pii*gg*tmass/4./xl/yl fac=(1.+(enclim/z0p)**2)**(-.5) fzapu=fzapu-finf*fac*z0p/abs(z0p) endif fxapu31(i)=fxapu/omega2 fyapu31(i)=fyapu/omega2 fzapu31(i)=fzapu/omega2 63 continue write(31) (xtest(i),i=1,n30) write(31) (ytest(i),i=1,n30) write(31) (ztest(i),i=1,n30) write(31) (fxapu31(i),i=1,n30) write(31) (fyapu31(i),i=1,n30) write(31) (fzapu31(i),i=1,n30) endif endif CALL CPU(TTIME2) cpugrav4=cpugrav4+ttime2-ttime1 c encmin = limit for nearby gravity c enclim = limit for near+far gravity c encmin < enclim -> use fft for the zone from encmin to enclim c encmin no gravity due to infinite sheet c gsheet=1 -> add gravity due to infinite sheet c***************************************************************** c***************************************************************** apuapu=0.d0 apuapu2=0.d0 if(gsheet.ge.1) then c normal particles finf=2.d0*pii*gg*tmass/4.d0/xl/yl do 80 i=1,lkm agrav_sheet(i)=0.d0 if(r(i,3).ne.0) then rajainf=enclim if(abs(r(i,3)).gt.zmax) rajainf=encmin fac=1.d0/sqrt((1.d0+(rajainf/r(i,3))**2)) agrav(i,3)=agrav(i,3)-finf*fac*r(i,3)/dabs(r(i,3)) agrav_sheet(i)=-finf*fac*r(i,3)/dabs(r(i,3)) apuapu=apuapu-finf*fac*r(i,3)/dabs(r(i,3))*xmass(i) apuapu2=apuapu2+agrav(i,3)*xmass(i) endif 80 continue c tracer particles c for particles beyond grid: add forces from > encmin c for other particles: add forces from > enclim if(n_tracer.gt.0) then do 180 i=1,n_tracer if(r_tracer(i,3).ne.0) then rajainf=enclim if(abs(r_tracer(i,3)).gt.zmax) rajainf=encmin fac=(1.+(rajainf/r_tracer(i,3))**2)**(-.5) agrav_tracer(i,3)=agrav_tracer(i,3)- + finf*fac*r_tracer(i,3)/abs(r_tracer(i,3)) endif 180 continue endif c add checks: sum of vertical gravity should be zero c contribution from near+far particles+distgrav: fz_check3 if(mod(simustep,i200).eq.0) then fz_check3=0. azmin=1d8 azmax=-1d8 do 103 i=1,lkm fz_check3=fz_check3+xmass(i)*agrav(i,3) if(agrav(i,3).lt.azmin) azmin=agrav(i,3) if(agrav(i,3).gt.azmax) azmax=agrav(i,3) 103 continue fz_check3=fz_check3/tmass/omega2 azmin3=azmin/omega2 azmax3=azmax/omega2 write(6,1104) 'f_check3:',orb,fz_check3 write(6,1104) 'f_amin3: ',orb,azmin3 write(6,1104) 'f_amax3: ',orb,azmax3 write(6,*) 'apuapu',apuapu/omega2/tmass,apuapu2/tmass endif c gsheet=1 ? endif c***************************************************************** c***************************************************************** c CORRECT MEAN Az to zero? c addition of distgrav can give a slightly non-zero mean acceleration, c which can gradually accumulate to give large non-zero mean z! c correct this c unless it is caused by the moonlet c if(n_moon.eq.0) then c do 81 i=1,lkm c agrav(i,3)=agrav(i,3)-apuapu2/tmass c 81 continue c endif c the correction is now made optional c gcorr=0 -> no correction c gcorr=1 -> correction if(gcorr.ge.1) then do 81 i=1,lkm agrav(i,3)=agrav(i,3)-apuapu2/tmass 81 continue if(mod(simustep,i200).eq.0) then fz_check4=0. azmin=1d8 azmax=-1d8 do 104 i=1,lkm fz_check4=fz_check4+xmass(i)*agrav(i,3) if(agrav(i,3).lt.azmin) azmin=agrav(i,3) if(agrav(i,3).gt.azmax) azmax=agrav(i,3) 104 continue fz_check4=fz_check4/tmass/omega2 azmin4=azmin/omega2 azmax4=azmax/omega2 write(6,1104) 'f_check4:',orb,fz_check4 write(6,1104) 'f_amin4: ',orb,azmin4 write(6,1104) 'f_amax4: ',orb,azmax4 write(6,*) '-------------------------------------------' endif endif 1104 format(a10,f12.3,3e20.6) 1105 format(a35,g14.3,2f12.3,i6) c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c**************************************************** c**************************************************** c OUTPUT TO GIDL: 30: dgravi1_read.pro c includes near+grid+sheet c add separately that due to sheet c**************************************************** c**************************************************** if(ival.eq.2) then write(30) (agrav(i,1)/omega2,i=1,lkm) write(30) (agrav(i,2)/omega2,i=1,lkm) write(30) (agrav(i,3)/omega2,i=1,lkm) c added with version 3 write(30) (agrav_sheet(i)/omega2,i=1,lkm) write(30) tgrid,igrav nnxg=nxg nnyg=nyg nnzg=nzg write(30) + (((fxgrid3(i,j,k)/omega2,k=1,nnzg),j=1,nnyg),i=1,nnxg) write(30) + (((fygrid3(i,j,k)/omega2,k=1,nnzg),j=1,nnyg),i=1,nnxg) write(30) + (((fzgrid3(i,j,k)/omega2,k=1,nnzg),j=1,nnyg),i=1,nnxg) write(30) cpugrav1,cpugrav2,cpugrav3,cpugrav4 endif c ************************************************** c ************************************************** c data collection (moved from lim1) c ************************************************** c ************************************************** C ******* TABULATE IXMIN ETC IF(IX.LT.IXMIN) IXMIN=IX IF(IX.GT.IXMAX) IXMAX=IX IXL=IXL+1 IXS=IXS+IX C ******* WARNING MESSAGE IF TOO MANY PAIRS IF(IX.GT..9*NIX) WRITE(IDEV6,*)'WARN: IX > ',0.9*NIX,IX CALL CPU(TIME2) CPUenc=CPUenc+TIME2-TIME1 c cpugrav1 = CPU spent in sorting c cpugrav2 = CPU spent in tabulating distant gravity c cpugrav3 = CPU spent in calculation of nearby gravity c cpugrav4 = CPU spent in interpolation of distant gravity if(mod(simustep,20).eq.0.or.simustep.le.5) then write(6,9049) 'grav: sort,tabu,near,far ', + cpugrav1,cpugrav2,cpugrav3,cpugrav4 write(49,4949) orb,cpugrav1,cpugrav2,cpugrav3,cpugrav4, + cpupota,cpuint,cputot WRITE(*,*) 'POT1,POT2,POT3,POT4', + npotsum1,npotsum2,npotsum3,npotsum4 endif 4949 format(f8.3,f6.1,2f10.1,f6.1,3f10.1) 9049 format(A30, 4g10.3) c---------------------------------------------------------- c GRAVITY FELT BY THE MOONLET IS SET TO ZERO c added 15/01/03 if(n_moon.ne.0) then agrav(1,1)=0. agrav(1,2)=0. agrav(1,3)=0. aagrav(1,1)=0. aagrav(1,2)=0. aagrav(1,3)=0. endif RETURN END c--------------------------------------------------------------------- SUBROUTINE FOURN2007(DATA,NN,NDIM,ISIGN) double precision WR,WI,WPR,WPI,WTEMP,THETA dimension NN(NDIM),DATA(*) NTOT=1 DO 11 IDIM=1,NDIM NTOT=NTOT*NN(IDIM) 11 CONTINUE NPREV=1 DO 18 IDIM=1,NDIM N=NN(IDIM) 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 c*********************************** subroutine tab_init c*********************************** include 'ssbox.inc' dimension ix_mc(nstar),iy_mc(nstar) c construct tables used for limiting the number of c particles checked for intersections: c divides particles into mx*my subcells c cxtab_mc(ic),cytab_mc(ic) = centers of cell ic c dtab_mc(ic) = maximum extent of particles in cell ic c dtab2=dtab**2 c indtab_mc contain indices of all particles, c sorted according to subcells: c indtab_mc(mintab_mc(ic):mintab_mc(ic)+numtab_mc(ic)) c --> particles in cell ic c---------------------------------------------- c 1) search subcell-index for each particle c there are particles whose centers are outside the actual region c (= dublicate particles in pos_read) c -> sub-cells must cover also these particles c (this is why for example xmax is not equal to xl) xmin=-xl ymin=-yl xmax=xl ymax=yl npart=lkm xmin=xmin*1.01d0 ymax=ymax*1.01d0 ymin=ymin*1.01d0 dxcell=(xmax-xmin)/mx_mc dycell=(ymax-ymin)/my_mc do 20 i=1,npart ix_mc(i)=int((r(i,1)-xmin)/dxcell)+1 iy_mc(i)=int((r(i,2)-ymin)/dycell)+1 20 continue c---------------------------------------------- c 2) place indices etc. lask=0 indmui=0 do 30 i=1,mx_mc do 31 j=1,my_mc lask=lask+1 cxtab_mc(lask)=xmin+(i-0.5)*dxcell cytab_mc(lask)=ymin+(j-0.5)*dycell nummui=0 do 32 k=1,npart if(ix_mc(k).eq.i.and.iy_mc(k).eq.j) then nummui=nummui+1 indmui=indmui+1 indtab_mc(indmui)=k endif 32 continue numtab_mc(lask)=nummui mintab_mc(lask)=indmui-nummui+1 dtab_mc=sqrt((dxcell/2.)**2+(dycell/2.)**2) dtab2_mc(lask)=(dtab_mc+encmin)**2 31 continue 30 continue return end