from CoolProp.Plots import Ph,Ts
from CoolProp.CoolProp import Props
from scipy.optimize import newton
import pylab,numpy as np

params = {'axes.labelsize': 12,
          'text.fontsize': 10,
          'legend.fontsize': 8,
          'legend.labelspacing':0.2,
          'xtick.labelsize': 10,
          'ytick.labelsize': 10,
          'lines.linewidth': 0.5,
          'text.usetex': False,
          'font.family':'Times New Roman'}
pylab.rcParams.update(params)

def SimpleCycle(Ref,Te,Tc,DTsh,DTsc,eta_a,Ts_Ph='Ph',skipPlot=False,axis=None):
    """
    This function plots a simple four-component cycle, on the current axis, or that given by the optional parameter *axis*
    
    Required parameters:
    
    * Ref : A string for the refrigerant
    * Te : Evap Temperature in K
    * Tc : Condensing Temperature in K
    * DTsh : Evaporator outlet superheat in K
    * DTsc : Condenser outlet subcooling in K
    * eta_a : Adiabatic efficiency of compressor (no units) in range [0,1]
    
    Optional parameters:
    
    * Ts_Ph : 'Ts' for a Temperature-Entropy plot, 'Ph' for a Pressure-Enthalpy
    * axis : An axis to use instead of the active axis
    * skipPlot : If True, won't actually plot anything, just print COP
    
    """
    T=np.zeros((6))
    h=np.zeros_like(T)
    p=np.zeros_like(T)
    s=np.zeros_like(T)
    T[1]=Te+DTsh
    pe=Props('P','T',Te,'Q',1.0,Ref)
    pc=Props('P','T',Tc,'Q',1.0,Ref)
    h[1]=Props('H','T',T[1],'P',pe,Ref)
    s[1]=Props('S','T',T[1],'P',pe,Ref)
    T2s=newton(lambda T: Props('S','T',T,'P',pc,Ref)-s[1],T[1]+30)
    h2s=Props('H','T',T2s,'P',pc,Ref)
    h[2]=(h2s-h[1])/eta_a+h[1]
    T[2]=Props('T','H',h[2],'P',pc,Ref)
    s[2]=Props('S','T',T[2],'P',pc,Ref)
    
    sbubble_c=Props('S','P',pc,'Q',0,Ref)
    sdew_c=Props('S','P',pc,'Q',1,Ref)
    sbubble_e=Props('S','P',pe,'Q',0,Ref)
    sdew_e=Props('S','P',pe,'Q',1,Ref)
    T[3]=Tc-DTsc
    h[3]=Props('H','T',T[3],'P',pc,Ref)
    s[3]=Props('S','T',T[3],'P',pc,Ref)
    h[4]=h[3]
    h[5]=h[1]
    s[5]=s[1]
    T[5]=T[1]
    p=[np.nan,pe,pc,pc,pe,pe]
    COP=(h[1]-h[4])/(h[2]-h[1])
    COPH=(h[2]-h[3])/(h[2]-h[1])
    
    hsatL=Props('H','T',Te,'Q',0,Ref)
    hsatV=Props('H','T',Te,'Q',1,Ref)
    ssatL=Props('S','T',Te,'Q',0,Ref)
    ssatV=Props('S','T',Te,'Q',1,Ref)
    vsatL=1/Props('D','T',Te,'Q',0,Ref)
    vsatV=1/Props('D','T',Te,'Q',1,Ref)
    x=(h[4]-hsatL)/(hsatV-hsatL)
    s[4]=x*ssatV+(1-x)*ssatL
    T[4]=x*Te+(1-x)*Te
    
    print COP,COPH
    if skipPlot==False:
        if axis==None:
            ax=pylab.gca()
        if Ts_Ph in ['ph','Ph']:
            ax.plot(h,p)
        elif Ts_Ph in ['Ts','ts']:
            s=list(s)
            T=list(T)
            s.insert(5,sdew_e)
            T.insert(5,Te)
            s.insert(3,sbubble_c)
            T.insert(3,Tc)
            s.insert(3,sdew_c)
            T.insert(3,Tc)
            ax.plot(s[1::],T[1::],'b')
        else:
            raise TypeError('Type of Ts_Ph invalid')

