;
; Plot the particle density averaged over the xy plane
;
!p.charsize=2
!p.multi=[0,1,1]
device,decompose=0
loadct,5
;
; Read in required data 
;
pc_read_param,obj=start
pc_read_pstalk,obj=p
pc_read_dim,obj=dim
;
; Set some defaults
;
default,nbins,20
default,wait,0.01
default,mintime,0.0
default,showplot,1
;
; Read other values than the default ones from file
;
@distfile
nbins=dim.nz
;
; Print some defaults
;
print,'nbins=',nbins
print,'wait=',wait
print,'mintime=',mintime
print,'showplot=',showplot
;
; Find the length of the domain and width of each bin
;
Lz=start.xyz1[2]-start.xyz0[2]
dz=Lz/nbins
;
; Define some arrays
;
bin=indgen(nbins)*0.
zbin=indgen(nbins)*Lz/nbins-Lz/2.
minbin=bin
maxbin=bin
;
; Find the size of the arrays we are handling
;
npart=n_elements(p.zp[*,0])*1l
itime_max=n_elements(p.zp[0,*])*1l
nsliceaver=10
print,'npart=',npart
print,'Number of times=',itime_max
;
; Perform timeloop
;
itimemin=0
for itime=0,itime_max-1 do begin
;
; Do loop over all particles and put them into bins
;
    bin=bin*0.
    for i=0l,npart-1 do begin
        ibin=floor(nbins*(p.zp[i,itime]+Lz/2.)/Lz)
        bin[ibin]=bin[ibin]+1
    end
;
; Plot results for each timestep
;
    title='time='+str(p.t[itime])
    if (itime gt 0) then begin
        rat=max(bin_aver)/min(bin_aver)
        rat3=max(bin)/min(bin)
        bin_aver=(bin_aver*itime+bin)/(itime+1)
        oplot,zbin,bin_aver,li=2
        ratio=[ratio,rat]
        ratio2=[ratio2,rat2]
        ratio3=[ratio3,rat3]
        time=[time,p.t[itime]]
        minbin=[minbin,min(bin)]
        maxbin=[maxbin,max(bin)]
    endif else begin
        bin_aver=bin
        rat=max(bin_aver)/min(bin_aver)
        ratio=rat
        ratio2=rat
        ratio3=rat
        time=p.t[itime]
        minbin=min(bin)
        maxbin=max(bin)
    endelse
    amod=itime MOD nsliceaver
    if (amod eq 0) then begin
        bin_aver2=bin
        rat2=max(bin_aver2)/min(bin_aver2)
        iitime=0
    endif else begin
        iitime=iitime+1
        bin_aver2=(bin_aver2*iitime+bin)/(iitime+1)
        rat2=max(bin_aver2)/min(bin_aver2)
;        print,'numberdensity ratios=',rat,rat2,rat3
    endelse
    if (showplot) then begin
        plot,zbin,bin,yrange=[0,npart*2.5/nbins],title=title,$
          xtitle='!8z!6',ytit='!8N!6',ps=-1
        oplot,zbin,bin_aver2,li=1
        wait,wait
    endif
end
;
; Find total average for time larger than mintime
;
index=where(time gt mintime)
average=mean(ratio3(index))
print,'average=',average
;
; Plot final results
;
!p.multi=[0,1,3]
plot,zbin,bin_aver,yrange=[min(bin_aver)*0.9,max(bin_aver)*1.1],title=title,$
  xtitle='!8z!6',ytit='!8N!6',ps=-1
plot,time,ratio,yrange=[1.,max([max(ratio),average])*1.1],ytit='max(N)/min(N)',xtit='time'
oplot,time,ratio3,li=2
plot_io,time,minbin,ytit='max(N),min(N)',xtit='time',yrange=[max([1,min(minbin)]),max(maxbin)]
oplot,time,maxbin
oplot,time,maxbin*0.+maxbin[n_elements(maxbin)-1],li=2,col=122
oplot,time,maxbin*0.+minbin[n_elements(maxbin)-1],li=2,col=122
!p.multi=[0,1,1]
oplot,time,ratio2,li=2
oplot,time(index),time(index)*0.+average,li=1,col=122
;
; Save data
;
save,file='./data/N_N0.dat',average,zbin,bin_aver
;
END