post processing of the results of SRIM COLLISON.txt files

srim-pp.py

import re
import string
import math

IN = 'COLLISON1.txt'
DELIMER = '\xb3'
KeV_in_Mev=1000


data_line_re = re.compile(".+\d+\.E\+0.*")
dead_line_re = re.compile(".*(Target|Replacements).*")

def get_data_lines(raw_line):
  return data_line_re.match(raw_line) and not dead_line_re.match(raw_line)

def line_to_record(data_line):

	rf = string.split(data_line, DELIMER)
	record = {}
	record['ion'] = rf[1]
	record['energy'] = float(rf[2])/KeV_in_Mev
	record['depth'] = float(rf[3])
	record['y'] = float(rf[4])
	record['z'] = float(rf[5])
	record['se'] = float(rf[6])
	record['atom'] = string.strip(rf[7])
	record['recoil_energy'] = float(rf[8])
	record['target_disp'] = float(rf[9])

	return record

def read_records(data_file):
	f = open(data_file)
	data_lines_only = filter(get_data_lines, f.readlines())
	records = map(line_to_record, data_lines_only)
        
	f.close()

	return sorted (records, key = lambda record: record['depth'] )


if __name__ == '__main__':

	last_level = None
	level_index = -1
	levels = []
	for r in read_records(IN):
		atom = r['atom']
		if atom != last_level:
			level_index += 1
			last_level = atom
			levels.append({})
			levels[level_index]['atom'] = atom
			levels[level_index]['avg'] = 0.0
			levels[level_index]['sigma'] = 0.0
			levels[level_index]['count'] = 0.0
			levels[level_index]['layer'] = level_index
			levels[level_index]['depth0'] = 1 << 63
			levels[level_index]['depth1'] = 0

		run_info = levels[level_index]
		current_avg = run_info['avg']
		current_sigma = run_info['sigma']
		current_count = run_info['count']

		if run_info['depth0'] > r['depth']:
			run_info['depth0'] = r['depth']
		if run_info['depth1'] < r['depth']:
			run_info['depth1'] = r['depth']

		run_info['avg'] = (current_avg*current_count + r['energy']) / (current_count + 1)
		run_info['sigma'] = (current_sigma*current_count + r['energy']*r['energy']) / (current_count + 1)
		run_info['count'] = current_count + 1

	final_layer = 0
	for l in levels:
		if l['depth1'] - l['depth0'] == 0:
			continue
		final_layer += 1
		sigma = math.sqrt(l['sigma'] - l['avg']*l['avg'])
		print "layer %2d %3s %.2f %.2f" % (l['layer'], l['atom'], l['avg'], sigma)