def TwoStage(Ref,Q,Te,Tc,DTsh,DTsc,eta_oi,f_p,Tsat_ic,DTsh_ic,Ts_Ph='Ph',prints=False,skipPlot=False,axis=None,**kwargs):
    """
    This function plots a two-stage cycle, on the current axis, or that given by the optional parameter *axis*
    
    Required parameters:
    
    * Ref : Refrigerant [string]
    * Q : Cooling capacity [W]
    * Te : Evap Temperature [K]
    * Tc : Condensing Temperature [K]
    * DTsh : Evaporator outlet superheat [K]
    * DTsc : Condenser outlet subcooling [K]
    * eta_oi : Adiabatic efficiency of compressor (no units) in range [0,1]
    * f_p : fraction of compressor power lost as ambient heat transfer in range [0,1]
    * Tsat_ic : Saturation temperature corresponding to intermediate pressure [K]
    * DTsh_ic : Superheating at outlet of intermediate stage [K]
    
    Optional parameters:
    
    * Ts_Ph : 'Ts' for a Temperature-Entropy plot, 'Ph' for a Pressure-Enthalpy
    * prints : True to print out some values
    * axis : An axis to use instead of the active axis
    * skipPlot : If True, won't actually plot anything, just print COP
    
    """
    
    T=np.zeros((8))
    h=np.zeros_like(T)
    p=np.zeros_like(T)
    s=np.zeros_like(T)
    rho=np.zeros_like(T)
    T[0]=np.NAN
    s[0]=np.NAN
    T[1]=Te+DTsh
    pe=Props('P','T',Te,'Q',1.0,Ref)
    pc=Props('P','T',Tc,'Q',1.0,Ref)
    pic=Props('P','T',Tsat_ic,'Q',1.0,Ref)
    Tbubble_c=Props('T','P',pc,'Q',Ref)
    Tbubble_e=Props('T','P',pe,'Q',Ref)

    h[1]=Props('H','T',T[1],'P',pe,Ref)
    s[1]=Props('S','T',T[1],'P',pe,Ref)
    rho[1]=Props('D','T',T[1],'P',pe,Ref)
    T[5]=Tbubble_c-DTsc
    h[5]=Props('H','T',T[5],'P',pc,Ref)
    s[5]=Props('S','T',T[5],'P',pc,Ref)
    rho[5]=Props('D','T',T[5],'P',pc,Ref)
    mdot=Q/(h[1]-h[5])
    
    rho1=Props('D','T',T[1],'P',pe,Ref)
    T2s=newton(lambda T: Props('S','T',T,'P',pic,Ref)-s[1],T[1]+30)
    h2s=Props('H','T',T2s,'P',pic,Ref)
    Wdot1=mdot*(h2s-h[1])/eta_oi
    h[2]=h[1]+(1-f_p)*Wdot1/mdot
    T[2]=Props('T','H',h[2],'P',pic,Ref)
    s[2]=Props('S','T',T[2],'P',pic,Ref)
    rho[2]=Props('D','T',T[2],'P',pic,Ref)
    T[3]=288
    p[3]=pic
    h[3]=Props('H','T',T[3],'P',pic,Ref)
    s[3]=Props('S','T',T[3],'P',pic,Ref)
    rho[3]=Props('D','T',T[3],'P',pic,Ref)
    rho3=Props('D','T',T[3],'P',pic,Ref)
    T4s=newton(lambda T: Props('S','T',T,'P',pc,Ref)-s[3],T[3]+30)
    h4s=Props('H','T',T4s,'P',pc,Ref)
    Wdot2=mdot*(h4s-h[3])/eta_oi
    h[4]=h[3]+(1-f_p)*Wdot2/mdot
    T[4]=Props('T','H',h[4],'P',pc,Ref)
    s[4]=Props('S','T',T[4],'P',pc,Ref)
    rho[4]=Props('D','T',T[4],'P',pc,Ref)
    
    sbubble_e=Props('S','T',Tbubble_e,'Q',0,Ref)
    sbubble_c=Props('S','T',Tbubble_c,'Q',0,Ref)
    sdew_e=Props('S','T',Te,'Q',1,Ref)
    sdew_c=Props('S','T',Tc,'Q',1,Ref)
    
    hsatL=Props('H','T',Tbubble_e,'Q',0,Ref)
    hsatV=Props('H','T',Te,'Q',1,Ref)
    ssatL=Props('S','T',Tbubble_e,'Q',0,Ref)
    ssatV=Props('S','T',Te,'Q',1,Ref)
    vsatL=1/Props('D','T',Tbubble_e,'Q',0,Ref)
    vsatV=1/Props('D','T',Te,'Q',1,Ref)
    x=(h[5]-hsatL)/(hsatV-hsatL)
    s[6]=x*ssatV+(1-x)*ssatL
    T[6]=x*Te+(1-x)*Tbubble_e
    rho[6]=1.0/(x*vsatV+(1-x)*vsatL)
    
    h[6]=h[5]
    h[7]=h[1]
    s[7]=s[1]
    T[7]=T[1]
    p=[np.nan,pe,pic,pic,pc,pc,pe,pe]
    COP=Q/(Wdot1+Wdot2)
    RE=h[1]-h[6]
    
    if prints==True:
        print 'x5:',x
        print 'COP:',COP
        print T[2]-273.15,T[4]-273.15,p[2]/p[1],p[4]/p[3]
        print mdot,mdot*(h[4]-h[5]),pic
        print 'Vdot1',mdot/rho1,'Vdisp',mdot/rho1/(3500/60.)*1e6/0.7
        print 'Vdot2',mdot/rho3,'Vdisp',mdot/rho3/(3500/60.)*1e6/0.7
        print mdot*(h[4]-h[5]),Tc-273.15
        for i in range(1,len(T)-1):
            print '%d & %g & %g & %g & %g & %g \\\\' %(i,T[i]-273.15,p[i],h[i],s[i],rho[i])
            
    else:
        print Tsat_ic,COP
    
    if skipPlot==False:
        if axis==None:
            ax=pylab.gca()
        else:
            ax=axis
        if Ts_Ph in ['ph','Ph']:
            ax.plot(h,p)
        elif Ts_Ph in ['Ts','ts']:
            s_copy=s.copy()
            T_copy=T.copy()
            for i in range(1,len(s)-1):
                ax.plot(s[i],T[i],'bo',mfc='b',mec='b')
                dT=[0,-5,5,-20,5,5,5]
                ds=[0,0.05,0,0,0,0,0]
                ax.text(s[i]+ds[i],T[i]+dT[i],str(i))
                
            s=list(s)
            T=list(T)
            s.insert(7,sdew_e)
            T.insert(7,Te)
            s.insert(5,sbubble_c)
            T.insert(5,Tbubble_c)
            s.insert(5,sdew_c)
            T.insert(5,Tc)
            
            ax.plot(s,T)
            s=s_copy
            T=T_copy
        else:
            raise TypeError('Type of Ts_Ph invalid')
    return COP
    
