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 
    """ 
    if far< 0.: 
        return NAN 
    elif far < 0.005:
        ag= air_gamma_1(t)
        ag[np.logical_or(t< 379., t> 4731.)]= NAN
        return ag
    elif far< 0.069:
        ag= air_gamma_2(t, far)
        ag[np.logical_or(t< 699., t> 4731.)]= NAN
        return ag
    else: 
        return NAN 

def air_gamma_1(t):
    return -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

def air_gamma_2(t, far):
    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 
    return a6* t**6.+ \
           a5* t**5.+ \
           a4* t**4.+ \
           a3* t**3.+ \
           a2* t**2.+ \
           a1* t+ a0 

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)