import os, json


Simon_curves = {
    "n-Propane" : {
        "BibTeX" : "Reeves-JCP-1964", "T_m": -187.75 + 273.15, "parts": [{"T_0" : 85.3, "a" : 7.180e8, "c" : 1.283, "p_0" : 0.0, "T_max" : 168.63}]
    },
    "n-Pentane" : {
        "BibTeX" : "Reeves-JCP-1964", "T_m": -129.89 + 273.15, "parts": [{"T_0" : 143.5, "a" : 6.600e8, "c" : 1.649, "p_0" : 0.0, "T_max" : 156.2}]
    },
    "Isopentane" : {
        "BibTeX" : "Reeves-JCP-1964", "T_m": -159.92 + 273.15, "parts": [{"T_0" : 112.5, "a" : 5.916e8, "c" : 1.563, "p_0" : 0, "T_max" : 212.16}]
    },
    "Propylene" : {
        "BibTeX" : "Reeves-JCP-1964", "T_m": -185.09 + 273.15, "parts": [{"T_0" : 86.0, "a" : 3.196e8, "c" : 2.821, "p_0" : 0, "T_min": 86.0, "T_max" : 129},
                                                                         {"T_0" : 109.6, "a" : 3.064e8, "c" : 3.871, "p_0" : 4.450e8, "T_min": 129, "T_max" : 145.3}]
    },
    "Cyclohexane" : {
        "BibTeX" : "Penoncello-IJT-1995", "T_m": 6.81 + 273.15, "parts": [{"T_0" : 279.7, "a" : 383.4e6, "c" : 1.41, "p_0" : 0, "T_max" : 401.7}]
    },
    "Krypton" : {
        "BibTeX" : "Michels-PHYSICA-1962", "T_m": 115.95, "parts": [{"T_0" : 1, "a" : 109479.2307, "c" : 1.6169841, "p_0" : -237497645.7, "T_max" : 168.7}]
    },
    "Xenon" : {
        "BibTeX" : "Michels-PHYSICA-1962", "T_m": 165.02, "parts": [{"T_0" : 1, "a" : 80890.5544859, "c" : 1.5891650, "p_0" : -260932309.446, "T_max" : 366.4}]
    },
    "CarbonMonoxide" : {
        "BibTeX" : "Barreiros-JCT-1982", "T_m": 68.3, "parts": [{"T_0" : 1, "a" : 19560.8, "c" : 2.10747, "p_0" : -142921439.2, "T_max" : 87.5}]
    },
    "Oxygen": {
        "BibTeX" : "Younglove-NIST-1982", "T_m": 54.75, "parts": [{"T_0" : 1, "a" : 227606.348, "c" : 1.769, "p_0" : -266999247.652, "T_max" : 63.1}]
    },
    "ParaHydrogen": {
        "BibTeX" : "Younglove-NIST-1982", "T_m": 18.9, "parts": [{"T_0" : 1, "a" : 125746.643, "c" : 1.955, "p_0" : -21155737.752, "T_min" : 13.8033, "T_max" : 22},
                                                    {"T_0" : 1, "a" : 248578.596, "c" : 1.764739, "p_0" : -26280332.904, "T_min" : 22, "T_max" : 164.5}]
    },
    "Methane": {
        "BibTeX" : "Abramson-HPR-2011", "T_m": 90.7, "parts": [{"T_0" : 90.6941, "a" : 0.208e9, "c" : 1.698, "p_0" : 1.17e4, "T_max" : 600}]
    },
    "Helium": {
        "BibTeX" : "Datchi-PRB-2000", "T_m": 1.15, "parts": [{"T_0" : 1, "a" : 1.6067e6, "c" : 1.565, "p_0" : -1.6067e6, "T_max" : 700}]
    },
    "Neon": {
        "BibTeX" : "SantamariaPerez-PRB-2010", "T_m": -1, "parts": [{"T_0" : 24.4, "a" : 1.7e9, "c" : 1/0.77, "p_0" : 101325, "T_max" : 700}]
    },
    "Hydrogen": {
        "BibTeX" : "Datchi-PRB-2000", "T_m": 14.009985, "parts": [{"T_0" : 1, "a" : 2.31e5, "c" : 1.7627, "p_0" : -0.0052e6-2.31e5, "T_max" : 700}]
    }
}

