;*********************************
;slice.pro
;*********************************
; plot slice of the simulation system at given z

pro slice,xx1,yy1,zz1,ss1,zslice=zslice,xcen=xcen,ycen=ycen,wid=wid,$
plot=plot,warn=warn,pos=pos,title=title,index=index,vx=vx,vy=vy,dt=dt,$
indexr=indexr,radi=radi,ff=ff,oplot=oplot,noplot=noplot,lim=lim,col=col,xl=xl,$
suppress=suppress,n50_0=n50_0,top=top,meanslice=meanslice,meandens=meandens


if(n_params() le 0) then begin
print,'pro slice,xx,yy,zz,ss,zslice=zslice,xcen=xcen,ycen=ycen,wid=wid'
print,'slice of the simulation systen at given z'
print,'/plot -> plot all centers'
print,'/oplot -> plot all centers'
print,'ff -> returns total area cut by the plane z=zslice'
print,'/warn -> plot overlaps'
print,'vx,vy,dt'
print,'pos -> define position of plot'
print,'title'
print,'lim=index only these particles'
print,'xcen,ycen'
print,'/appl1/heikki/UNAMIDL/slice.pro'
return
end


if(keyword_set(title)) then begin
nwin
!p.title=title
endif

if(not keyword_set(zslice)) then zslice=0
if(not keyword_set(xcen)) then xcen=0
if(not keyword_set(ycen)) then ycen=0
if(not keyword_set(dt)) then dt=314.
if(keyword_set(wid)) then xylim,wid

xx=xx1-xcen
yy=yy1-ycen
zz=zz1
ss=ss1

if(keyword_set(wid)) then begin
tindex=where(abs(xx) le wid*1.1 and abs(yy) le wid*1.1)
xx=xx(tindex)
yy=yy(tindex) 
zz=zz(tindex)
ss=ss(tindex)
endif

if(keyword_set(lim)) then begin
xx=xx(lim)
yy=yy(lim) 
zz=zz(lim)
ss=ss(lim)
endif

pcol=!p.color
if(keyword_set(col)) then !p.color=col

n50=50
if(keyword_set(n50_0)) then n50=n50_0
n51=n50+1

fii=findgen(n51)/n50*2.*!pi
sinf=sin(fii)
cosf=cos(fii)

index=where(abs(zz-zslice) le ss, count)
if(count lt 1) then begin
return
endif

indexr=index
xx=xx(index)
yy=yy(index) 
zz=zz(index)
ss=ss(index)
radi=sqrt(ss^2-(zz-zslice)^2)
ff=total(radi^2*!pi)
if(not keyword_set(noplot)) then begin

if(not keyword_set(top)) then begin
if(not keyword_set(oplot)) then begin
if(keyword_set(pos)) then plot,xx,yy,psym=3,pos=pos
if(not keyword_set(pos)) then plot,xx,yy,psym=3
endif
if(keyword_set(oplot)) then oplot,xx,yy,psym=3
endif

if(keyword_set(vx)) then begin
vxmean=mean(vx(index))
vymean=mean(vy(index))
endif

for i=0l,count-1 do begin
plots,xx(i)+radi(i)*cosf,yy(i)+radi(i)*sinf

if(keyword_set(vx)) then begin
oplot,[xx(i),xx(i)+(vx(i)-vxmean)*dt],[yy(i),yy(i)+(vy(i)-vymean)*dt],col=3
if(vx(i)-vxmean gt 0) then begin
oplot,[xx(i),xx(i)+(vx(i)-vxmean)*dt],[yy(i),yy(i)+(vy(i)-vymean)*dt],col=2
endif
endif

if(not keyword_set(suppress)) then begin

;********** warn *************
if(keyword_set(warn)) then begin
dist2=sqrt((xx-xx(i))^2+(yy-yy(i))^2+(zz-zz(i))^2)
dist2(i)=1000.
index2=where(dist2 le warn*(ss+ss(i)),count2)
if(count2 ge 1) then begin
plots,xx(i)+radi(i)*cosf,yy(i)+radi(i)*sinf,col=2,thick=3
for j=0l,count2-1 do begin
plots,[xx(i),xx(index2(j))],[yy(i),yy(index2(j))],col=3,thick=1
endfor
endif
endif

if(not keyword_set(warn)) then begin
dist2=sqrt((xx-xx(i))^2+(yy-yy(i))^2+(zz-zz(i))^2)
dist2(i)=1000.

index2=where(dist2 le 0.99*(ss+ss(i)),count2)
if(count2 ge 1) then begin
plots,xx(i)+radi(i)*cosf,yy(i)+radi(i)*sinf,col=3,thick=3
print,'dist <0.99',1.*count2
endif

index2=where(dist2 le 0.95*(ss+ss(i)),count2)
if(count2 ge 1) then begin
plots,xx(i)+radi(i)*cosf,yy(i)+radi(i)*sinf,col=8,thick=3
print,'dist <0.95',1.*count2
endif

index2=where(dist2 le 0.90*(ss+ss(i)),count2)
if(count2 ge 1) then begin
plots,xx(i)+radi(i)*cosf,yy(i)+radi(i)*sinf,col=2,thick=3
print,'dist <0.90',1.*count2
endif

index2=where(dist2 le 0.8*(ss+ss(i)),count2)
if(count2 ge 1) then begin
plots,xx(i)+radi(i)*cosf,yy(i)+radi(i)*sinf,col=8,thick=3
print,'dist <0.8',1.*count2
endif
index2=where(dist2 le 0.5*(ss+ss(i)),count2)
if(count2 ge 1) then begin
plots,xx(i)+radi(i)*cosf,yy(i)+radi(i)*sinf,col=7,thick=3
print,'dist <0.5',1.*count2
endif

endif
endif

endfor

if(keyword_set(plot)) then oplot,xx,yy,psym=3


endif	;noplot
;print,'total area: ',ff
if(keyword_set(xl)) then begin
print,' TAU   = ',total(!pi*ss1^2)/xl/xl/4.
print,' FF(0) = ',ff/xl/xl/4.
print,'STD(Z) = ',stdev(zz1)
endif

if(keyword_set(meanslice)) then begin
H=0.5*sqrt(6.)*stdev(zz1)
ztab=-h+(findgen(abs(meanslice))+.5)/abs(meanslice)*2*h
;print,ztab
ffsum=0.
for iz=0,abs(meanslice)-1 do begin
index=where(abs(zz1-ztab(iz)) le ss1, count)
if(count gt 1) then begin
indexr=index
zapu=zz1(index)-ztab(iz)
sapu=ss1(index)
radi2=sapu^2-zapu^2
ff=total(radi2*!pi)
if(meanslice lt 0) then print,ztab(iz),ff/xl/xl/4.

ffsum=ffsum+ff
endif
endfor
ffsum=ffsum/abs(meanslice)
meandens=ffsum/xl/xl/4.
print,'meandens=',meandens
endif



!p.color=pcol
if(keyword_set(wid)) then xylim,0
if(keyword_set(title)) then !p.title=''
end