%matplotlib inline

import math
from numpy import *
import matplotlib.pyplot as plt
from matplotlib import cm
from pylab import *
from operator import itemgetter

import hygrometry

def frange(x, y, jump):
    """Like range(), but works with floats."""
    x=float(x)
    y=float(y)
    jump = float(jump)
    while x < y:
        yield x
        x += jump


def plot2derrors(all_maxerrors):
    mmap=[{}, {}, {}]
    
    for i in range(0, len(all_maxerrors[0])):
        for j in [0, 1, 2]:
            key="%.1f" % (all_maxerrors[j][i],)
            if key in mmap[j]:
                if mmap[j][key][1] < all_maxerrors[3][i]:
                    mmap[j][key] = [all_maxerrors[j][i], all_maxerrors[3][i]]
            else:
                mmap[j][key] = [all_maxerrors[j][i], all_maxerrors[3][i]]
    
    xl=['T', 'H', 'P']        
    for j in [0, 1, 2]:
        X=[]
        Y=[]
        for key in mmap[j]:
            X.append(mmap[j][key][0])
            Y.append(mmap[j][key][1])
        X, Y = [list(x) for x in zip(*sorted(zip(X, Y), key=lambda pair: pair[0]))]    
        figure()
        
        plot(X, Y, 'r')
        xlabel(xl[j])
        ylabel('max(E)')
        show()


def calc_maxerrors(wet_bulb_calc_func, sensor_range, sensor_errors, steps, pressure_sensor=True, pressure_constant=977):
    tidx=0
    rhidx=1
    pidx=2
        
    merr = 0
    minfo = []

    all_maxerrors = [[], [], [], []]
    count = 0.0
    mserr = 0.0
    for Tc in frange(sensor_range[tidx][0], sensor_range[tidx][1], steps[tidx]):
        for P in frange(sensor_range[pidx][0], sensor_range[pidx][1], steps[pidx]):
            for RH in frange(sensor_range[rhidx][0], sensor_range[rhidx][1], steps[rhidx]):
                sterr = 0
                Ptruth=P
                if not pressure_sensor:
                    Ptruth=pressure_constant
                st = [Tc, RH, P]
                truth = wet_bulb_calc_func(Tc, RH, Ptruth)
                
                eTc=sensor_errors[tidx](Tc)
                eRH=sensor_errors[rhidx](RH)
                eP=sensor_errors[pidx](P)

                
                combos = [[eTc, 0, 0], [-eTc, 0, 0], \
                          [0, eRH, 0], [0, -eRH, 0], \
                          [0, 0, eP], [0, 0, -eP], \
                          [eTc, eRH, 0], [eTc, -eRH, 0], [-eTc, eRH, 0], [-eTc, -eRH, 0], \
                          [eTc, 0, eP], [eTc, 0, -eP], [-eTc, 0, eP], [-eTc, 0, -eP], \
                          [0, eRH, eP], [0, eRH, -eP], [0, -eRH, eP], [0, -eRH, -eP], \
                          [eTc, eRH, eP], [eTc, -eRH, eP], [eTc, -eRH, -eP],\
                          [-eTc, eRH, eP], [-eTc, -eRH, eP], [-eTc, -eRH, -eP]]
            
                for c in combos:
                    tmp = wet_bulb_calc_func(Tc+c[0], RH+c[1], P+c[2])
                    nerr = math.fabs(truth-tmp)
                    mserr += nerr
                    count += 1.0
                    if nerr > merr:
                        merr = nerr
                        minfo = [merr, [Tc, RH, P], c]
                    if nerr > sterr:
                        sterr = nerr
                all_maxerrors[0].append(Tc)
                all_maxerrors[1].append(RH)
                all_maxerrors[2].append(P)
                all_maxerrors[3].append(sterr)
    mserr=mserr/count
    minfo.append(mserr)

    return [minfo, all_maxerrors]
    

def calc_uncertainty(wet_bulb_calc_func, sensor_range, sensor_errors, steps, pressure_sensor=True, pressure_constant=977):
    tidx=0
    rhidx=1
    pidx=2
    
    
    merr = 0
    minfo = []
    h=0.01
    h_T=h
    h_RH=h
    h_P=h
    all_maxerrors = [[], [], [], []]
    for Tc in frange(sensor_range[tidx][0], sensor_range[tidx][1], steps[tidx]):
        for RH in frange(sensor_range[rhidx][0], sensor_range[rhidx][1], steps[rhidx]):
            for P in frange(sensor_range[pidx][0], sensor_range[pidx][1], steps[pidx]):
                if not pressure_sensor:
                    eP=math.fabs(P-pressure_constant)
                #eP=0
                eTc=sensor_errors[tidx](Tc)
                eRH=sensor_errors[rhidx](RH)
                eP=sensor_errors[pidx](P)

                
                dTc_dF=(wet_bulb_calc_func(Tc+h_T, RH, P) - wet_bulb_calc_func(Tc-h_T, RH, P))/(2.0*h_T)
                dRH_dF=(wet_bulb_calc_func(Tc, RH+h_RH, P) - wet_bulb_calc_func(Tc, RH-h_RH, P))/(2.0*h_RH)
                dP_dF=(wet_bulb_calc_func(Tc, RH, P+h_P) - wet_bulb_calc_func(Tc, RH, P-h_P))/(2.0*h_P)
                dF = math.sqrt( (dTc_dF*eTc)**2.0 + (dRH_dF*eRH)**2.0 + (dP_dF*eP)**2.0 )
                
                if dF > merr:
                    merr = dF
                    minfo = [merr, [Tc, RH, P]]
                all_maxerrors[0].append(Tc)
                all_maxerrors[1].append(RH)
                all_maxerrors[2].append(P)
                all_maxerrors[3].append(dF)
    return [minfo, all_maxerrors]
    