polynomial_in_Tr = {
    "Argon" : {
        "BibTeX" : "Tegeler-JPCRD-1999", "T_m": 87.28, "parts": [{"T_0" : 83.8058, "a" : [-7476.2665, 9959.0613], "t" : [1.05,1.275], "p_0" : 68891, "T_max" : 254.0}]
    },
    "Fluorine" : {
        "BibTeX" : "deReuck-BOOK-1990", "T_m": 53.15, "parts": [{"T_0" : 53.4811, "a" : [988043.478261], "t" : [2.1845], "p_0" : 252, "T_max" : 55.4}]
    },
    "Nitrogen" : {
        "BibTeX" : "Span-JPCRD-2000", "T_m": 77.34, "parts": [{"T_0" : 63.151, "a" : [12798.61], "t" : [1.78963], "p_0" : 12523, "T_max" : 283.8}] 
    },
    "Ethane" : {
        "BibTeX" : "Buecker-JCRD-2006", "T_m": 90.4, "parts": [{"T_0" : 90.368, "a" : [2.23626315e8, 1.05262374e8], "t" : [1.0, 2.55], "p_0" : 1.14, "T_max" : 110.2}]
    },
    "Isobutane" : {
        "BibTeX" : "Buecker-JPCRD-2006B", "T_m": 113.55, "parts": [{"T_0" : 113.73, "a" : [1.9536371309e9], "t" : [6.12], "p_0" : 0.0219, "T_max" : 124.9}]
    },
    "Ethylene" : {
        "BibTeX" : "Smukala-JPCRD-2000", "T_m": 169, "parts": [{"T_0" : 103.989, "a" : [2947001.84], "t" : [2.045], "p_0" : 122.65, "T_min" : 103.989, "T_max" : 110.369},
                                                               {"T_0" : 110.369, "a" : [6.82693421], "t" : [1.089], "p_0" : 46.8e6, "T_min" : 110.369, "T_max" : 188}]
    },
    "n-Butane" : {
        "BibTeX" : "Buecker-JPCRD-2006B", "T_m": -137.92 + 273.15, "parts": [{"T_0" : 134.895, "a" : [5.585582364e8], "t" : [2.206], "p_0" : 0.653, "T_max" : 163.9}]
    },
    "Water" : {
        "BibTeX" : "IAPWS", "T_m": -1, "parts": [{"T_0" : 273.16, "a" : [-0.119539337e7,-0.808183159e5,-0.333826860e4], "t" : [0.3000000e1, 0.257500e2, 0.103750e3], "p_0" : 611.657, "T_min": 273.16, "T_max" : 251.165},
                                                 {"T_0" : 251.165, "a" : [0.299948], "t" : [60], "p_0" : 208.566e6, "T_min": 251.165, "T_max" : 256.164},
                                                 {"T_0" : 256.164, "a" : [1.18721], "t" : [8], "p_0" : 350.1e6, "T_min": 256.164, "T_max" : 273.31},
                                                 {"T_0" : 273.31, "a" : [1.07476], "t" : [4.6], "p_0" : 623.4e6, "T_min": 273.31, "T_max" : 355}
        ]
    }
}

polynomial_in_theta = {
    "Methanol" : {
        "BibTeX" : "deReuck-BOOK-1993", "T_m": 337.8, "parts": [{"T_0" : 175.61, "a" : [5.330770e9, 4.524780e9, 3.888861e10], "t" : [1, 1.5, 4], "p_0" : 0.187, "T_max" : 245.9}]
    },
    "CarbonDioxide" : {
        "BibTeX" : "Span-JPCRD-1996", "T_m": 216.58, "parts": [{"T_0" : 216.592, "a" : [1955.5390, 2055.4593], "t" : [1, 2], "p_0" : 517950, "T_max" : 327.6}]
    }
}

import CoolProp
__ = 0
for fluid in CoolProp.__fluids__:
    if fluid not in Simon_curves and fluid not in polynomial_in_Tr and fluid not in polynomial_in_theta:
        print fluid
        __ += 1
    else:
        print ' '*30, fluid
print __

import CoolProp.CoolProp as CP
import json, numpy as np, matplotlib.pyplot as plt, pandas

ip = 1
irho = 1
Nrow,Ncol = 5,5

figp = plt.figure(figsize = (20,20))
figrho = plt.figure(figsize = (20,20))

def plot_rho(T, rho, fit = False):
    x, y = (T-T[0])/(T[len(T)-1]-T[0]), (rho-rho[0])/(rho[len(rho)-1]-rho[0])
    
    c = np.polyfit(x, y, 3)
    yfit = np.polyval(c, x)
    err = yfit - y
    rms = np.sqrt(np.mean(np.power(err,2)))

    rhofit = yfit*(rho[len(rho)-1]-rho[0])+rho[0]
    if fit:
        return T, (rhofit/rho-1)*100
    else:
        return x, y
    