def EconomizedCycle(Ref,Qin,Te,Tc,DTsh,DTsc,eta_oi,f_p,Ti,Ts_Ph='Ts',skipPlot=False,axis=None,**kwargs):
    """
    This function plots an economized cycle, on the current axis, or that given by the optional parameter *axis*
    
    Required parameters:
    
    * Ref : Refrigerant [string]
    * Qin : Cooling capacity [W]
    * Te : Evap Temperature [K]
    * Tc : Condensing Temperature [K]
    * DTsh : Evaporator outlet superheat [K]
    * DTsc : Condenser outlet subcooling [K]
    * eta_oi : Adiabatic efficiency of compressor (no units) in range [0,1]
    * f_p : fraction of compressor power lost as ambient heat transfer in range [0,1]
    * Ti : Saturation temperature corresponding to intermediate pressure [K]
    
    Optional parameters:
    
    * Ts_Ph : 'Ts' for a Temperature-Entropy plot, 'Ph' for a Pressure-Enthalpy
    * axis : An axis to use instead of the active axis
    * skipPlot : If True, won't actually plot anything, just print COP
    
    """
    
    m=1
    
    T=np.zeros((11))
    h=np.zeros_like(T)
    p=np.zeros_like(T)
    s=np.zeros_like(T)
    rho=np.zeros_like(T)
    
    T[0]=np.NAN
    s[0]=np.NAN
    T[1]=Te+DTsh
    pe=Props('P','T',Te,'Q',1.0,Ref)
    pc=Props('P','T',Tc,'Q',1.0,Ref)
    pi=Props('P','T',Ti,'Q',1.0,Ref)
    p[1]=pe
    h[1]=Props('H','T',T[1],'P',pe,Ref)
    s[1]=Props('S','T',T[1],'P',pe,Ref)
    rho[1]=Props('D','T',T[1],'P',pe,Ref)
    T2s=newton(lambda T: Props('S','T',T,'P',pi,Ref)-s[1],T[1]+30)
    h2s=Props('H','T',T2s,'P',pi,Ref)
    wdot1=(h2s-h[1])/eta_oi
    h[2]=h[1]+(1-f_p[0])*wdot1
    p[2]=pi
    T[2]=T_hp(Ref,h[2],pi,T2s)
    s[2]=Props('S','T',T[2],'P',pi,Ref)
    rho[2]=Props('D','T',T[2],'P',pi,Ref)
    
    T[5]=Tc-DTsc
    h[5]=Props('H','T',T[5],'P',pc,Ref)
    s[5]=Props('S','T',T[5],'P',pc,Ref)
    rho[5]=Props('D','T',T[5],'P',pc,Ref)
    
    p[5]=pc
    p[6]=pi
    h[6]=h[5]
    
    p[7]=pi
    p[8]=pi
    p[6]=pi
    T[7]=Ti
    h[7]=Props('H','T',Ti,'Q',1,Ref)
    s[7]=Props('S','T',Ti,'Q',1,Ref)
    rho[7]=Props('D','T',Ti,'Q',1,Ref)
    T[8]=Ti
    h[8]=Props('H','T',Ti,'Q',0,Ref)
    s[8]=Props('S','T',Ti,'Q',0,Ref)
    rho[8]=Props('D','T',Ti,'Q',0,Ref)
    x6=(h[6]-h[8])/(h[7]-h[8]) #Vapor Quality
    s[6]=s[7]*x6+s[8]*(1-x6)
    rho[6]=1.0/(x6/rho[7]+(1-x6)/rho[8])
    T[6]=Ti
    
    #Injection mass flow rate
    x=m*(h[6]-h[8])/(h[7]-h[6])
    
    
    p[3]=pi
    h[3]=(m*h[2]+x*h[7])/(m+x)
    T[3]=T_hp(Ref,h[3],pi,T[2])
    s[3]=Props('S','T',T[3],'P',pi,Ref)
    rho[3]=Props('D','T',T[3],'P',pi,Ref)
    T4s=newton(lambda T: Props('S','T',T,'P',pc,Ref)-s[3],T[2]+30)
    h4s=Props('H','T',T4s,'P',pc,Ref)
    p[4]=pc
    wdot2=(h4s-h[3])/eta_oi
    h[4]=h[3]+(1-f_p[1])*wdot2
    T[4]=T_hp(Ref,h[4],pc,T4s)
    s[4]=Props('S','T',T[4],'P',pc,Ref)
    rho[4]=Props('D','T',T[4],'P',pc,Ref)
    
    p[9]=pe
    h[9]=h[8]
    T[9]=Te
    hsatL_e=Props('H','T',Te,'Q',0,Ref)
    hsatV_e=Props('H','T',Te,'Q',1,Ref)
    ssatL_e=Props('S','T',Te,'Q',0,Ref)
    ssatV_e=Props('S','T',Te,'Q',1,Ref)
    vsatL_e=1/Props('D','T',Te,'Q',0,Ref)
    vsatV_e=1/Props('D','T',Te,'Q',1,Ref)
    x9=(h[9]-hsatL_e)/(hsatV_e-hsatL_e) #Vapor Quality
    s[9]=ssatV_e*x9+ssatL_e*(1-x9)
    rho[9]=1.0/(x9*vsatV_e+(1-x9)*vsatL_e)
    
    s[10]=s[1]
    T[10]=T[1]
    h[10]=h[1]
    p[10]=p[1]
    
    Tbubble_e=Te
    Tbubble_c=Tc
    sbubble_e=Props('S','T',Tbubble_e,'Q',0,Ref)
    sbubble_c=Props('S','T',Tbubble_c,'Q',0,Ref)
    sdew_e=Props('S','T',Te,'Q',1,Ref)
    sdew_c=Props('S','T',Tc,'Q',1,Ref)
    
    Wdot1=m*wdot1
    Wdot2=(m+x)*wdot2
    if skipPlot==False:
        if axis==None:
            ax=pylab.gca()
        else:
            ax=axis
        if Ts_Ph in ['ph','Ph']:
            ax.plot(h,p)
            ax.set_yscale('log')
        elif Ts_Ph in ['Ts','ts']:
            ax.plot(np.r_[s[7],s[3]],np.r_[T[7],T[3]],'b')
            s_copy=s.copy()
            T_copy=T.copy()
            dT=[0,-5,5,-12,5,12,-12,0,0,0]
            ds=[0,0.05,0.05,0,0.05,0,0.0,0.05,-0.05,-0.05]
            for i in range(1,len(s)-1):
                ax.plot(s[i],T[i],'bo',mfc='b',mec='b')
                ax.text(s[i]+ds[i],T[i]+dT[i],str(i),ha='center',va='center')
                
            s=list(s)
            T=list(T)
            s.insert(10,sdew_e)
            T.insert(10,Te)
            s.insert(5,sbubble_c)
            T.insert(5,Tbubble_c)
            s.insert(5,sdew_c)
            T.insert(5,Tc)
            ax.plot(s,T,'b')
            
            s=s_copy
            T=T_copy
        else:
            raise TypeError('Type of Ts_Ph invalid')
            
    COP=m*(h[1]-h[9])/(m*(h[2]-h[1])+(m+x)*(h[4]-h[3]))
    for i in range(1,len(T)-1):
            print '%d & %g & %g & %g & %g & %g \\\\' %(i,T[i]-273.15,p[i],h[i],s[i],rho[i])
    print x,m*(h[1]-h[9]),(m*(h[2]-h[1])+(m+x)*(h[4]-h[3])),COP
    mdot=Qin/(h[1]-h[9])
    mdot_inj=x*mdot
    print 'x9',x9,
    print 'Qcond',(mdot+mdot_inj)*(h[4]-h[5]),'T4',T[4]-273.15
    print mdot,mdot+mdot_inj
    f=3500/60.
    eta_v=0.7
    print 'Vdisp1: ',mdot/(rho[1]*f*eta_v)*1e6,'cm^3'
    print 'Vdisp2: ',(mdot+mdot_inj)/(rho[1]*f*eta_v)*1e6,'cm^3'
    return COP

if __name__=='__main__':
    
    Ref='R290'
    fig=pylab.figure(figsize=(4,3))
    ax=fig.add_axes((0.15,0.15,0.8,0.8))
    Ph(Ref,Tmin=273.15-100,hbounds=[0,600],axis=ax)
    COP=SimpleCycle(Ref,273.15-5,273.15+45,5,7,0.7,Ts_Ph='Ph')
    pylab.show()
    
    Ref='R410A'
    fig=pylab.figure(figsize=(4,3))
    ax=fig.add_axes((0.15,0.15,0.8,0.8))
    Ts(Ref,Tmin=273.15-100,sbounds=[0,600],axis=ax)
    COP=SimpleCycle(Ref,273.15-5,273.15+45,5,7,0.7,Ts_Ph='Ts')
    pylab.show()    
    
##     for x in np.linspace(0,1):
##         Ref='REFPROP-MIX:R152A[%g]&R32[%g]' %(x,1-x)
##         COP=SimpleCycle(273.15+8,273.15+44,5,7,0.7,skipPlot=True,Ts_Ph='Ph')
##     pylab.show()
