#! /usr/bin/python

# MeasureTau.py
# program to compute optical depth (tau) and Volume filling density 
# for a particle file generated by BEAM setup code.

import numpy as np
import matplotlib.pyplot as plt
import sys
import os
from math import *


inputFile = sys.argv[1]
print inputFile

# make some lists

radius = []
ParticleArea = []
ParticleVolume = []
xcoord = []
ycoord = []
zcoord = []

pi = 4*atan(1.0)

f = open(inputFile,'r')

# skip header but keep in buffers 
header1 = f.readline()
header2 = f.readline()
header3 = f.readline()

# extract length of system from second line
str1 = header2.lstrip() # remove leading spaces
ispace = str1.index(' ') # find next space to isolate L/2
for line in f:
     line = line.strip()
     columns = line.split()
     radius.append(float(columns[0]))
     ParticleArea.append(pi*float(columns[0])**2)  
     ParticleVolume.append((4.0*pi/3.0)*float(columns[0])**3)  
     xcoord.append(float(columns[1]))
     ycoord.append(float(columns[2]))
     zcoord.append(float(columns[3]))


L = max(xcoord)-min(xcoord)
print 'L = ', L
print 'zmin, zmax = ',min(zcoord), max(zcoord)
h = min(abs(min(zcoord)), max(zcoord))
print 'h = ', h
hmin = np.mean(np.array(zcoord)) - 3*np.std(np.array(zcoord))
hmax = np.mean(np.array(zcoord)) + 3*np.std(np.array(zcoord))

print hmin,hmax



tau = sum(ParticleArea)/(L**2)
D = sum(ParticleVolume)/(h*L**2)

print 'tau = ', tau
print 'D = ', D

plt.hist(zcoord,bins=50)
plt.show()



