import numpy as np

# Calculations assume IP units are: F, psi, gpm, cfm
# Calcuations assume SI units are: C, Pa, m3/s, m3/s; in some cases K and kPa are used internal to a function
# t = dry bulb temperature
class Psych():
    '''Calculate psychrometric parameters.'''
    def pSat(self,t,units = 'IP'): # units are Pa
        '''Calculate the saturation pressure given the dry bulb temperature.'''
        # ASHRAE Fundamentals, Chapter 1, 2017, Eqn. 6
        # Saturation pressure over liquid water for 0 °C to 200 °C (32 °F to 392 °F)
        # Convert to SI units if units is IP
        if units == 'IP':
            tC = 5/9*(t - 32) # [C]
        else:
            tC = t

        tK = tC + 273.15 # [K]

        C8 = -5.8002206e3
        C9 = 1.3914993
        C10 = -4.8640239e-2
        C11 = 4.1764768e-5
        C12 = -1.4452093e-8
        C13 = 6.5459673 

        pWS = np.exp(C8/tK + C9 + C10*tK + C11*tK**2 + C12*tK**3 + C13 * np.log(tK)) # [Pa]
        return pWS # Pa

    def pVap(self,t,RH,units = 'IP'):
        '''Calculate the vapor pressure given dry bulb temperature and relative humidity.'''
        if np.size(RH) > 1:            
            if (max(RH) > 1):
                RH = [x/100 for x in RH]
        else:
            if RH > 1:
                RH = RH/100
        pV = RH*self.pSat(t,units)
        return pV # Pa

    def relHum(self,t,w,pATM,units = 'IP'):
        '''Calculate the relative humidity given the dry bulb temperature, humidity ratio, and pressure.'''
        # ASHRAE Fundamentals, Chapter 1, 2017, Eqn. 23 
        # Conver to SI units if in IP
        if units == 'IP':
            tC = 5/9*(t - 32) # [C] 
            pATM = pATM*6895 # Pa
        else:
            tC = t # C
        # Calculate the saturation parameters
        pWS = self.pSat(tC,'C') # Pa
        wSat = self.omega(tC,np.array([]),pATM,pWS,'C') # kgWater/kgDryAir
        mu = w/wSat
        RH = mu/(1-(1-mu)*(pWS/pATM))
        if np.size(RH) > 1:
            RH = [max(0,x) for x in RH]
            RH = [min(1,x) for x in RH]
        else:
            RH = max(0,RH)
            RH = min(1,RH)
        return RH # 0 to 1

    def enthalpy(self,t,omega,units = 'IP'):
        '''Calculate the enthalpy of moist air given dry bulb temperature and humidity ratio.'''
        # ASHRAE Fundamentals, Chapter 1, 2017, Eqn. 30.
        if units == 'IP':
            tC = 5/9*(t - 32) # [C]   
        else:
            tC = t # C
        h = 1.006*tC+omega*(2501+1.86*tC)
        return h # kJ/kgDA

    def enthalpyDA(self,t,RH = [50],units = 'IP'):
        '''Calculate the enthalpy of dry air given dry bulb temperature.'''
        # ASHRAE Fundamentals, Chapter 1, 2017, Eqn. 27
        # specific enthalpy of dry air
        if units == 'IP':
            tC = 5/9*(t - 32) # [C]   
        else: 
            tC = t # C
        return 1.006*tC # kJ/kgDA

    def enthalpyG(self,t,units = 'IP'):
        '''Caclulate the enthalpy of the moist air given the dry bulb temperature.'''
        # ASHRAE Fundamentals, Chapter 1, 2017, Eqn. 29.
        if units == 'IP':
            tC = 5/9*(t - 32) # [C]    
        else: 
            tC = t # C
        return 2501+1.86*tC # kJ/kgW

    def omega(self,t,RH,pATM,p = np.array([]),units = 'IP'):
        '''Calculate the humidity ratio given dry bulb temperautre, relative humidity, and pressure.'''
        # ASHRAE Fundamentals, Chapter 1, 2017, Eqn. 20
        # humidity ratio
        if np.size(RH) >= 1:
            pV = self.pVap(t,RH,units)
        else:
            pV = p
        if units == 'IP':
            pATM = pATM*6895 # converts from psi to Pa
        omega = 0.622*pV/(pATM-pV+.01)
        return omega # kgWater/kgDA

    def cpTDA(self,omega,t,units = 'IP'):
        '''Calculate cp*T of dry air, which is roughly the sensible enthalpy, given humidity ratio and dry bulb temperature.'''
        # ASHRAE Fundamentals, Chapter 1, 2017, Eqn. 30
        if units == 'IP':
            tC = 5/9*(t - 32) 
        else:
            tC = t
        return (1.006+1.86*omega)*tC # kJ/kgDA

    def cp(self,omega,t = 70,units = 'IP'):
        '''Calculate cp, the specific heat, given humidity ratio.'''
        # ASHRAE Fundamentals, Chapter 1, 2017, Eqn. 30
        # Note: may want to double check that the T in the equation is °C and not K
        if units == 'IP':
            tC = 5/9*(t - 32) 
        else:
            tC = t
        return (1.006+1.86*omega) # kJ/kgDA

    def nuDA(self,omega,t,p,units = 'IP'):
        '''Calculate the specific volume of dry air, given humidity ratio, dry bulb temperature, and pressure.'''
        # ASHRAE Fundamentals, Chapter 1, 2017, Eqn. 26
        if units == 'IP':
            tC = 5/9*(t - 32) 
            p = p*6895/1000 # convert from psi to kPa
        else:
            tC = t  
            p = p/1000 # convert from Pa to kPa
        return 0.287042*(tC+273.15)*(1+1.607858*omega)/(p+.01) # m3/kgDA

    def mDA(self,flow,t,omega,p,units = 'IP'):
        '''Calculate the dry air mass flow rate for use in heat transfer equations given the volumetric airflow rate, 
           the dry bulb temperature, the humidity ratio, and pressure.''' 
        if units == 'IP':
            tC = 5/9*(t - 32)  
            p = p*6895 # convert from psi to Pa
            if np.size(flow) == 1:
                flow = flow*0.0004719 # cfm to m3/s
            else:
                flow = [x*0.0004719 for x in flow]
        else:
            tC = t  
        nu = self.nuDA(omega,tC,p,'SI') # m3/kgDA
        return flow/nu,nu # kgDA/s

    def dewPoint(self,t,RH,pATM,omega = np.array([]),units = 'IP'):
        '''calculate the dew point temperature given dry bulb temperature, relative humidity or humidity ratio, and pressure. If not using 
           relative humidity, enter [] for RH.'''
        # ASHRAE Fundamentals, Chapter 1, 2017, Eqn. 37, 0 °C to 93 °C
        # dewpoint temperature
        if np.size(RH) < 1:
            RH = self.relHum(t,omega,pATM,units)
        pW = self.pVap(t,RH,units)/1000 #convert from Pa to kPa
        alpha = np.log(pW)
        C14 = 6.54
        C15 = 14.526
        C16 = 0.7389
        C17 = 0.09486
        C18 = 0.4569
        tDP = C14 + C15*alpha + C16*alpha**2 + C17*alpha**3 + C18*(pW)**(0.1984)
        return tDP # C

    def specHeatMoist(self,t,RH,units = 'IP'):
        '''Calculate the moist air specific heat as calculated from a curve fit to EES data given dry bulb temperature and relative humidity.'''
        coefs = [0.00000000e+00,5.93474555e-05,5.16117690e-03,-1.28525163e-05,
                 3.70172025e-04,-1.41464346e-03,3.25406975e-07,3.93153415e-05,
                 7.16297464e-05,7.82887287e-04]
        if units == 'IP':
            tC = 5/9*(t-32) # [C]
        else:
            tC = t
        a = tC
        if np.size(RH) > 1:
            if max(RH) > 1:
                b = [x/100 for x in RH]
        else:
            if max(RH) > 1:
                b = RH/100

        if a.size > 1:
            cP = []
            for a1,b1 in zip(a,b):
                cP.append(coefs[0]*1+coefs[1]*a1+coefs[2]*b1+coefs[3]*a1**2+coefs[4]*a1*b1+coefs[5]*b1**2+coefs[6]*a1**3+ \
                     coefs[7]*a1**2*b1+coefs[8]*a1*b1**2+coefs[9]*b1**3+1.00487454)    
        else:
            cP = coefs[0]*1+coefs[1]*a+coefs[2]*b+coefs[3]*a**2+coefs[4]*a*b+coefs[5]*b**2+coefs[6]*a**3+ \
                 coefs[7]*a**2*b+coefs[8]*a*b**2+coefs[9]*b**3+1.00487454  
        return cP # kJ/kgDA-K

    def rhoMoist(self,t,RH,units = 'IP'):
        '''Calculate the moist air density as calculated from a curve fit to EES data given dry bulb temperature and relative humidity.'''
        coefs = [0.00000000e+00,-4.57636459e-03,-4.10246316e-03,9.08477714e-06,
                 -4.97255831e-04,6.44324044e-04]
        if units == 'IP':
            tC = 5/9*(t - 32) # [C]   
        else:
            tC = t
        a = tC
        if np.size(RH) > 1:
            if max(RH) > 1:
                b = [x/100 for x in RH]
        else:
            if max(RH) > 1:
                b = RH/100
        if np.size(b) <= 1:
            rho = coefs[0]*1 + coefs[1]*a + coefs[2]*b + coefs[3]*a**2 + coefs[4]*a*b + coefs[5]*b**2+1.29321441
        else:
            rho = [coefs[0]*1 + coefs[1]*x + coefs[2]*y + coefs[3]*x**2 + coefs[4]*x*y + coefs[5]*y**2+1.29321441 for x,y in zip(a,b)]
        return rho # kgDA/m3

    def rhoAir(self,t,units = 'IP'):
        '''Calculate the density and specific heat of dry air based on curve fits in EES.'''
        if units != 'IP':
            t = t*9/5+32
        # Based on curve fits in EES
        rho = 1.33822-0.00242228*t
        cp = 0.981398+0.000559998*t
        return rho,cp # kgDA/m3, kJ/kgDA-K