from ClimateUtilities import *
from string import atof
from glob import glob
import Numeric
import phys
a = phys.CO2.R/phys.CO2.cp

dataPrefix = '/Users/rtp1/bigdata/MarsTPS/'

MarsRadioHeader = ['radius','lat','lon','geopotential','p','sigp','T','sigT','n','sign']
#Dictionary of indices of useful stuff in PDS header record
hdict = {}
hdict['lat'] = 7
hdict['lon'] = 9
hdict['latsun'] = 11
hdict['Ls'] = 13 #Longitude of the sun. A measure of the season
hdict['psurf'] = 16 #Surface pressure
hdict['solarTime'] = 20 #Local true solar time, hours
hdict['sza'] = 21 #Solar zenith angle

filenames = glob(dataPrefix+'*.TPS')
for filename in filenames:
    #Get information from PDS header table
    f = open(filename)
    h = f.readline().split(',')
    lat = atof(h[hdict['lat']])
    lon = atof(h[hdict['lon']])
    latsun = atof(h[hdict['latsun']])
    solarTime = atof(h[hdict['solarTime']])
    sza = atof(h[hdict['sza']])
    psurf = atof(h[hdict['psurf']])
    Ls = atof(h[hdict['Ls']])
    #print filename,lat,lon,latsun,solarTime,sza
    c = readTable(filename,MarsRadioHeader,',')
    #Subset data for output
    c1 = Curve()
    c1.description =\
        'lat = %f,lon=%f,latsun = %f,solarTime = %f,sza = %f,ps=%f,Ls=%f'%(lat,lon,latsun,solarTime,sza,psurf,Ls)
    c1.addCurve(c['p'],'p','Pressure (Pa)')
    c1.addCurve(c['T'],'T','Temperature (K)')
    Tmax = max(c1['T'])
    #Modify the if clause to select the kind of data you want
    if (latsun > 15) and Tmax > 220. :
        name = filename.split('/')[-1]
        print name,psurf,Tmax,solarTime,lat,latsun,sza
        c1.dump(name.split('.')[0]+'.txt')