''' This is an example script creating subplots
of temperature and its standard deviation'''

__author__ = "Maarten J. Waterloo <maarten.waterloo@falw.vu.nl>"
__version__ = "1.0"

# Import some useful python modules
from scipy import array              # To use arrays
from pylab import *                  # Use matlab style commands
from os import getcwd,listdir,chdir  # For system commands
from meteolib import *               # Use meteorology functions

# Set working directory variable
#
# For UNIX
#workdir = '/mnt/master/courses/450060_soil_vegetation_atmosphere_exchange/svae_meteo/pyscripts'
workdir = '/media/usb0/svae_meteo/pyscripts'
# For MS Windows
#workdir = 'e:/courses/450060_soil_vegetation_atmosphere_exchange/svae_meteo/pyscripts'

# Change to workdir
chdir(workdir)
# Print working directory an give a newline (\n)
print 'Working directory set to: ',getcwd(),'\n'
print 'Files in workdir:\n',listdir(getcwd()),'\n'

# Activate LaTeX interpretation for figure texts
rc('text', usetex=True)

# Change default settings for figures
newdefaults = {'fontname':    'Arial',  # Use Arial font
               'backend':       'eps',  # Save figure as EPS file   
               'axes.labelsize':   14,  # Axis label size in points
               'text.fontsize':    14,  # Text size in points
               'legend.fontsize':  14,  # Legend label size in points
               'xtick.labelsize':  14,  # x-tick label size in points
               'ytick.labelsize':  14,  # y-tick label size in points
               'lines.markersize': 10,  # markersize, in points
               'lines.linewidth':   1.5 # line width in points
               }
rcParams.update(newdefaults)

# Create array with temperature data
t_data = array([[2007,121,1435,22.95,.287,.536],\
                [2007,121,1440,22.45,.188,.434],\
                [2007,121,1445,21.84,.207,.455],\
                [2007,121,1450,21.90,.265,.514],\
                [2007,121,1455,22.34,.476,.69]])

# Create some variables from t_data columns
T = t_data[:,3]    # Put temperature values in variable T
stdT = t_data[:,5] # Put st.dev. temperature values in stdT
time = t_data[:,2] # Put time in variable time

# Get coefficients for line through T and stdT (stdT = aT + b)
a, b = polyfit(T, stdT, 1)
print 'Linear fit of T versus stdT'
print 'Multiplier a= ',a   # multiplier 
print 'Intercept b= ',b,'\n'  # intercept

# Calculate some statistics of temperature
Tavg = average(T) # Average
Tmed = median(T)  # Median
std = std(T)      # Standard deviation
Tsum = sum(T)     # Sum

# Calculate the latent heat of vapourisation (T must be in Celsius)
L = L_calc(T)
print 'Latent heat of vapourisation from L_calc()\n',L

# Use the find command to find data that satisfies a certain condition
#year = find(t_data==2007)   # Find the indices for the year
tmax = find(T==max(T))      # Find the index of the maximum temperature
tmin = find(T==min(T))      # Find the index of the minimum temperature

# Now create figure 1
figure(1)     # Define figure 1
plot(T,stdT,'o',label='data') # Plot stdT versus T, use round marker
plot(T,a*T+b,label='fitted')  # plot line using fitting coefficients
xlabel(r'${\rm Temperature~[}^\circ {\rm C]}$')
ylabel(r'${\rm St.~dev.~Temperature~[}^\circ {\rm C]}$')
legend()

# Now create figure 2
figure(2)     # Define figure 2
subplot(211)  # upper plot from 2 plots
plot(time,T, label='Air Temperature')
ylabel(r'${\rm Temperature~[}^\circ {\rm C]}$')
legend(loc=1)
subplot(212)  # lower plot from 2 plots
# Plot line (-), red colour (r) and + marker (+)
plot(time,stdT, '-r+', label='St. dev. T')
ylabel(r'${\rm St.~dev.~Temperature~[}^\circ {\rm C]}$')
xlabel('Time')
legend(loc=2)
savefig('duotemp.eps',dpi=300)

# Create figure 3, now with two y-axes
figure(3)
xlabel(r'${\rm Time}$')
ylabel(r'${\rm Temperature~[}^\circ {\rm C]}$')
plot(time,T,'ro',label='Temperature') # Plot stdT versus T, use round marker
legend(loc=1)
twinx()
plot(time,stdT,label='St.dev')  # plot line using fitting coefficients
ylabel(r'${\rm St.~dev.~Temperature~[}^\circ {\rm C]}$')
legend(loc=2)
savefig('twoaxes.eps',dpi=300)

print '\nType show() in command window to see figures...'