;
; 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,mintime,0.0
default,wait,0.0
;
; Read other values than the default ones from file
;
nbins=dim.nz
@distfile
;
; Print some defaults
;
print,'nbins=',nbins
print,'wait=',wait
print,'mintime=',mintime
;
; Find the length of the domain and width of each bin
;
Lz=start.xyz1[2]-start.xyz0[2]
dz=Lz/nbins
;
; Find how many different particle radii we have
;
npart_radii=0
dims=size(start.ap0)
if (dims[0] gt 0) then begin
    ninit=dims[1]
    for i=0,ninit-1 do begin
        if (start.ap0[i] ne 0) then begin
            npart_radii=npart_radii+1
        endif
    end
endif else begin
    npart_radii=1
end
print,'npart_radii=',npart_radii
;
; Define some arrays
;
bin=indgen(nbins,npart_radii)*0.
zbin=indgen(nbins)*Lz/nbins-Lz/2.
minbin=bin
maxbin=bin
;
; Find the size of the arrays we are handling
;
itime_max=n_elements(p.zp[0,*])*1l
ratio =fltarr(npart_radii,itime_max)
ratio3=fltarr(npart_radii,itime_max)
npart=n_elements(p.zp[*,0])*1l
print,'npart=',npart
print,'Number of times=',itime_max

radius=p.ap[*,0]
iradius=intarr(npart)*0
small_number=1e-8
for iap=0,npart_radii-1 do begin
    iradius(where((radius lt (start.ap0(iap)+small_number)) and $
                  (radius gt (start.ap0(iap)-small_number))))=iap
end
;
; Check if there is gravity
;
lgrav=0
grav=0.
if start.gravz gt 0 then begin
    lgrav=1
    grav=start.gravz
endif
;
; Define which indeces are of interest for plotting
;
ind2=nbins/2
if (lgrav) then begin
    ind1=nbins/4
endif else begin
    ind1=0
end
;
; Perform timeloop
;
itimemin=0
for itime=0,itime_max-1 do begin
    if (p.t[itime] gt mintime) then begin
        bin=bin*0.
;
; Find the index in the bin directory where all particles belong
;
        ibin=floor(nbins*(p.zp[*,itime]+Lz/2.)/Lz)
;
; Loop over all particle sizes and bins and count particles
;
        for iap=0,npart_radii-1 do begin
            for i=0,nbins-1 do begin
                bin[i,iap]=bin[i,iap]+$
                  n_elements(where((ibin eq i)and(iradius eq iap)))
                if (n_elements(where((ibin eq i)and(iradius eq iap)))eq 1) then begin
                    if (where((ibin eq i)and(iradius eq iap))eq -1) then begin
                        bin[i,iap]=bin[i,iap]-1
                    endif
                endif
            end
        end
;
; Store results for each timestep in arrays and find average over all times
;
        if (itimemin gt 0) then begin
            for i=0,npart_radii-1 do begin
                rat=bin_aver[ind2,i]/bin_aver[ind1,i]
                rat3=bin[ind2,i]/bin[ind1,i]
                bin_aver[*,i]=(bin_aver[*,i]*itimemin+bin[*,i])/(itimemin+1)
                ratio[i,itimemin]=rat
                ratio3[i,itimemin]=rat3
            end
            time=[time,p.t[itime]]
            minbin=[minbin,min(bin)]
            maxbin=[maxbin,max(bin)]
        endif else begin
            bin_aver=bin
            for i=0,npart_radii-1 do begin
                rat=max(bin_aver[*,i])/min(bin_aver[*,i])
                ratio[i,0]=rat
                ratio3[i,0]=rat
            end
            time=p.t[itime]
            minbin=min(bin)
            maxbin=max(bin)
        endelse
;
; Plot results for each timestep
;
        title='time='+str(p.t[itime])
        plot,zbin,bin[*,0],title=title
        oplot,zbin,bin_aver[*,0],col=122
        for iap=1,npart_radii-1 do begin
            oplot,zbin,bin[*,iap],li=iap
            oplot,zbin,bin_aver[*,iap],li=iap,col=122
        end
        wait,wait
        itimemin=itimemin+1
    end
end
;
; Find total average for time larger than mintime
;
index=where(time gt mintime)
average=fltarr(npart_radii)
for i=0,npart_radii -1 do begin
    average[i]=mean(ratio3(i,index))
    print,'average(',start.ap0[i],')=',average[i]
end
;
; Find total Ntilde for time larger than mintime
;
Ntilde=fltarr(npart_radii)
for i=0,npart_radii -1 do begin
    Ntilde[i]=bin_aver[nbins/2,i]/bin_aver[0,i]
    print,'Ntilde(',start.ap0[i],')=',Ntilde[i]
end
;
; Plot particle number density profiles averaged over all times
;
!p.multi=[0,1,2]
plot,zbin,bin_aver[*,0],$
  yrange=[min(bin_aver)*0.9,max(bin_aver)*1.1],title=title,$
  xtitle='!8z!6',ytit='!8N!6'
for iap=1,npart_radii-1 do begin
    oplot,zbin,bin_aver[*,iap],li=iap
end
;
; Find yrange
;
if (lgrav) then begin
    ymax=max([max(ratio[5:*]),max(ratio3[5:*])])*1.1
    ymin=min([min(ratio[5:*]),min(ratio3[5:*])])*1.1
endif else begin
    ymin=1.
    ymax=max([max(ratio),max(ratio3)])*1.1
end
;
; Plot time evolution of min and max ratios for the different particle sizes
;
plot,time,ratio[0,*],yrange=[ymin,ymax],$
  ytit='max(N)/min(N)',xtit='time'
oplot,time,ratio3[0,*],col=122
for iap=1,npart_radii-1 do begin
    oplot,time,ratio[iap,*],li=iap
    oplot,time,ratio3[iap,*],li=iap,col=122
end
;
; Save data
;
save,file='./data/N_N0.dat',average,zbin,bin_aver,Ntilde
;
!p.multi=[0,1,1]

END