def batch_WB(wet_bulb_func, range_t, range_rh, range_p):
    ret = []
    for Tc in frange(range_t[0], range_t[1], range_t[2]):
        for RH in frange(range_rh[0], range_rh[1], range_rh[2]):
            for P in frange(range_p[0], range_p[1], range_p[2]):
                value = wet_bulb_func(Tc, RH, P)
                ret.append([Tc, RH, P, value])
    return ret


"""
Explanation: Wet Bulb Calculation Analysis
Common/Helper Methods
End of explanation
"""


range_Tc = [hygrometry.conv_f2c(85.), hygrometry.conv_f2c(110.)]
range_P = [904., 1050.]
range_RH = [20., 100.]
srange=[range_Tc, range_RH, range_P]

"""
Explanation: Ranges and Sensor Accuracy
Analysis Assumptions:
* Temperatures range from 85F - 110F
* Pressures range in the from 904 hPa - 1050 hPa
* RH in from 20% - 100%
Example sensor to show error propgation will be the SHT-75.
End of explanation
"""



#These errors are from SHT75 datasheet, Figure 2 and Figure 3
def error_Tc(Tc):
    #-40C - 10
    if Tc <= 10:
        return (-1.1/50.)*(Tc-10) + 0.4
    #10C - 25C
    elif Tc <= 25:
        return (-0.1/15.)*(Tc-25) + 0.3
    #25C - 40C
    elif Tc <= 40:
        return (0.1/15.)*(Tc-25) + 0.3
    #40C - 100C
    elif Tc <= 100:
        return (1.1/50.)*(Tc-40) + 0.4
    else:
        return Inf
#These errors are from SHT75 datasheet
def error_RH(RH):
    if RH >= 10.0 and RH <= 90:
        return 1.8
    elif RH < 10.0:
        return (-2.2/10.0)*(RH-10.0) + 1.8
    else:
        return (2.2/10.0)*(RH-90.0) + 1.8

#TODO: Look at pressure datasheet and our calculations some more.        
def error_P(P):
    return 2.7

#Plot sensor errors
figure()
T = []
DT = []
for t in frange(-40., 100., 1):
    T.append(t)
    DT.append(error_Tc(t))
plot(T, DT, 'r')
xlabel('T (C)')
ylabel('DT (C)')
title('SHT-75 Temperature Error')
show()

figure()
RH = []
DRH = []
for rh in frange(0., 101., 1):
    RH.append(rh)
    DRH.append(error_RH(rh))
plot(RH, DRH, 'r')
xlabel('RH (%RH)')
ylabel('DRH (%RH)')
title('SHT-75 Humidity Error')
show()


figure()
P = []
DP = []
for p in frange(range_P[0], range_P[1], 1):
    P.append(p)
    DP.append(error_P(p))
plot(P, DP, 'r')
xlabel('P (hPa)')
ylabel('DP (hPa)')
title('Pressure Module Error')
show()


"""
Explanation: Sensor Errors
End of explanation
"""


print (hygrometry.wetbulb(20, 0, 985) - 5.85) < 0.01
print (hygrometry.wetbulb(20, 0, 985) - 20.) < 0.01
print math.fabs(hygrometry.wetbulb(20, 50, 985) - 13.8) < 0.01
print (hygrometry.wetbulb(30, 25, 985) - 16.92) < 0.01
print (hygrometry.wetbulb(40, 50, 3) - 27.67) < 0.01
print (hygrometry.wetbulb(41.3, 24.2, 351.6) - 20.5) < 0.01
print (hygrometry.wetbulb(22.9, 81.2, 476.1) - 20.12) < 0.01
print (hygrometry.wetbulb(41.2, 24, 900) - 23.83) < 0.01
print math.fabs(hygrometry.wetbulb(41, 12, 1000) - 20.08) < 0.01 

"""
Explanation: Tim Brice and Todd Hall Wet Bulb calculation

Original Javascript version: http://www.srh.noaa.gov/epz/?n=wxcalc_rh
Easy Calculation tries to explian it: https://www.easycalculation.com/weather/learn-dewpoint-wetbulb.php

vapor pressure, saturation vapor pressure, dew point:
$$ es = 6.112 \cdot e^{\frac{17.67 \cdot T_c}{T_c+243.5}} $$
$$ V = \frac{es \cdot RH}{100} $$
$$ T_d = \frac{243.5 \cdot \ln(V/6.112)}{17.67 - \ln(V/6.112)} $$
Then sent to some incremental loop function, calc_wb()
Error and Stability NOAA
Forward Stability
End of explanation
"""


