I want to get a contour plot of number density of atoms in r-z plane. I have xyz trajectory of system which can be converted to (r,z,theta) trajectory where r=sqrt(x^2 + y^2), z=z, and theta = atan2(y/x). I am trying to write a FORTRAN code for this. I am facing difficulty in creating 2D histogram for this. Here, r will be always positive, but in the figure which I attach here from literature has negative value of r also.  Please suggest the way to construct this figure or binning method.