import numpy as np
NAN= float('nan') 

def air_gamma(t, far=0.0): 
    """ 
    Specific heat ratio (gamma) of Air/JP8 
    t - static temperature, Rankine 
    [far] - fuel air ratio [- defaults to 0.0 (dry air)] 
    air_gamma - specific heat ratio 
    """ 
    if far < 0.: 
        return NAN 
    elif far < 0.005: 
        if t < 379. or t > 4731.: 
            return NAN 
        else: 
            air_gamma = -3.472487e-22 * t ** 6. + 6.218811e-18 * t ** 5. - 4.428098e-14 * t ** 4. + 1.569889e-10 * t ** 3. - 0.0000002753524 * t ** 2. + 0.0001684666 * t + 1.368652 
    elif far < 0.069: 
        if t < 699. or t > 4731.: 
            return NAN 
        else: 
            a6 = 4.114808e-20 * far ** 3. - 1.644588e-20 * far ** 2. + 3.103507e-21 * far - 3.391308e-22 
            a5 = -6.819015e-16 * far ** 3. + 2.773945e-16 * far ** 2. - 5.469399e-17 * far + 6.058125e-18 
            a4 = 4.684637e-12 * far ** 3. - 1.887227e-12 * far ** 2. + 3.865306e-13 * far - 4.302534e-14 
            a3 = -0.00000001700602 * far ** 3. + 0.000000006593809 * far ** 2. - 0.000000001392629 * far + 1.520583e-10 
            a2 = 0.00003431136 * far ** 3. - 0.00001248285 * far ** 2. + 0.000002688007 * far - 0.0000002651616 
            a1 = -0.03792449 * far ** 3. + 0.01261025 * far ** 2. - 0.002676877 * far + 0.0001580424 
            a0 = 13.65379 * far ** 3. - 3.311225 * far ** 2. + 0.3573201 * far + 1.372714 
            air_gamma = a6 * t ** 6. + a5 * t ** 5. + a4 * t ** 4. + a3 * t ** 3. + a2 * t ** 2. + a1 * t + a0 
    elif far >= 0.069: 
        return NAN 
    else: 
        return NAN 
    return air_gamma
    
def air_gamma_0(t, far=0.0): 
    """ 
    Specific heat ratio (gamma) of Air/JP8 
    t - static temperature, Rankine 
    [far] - fuel air ratio [- defaults to 0.0 (dry air)] 
    air_gamma - specific heat ratio 
    """
    def ag(t, g, f, c1, c2):
        ag= g(t, f)
        ag[np.logical_or(c1, c2)]= NAN
        return ag
        
    if far< 0.: 
        return NAN 
    elif far< 0.005:
        return ag(t, ag_1, far, t< 379., t> 4731.)
    elif far< 0.069:
        return ag(t, ag_2, far, t< 699., t> 4731.)
    else: 
        return NAN 

def ag_1(t, far):
    a= [-3.472487e-22, 6.218811e-18, -4.428098e-14, \
         1.569889e-10, -0.0000002753524, 0.0001684666, 1.368652]
    return np.dot(np.vander(t, 7), np.array(a))

def ag_2(t, far):
    a= [[ 4.114808e-20, -1.644588e-20, 3.103507e-21, -3.391308e-22], \
        [-6.819015e-16, 2.773945e-16, -5.469399e-17, 6.058125e-18], \
        [ 4.684637e-12, -1.887227e-12, 3.865306e-13, -4.302534e-14], \
        [-0.00000001700602, 0.000000006593809, -0.000000001392629, 1.520583e-10], \
        [ 0.00003431136, -0.00001248285, 0.000002688007, -0.0000002651616], \
        [-0.03792449, 0.01261025, -0.002676877, 0.0001580424], \
        [13.65379, -3.311225, 0.3573201, 1.372714]]
    a= np.dot(a, np.vander([far], 4).T) 
    return np.dot(np.vander(t, 7), np.array(a)).flatten()

if __name__ == '__main__':
    """Some simple tests."""
    T= [400., 500., 600., 700., 800.]
    print [air_gamma(t) for t in T]
    print air_gamma_0(np.array(T))
    print [air_gamma(t, .06) for t in T]
    print air_gamma_0(np.array(T), .06)