% create circular sampling rings
%t=ones(Nx,Ny,Nz);
Mp1=Mp;
Dp1=Dp;

  for kk=1:4
  pv1=smooth3(pv1);
  Mp1=smooth3(Mp1);
  Dp1=smooth3(Dp1);
  end


Nn=Nx/2;

rad=[1:1:(Nn-1)];
ring=ones(Nx,Ny,Nz);
weight=zeros(Nn-1,1);
Dring=zeros(Nx,Ny,Nz);
Mring=zeros(Nx,Ny,Nz);
radialmean=zeros(1,Nn-1);
XCenter = (Nx+1)/2 ;
YCenter = (Ny+1)/2 ;


for k=1:Nn-1
  for i=1:Nx
    for j=1:Ny
      radius(i,j) = sqrt((i-XCenter)^2+(j-YCenter)^2);
    end
  end
  radius(find( round(radius) > (k+1) ))=radius(find (round(radius) > (k+1) ))*0;
  radius(find( round(radius) < (k)   ))=radius(find (round(radius) < (k)   ))*0;
  radius(find(radius>0))=1.0;
  ring(:,:,k)=radius;
  weight(k)=length(find(ring(:,:,k) > 0));
end


for k=1:Nz
for i=1:Nn-1
Tring(:,:,k)=meantheta(:,:,k).*ring(:,:,i);
Dring(:,:,k)=Dp1(:,:,k).*ring(:,:,i);
Mring(:,:,k)=Mp1(:,:,k).*ring(:,:,i);
meanDp(i,k)=Nx*Ny*mean(mean(squeeze(Dring(:,:,k))))/weight(i);
meanMp(i,k)=Nx*Ny*mean(mean(squeeze(Mring(:,:,k))))/weight(i);
meanT(i,k)=Nx*Ny*mean(mean(squeeze(Tring(:,:,k))))/weight(i);
end
end

rad=dx:dx:(Nn-1)*dx;

for k=1:Nz
radd(:,k)=rad(:);
end

c=-0.05:0.005:0.005;
figure
A=meanDp(:,5:26).*radd(:,5:26);

i=find(A>0.00);
A(i)=A(i)/2;
%contourf(flipud(A'),c)
imagesc(flipud((A(1:Nn-5,:)')));
set(gca,'DataAspectRatio',[2,2,2])
        set(gca,'ydir','norm')
xlabel('r ','FontSize',20)
ylabel('z ','Rotation',[0],'FontSize',20)
set(gca,'XtickLabel','|||')
set(gca,'YtickLabel','|||')
% title('D_p, azimuthal average ','Fontsize',20);
text(20,Nz-5,'D_p, azimuthal average ','Fontsize',20); 

hold on
contour(flipud(meanT(1:Nn-5,5:26)')/1000,10);

%subplot(2,1,2)
%contourf(flipud((meanMp(:,5:26).*radd(:,5:26))'),-c);colorbar;

%--------TEST---------------------

av=6.28*mean(mean(meanDp(:,5:26).*radd(:,5:26)))*(Nn-1)*Nz*dx*dz
av1=mean(mean(mean(Dp1(:,:,5:26))))*Nx*Ny*(26-5)*dx*dy*dz
% vol=Nz*dz*3.14*((Nn-1)*dx)*((Nn-1)*dx)

avm=6.28*mean(mean(meanMp(:,5:26).*radd(:,5:26)))*(Nn-1)*Nz*dx*dz
avm1=mean(mean(mean(Mp1(:,:,5:26))))*Nx*Ny*(26-5)*dx*dy*dz