def simon():
    global ip, irho
    
    for fluid, values in Simon_curves.iteritems():
        axp = figp.add_subplot(Nrow, Ncol, ip); ip += 1
        axrho = figrho.add_subplot(Nrow, Ncol, irho); irho += 1
        axp.set_xlabel('T [K]')
        axp.set_ylabel('p [Pa]')
        axrho.set_xlabel('T [K]')
        axrho.set_ylabel('rho [mol/m$^3$]')
        axp.set_title(fluid+' - '+str(round(CP.Props(fluid,"molemass"),2)))
        axrho.set_title(fluid)
        
        fname = os.path.join('fluids',fluid+'.json')
        j = json.load(open(fname,'r'))
        for part in values['parts']:
            if 'T_min' not in part:
                part['T_min'] = round(CP.Props(fluid,"Tmin"),4)
        values['type'] = 'Simon'
                
        j['ANCILLARIES']['melting_line'] = values
        
        fp = open(fname,'w')
        from package_json import json_options
        fp.write(json.dumps(j,**json_options))
        fp.close()
        
#         if not isinstance(values, list):
#             values = [values]
#             df = pandas.read_csv('melting_curves/'+fluid+'.mlt',names=['T','p','rho'])
#             axp.plot(df['T'], df['p'], 'o', mfc='none')
#             x,y = plot_rho(df['T'],df['rho'],fit = True)
#             axrho.plot(x,y, 'o', mfc='none')
#         else:
#             for i in ['I','II']:
#                 df = pandas.read_csv('melting_curves/'+fluid+'-'+i+'.mlt',names=['T','p','rho'])
#                 axp.plot(df['T'], df['p'], 'o', mfc='none')
#                 x,y = plot_rho(df['T'],df['rho'],fit = True)
#                 axrho.plot(x,y, 'o', mfc='none')
        
        T_m = values['T_m']
        for i, value in enumerate(values['parts']):
            
            Tmin = value.get('T_min',CP.Props(fluid,"Tmin"))
            Tmax = value['T_max']
            
            T = np.linspace(Tmin, Tmax, 200)
            T_0 = value['T_0']
            p_0 = value['p_0']
            a = value['a']
            c = value['c']
            
            p = p_0 + a*((T/T_0)**c - 1)
            
            axp.plot(T, p)
            
            cc = 1.75
            aa = 3e8#(101325-p_0)/((T_m/T_0)**cc-1)
            pt = CP.Props(fluid,'ptriple')
            pp = pt + aa*((T/Tmin)**cc - 1)
            axp.plot(T_m,101325,'*')
            axp.plot(T,pp,'--')
            
            print fluid, CP.Props(fluid,"molemass"), CP.Props(fluid, 'accentric'), pp[-1]/p[-1]-1
            
#             if fluid == 'Helium':
#                 T = np.array([326.2,345.1,362.8,385.1,419.4,459,499,535.7,570,608])
#                 p = p_0 + a*((T/T_0)**c - 1)
#                 print p

def Tr():
    global ip, irho
    for fluid, values in polynomial_in_Tr.iteritems():
        
        axp = figp.add_subplot(Nrow, Ncol, ip); ip += 1
        axrho = figrho.add_subplot(Nrow, Ncol, irho); irho += 1
        axp.set_xlabel('T [K]')
        axp.set_ylabel('p [Pa]')
        axrho.set_xlabel('T [K]')
        axrho.set_ylabel('rho [mol/m$^3$]')
        axp.set_title(fluid+' - '+str(round(CP.Props(fluid,"molemass"),2)))
        axrho.set_title(fluid)
        
        fname = os.path.join('fluids',fluid+'.json')
        j = json.load(open(fname,'r'))
        for part in values['parts']:
            if 'T_min' not in part:
                part['T_min'] = round(CP.Props(fluid,"Tmin"),4)
        values['type'] = 'polynomial_in_Tr'
                
        j['ANCILLARIES']['melting_line'] = values
        
        fp = open(fname,'w')
        from package_json import json_options
        fp.write(json.dumps(j,**json_options))
        fp.close()
        
        if fluid == 'Ethylene':
            T = [104.003, 104.059, 104.13, 104.2, 104.27, 104.41, 104.55, 104.69, 104.83, 104.969, 105.108, 105.386, 106.077, 106.764, 107.446, 111.384, 119.283, 127.136, 158.146, 188.621]
            p = np.array([0.1, 0.5, 1, 1.5, 2, 3, 4, 5, 6, 7, 8, 10, 15, 20, 25, 50, 75, 100, 200, 300])*1e6
            
            axp.plot(T,p,'*')
            
