{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "\n",
    "# symbolic variables for Helmholtz energy\n",
    "fi, fr, delta, tau = symbols(\"fi fr delta tau\")\n",
    "\n",
    "# symbolic variables for parameters in bank of terms\n",
    "li1, li2, ei1, ei2, ci1, ci2, si1, si2 = symbols(\"li1 li2 ei1 ei2 ci1 ci2 si1 si2\")\n",
    "pi1, pi2, pi3 = symbols(\"pi1 pi2 pi3\")\n",
    "bi1, bi2, bi3, bi4 = symbols(\"bi1 bi2 bi3 bi4\")\n",
    "gi1, gi2, gi3, gi4, gi5, gi6, gi7, gi8, gi9 = symbols(\"gi1 gi2 gi3 gi4 gi5 gi6 gi7 gi8 gi9\")\n",
    "ni1, ni2, ni3, ni4, ni5, ni6, ni7, ni8, ni9, ni10, ni11, ni12 = symbols(\"ni1 ni2 ni3 ni4 ni5 ni6 ni7 ni8 ni9 ni10 ni11 ni12\")\n",
    "Distance, Phi = symbols(\"Distance Phi\")\n",
    "\n",
    "# terms for ideal Helmholtz energy\n",
    "logTerms = + li1*log(tau**li2);\n",
    "idealPolyTerms = + pi1*tau**pi2;\n",
    "EinsteinTerms = + ei1*log(1 - exp(ei2*tau));\n",
    "coshTerms = - ci1*log((cosh(ci2*tau)));\n",
    "sinhTerms = + si1*log((sinh(si2*tau)));\n",
    "\n",
    "# terms for residual Helmholtz energy\n",
    "residualPolyTerms = pi1*tau**pi2*delta**pi3;\n",
    "BWRTerms = bi1*tau**bi2*delta**bi3*exp(-delta**bi4);\n",
    "GaussTerms = gi1*tau**gi2*delta**gi3*exp(gi6*(delta-gi9)**2 + gi7*(tau-gi8)**2);\n",
    "Distance = ((1-tau)+ni8*((delta-1)**2)**(1/(2*ni7)))**2 + ni11*((delta-1)**2)**ni12\n",
    "Phi = exp(-ni9*(delta-1)**2-ni10*(tau-1)**2)\n",
    "nonAnalyticalTerms = ni1*delta*Distance**ni6*Phi;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-ci1*log(cosh(ci2*tau)) + li1*log(tau**li2) + si1*log(sinh(si2*tau))\n",
      "delta*ni1*(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)**ni6*exp(-ni10*(tau - 1)**2 - ni9*(delta - 1)**2)\n"
     ]
    }
   ],
   "source": [
    "# Set the Helmholtz energy we want to look at\n",
    "# fi = +log(delta) +logTerms +idealPolyTerms +EinsteinTerms +coshTerms +sinhTerms\n",
    "fi = logTerms + coshTerms + sinhTerms\n",
    "# fr = +residualPolyTerms +BWRTerms +GaussTerms +nonAnalyticalTerms\n",
    "fr = nonAnalyticalTerms\n",
    "\n",
    "# and print terms, just to make sure\n",
    "print(fi)\n",
    "print(fr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-ci1*ci2*tanh(ci2*tau) + li1*li2/tau + si1*si2/tanh(si2*tau)\n",
      "-2*delta*ni1*(ni10*(tau - 1)*(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)**(ni6 + 1) + ni6*(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)**ni6*(ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1))*exp(-ni10*(tau - 1)**2 - ni9*(delta - 1)**2)/(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)\n",
      "ni1*(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)**(ni6 - 1)*(2*delta*ni6*(ni11*ni12*ni7*((delta - 1)**2)**ni12 + ni8*(ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)*((delta - 1)**2)**(1/(2*ni7))) + ni7*(delta - 1)*(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)*(-2*delta*ni9*(delta - 1) + 1))*exp(-ni10*(tau - 1)**2 - ni9*(delta - 1)**2)/(ni7*(delta - 1))\n"
     ]
    }
   ],
   "source": [
    "## First derivs\n",
    "# Calculate\n",
    "fit = simplify(diff(fi,tau))\n",
    "frt = simplify(diff(fr,tau))\n",
    "frd = simplify(diff(fr,delta))\n",
    "\n",
    "## Print results\n",
    "print(fit)\n",
    "print(frt)\n",
    "print(frd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Second derivs\n",
    "# Calculate\n",
    "fitt = simplify(diff(fit,tau))\n",
    "frtt = simplify(diff(frt,tau))\n",
    "frtd = simplify(diff(frt,delta))\n",
    "frdd = simplify(diff(frd,delta))\n",
    "\n",
    "# Print results\n",
    "print(fitt)\n",
    "print(frtt)\n",
    "print(frtd)\n",
    "print(frdd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Third derivs\n",
    "#Calculate\n",
    "fittt = simplify(diff(fitt,tau))\n",
    "frttt = simplify(diff(frtt,tau))\n",
    "frttd = simplify(diff(frtd,tau))\n",
    "frtdd = simplify(diff(frtd,delta))\n",
    "frddd = simplify(diff(frdd,delta))\n",
    "\n",
    "# Print results\n",
    "print(fittt)\n",
    "print(frttt)\n",
    "print(frttd)\n",
    "print(frtdd\n",
    "print(frddd)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