step_Tc = 0.3
step_RH = 1.8
step_P = 2.7
ssteps=[step_Tc, step_RH, step_P]
info, errs = calc_maxerrors(hygrometry.wetbulb, srange, [error_Tc, error_RH, error_P], ssteps, True)
print info
print "Max error: %f" % (info[0],)
plot2derrors(errs)

"""
Explanation: Max Error Propagation
End of explanation
"""



print "      RH/Temperature Wet Bulb Error Table in Celsius"
print "T\RH| ",
for xh in frange(0., 101., 10.):
    print " %-3s " % ("%d" % (xh,),),
print
print "-"*(11*6)+"------",

frow = True
for y in frange(0., 51., 5.):
    print
    print "%-3s | " % ("%d" % (y,),),    
    for x in frange(0., 101., 10.):
        info, errs = calc_maxerrors(hygrometry.wetbulb, [[y, y+1], [x, x+1], [904.0, 1050.0]], [error_Tc, error_RH, error_P], [5, 10, 2.7], True)
        merr = array(errs[3]).max()
        print "%.2f " % (merr,),
print



print "      RH/Temperature Wet Bulb Error Table in Fahrenheit"
print "T\RH| ",
for xh in frange(0., 101., 10.):
    print " %-3s " % ("%d" % (xh,),),
print
print "-"*(11*6)+"------",

frow = True
for y in frange(0., 51., hygrometry.conv_f2c(37.)):
    print
    print "%-3s | " % ("%d" % (hygrometry.conv_c2f(y),),),  
    for x in frange(0., 101., 10.):
        info, errs = calc_maxerrors(hygrometry.wetbulb, [[y, y+1], [x, x+1], [904.0, 1050.0]], [error_Tc, error_RH, error_P], [5, 10, 2.7], True)
        merr = array(errs[3]).max()
        print "%.2f " % (hygrometry.conv_c2f(merr)-32.0,),
print


"""
Explanation: Max Error Table
End of explanation
"""


values = batch_WB(hygrometry.wetbulb, [hygrometry.conv_f2c(85.), hygrometry.conv_f2c(110.), 0.01], [65., 66., 65.], [977., 978., 977.])
Tv, RHv, Pv, Vv = zip(*values)
figure()
        
plot(Tv, Vv, 'r')
xlabel('T')
ylabel('WB')
show()

"""
Explanation: Backward Stability
If we use the same RH and P and only change T, it should be smooth or that might indicate not backward stable.
End of explanation
"""


values = batch_WB(hygrometry.wetbulb, [37., 38., 37.], [20, 100, 0.1], [977., 978., 977.])
Tv, RHv, Pv, Vv = zip(*values)
figure()
plot(RHv, Vv, 'r')
xlabel('RH')
ylabel('WB')
show()

"""
Explanation: Lets do the same check with RH
End of explanation
"""


values = batch_WB(hygrometry.wetbulb, [37., 38., 37.], [65., 66., 65.], [904., 1050., 0.1])
Tv, RHv, Pv, Vv = zip(*values)
figure()
        
plot(Pv, Vv, 'r')
xlabel('P')
ylabel('WB')
show()

"""
Explanation: Backward stability check with P
End of explanation
"""


step_Tc = 0.3
step_RH = 1.8
step_P = 2.7
ssteps=[step_Tc, step_RH, step_P]
info, errs = calc_uncertainty(hygrometry.wetbulb, srange, [error_Tc, error_RH, error_P], ssteps, True)
print info
print "Uncertainty: %f" % (info[0],)
plot2derrors(errs)

"""
Explanation: Uncertainty
The algorithm I use to calculate uncertainty is also a good test for backward stability.
End of explanation
"""


def error_ambientP(P):
    return 17.0

print "Change of Wet Bulb Error in Fahrenheit caused from lack of pressure sensor"
print "T\RH| ",
for xh in frange(0., 101., 10.):
    print " %-3s " % ("%d" % (xh,),),
print
print "-"*(11*6)+"------",

frow = True
for y in frange(0., 51., hygrometry.conv_f2c(37.)):
    print
    print "%-3s | " % ("%d" % (hygrometry.conv_c2f(y),),),  
    for x in frange(0., 101., 10.):
        info, errs = calc_maxerrors(hygrometry.wetbulb, [[y, y+1], [x, x+1], [904.0, 1050.0]], [error_Tc, error_RH, error_P], [5., 10., 2.7], True)
        merr = array(errs[3]).max()
        info, errs = calc_maxerrors(hygrometry.wetbulb, [[y, y+1], [x, x+1], [904.0, 1050.0]], [error_Tc, error_RH, error_ambientP], [5., 10., 2.7], True)
        merr2 = array(errs[3]).max()
        print "%.2f " % (hygrometry.conv_c2f(math.fabs(merr-merr2))-32.0,),
print

"""
Explanation: Pressure Sensor Configuration
Does lack of having a pressure sensor affect our errors if pressure can vary from 17hPa from some known value?
End of explanation
"""


