In [1]:

    

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (12.0, 9.0)
%matplotlib inline
import csv


    

In [2]:

    

filenames=[]

for f in range(1,4):
    filename="1E9inserts_1s_{}.csv".format(f)
    filenames.append(filename)
    print filename
#filenames=["2E6inserts_1s.csv"]


    

    




1E9inserts_1s_1.csv
1E9inserts_1s_2.csv
1E9inserts_1s_3.csv

In [3]:

    

arrays=[]
for f in filenames:
     with open(filename,'r') as csvfile:
        for line in csvfile:
            if line[0]=="#":
                continue
            l=line.split()
            arrays.append(np.array([float(l[-3]),float(l[-2]),float(l[-1])]))


    

In [11]:

    

#array=np.loadtxt(filenames[0])


    

In [17]:

    

arrays=np.array(arrays)
max=np.max(pos)
min=np.min(pos)
print min,max


    

    




-8.9514e-07 9.14442e-07

In [46]:

    

hist,edges=np.histogramdd(pos, bins=200,range=[[min,max],[min,max],[min,max]])
print np.sum(hist)
heatmap=hist[:,:,100]
xedges=edges[0]
yedges=edges[1]
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]


    

    




43491259.0

In [57]:

    

print np.sum(hist>0)
print np.sum(hist==0)


    

    




122002
7877998

In [4]:

    

plt.clf()
plt.imshow(np.log(heatmap),extent=extent,interpolation='nearest')
plt.colorbar()
plt.grid()
plt.ticklabel_format(style='sci',scilimits=(0,0))
plt.xticks(rotation='vertical')
plt.gca().set_aspect('equal', adjustable='box')
plt.savefig('mobility.svg',dpi=300)
plt.show()


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-fbdc8f56cabc> in <module>()
      1 plt.clf()
----> 2 plt.imshow(np.log(heatmap),extent=extent,interpolation='nearest')
      3 plt.colorbar()
      4 plt.grid()
      5 plt.ticklabel_format(style='sci',scilimits=(0,0))

NameError: name 'heatmap' is not defined





    






<matplotlib.figure.Figure at 0x7ff742c24e50>

In [ ]:

    

dist=np.sqrt(np.sum(pos**2,axis=1))
print np.mean(dist)
print np.std(dist)


    

In [ ]:

    

x2=np.sqrt(np.mean(pos[:,0]**2))
y2=np.sqrt(np.mean(pos[:,1]**2))
z2=np.sqrt(np.mean(pos[:,2]**2))
print x2,y2,z2


    

In [ ]:

    

plt.hist(dist,1000)
plt.ylabel('Counts')
plt.xlabel('diffusionlength [m]')
plt.savefig("distdiff.svg")
plt.show()


    

In [ ]:

    

def fakecubefile(xyzarray,edges,filename,norm=True,log=False):
    xav=0.5*(edges[0][0]+edges[0][-1])
    yav=0.5*(edges[1][0]+edges[1][-1])
    zav=0.5*(edges[2][0]+edges[2][-1])
    xmin=0.5*(edges[0][0]-edges[0][1])+edges[0][0]
    xstep=np.absolute(edges[0][0]-edges[0][1])
    xsteps=len(edges[0])-1
    ymin=0.5*(edges[1][0]-edges[1][1])+edges[1][0]
    ystep=np.absolute(edges[1][0]-edges[1][1])
    ysteps=len(edges[1])-1
    zmin=0.5*(edges[2][0]-edges[2][1])+edges[2][0]
    zstep=np.absolute(edges[2][0]-edges[2][1])
    zsteps=len(edges[2])-1
    if norm:
        normfactor=np.min([xstep,ystep,zstep])
        print normfactor
        xmin=xmin/normfactor
        ymin=ymin/normfactor
        zmin=zmin/normfactor
        xstep=xstep/normfactor
        ystep=ystep/normfactor
        zstep=zstep/normfactor
        xav=xav/normfactor
        yav=yav/normfactor
        zav=zav/normfactor
    with open(filename,'w') as f:
        f.write("Mobility plot faked by python\n")
        f.write("Created by Jens\n")
        f.write("1 {} {} {}\n".format(xmin,ymin,zmin))
        f.write("{} {} 0 0\n".format(xsteps,xstep))
        f.write("{} 0 {} 0 \n".format(ysteps,ystep))
        f.write("{} 0 0 {}\n".format(zsteps,zstep))
        f.write("55 0 {} {} {}\n".format(xav,yav,zav))
        N=0
        for i in range(xyzarray.shape[0]):
            for j in range(xyzarray.shape[1]):
                 for k in range(xyzarray.shape[2]):
                    value=xyzarray[i,j,k]
                    if log and value>0:
                        value=np.log(value)
                    N+=1  
                    if N==6:
                        f.write("{:E}\n".format(value))
                        N=0
                    else:
                        f.write("{:E} ".format(value))


    

In [75]:

    

fakecubefile(hist,edges,'test.cub',log=True)


    

    




9.04791e-09

In [ ]: