# calculates Preservation Index from t and rh records # requires python interpreter (www.python.org) # # command line example: # $>python twpi.py 5 6 datafile.txt # twpi.py is this script # temperature is in column 5 and rh in column 6 # NOTE: date-time formats count as several columns because: # a column is delimited by / \ : ; , space, tab or newline # you may need to experiment because # dates with backslash cause parsing problems # temperature must be in celsius, rh in %, decimal point . (not comma) # comment lines in data file must start with # # the data file must not contain empty rows # the file name must not contain spaces # the file path should have forward slashes: path/to/file/name # this works from windows98 and later # Output format: # output file has name datafile.ipipi # it contains the original data with PI and 1/PI added to the end of each row # note that the calculation assumes equal time intervals between each # row in the data set. It does not use the date-time recorded in the data. # # programmer: Tim Padfield # this program is open source and may freely be copied and altered # but please describe your alterations in this comment section. import sys, os, re, math # rewrite calculations as you wish! def calcpi(t, rh, method): try: t += 273 # celsius to kelvin # IPI formula for preservation index if method == 'IPI': years = math.exp((95220 - 134.9 * rh)/(8.314*t) + (0.0284 * rh) - 28.023)/360 # Tim Padfield's arrhenius/mass action approximation if method == 'TP': years = 1.0/(rh * 5.9 * (10**12) * math.exp(-90300/(8.314 * t))) rate = 1.0/years years = str(int(years)) return rate, years except: print 'calculation out of range' sys.exit(1) # Regular expression to slice data into number tokens dataslicer = re.compile(r'[\\/:;, \t\n]*') print '\nProgram calculates PI, TWPI from climate data file\n' print 'type -h on command line for help\n' if len(sys.argv) > 1: if sys.argv[1] == '-h': # write help and stop print 'COMMAND LINE:' print 'python twpi.py ' print 'temperature column number ' print 'RH column number ' print 'data file name \n' print 'EXAMPLE:' print 'prompt>python twpi.py 5 6 climate.txt\n' print 'The command line can be composed interactively:' print 'just type:' print 'prompt> python twpi.py\n' print 'NOTES:' print 'date-time formats count as several columns because:' print 'a column is delimited by / \ : ; , space, tab or newline' print 'temperature must be in celsius, rh in %' print 'decimal point must be . (not comma)' print "filenames with spaces must be quoted: 'file with space.dat'" print 'comment lines in the data file should start with #' print 'otherwise, by chance, data columns may be recognised' print 'in the comment text and processed,' print 'causing obvious, or subtle errors in the result' print '\n' print 'OUTPUT FORMAT:' print 'output file has data file root with extension .ipipi' print 'it contains the original data rows with these numbers added:' print 'mass-action/Arrh rate, PI as rate, ' print ' m-a/A as lifetime, PI' print 'At the end of the file are two values of the TWPI\n' sys.exit(1) if len(sys.argv) < 4: print 'needs two column numbers and one input file' success = 0 while success == 0: try: print 'temperature column: ', tcol = raw_input() success = 1 except: print 'Something wrong, try again' while success == 1: try: print '\nRH column: ', rhcol = raw_input() success = 2 except: print 'Something wrong, try again' while success == 2: try: print '\ndata file: ', dfilename = raw_input() dfile = open(dfilename, 'r') success = 3 except: print 'something wrong, try again' dfilebits = dfilename.split(".") dfileplus = dfilebits[0] + ".ipipi" ofile = open(dfileplus,'w') print 'Writing results to ' , print dfileplus else: try: tcol = sys.argv[1] rhcol = sys.argv[2] dfilename = sys.argv[3] dfile = open(dfilename, 'r') dfilebits = dfilename.split(".") dfileplus = dfilebits[0] + ".ipipi" ofile = open(dfileplus,'w') print '\nWriting arrh-rate, 1/PI and PI to end of each line in file: ', print dfileplus ofile.write('added to end of each row:\n') ofile.write('arrh-rate, IPI-rate, arrh-PI, IPI-PI\n') ofile.write('calculated from t in column '+sys.argv[1]) ofile.write(' and RH in column '+sys.argv[2]+'\n') except: print 'Unable to read files or make output file' print tcol,rhcol,dfilename,dfileplus sys.exit(1) #start reading and processing # subtract one from column count to make first col 0 tcol = int(tcol)-1 rhcol = int(rhcol)-1 linecount = 0 rate_accum_tp = 0 rate_accum_ipi = 0 while 1: try: oneline = dfile.readline() if oneline == '': # EOF break if oneline[0] == '#': #comment ofile.write(oneline) continue if len(oneline) < 2: # blank line ofile.write(oneline) continue splitline = dataslicer.split(oneline) try: t = float(splitline[tcol]) rh = float(splitline[rhcol]) except(IndexError, ArithmeticError,ValueError), e: ofile.write(oneline) print 'error at line: ',linecount continue result_tp = calcpi(t,rh,'TP') rate_tp = result_tp[0] rate_accum_tp += rate_tp rate_tp = str(rate_tp) rate_tp = rate_tp[0:6] result_ipi = calcpi(t,rh,'IPI') rate_ipi = result_ipi[0] rate_accum_ipi += rate_ipi rate_ipi = str(rate_ipi) rate_ipi = rate_ipi[0:6] oneline = oneline[:-1]+' '+rate_tp+' '+rate_ipi+' '\ + result_tp[1]+' '+result_ipi[1]+'\n' ofile.write(oneline) linecount += 1 except EOFError, e: print 'end of file' break except (IndexError, ArithmeticError), e: print 'Arithmetic error: ', e, print 'at line: ', linecount #sys.exit(1) ofile.write('#\n') if rate_accum_tp > 0: lifetime_tp = str(linecount/rate_accum_tp) lifetime_tp = lifetime_tp.split('.') lifetime_tp = lifetime_tp[0] ofile.write('#lifetime from mass-action/Arrh: '\ + lifetime_tp + ' years from ' + str(linecount)\ + ' measuring periods\n' ) print '#lifetime from mass-action/Arrh: '\ + lifetime_tp + ' years from ' + str(linecount)\ + ' measuring periods' if rate_accum_ipi > 0: lifetime_ipi = str(linecount/rate_accum_ipi) lifetime_ipi = lifetime_ipi.split('.') lifetime_ipi = lifetime_ipi[0] ofile.write('#TWPI: ' + lifetime_ipi + ' years from '\ + str(linecount) + ' measuring periods\n' ) print '#TWPI: ' + lifetime_ipi + ' years from '\ + str(linecount) + ' measuring periods' ofile.close() dfile.close() sys.exit(1)