#------------------------------------------------------
#Processes data from the Zachos2001 compilation of
#18-O a d 13C, and produces output file suitable for
#plotting.  The original dataset contains data from
#a variety of different forams, and the function extract(...)
#gives the user the choice of which ones to aggregate together
#into the dataset.  It also does smoothing of the data
#(controlled by nSmooth). The argument index controls
#which data column (O18 or C13) is fetched.
#This script is meant to be used interactively.
#
#Available foram genera are:
#    CIB , PWUELL, CIBKULL, PWUELCIB, NUT, ARA
#
#Typical usage:
#       c = extract(['CIB','PWUELL','NUT'],O18index)
#then plot or save
#------------------------------------------------------

from ClimateUtilities import *
import string

fileName = 'Zachos2001.txt'

#Because this file contains text data in one of the columns,
#we can't just do the input using readTable, which is only
#designed for numerical data. We have to use low level input

#Open the file and read it in
infile = open(fileName)
lines = infile.readlines()

#Function to process the lines, and extract just the 18-O data
#for the kind of foram you want

O18index = 3
C13index = 4
def extract(foramList,index = 3,nSmooth = 5):
    ageList = []
    datList = []
    datListSmooth = []
    start = True
    for line in lines:
        data = line.split()
        #This try/except block finds the first line of data
        try:
            age = string.atof(data[1])
            genus = data[2]
            dat = string.atof(data[index])
        except:
            continue
        
        if genus in foramList:
            if start:
                start = False
                datSmooth = dat
            else:
                datSmooth = (1.-1./nSmooth)*datSmooth + dat/nSmooth
            ageList.append(age)
            datList.append(dat)
            datListSmooth.append(datSmooth)
    c = Curve()
    c.addCurve(ageList,'age')
    c.addCurve(datList,'data')
    c.addCurve(datListSmooth,'Smoothed')
    c.reverseX = True
    return(c)