c simple1707_gravi.f c*********************************************************** SUBROUTINE GRAVIGRID_mc(ival) c*********************************************************** c calculate separately FFT-grid gravity and c PP gravity due to nearby particles: HS 11/08/98 --> c in every step: c 1) calculate direct PP 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 an infinite sheet c from ENCLIM < DIST c----------------------------------------------------------------------- c GRAVITY-check output (subroutines gidl30 and gidl31) 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 -> gravity calculation of gravity can be checked in detail 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----------------------------------------------------------------------- c CPU-counters: c cpugrav1 = CPU spent in tab_init (2D-tabulation) c cpugrav2 = CPU spent in constructing fft gravity grid c cpugrav3 = CPU spent in calculation of pp gravity c cpugrav4 = CPU spent in interpolation of fft gravity forces c PARTICLE_PAIR counters: c npotsum1 = checked for distance with tab c npotsum2 = accepted to contribute to gravity c npotsum3 = ero < raja1 c npotsum4 = ero < dero --> makes to impact lists (same as ix) C*********************************************************** INCLUDE 'ssbox.inc' C SIIRRA INCLUDEEN: c gravity from infite sheet c dimension agrav_sheet(nstar) c for FFT c 3D example : maximum true grid-size 512*256*8 -> 176 megas c parameter (nxmax=512,nymax=256,nzmax=8) c parameter (nzmax2=2*nzmax,ndatamax=4*nxmax*nymax*nzmax) c parameter declaration moved to include-file 310503 C SIIRRA INCLUDEEN c real rhogrid3(nxmax,nymax,nzmax) c real fxgrid3(nxmax,nymax,nzmax) c real fygrid3(nxmax,nymax,nzmax) c real fzgrid3(nxmax,nymax,nzmax) c real frrho3(nxmax,nymax,nzmax2) c real firho3(nxmax,nymax,nzmax2) c real distx3(nxmax,nymax,nzmax2) c real disty3(nxmax,nymax,nzmax2) c real distz3(nxmax,nymax,nzmax2) c real frdist3(nxmax,nymax,nzmax2) c real fidist3(nxmax,nymax,nzmax2) real data(ndatamax) dimension nnn(3) c 2d-partitioning: local variables dimension pottab_mc(nstar2) dimension n_pottab_mc(nstar2) dimension m_pottab_mc(nstar2) CALL CPU(TIME1) c------------------------------------------------------------ c MISCELLENEOUS VARIABLES c interval for f_checks (in steps) i200=100 c relative tangential shift of image boxes svak=bvak/2.d0/xl encmin2=encmin**2 c FFT-grid related 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 of FFT grid defined by fzrange c fzrange > 0 use zmax=fzrange*stdev(zz) c fzrange < 0 use zmax=abs(fzrange) c----------------------------------------------- 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 apu_zmin=zmin apu_zmax=zmax c apu_zmin,apu_zmax added to common c*********************************************** c*********************************************** c TABULATE FFT GRAVITY c forces from xp,yp,zp density with FFT c encmin = limit for pp gravity c enclim = limit for pp+fft gravity c encmin < enclim -> use fft for the zone from encmin < delta < enclim if(encmin.lt.enclim) then CALL CPU(TTIME1) c---------------------------------------------------------------- c calculate density rhogrid3 c---------------------------------------------------------------- c in the beginning of the run check the time consumptions if(time.lt.period/20.) CALL CPU(TTTIME1) c clear density arrays 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 c use cloud-in-cell (cic)-assignment to 8 nearest grid-points 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 if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET, CPU-RHOGRID',tttime2-tttime1 write(6,*) rhogrid3(1,1,nzg/2),rhogrid3(1,1,nzg/2+1) CALL CPU(TTTIME1) endif c------------------------------------------------------- c construct distx3,disty3,distz3 c------------------------------------------------------- c = dx/(dx**2+dy**2+dz**2+eps**2) etc c make sure that forces are symmetric i_yla=n64x-1 j_yla=n64y-1 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 422 continue 421 continue 420 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-DIST3',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 END CALCULATION OF FFT GRID-FORCES c encmin makes to impact lists (same as ix) npotsum1=0 npotsum2=0 npotsum3=0 npotsum4=0 c go through all target particles ii to which gravity is calculated do 20 ii=1,lkm-1 xt=r(ii,1) yt=r(ii,2) c choose particles jj potentially contributing to gravity at particle ii c for particles ii near the borders of the box take c into account forces exerted on their image particles c at xt + n*avak c yt + n*bvak + m*cvak c to find jj's use 2d-partititoning: c npot=temporary counter of number of potential impactors jj c these will be stored to c pottab_mc(npot) = index of jj c n_pottab_mc(npot)= radial image shift n c m_pottab_mc(npot)= tangential image shift m npot=0 c limits for n and m c defauls n_mc1=0 n_mc2=0 m_mc1=0 m_mc2=0 c near borders 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 c loops over potential images of ii c (usually just the particle itself) do 9008 n_mc=n_mc1,n_mc2 do 9009 m_mc=m_mc1,m_mc2 c loop over 2D-partitioning table: c indtab_mc(mintab_mc(imc):mintab_mc(imc)+numtab_mc(imc)) c --> indexes of particles in cell imc c the particles in cell imc can contribute at c squared distance dtab2(imc) from its center [cxtab(imc),cytab(imc)] 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 c can this cell contribute? if(d2.gt.dtab2_mc(imc)) goto 9010 c yes -> add the particles of the cell imc top 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 c-------------------------------------------------------------- c have found npot potential contributors for dist < encmin fron 2D-tabulation c next check which of them are really inside encmin: do 30 j=1,npot jj=pottab_mc(j) c make sure each pair is included only once 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 pp gravitational forces for dist 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 call gravi_out_terminal(4) endif endif c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c**************************************************** c**************************************************** c OUTPUT TO GIDL: 30: c total near+grid+sheet if(ival.eq.2) then call gravi_out_gidl30(2) 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 2D-partitioning c cpugrav2 = CPU spent in tabulating FFT gravity c cpugrav3 = CPU spent in calculation of nearby PP gravity c cpugrav4 = CPU spent in interpolation of FFT gravity if(mod(simustep,100).eq.0.or.simustep.le.5) then WRITE(*,*) 'GRAV: orb= ',orb write(6,9049) ' GRAV: 2dtab, pp, grid, interp ', + cpugrav1,cpugrav3,cpugrav2,cpugrav4 write(*,*) 'GRAV: npot1,2,3,4 = ', + npotsum1,npotsum2,npotsum3,npotsum4 endif 4949 format(f8.3,f6.1,2f10.1,f6.1,3f10.1) 9049 format(A32, 4g10.3) 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