#-------------------------------------------------------------------------
#This script uses the "vanilla" executable to extract data from
#the Mars Global Surveyor TES raw data files. It subsets according
#to estimated target surface temperature, and writes out the
#results in a convenient form for further analysis. The raw data
#set from which the processed data in the parent directory was
#created is not included here, but can be obtained from the NASA
#Planetary Data System.
#
#With thanks to Kat Scanlon for helping to figure out how to use vanilla
#to extract data.
#--------------------------------------------------------------------------

#Example of use of vanilla to extract spectral
#data from sensor 1, with target planet surface
#temperature between 240K and 241K:
#
#vanilla /Volumes/MGST_0156/data/mars -fields "target_temperature detector calibrated_radiance[]" -select "target_temperature 260 261 detector 1 1" > slask.out
#
#
#ToDo:
#         *How do we know if data is single scan or double scan?
#          This affects the wavenumbers. wavnumss.tab is the double-scan file.
#          The double scan wavenumbers differ a lot, and these are mixed
#          in with the single-scan values
#
#          Also, why is length of waveNum file longer than spectral data?

from ClimateUtilities import *
import phys

#This allows us to execute system commands. It is
#how we call the executable "vanilla" 
from os import system

#The following extracts data from a NASA PDS CD-ROM of the TES data
#Mount the CD-rom, and put in the mount directory and name of
#the CDROM below. To use a different CD, just change the CD name
#This is the mountpoint on Mac OSX. If you have downloaded a data
#volumen from the NASA PDS, just replace the mountpoint with the
#directory containing the data.
mountPoint = '/Volumes/MGST_0156/'
dataDir = mountPoint+'data/mars'
#Specify the target temperature range. You can edit
#the command string to select on different criteria
#You can make different subsets if you want, as long as
#you write out the same selection of data variables as used here.
#The getData function takes the desired surface temperature as an argument
#and extracts all data meeting the criterion. For example, type
#getData(220.) to extract all spectra with a surface temperature
#between 220K and 221K.  You can easily rewrite the filter in
#getData to do other forms of subsetting if you want. You can
#follow this example to extract types of data other than spectra
#(e.g. if you want to make plots of surface temperature vs. latitude
#and longitude).

#vanilla_options =' -fields  -select ')
def getData(Tsurf):
    outname = 'TES%d.out'%Tsurf
    #
    #Select the output data to extract. Note that if you edit this line
    #to output different metadata, you'll have to change the file
    #read and line parsing in the the extract(...) function below
    #
    output_selection = ' "target_temperature latitude longitude sub_solar_latitude sub_solar_longitude scan_length detector calibrated_radiance[]" '
    #
    #Subset the data according to criteria you choose
    data_filter = ' "target_temperature %d %d detector 1 1 scan_length 1 1" > %s'%(Tsurf,Tsurf+1,outname)
    system('./vanilla '+dataDir+ ' -fields '+ output_selection + ' -select '+ data_filter)



##def extract(line):
##	line = line.split('\t')[7:]
##	data = Numeric.array([float(item) for item in line])
##	return data
##
##def Planck(wave,T):
##    return 100.*phys.c*phys.B(100.*phys.c*wave,T)
##
##detector = 1 #Assumes we read single scan, first detector
##waveFile = open("wavnumss.tab")
##lines = waveFile.readlines()
##array = [float(line.split(',')[detector]) for line in lines]
##waveNums = Numeric.array(array)[:-5]
##
##print array
##waveFile.close()
##
##
##getData()#Get the data from the CD-ROM. Comment out if you already
##         #have the data and just specify outname
##
##specFile = open(outname)
##lines = specFile.readlines()[1:] #First record is the header
##
###Search the spectra and find the one with the most net irradiance.
###This is the clearest sky case. Also check for missing data
###at end of spectrum.
##\
###Get metadata:
##TargetTList = [float(line.split('\t')[0]) for line in lines]
##Lats = [float(line.split('\t')[1]) for line in lines]
##Lons = [float(line.split('\t')[2]) for line in lines]
##SunLats = [float(line.split('\t')[3]) for line in lines]
##SunLons = [float(line.split('\t')[4]) for line in lines]
##scanType =  [int(line.split('\t')[5]) for line in lines]
###
##BadDatList =  [float(line.split('\t')[-1]) for line in lines]
##TotRadList = []
##for line in lines:
##    TotRadList.append(sum(extract(line)))
##print max(TotRadList),min(TotRadList)
##imax = Numeric.argmax(TotRadList)
##
###Function to show deviation from blackbody
##def BBfit(i):
##    T = TargetTList[i]
##    data = extract(lines[i])
##    P = Numeric.array([Planck(w,T) for w in waveNums])
##    norm = sum(data[90:])/sum(P[90:])
##    P = P*norm #Normalize to high-freq tail
##    return (sum((P-data)**2)/len(P))**.5/sum(P/len(P))
##
##BBfitL = [BBfit(i) for i in range(len(lines))]
##
##def doPlot(i):
##    print 'Lat Lon',Lats[i],Lons[i]
##    print 'Sunlat,Sunlon',SunLats[i],SunLons[i]
##    data = extract(lines[i])
##    T = TargetTList[i]
##    P = Numeric.array([Planck(w,T) for w in waveNums])
##    norm = sum(data[90:])/sum(P[90:])
##    P = P*norm #Normalize to high-freq tail
##    c = Curve()
##    c.addCurve(waveNums)
##    c.addCurve(data)
##    c.addCurve(P)
##    return c
##
##