; ericgr3d 9-20-99 Eric Weeks ; ; this function has no relationship with ericgr2d ; modified slightly ERW 6-7-04 ; ; see: http://www.physics.emory.edu/~weeks/idl/gofr0.html ; ;------------------------------------------------------------ function subgr,data,rmin,rmax,deltar,lilnbar=lilnbar,enone=enone dim=n_elements(data(*,0)) nr=(rmax-rmin)/deltar+1 result=fltarr(nr) xmin=min(data(0,*),max=xmax) ymin=min(data(1,*),max=ymax) if (dim eq 3) then begin zmin=min(data(2,*),max=zmax) endif x0=xmin+rmax & x1=xmax-rmax y0=ymin+rmax & y1=ymax-rmax w1=where(data(0,*) gt x0 and data(0,*) lt x1 and \$ data(1,*) gt y0 and data(1,*) lt y1,nw1) if (nw1 eq 0) then message,'no particles ???',/inf rmin2=rmin*rmin rmax2=rmax*rmax npts=n_elements(data(0,*)) one=fltarr(npts)+1.0 if (dim ne 2) then begin vol=(zmax-zmin)*(xmax-xmin)*(ymax-ymin) endif else begin vol=(xmax-xmin)*(ymax-ymin) endelse lilnbar=lilnbar+npts/vol enone=enone+nw1 ; ======================= ; sphere, center at origin, slice off top at z=z0: the surface ; area that is removed is 2piR^2(1-z0/R) so thus the remaining ; surface area (from origin up to z0) is 2 pi R z0, a nice ; simple formula. Makes sense by dimensional analysis... ; ======================= twopi=2*3.14159265358 rvec=findgen(nr)*deltar+rmin normz=1.0/(twopi*rvec*deltar) for i=0L,nw1-1L do begin d0=data(*,w1(i)) dd=one##d0-data dis=total(dd*dd,1) w2=where(dis gt rmin2 and dis lt rmax2,nw2) if (nw2 gt 0) then begin ; if (min(dis(w2)) lt 0.1) then message,'foo',/inf newdis=sqrt(dis(w2)) thehisto=histogram(newdis,max=rmax,min=rmin,binsize=deltar) if (dim eq 3) then begin z0=zmax-d0(2) z1=d0(2)-zmin normz=twopi*rvec*( (rvec