;
; 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=dims
;
; Set some defaults
;
default,wait,0.0
default,mintime,0.0
default,showplot,1
;
; Read other values than the default ones from file
;
@distfile
nbins=dims.nz
;
; Print some defaults
;
print,'nbins=',nbins
print,'wait=',wait
print,'mintime=',mintime
print,'showplot=',showplot
;
; Find length of domain and width of each bin
;
Lz=start.xyz1[2]-start.xyz0[2]
dz=Lz/nbins
;
; Define some arrays
;
bin=indgen(nbins)*0.
bin=dblarr(nbins)*0.D0
zbin=indgen(nbins)*Lz/nbins-Lz/2.
;
; Find the size of the arrays we are handling
;
npart=n_elements(p.zp[*,0])*1l
itime_max=n_elements(p.zp[0,*])*1l
print,'npart=',npart
print,'Number of times=',itime_max
;
; Perform timeloop
;
itimemin=0
for itime=0,itime_max-1 do begin
    if (p.t[itime] gt mintime) then 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
;
; Find average profile
;
        title='time='+str(p.t[itime])
        if (itimemin gt 0) then begin
            bin_aver=(bin_aver*itimemin+bin)/(itimemin+1)
            oplot,zbin,bin_aver,li=2
        endif else begin
            bin_aver=bin*1.D0
        endelse
;
; Plot results for each timestep
;
        if (showplot) then begin
            plot,zbin,bin,yrange=[0,npart*2.5/nbins],title=title,$
              xtitle='!8z!6',ytit='!8N!6',li=1,ps=-1
            oplot,zbin,bin_aver,li=0
            wait,wait
        endif
        itimemin=itimemin+1
    endif
end
;
; Find total average for time larger than mintime
;
average=max(bin_aver)/min(bin_aver)
print,'ratio=',average
;
; Save data
;
if (nbins eq dims.nx) then begin
    save,file='./data/N_N0_real.dat',average,zbin,bin_aver
endif else begin
    file='./data/N_N0'+str(nbins)+'.dat'
    save,file=file,average,zbin,bin_aver
end
;
END