glb.py

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random
import itertools


colors = ('b', 'g', 'r')

A = 208
rmax = 15
a = 1.
r0 = 1.25
R = r0*A**.3333

b = 8



N1 = np.empty((A,2))
N2 = np.empty((A,2))

for k,N in ((-1,N1),(1,N2)):
    for i in range(A):
        while True:
            r = rmax*random()
            if random() < 1/(1 + np.exp(r-R/a)):
                phi = 2*np.pi*random()
                N[i] = [k*b/2 + r*np.cos(phi),r*np.sin(phi)]
                break
            else:
                continue

BC = []
for n1,n2 in itertools.product(N1,N2):
    if np.sum(np.square(n1-n2)) < 1.:
        BC.append(np.mean([n1,n2],axis=0))
BC = np.array(BC)

fig = plt.figure(frameon=False)
ax = fig.add_axes([0,0,1,1])
ax.axis('off')

ax.scatter(*N1.T,s=30,facecolors='none',edgecolors=colors[0],linewidths=.6,label='Nucleus 1')
ax.scatter(*N2.T,s=30,facecolors='none',edgecolors=colors[1],linewidths=.6,label='Nucleus 2')
ax.scatter(*BC.T,s=30,facecolors=colors[2],alpha=.5,linewidths=.4,label='Overlap')

plt.xlim(-rmax,rmax)
plt.ylim(-rmax,rmax)
#plt.gca().xaxis.set_minor_locator(ticker.AutoMinorLocator(2))
#plt.gca().yaxis.set_minor_locator(ticker.AutoMinorLocator(2))

#plt.xlabel('x [fm]')
#plt.ylabel('y [fm]')
#plt.title('$b = {}$ fm'.format(b))
#plt.legend()

#plt.tight_layout(pad=0)
#plt.savefig('toyglauber{}.pdf'.format(b))
plt.show()