#! /usr/bin/python
# input letter name ph lowest and highest values, step and pkafile  as command line args
# calculates what fraction of molecules are protonated
# based on input file

import sys

lettername = sys.argv[1]
ph1 = float(sys.argv[2])
ph2 = float(sys.argv[3])
phstep = float(sys.argv[4])

try: 
	filename = open(sys.argv[5])	
except:
	print "Could not open file", sys.argv[5]
	sys.exit()

pkadictionary = {}
for line in filename.readlines():
	fields = line.split()
	pkadictionary[fields[0]]= float(fields[1])

try:
	pka = pkadictionary[lettername]
except:
	print "No pKa found for", lettername
	sys.exit()

def function(ph):
	return (1/((10**(ph-pka))+1))

testhalf = 0.5 - function(ph1)
while ph1 < ph2:
	fprot = function(ph1)
	print ph1, fprot
	if testhalf * (0.5-fprot) < 0:
		smallstep = phstep/10
		r = ph1 - 2*phstep
		m = function(r)
		e = 0.5 - m
		while r < ph1:
			m = function(r)
			if e * (0.5 - m) < 0:
				y = r	
				break
			r = r + smallstep
			e = 0.5 - m		
	ph1 = ph1 + phstep	
	testhalf = 0.5-fprot

#halfprotonated = 1/((10**(y-pka))+1)
# print y, halfprotonated