{
 "metadata": {
  "name": "",
  "signature": "sha256:4f22d42e96a1a086f159a429a1263cf7425ca5408c840ec424520d88bfd35839"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "import numpy as np\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from math import exp, fsum\n",
      "b = np.zeros((33,), dtype = np.double)\n",
      "\n",
      "fluid = 'R152a'\n",
      "if fluid == 'R123':\n",
      "    # Start with the coefficients from Younglove and McLinden as written\n",
      "\n",
      "    b = [0,-6.57453133659e-3,2.93479845842,-9.89140469845e1,2.01029776013e4,-3.83566527886e6,2.27587641969e-3,-9.08726819450,4.34181417995e3,3.54116464954e6,-6.35394849670e-4,3.20786715274,-1.31276484299e3,-1.16360713718e-1,-1.13354409016e1,-5.37543457327e3,2.58112416120,\n",
      "    -1.06148632128e-1,5.00026133667e1,-2.04326706346,-2.49438345685e6,-4.63962781113e8,-2.84903429588e5,9.74392239902e9,-6.37314379308e3,3.14121189813e5,-1.45747968225e2,-8.43830261449e6,-2.41138441593,1.08508031257e3,-1.06653193965e-2,-1.21343571084e1,-2.57510383240e2]\n",
      "\n",
      "    Tc = 456.831\n",
      "    rhoc = 3.596417\n",
      "    R = 0.08314510\n",
      "    M = 152.930\n",
      "    gamma = 1\n",
      "\n",
      "elif fluid == 'R152a':\n",
      "    \n",
      "    b = [0,-0.101623317e-1,+0.215677130e1,-0.648581254e2,+0.122535596e5,-0.206805988e7,-0.379836507e-3,-0.441333233e0,+0.158248875e3,+0.564062216e6,-0.124115350e-3,+0.494972179e0,-0.208058040e3,-0.131403187e-1,+0.212083849e0,-0.151263785e3,+0.311108025e-1,-0.115280980e-2,+0.437040026e0,-0.965596535e-2,-0.242705525e6,-0.518042519e8,-0.119070546e5,+0.459333195e9,-0.719317287e2,-0.840102861e4,-0.102910957e1,-0.325913881e5,-0.412362182e-2,+0.175102808e1,-0.198636625e-4,-0.421363036e-2,-0.198696761e1]\n",
      "    Tc = 386.411\n",
      "    rhoc = 5.57145\n",
      "    R = 0.08314471\n",
      "    M = 66.051\n",
      "    gamma = 1\n",
      "else:\n",
      "    raise ValueError(fluid)\n",
      "    \n",
      "# Coefficients for Equation 3.28 from Span, 2000\n",
      "n_2i = b[:]\n",
      "t_2i = np.array([0,  1,0.5,0,-1,-2,  1,0,-1,-2, 1,0,-1, 0, -1,-2, -1, -1,-2, -2, -2,-3, -2,-4, -2,-3, -2,-4, -2,-3, -2,-3,-4],dtype = np.longdouble)\n",
      "d_2i = np.array([0,  2,2,2,2,2, 3,3,3,3, 4,4,4, 5, 6,6, 7, 8,8, 9, 3,3, 5,5, 7,7, 9,9, 11,11, 13,13,13],dtype = np.longdouble)\n",
      "\n",
      "### Translate to the common Equation 3.26 formulation\n",
      "n_3_26 = [0] + [n_2i[i]*rhoc**(d_2i[i]-1)*Tc**(t_2i[i]-1)/R for i in range(1,33)]\n",
      "t_3_26 = [0] + [-t_2i[i]+1 for i in range(1,33)]\n",
      "d_3_26 = [0] + [d_2i[i]-1 for i in range(1,33)]\n",
      "\n",
      "# First polynomial part is easy (see Span, 2000, page 26)\n",
      "n = [0] + [n_3_26[i]/d_3_26[i] for i in range(1,20)]\n",
      "t = [0] + t_3_26[1:20]\n",
      "d = [0] + d_3_26[1:20]\n",
      "\n",
      "# Next part is harder\n",
      "n += [0]*21\n",
      "d += [0,0,0,0,0,0,2,2,2,4,4,4,6,6,6,8,8,8,10,10,10]\n",
      "t += [3,4,5,3,4,5,3,4,5,3,4,5,3,4,5,3,4,5,3,4,5]\n",
      "\n",
      "nn = n_3_26[:]\n",
      "n[20] = nn[20]/(2.0*gamma)+nn[22]/(2.0*gamma**2)+nn[24]/(gamma**3)+3*nn[26]/(gamma**4)+12*nn[28]/(gamma**5)+60*nn[30]/(gamma**6)\n",
      "n[21] = nn[21]/(2*gamma**2)+nn[25]/(gamma**3)+12*nn[29]/(gamma**5)+60*nn[31]/(gamma**6)\n",
      "n[22] = nn[23]/(2*gamma**2)+3*nn[27]/(gamma**4)+60*nn[32]/(gamma**6)\n",
      "n[23] = -n[20]\n",
      "n[24] = -n[21]\n",
      "n[25] = -n[22]\n",
      "\n",
      "n[26] = -nn[22]/(2*gamma)-nn[24]/(gamma**2)-3*nn[26]/(gamma**3)-12*nn[28]/(gamma**4)-60*nn[30]/(gamma**5)\n",
      "n[27] = -nn[25]/(gamma**2)-12*nn[29]/(gamma**4)-60*nn[31]/(gamma**5)\n",
      "n[28] = -nn[23]/(2*gamma)-3*nn[27]/(gamma**3)-60*nn[32]/(gamma**5)\n",
      "n[29] = -nn[24]/(2*gamma)-3*nn[26]/(2*gamma**2)-6*nn[28]/(gamma**3)-30*nn[30]/(gamma**4)\n",
      "n[30] = -nn[25]/(2*gamma)-6*nn[29]/(gamma**3)-30*nn[31]/(gamma**4)\n",
      "n[31] = -3*nn[27]/(2*gamma**2)-30*nn[32]/(gamma**4)\n",
      "n[32] = -nn[26]/(2*gamma)-2*nn[28]/(gamma**2)-10*nn[30]/(gamma**3)\n",
      "\n",
      "n[33] = -2*nn[29]/(gamma**2)-10*nn[31]/(gamma**3)\n",
      "n[34] = -nn[27]/(2*gamma)-10*nn[32]/(gamma**3)\n",
      "n[35] = -nn[28]/(2*gamma)-5*nn[30]/(2*gamma**2)\n",
      "n[36] = -nn[29]/(2*gamma)-5*nn[31]/(2*gamma**2)\n",
      "\n",
      "n[37] = -5*nn[32]/(2.0*gamma**2)\n",
      "n[38] = -nn[30]/(2.0*gamma)\n",
      "n[39] = -nn[31]/(2.0*gamma)\n",
      "n[40] = -nn[32]/(2.0*gamma)\n",
      "\n",
      "l = [0] + [0 for i in range(1,23)] + [2 for i in range(23,41)]\n",
      "\n",
      "# Reorder the terms in REFPROP ordering\n",
      "n = [0] + n[20:23] + n[1:20] + n[23::]\n",
      "d = [0] + d[20:23] + d[1:20] + d[23::]\n",
      "t = [0] + t[20:23] + t[1:20] + t[23::]\n",
      "l = [0] + l[20:23] + l[1:20] + l[23::]\n",
      "\n",
      "print 'n:', n[1:41]\n",
      "print 't:', t[1:41]\n",
      "print 'd:', d[1:41]\n",
      "print 'l:', l[1:41]\n",
      "\n",
      "for i in range(1,41):\n",
      "    print n[i],t[i],d[i],l[i]\n",
      "\n",
      "if fluid =='R123':\n",
      "    def p_328(T_K, rho_molL):\n",
      "        T = T_K\n",
      "        rho = rho_molL\n",
      "\n",
      "        s = n_2i*T**t_2i*rho**d_2i\n",
      "        sum1 = fsum(s[1:20])\n",
      "        sum2 = fsum(s[20:32])\n",
      "\n",
      "        # Form of equation 3.28 from Span, 2000\n",
      "        return rho*R*T + sum1 + exp(-(rho/rhoc)**2)*sum2\n",
      "\n",
      "    print p_328(146+273.15, 129.92/M), 'should yield', 19.5614, 'to within rounding error'\n",
      "\n",
      "    def p_326(T_K, rho_molL):\n",
      "        T = T_K\n",
      "        rho = rho_molL\n",
      "        tau = Tc/T\n",
      "        delta = rho/rhoc\n",
      "\n",
      "        sum1 = fsum([n_3_26[i]*tau**t_3_26[i]*delta**d_3_26[i] for i in range(1,20)])\n",
      "        sum2 = fsum([n_3_26[i]*tau**t_3_26[i]*delta**d_3_26[i] for i in range(20,33)])\n",
      "\n",
      "        # Form of equation 3.26 from Span, 2000\n",
      "        return rho*R*T*(1+sum1 + exp(-(rho/rhoc)**2)*sum2)\n",
      "\n",
      "    print p_326(146+273.15, 129.92/M), 'should yield', 19.5614, 'to within rounding error'\n",
      "\n",
      "    def p_325(T_K, rho_molL):\n",
      "        T = T_K\n",
      "        rho = rho_molL\n",
      "        tau = Tc/T\n",
      "        delta = rho/rhoc\n",
      "\n",
      "        sum1 = fsum([n[i]*d[i]*tau**t[i]*delta**(d[i]-1) for i in range(1,23)])\n",
      "        sum2 = fsum([n[i]*delta**(d[i]-1)*(d[i]-2*delta**2)*exp(-delta**2)*tau**t[i] for i in range(23,41)])\n",
      "\n",
      "        daddelta = sum1 + sum2\n",
      "        print tau, delta, daddelta\n",
      "        # Form of equation 3.26 from Span, 2000\n",
      "        return rho*R*T*(1+delta*daddelta)\n",
      "\n",
      "    print p_325(146+273.15, 129.92/M), 'should yield', 19.5614, 'to within rounding error'\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "n: [-3.5465795001851741, -0.36463127995636307, 0.033323333061128141, -0.6809684338301859, 7.3521264810280345, -11.247306378057127, 5.4991671429765772, -2.401863270213656, -0.070903644643947439, -0.21320088682136959, 0.19783973673284622, 1.824947698265976, -0.086054647670446016, 0.88813736685420042, -0.9661273471407773, -0.098522347852955536, 0.018341936863418316, -0.033855020406882257, 0.012492110085780536, -0.0022105670710377131, 0.0021687913327794162, -0.0002335976904702689, 3.5465795001851741, 0.36463127995636307, -0.033323333061128141, 2.7613383040168675, -0.069118571023672817, -0.033323333061128141, 0.7827613268561544, -0.034559285511836409, 0.13781353206688696, 0.18617312581521558, -0.03411193928992716, 0.045937844022295643, 0.021647001270605898, -0.0085279848224817899, 0.006203940397171849, 0.0018521029114874563, 0.00101674662708817, 0.0012407880794343697]\n",
        "t: [3, 4, 5, 0.0, 0.5, 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 2.0, 3.0, 3.0, 3, 4, 5, 3, 4, 5, 3, 4, 5, 3, 4, 5, 3, 4, 5, 3, 4, 5]\n",
        "d: [0, 0, 0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 4.0, 5.0, 5.0, 6.0, 7.0, 7.0, 8.0, 0, 0, 0, 2, 2, 2, 4, 4, 4, 6, 6, 6, 8, 8, 8, 10, 10, 10]\n",
        "l: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]\n",
        "-3.54657950019 3 0 0\n",
        "-0.364631279956 4 0 0\n",
        "0.0333233330611 5 0 0\n",
        "-0.68096843383 0.0 1.0 0\n",
        "7.35212648103 0.5 1.0 0\n",
        "-11.2473063781 1.0 1.0 0\n",
        "5.49916714298 2.0 1.0 0\n",
        "-2.40186327021 3.0 1.0 0\n",
        "-0.0709036446439 0.0 2.0 0\n",
        "-0.213200886821 1.0 2.0 0\n",
        "0.197839736733 2.0 2.0 0\n",
        "1.82494769827 3.0 2.0 0\n",
        "-0.0860546476704 0.0 3.0 0\n",
        "0.888137366854 1.0 3.0 0\n",
        "-0.966127347141 2.0 3.0 0\n",
        "-0.098522347853 1.0 4.0 0\n",
        "0.0183419368634 2.0 5.0 0\n",
        "-0.0338550204069 3.0 5.0 0\n",
        "0.0124921100858 2.0 6.0 0\n",
        "-0.00221056707104 2.0 7.0 0\n",
        "0.00216879133278 3.0 7.0 0\n",
        "-0.00023359769047 3.0 8.0 0\n",
        "3.54657950019 3 0 2\n",
        "0.364631279956 4 0 2\n",
        "-0.0333233330611 5 0 2\n",
        "2.76133830402 3 2 2\n",
        "-0.0691185710237 4 2 2\n",
        "-0.0333233330611 5 2 2\n",
        "0.782761326856 3 4 2\n",
        "-0.0345592855118 4 4 2\n",
        "0.137813532067 5 4 2\n",
        "0.186173125815 3 6 2\n",
        "-0.0341119392899 4 6 2\n",
        "0.0459378440223 5 6 2\n",
        "0.0216470012706 3 8 2\n",
        "-0.00852798482248 4 8 2\n",
        "0.00620394039717 5 8 2\n",
        "0.00185210291149 3 10 2\n",
        "0.00101674662709 4 10 2\n",
        "0.00124078807943 5 10 2\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}