#-------------------------------------------------------------------
#Reads and plots infrared spectra from the Mars Global Surveyor
#TES instrument. This script reads processed spectral data files
#that were extracted using the "vanilla" program.  A small subset
#of preprocessed files are included in the Chapter 4 data directory.
#These particular files are subsetted according to the TES estimate
#of the surface temperature of the footprint being looked at
#(e.g. TES220.out contains data with estimated surface temperature
#beween 220K and 221K) but this program will work with any other
#subset criteria, so long as the same set of variables (temperature,
#spectra, etc) are selected for output.
#
#For information on "vanilla" and how to use it to extract data of
#this type, and other Mars data, from TES files you can download
#from the NASA Planetary Data System, see the script processTES.py
#in the Original subdirectory of the MarsTES data directory.
#
#This script makes use of a simple trick for comparing the TES spectra
#with the blackbody Planck curve.  The spectra are the irradiances measured
#by the satellite in orbit, in W/m**2/steradian. To compare these with
#what the instrument would see if the surface blackbody emission were
#unimpeded by an atmosphere, one would have to really do a somewhat
#complicated navigational calculation to determine the distance from
#the satellite from the target point, the angular resolution of the
#instrument, and so forth.  That's what a real professional would do,
#but here, instead we use the radiatively estimated surface temperature (recorded
#in the metadata), and then normalize the Planck curve based on the
#near-infrared portion of the measured spectrum, which is relatively
#unaffected by the Mars atmosphere.  This process can be fooled sometimes
#when there is missing data on the shortwave end of the spectrum.
#
#The main function you will want to use is doPlot(i), which returns
#a Curve object giving the spectrum of the ith' observation in
#the file.  This Curve also contains the fitted Planck function for
#comparison. 
#-------------------------------------------------------------------

#Data on section of text which this script is associated with
Chapter = '4'
Section = 'section:RealGasOLROverview'
Figure = 'fig:MarsTES'

from ClimateUtilities import *
import phys
from glob import glob #Used to find what data files are available

#Path to Workbook datasets
datapath = '/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/'

def extract(line):
	line = line.split('\t')[7:]
	data = numpy.array([float(item) for item in line])
	return data

def Planck(wave,T):
    return 100.*phys.c*phys.B(100.*phys.c*wave,T)

#Function to show deviation from blackbody
def BBfit(i):
    T = TargetTList[i]
    data = extract(lines[i])
    P = numpy.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))


#Get the wavenumber data
detector = 1 #Assumes we read single scan, first detector
waveFile = open(datapath+'Chapter4Data/MarsTES/'+'wavnumss.tab')
lines = waveFile.readlines()
array = [float(line.split(',')[detector]) for line in lines]
waveNums = numpy.array(array)[:-5]
waveFile.close()

#Get spectral data from the preprocessed TES spectra files
#
#First get a list of what's available
specFileNames = glob(datapath+'Chapter4Data/MarsTES/TES*.out')

print "Available data files are:"
for i in range(len(specFileNames)):
    print i,specFileNames[i].split('/')[-1] #Prints out just filename

j = raw_input("Give the number of the file you want")
specFileName = specFileNames[int(j)]

#Open the file we want and read in the data. Each
#data line is one spectrum plus metadata
specFile = open(specFileName)
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]
#
#A few diagnostics to help us sort out what we want to plot
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 = numpy.argmax(TotRadList)
BBfitL = [BBfit(i) for i in range(len(lines))]

#Function to plot the i'th spectrum in the file.
#Returns a Curve object for plotting or writing out.
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 = numpy.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,'Wavenumber')
    c.addCurve(data,'TESIrradiance')
    c.addCurve(P,'Planck')
    return c