#         if not isinstance(values, list):
#             values = [values]
#             df = pandas.read_csv('melting_curves/'+fluid+'.mlt',names=['T','p','rho'])
#             axp.plot(df['T'], df['p'], 'o', mfc='none')
#             x,y = plot_rho(df['T'],df['rho'],fit = True)
#             axrho.plot(x,y, 'o', mfc='none')
#         
#         else:
#             for i in ['I','II']:
#                 df = pandas.read_csv('melting_curves/'+fluid+'-'+i+'.mlt',names=['T','p','rho'])
#                 axp.plot(df['T'], df['p'], 'o', mfc='none')
#                 x,y = plot_rho(df['T'],df['rho'],fit = True)
#                 axrho.plot(x,y, 'o', mfc='none')
                
        T_m = values['T_m']
        for i,value in enumerate(values['parts']):
            
            Tmin = value.get('T_min',CP.Props(fluid,"Tmin"))
            Tmax = value['T_max']
            T = np.linspace(Tmin, Tmax, 200)
                
            a = value['a']
            t = value['t']
            T_t = value['T_0']
            p_t = value['p_0']
            
            RHS = 0
            for i in range(len(a)):
                RHS += a[i]*((T/T_t)**t[i] - 1)
        
            p = p_t*(RHS + 1)
            
            axp.plot(T, p)
            
            cc = 1.75
            aa = 3e8#(101325-p_0)/((T_m/T_0)**cc-1)
            pt = CP.Props(fluid,'ptriple')
            pp = pt + aa*((T/Tmin)**cc - 1)
            axp.plot(T_m,101325,'*')
            axp.plot(T,pp,'--')
            
            print fluid, CP.Props(fluid,"molemass"), CP.Props(fluid, 'accentric'), pp[-1]/p[-1]-1
        
    
def theta():
    global ip, irho
    for fluid, values in polynomial_in_theta.iteritems():
        
        axp = figp.add_subplot(Nrow, Ncol, ip); ip += 1
        axrho = figrho.add_subplot(Nrow, Ncol, irho); irho += 1
        axp.set_xlabel('T [K]')
        axp.set_ylabel('p [Pa]')
        axrho.set_xlabel('T [K]')
        axrho.set_ylabel('rho [mol/m$^3$]')
        axp.set_title(fluid+' - '+str(round(CP.Props(fluid,"molemass"),2)))
        axrho.set_title(fluid)
        
        fname = os.path.join('fluids',fluid+'.json')
        j = json.load(open(fname,'r'))
        for part in values['parts']:
            if 'T_min' not in part:
                part['T_min'] = round(CP.Props(fluid,"Tmin"),4)
        values['type'] = 'polynomial_in_Theta'
                
        j['ANCILLARIES']['melting_line'] = values
        
        fp = open(fname,'w')
        from package_json import json_options
        fp.write(json.dumps(j,**json_options))
        fp.close()
        
        T_m = values['T_m']
        for value in values['parts']:
            
            a = value['a']
            t = value['t']
            T_t = value['T_0']
            p_t = value['p_0']
            
            Tmin = T_t
            Tmax = value['T_max']
            T = np.linspace(Tmin, Tmax, 200)
            
            RHS = 0
            for i in range(len(a)):
                RHS += a[i]*(T/T_t - 1)**t[i]
            
            p = p_t*(RHS + 1)
            
            #df = pandas.read_csv('melting_curves/' + fluid + '.mlt', names=['T','p','rho'])
            
            #axp.plot(df['T'], df['p'], 'o', mfc='none')
            
            axp.plot(T, p)
        
            #x,y = plot_rho(df['T'],df['rho'],fit = True)
            #axrho.plot(x,y, 'o', mfc='none')
            
            cc = 1.75
            aa = 3e8#(101325-p_0)/((T_m/T_0)**cc-1)
            pt = CP.Props(fluid,'ptriple')
            pp = pt + aa*((T/Tmin)**cc - 1)
            axp.plot(T_m,101325,'*')
            axp.plot(T,pp,'--')
            
            print fluid, CP.Props(fluid,"molemass"), CP.Props(fluid, 'accentric'), pp[-1]/p[-1]-1
        
if __name__=='__main__':
    simon()
    Tr()
    theta()

    figp.tight_layout()
    figrho.tight_layout()
    
    figp.savefig('p.pdf')
    figrho.savefig('rho.pdf')
    
    plt.close()