{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "
\n", "" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import jinja2\n", "import json\n", "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "import matplotlib.path as mpath\n", "import matplotlib.patches as mpatches\n", "\n", "import mpld3\n", "from mpld3 import plugins, utils\n", "\n", "class PPlanePlugin(plugins.PluginBase):\n", " JAVASCRIPT = r\"\"\"\n", " // init custom PPLane plugin\n", " mpld3.register_plugin(\"drag\", PPlanePlugin);\n", " PPlanePlugin.prototype = Object.create(mpld3.Plugin.prototype);\n", " PPlanePlugin.prototype.constructor = PPlanePlugin;\n", " PPlanePlugin.prototype.requiredProps = [\"id\", \"idmpoint\", \"idlinesl\", \"idlinesr\", \n", " \"idqlm\", \"idqmm\", \"idqrm\", \"idqone\", \"idqtwo\", \"idzz\"];\n", " PPlanePlugin.prototype.defaultProps = {}\n", " function PPlanePlugin(fig, props){\n", " mpld3.Plugin.call(this, fig, props);\n", " mpld3.insert_css(\"#\" + fig.figid + \" path.dragging\",\n", " {\"fill-opacity\": \"1.0 !important\",\n", " \"stroke-opacity\": \"1.0 !important\"});\n", " };\n", "\n", " // Call draw function, this function is being looped all the time\n", " PPlanePlugin.prototype.draw = function(){\n", " // Get elements into script variables\n", " var obj = mpld3.get_element(this.props.id);\n", " var midpoint = mpld3.get_element(this.props.idmpoint);\n", " var linesl = mpld3.get_element(this.props.idlinesl);\n", " var linesr = mpld3.get_element(this.props.idlinesr);\n", " var qlm = mpld3.get_element(this.props.idqlm);\n", " var qmm = mpld3.get_element(this.props.idqmm);\n", " var qrm = mpld3.get_element(this.props.idqrm);\n", " var qone = mpld3.get_element(this.props.idqone);\n", " var qtwo = mpld3.get_element(this.props.idqtwo);\n", " var zz = this.props.idzz;\n", "\n", " // Set initial conditions for javascript calculations\n", " var qleft = obj.offsets[0];\n", " var qright = obj.offsets[1];\n", " var qmid = midpoint.offsets[0];\n", " var off = 13;\n", " \n", " // Main d3 drag function\n", " var drag = d3.behavior.drag()\n", " .origin(function(d) { return {x:obj.ax.x(d[0]),\n", " y:obj.ax.y(d[1])}; }) \n", " .on(\"dragstart\", dragstarted)\n", " .on(\"drag\", dragged)\n", " .on(\"dragend\", dragended);\n", " \n", " // Set elements of ql and qr points and call main drag function\n", " obj.elements()\n", " .data(obj.offsets)\n", " .style(\"cursor\", \"default\")\n", " .call(drag); \n", "\n", " // Begin drag function\n", " function dragstarted(d) {\n", " d3.event.sourceEvent.stopPropagation();\n", " d3.select(this).classed(\"dragging\", true);\n", " }\n", " \n", " // The drag function called while dragging is happening (meat of code here)\n", " function dragged(d,i) {\n", " // Convert mouse coordinates in drag event (d3.event) to python coordinates d\n", " d[0] = obj.ax.x.invert(d3.event.x);\n", " d[1] = obj.ax.y.invert(d3.event.y); \n", " // Move ql and qr stored in obj (they have been selected in drag)\n", " d3.select(this)\n", " .attr(\"transform\", \"translate(\" + [d3.event.x, d3.event.y] + \")\");\n", " // If obj corresponds to ql, move all the other left elements \n", " if (i==0){\n", " // Move lines and text marker\n", " linesl.elements().transition().duration(5)\n", " .attr(\"transform\", \"translate(\" + [d3.event.x, d3.event.y] + \")\");\n", " qlm.elements().transition().duration(1)\n", " .attr(\"transform\", \"translate(\" + [d3.event.x + off, d3.event.y + off] + \")\");\n", " // In script calculations of middle state\n", " qleft = [d[0], d[1]];\n", " var alphal = (-(qright[0] - qleft[0]) + zz*(qright[1] - qleft[1]))/(2*zz);\n", " qmid[0] = qleft[0] - alphal*zz;\n", " qmid[1] = qleft[1] + alphal;\n", " var xx = obj.ax.x(qmid[0]);\n", " var yy = obj.ax.y(qmid[1]); \n", " \n", " }\n", " // if element corresponds to qr\n", " else {\n", " // Move lines and text marker\n", " linesr.elements().transition().duration(5)\n", " .attr(\"transform\", \"translate(\" + [d3.event.x, d3.event.y] + \")\");\n", " qrm.elements().transition().duration(1)\n", " .attr(\"transform\", \"translate(\" + [d3.event.x + off, d3.event.y + off] + \")\");\n", " // In script calculations of middle state \n", " qright = [d[0], d[1]];\n", " var alphal = (-(qright[0] - qleft[0]) + zz*(qright[1] - qleft[1]))/(2*zz);\n", " qmid[0] = qleft[0] - alphal*zz;\n", " qmid[1] = qleft[1] + alphal;\n", " var xx = obj.ax.x(qmid[0]);\n", " var yy = obj.ax.y(qmid[1]);\n", " }\n", " // Update middle state point and marker position\n", " midpoint.elements().transition().duration(5)\n", " .attr(\"transform\", \"translate(\" + [xx, yy] + \")\");\n", " qmm.elements().transition().duration(5)\n", " .attr(\"transform\", \"translate(\" + [xx + 0.7*off, yy + 0.7*off] + \")\"); \n", " \n", " // Update subplots of q1 and q2\n", " qone.data[0][1] = qleft[0];\n", " qone.data[1][1] = qleft[0];\n", " qone.data[2][1] = qmid[0];\n", " qone.data[3][1] = qmid[0];\n", " qtwo.data[0][2] = qleft[1];\n", " qtwo.data[1][2] = qleft[1];\n", " qtwo.data[2][2] = qmid[1];\n", " qtwo.data[3][2] = qmid[1];\n", " qone.elements().transition()\n", " .attr(\"d\", qone.datafunc(qone.data));\n", " qtwo.elements().transition()\n", " .attr(\"d\", qtwo.datafunc(qtwo.data));\n", " }\n", " // End dragging\n", " function dragended(d) {\n", " d3.select(this).classed(\"dragging\", false);\n", " }\n", " }\n", " mpld3.register_plugin(\"drag\", PPlanePlugin); \n", " \"\"\"\n", "\n", " def __init__(self, points, midpoint, linesl, linesr, qlmarker, qmmarker, qrmarker, qone, qtwo,zz):\n", " if isinstance(points, mpl.lines.Line2D):\n", " suffix = \"pts\"\n", " else:\n", " suffix = None\n", "\n", " self.dict_ = {\"type\": \"drag\",\n", " \"id\": utils.get_id(points, suffix),\n", " \"idmpoint\": utils.get_id(midpoint,suffix),\n", " \"idlinesl\": utils.get_id(linesl,suffix),\n", " \"idlinesr\": utils.get_id(linesr,suffix),\n", " \"idqlm\": utils.get_id(qlmarker,suffix),\n", " \"idqmm\": utils.get_id(qmmarker,suffix),\n", " \"idqrm\": utils.get_id(qrmarker,suffix),\n", " \"idqone\": utils.get_id(qone),\n", " \"idqtwo\": utils.get_id(qtwo),\n", " \"idzz\": zz}\n", "\n", "\n", "# Create figure\n", "# Create a figure\n", "fig, ax = plt.subplots(1,3, figsize=(13.5, 4.5))\n", "fig.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.95,\n", " hspace=0.3, wspace=0.15)\n", "\n", "#Create subfigure ax1 for phaseplane\n", "#ax[0] = fig.add_subplot(2,2,1)\n", "\n", "# Calculate ql and qr randomly\n", "qL = np.array([-2.0, 3.0])\n", "qR = np.array([0.0, -2.0])\n", "dq = np.array([qR[0]-qL[0], qR[1]-qL[1]])\n", "\n", "# Plot eigenlines using paths as markers\n", "zz = 2.0\n", "rl = np.array([-zz, 1])\n", "rr = np.array([zz, 1])\n", "# Create path following eigenvalues r1 and r2\n", "verts = [\n", " (-rl[0], -rl[1]), # left, bottom\n", " (rl[0], rl[1]), # left, top\n", " (-rr[0], -rr[1]), # left, bottom\n", " (rr[0], rr[1]), # left, top\n", " ]\n", "codes = [mpath.Path.MOVETO,\n", " mpath.Path.LINETO,\n", " mpath.Path.LINETO,\n", " mpath.Path.LINETO\n", " ]\n", "path = mpath.Path(verts, codes)\n", "\n", "#Plot lines across ql and qr using paths as markers\n", "linesl = ax[0].plot(qL[0],qL[1], 'k', marker=path, markersize=2550, fillstyle='none')\n", "linesr = ax[0].plot(qR[0],qR[1], 'k', marker=path, markersize=2550, fillstyle='none')\n", "\n", "# Plot ql and qr\n", "points = ax[0].plot([qL[0],qR[0]], [qL[1], qR[1]], 'or', alpha=0.7, markersize=15, markeredgewidth=1)\n", "data = [\"q_l\", \"q_r\"]\n", "offset = 0.4\n", "qlmarker = ax[0].plot(qL[0] + offset, qL[1] - offset, 'ok', marker=(r\"$ q_l $\"), markersize=15)\n", "qrmarker = ax[0].plot(qR[0] + offset, qR[1] - offset, 'ok', marker=(r\"$ q_r $\"), markersize=15)\n", "\n", "#Plot midpoint\n", "alL = (-dq[0] + dq[1]*zz)/(2*zz)\n", "qm = qL + alL*rl \n", "midpoint = ax[0].plot(qm[0],qm[1],'ob', alpha=0.9, markersize=8, markeredgewidth=1)\n", "qmmarker = ax[0].plot(qm[0]+offset,qm[1]-0.7*offset, 'k',marker=(r\"$ q_m $\"),markersize=12)\n", "\n", "# Set axis 1 properties\n", "ax[0].set_title(\"Phase Plane\", fontsize=18)\n", "ax[0].axis([-5,5,-5,5])\n", "ax[0].set_aspect('equal')\n", "ax[0].grid(alpha=0.1,color='k', linestyle='--')\n", "\n", "# Remove defult mpld3 plugins\n", "plugins.clear(fig)\n", "\n", "#Create subfigure ax2 for solution plane\n", "#ax[1] = fig.add_subplot(2,2,3)\n", "# Create solutionl line with six points\n", "xsol = np.array([-5.0,-2.0,-2.0,2.0,2.0,5.0])\n", "qsol1 = 1*xsol\n", "qsol1[0:2] = qL[0]\n", "qsol1[2:4] = qm[0]\n", "qsol1[4:6] = qR[0]\n", "\n", "# Set axis 2 properties\n", "ax[1].set_title(\"q1 at time = 3\", fontsize=18)\n", "ax[1].axis([-5,5,-5,5])\n", "ax[1].grid(alpha=0.1,color='k', linestyle='--')\n", "ax[1].set_aspect('equal')\n", "\n", "# Plot solution\n", "qone = ax[1].plot(xsol, qsol1, '-k', linewidth = 4, alpha = 1.0)\n", "\n", "#Create subfigure ax2 for solution plane\n", "#ax[2] = fig.add_subplot(2,2,2)\n", "# Create solutionl line with six points\n", "xsol2 = np.array([-5.0,-2.0,-2.0,2.0,2.0,5.0])\n", "qsol2 = 1*xsol2\n", "qsol2[0:2] = qL[1]\n", "qsol2[2:4] = qm[1]\n", "qsol2[4:6] = qR[1]\n", "\n", "# Set axis 2 properties\n", "ax[2].set_title(\"q2 at time = 3\", fontsize=18)\n", "ax[2].axis([-5,5,-5,5])\n", "ax[2].grid(alpha=0.1,color='k', linestyle='--')\n", "ax[2].set_aspect('equal')\n", "\n", "# Plot solution\n", "qtwo = ax[2].plot(xsol2, qsol2, '-k', linewidth = 4, alpha = 1.0)\n", "\n", "# Call mpld3 custom PPLane plugin to interact with plot\n", "plugins.connect(fig, PPlanePlugin(points[0],midpoint[0],linesl[0],linesr[0],qlmarker[0],\n", " qmmarker[0],qrmarker[0],qone[0],qtwo[0],zz))\n", "\n", "mpld3.display()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "
\n", "" ], "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import jinja2\n", "import json\n", "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "import matplotlib.path as mpath\n", "import matplotlib.patches as mpatches\n", "\n", "\n", "import mpld3\n", "from mpld3 import plugins, utils\n", "\n", "class PPlaneNLPlugin(plugins.PluginBase):\n", " JAVASCRIPT = r\"\"\"\n", " // init custom PPLaneNL plugin\n", " mpld3.register_plugin(\"drag\", PPlaneNLPlugin);\n", " PPlaneNLPlugin.prototype = Object.create(mpld3.Plugin.prototype);\n", " PPlaneNLPlugin.prototype.constructor = PPlaneNLPlugin;\n", " PPlaneNLPlugin.prototype.requiredProps = [\"id\", \"idmpoint\", \"idg\", \"iditer\", \"iditer_charac\", \"idtime\",\n", " \"idhugol\", \"idhugor\", \"idintcl\", \"idintcr\",\n", " \"idhugol2\", \"idhugor2\", \"idintcl2\", \"idintcr2\", \n", " \"idqlm\", \"idqmm\", \"idqrm\", \n", " \"idshock1\", \"idshock2\", \"idrar1\", \"idrar2\", \n", " \"idtimedot\", \"idtimeline\",\n", " \"idq1\", \"idq2\"];\n", " PPlaneNLPlugin.prototype.defaultProps = {}\n", " function PPlaneNLPlugin(fig, props){\n", " mpld3.Plugin.call(this, fig, props);\n", " mpld3.insert_css(\"#\" + fig.figid + \" path.dragging\",\n", " {\"fill-opacity\": \"1.0 !important\",\n", " \"stroke-opacity\": \"1.0 !important\"});\n", " };\n", "\n", " // Call draw function, this function is being looped all the time\n", " PPlaneNLPlugin.prototype.draw = function(){\n", " // Get elements into script variables\n", " var obj = mpld3.get_element(this.props.id);\n", " var midpoint = mpld3.get_element(this.props.idmpoint);\n", " var g = this.props.idg;\n", " var iter = this.props.iditer;\n", " var iter_charac = this.props.iditer_charac;\n", " var time = this.props.idtime;\n", " var hugol = mpld3.get_element(this.props.idhugol);\n", " var hugor = mpld3.get_element(this.props.idhugor);\n", " var intcl = mpld3.get_element(this.props.idintcl);\n", " var intcr = mpld3.get_element(this.props.idintcr);\n", " var hugol2 = mpld3.get_element(this.props.idhugol2);\n", " var hugor2 = mpld3.get_element(this.props.idhugor2);\n", " var intcl2 = mpld3.get_element(this.props.idintcl2);\n", " var intcr2 = mpld3.get_element(this.props.idintcr2);\n", " var qlm = mpld3.get_element(this.props.idqlm);\n", " var qmm = mpld3.get_element(this.props.idqmm);\n", " var qrm = mpld3.get_element(this.props.idqrm);\n", " var shock1 = mpld3.get_element(this.props.idshock1);\n", " var shock2 = mpld3.get_element(this.props.idshock2);\n", " var rar1 = mpld3.get_element(this.props.idrar1);\n", " var rar2 = mpld3.get_element(this.props.idrar2);\n", " var timedot = mpld3.get_element(this.props.idtimedot);\n", " var timeline = mpld3.get_element(this.props.idtimeline);\n", " var q1 = mpld3.get_element(this.props.idq1);\n", " var q2 = mpld3.get_element(this.props.idq2);\n", " \n", " // Set initial conditions for javascript calculations\n", " //var qleft = obj.offsets[0];\n", " //var qright = obj.offsets[1];\n", " //var qmid = midpoint.offsets[0];\n", " var offx = obj.ax.x(0.3);\n", " var offy = offx;\n", " // Define initial values of hl, hr, hul and hur\n", " var hl = obj.offsets[0][0];\n", " var hul = obj.offsets[0][1];\n", " var hr = obj.offsets[1][0];\n", " var hur = obj.offsets[1][1];\n", " hmfinal = 0;\n", " humfinal = 0;\n", " \n", " // Main d3 drag function\n", " var drag = d3.behavior.drag()\n", " .origin(function(d) { return {x:obj.ax.x(d[0]),\n", " y:obj.ax.y(d[1])}; }) \n", " .on(\"dragstart\", dragstarted)\n", " .on(\"drag\", dragged)\n", " .on(\"dragend\", dragended);\n", " \n", " // Dragtime function for draggable time-dot (analog to main)\n", " var dragtime = d3.behavior.drag()\n", " .origin(function(d) { return {x:timedot.ax.x(d[0]),\n", " y:timedot.ax.y(d[1])}; }) \n", " .on(\"dragstart\", dragstarted)\n", " .on(\"drag\", dragged)\n", " .on(\"dragend\", dragended);\n", " \n", " // Set elements of ql and qr draggable points and call main drag function\n", " obj.elements()\n", " .data(obj.offsets)\n", " .style(\"cursor\", \"default\")\n", " .call(drag);\n", " \n", " // Set elements for timedot draggable point and call dragtime function \n", " timedot.elements()\n", " .data(timedot.offsets)\n", " .style(\"cursor\", \"default\")\n", " .call(dragtime);\n", "\n", " // Begin phi and phi prime function\n", " function phi(h,hl,hr,hul,hur) {\n", " var ul = hul/hl;\n", " var ur = hur/hr;\n", " if (h < hl) {\n", " var termIC1 = 2.0*Math.sqrt(g*hl);\n", " var termIC2 = 2.0*Math.sqrt(g*h);\n", " var phil = ul + termIC1 - termIC2;\n", " var philp = -Math.sqrt(g/h);\n", " }\n", " if (h >= hl) {\n", " var termHL1 = (h - hl);\n", " var termHL2 = Math.sqrt(0.5*g*(1.0/h + 1.0/hl));\n", " var termHL3 = termHL1*termHL2;\n", " var phil = ul - termHL3 ;\n", " var nume = termHL2*4*h*h;\n", " var philp = -termHL2 + termHL1*g/nume;\n", " }\n", " if (h < hr) {\n", " var termIC1 = 2.0*Math.sqrt(g*hr);\n", " var termIC2 = 2.0*Math.sqrt(g*h);\n", " var phir = ur - termIC1 + termIC2;\n", " var phirp = Math.sqrt(g/h);\n", " }\n", " if (h >= hr) {\n", " var termHL1 = (h - hr);\n", " var termHL2 = Math.sqrt(0.5*g*(1.0/h + 1.0/hr));\n", " var termHL3 = termHL1*termHL2;\n", " var phir = ur + termHL3 ;\n", " var nume = termHL2*4*h*h;\n", " var phirp = termHL2 - termHL1*g/nume;\n", " }\n", " var phi = phil-phir;\n", " var phiprime = philp - phirp\n", " var um = 0.5*(phil + phir);\n", " return [phi,phiprime,um];\n", " }\n", " \n", " // Begin Newton iteration\n", " function newton(hstar,hl,hr,hul,hur) {\n", " var hn = hstar;\n", " var error = 1;\n", " while (error > 0.005) {\n", " var hold = hn;\n", " var phin = phi(hn,hl,hr,hul,hur);\n", " var hn = hold - phin[0]/phin[1];\n", " var um = phin[2];\n", " error = Math.abs(hold - hn);\n", " }\n", " return [hn,um];\n", " }\n", " \n", " // Calculate hugoniot loci\n", " function hugoloci(hm,hside,huside,sign) {\n", " var termHL1 = hm*(hm - hside);\n", " if (hm == 0) {\n", " var termHL2 = 1.0; }\n", " else {\n", " var termHL2 = Math.sqrt(0.5*g*(1.0/hm + 1.0/hside));}\n", " var termHL3 = termHL1*termHL2;\n", " var hloci = hm*huside/hside + sign*termHL3;\n", " return hloci;} \n", " \n", " // Calculate integral curve\n", " function integralcurve(hm,hside,huside,sign) {\n", " var termIC1 = 2.0*hm*Math.sqrt(g*hside);\n", " var termIC2 = 2.0*hm*Math.sqrt(g*hm);\n", " var intcurve = hm*huside/hside - sign*(termIC1 - termIC2);\n", " return intcurve;}\n", " \n", " // Calculate solution plot of h and hu as function of x\n", " function solplot(hl,hm,hr,hul,hum,hur) {\n", " var lam = [];\n", " lam.push([-1000000, 0,0])\n", " var ul = hul/hl;\n", " var um = hum/hm;\n", " var ur = hur/hr;\n", " var shock1 = false;\n", " var shock2 = false;\n", " var qsol1 = d3.range(iter);\n", " var qsol2 = d3.range(iter);\n", " if (hm >= hl) {\n", " var s1 = time*(hum - hul)/(hm - hl);\n", " shock1 = true;\n", " } else {\n", " var s1b = [0,0];\n", " s1b[0] = time*(ul - Math.sqrt(g*hl));\n", " s1b[1] = time*(um - Math.sqrt(g*hm));\n", " s1b.sort(d3.ascending);\n", " }\n", " if (hm >= hr) {\n", " var s2 = time*(hur - hum)/(hr - hm);\n", " shock2 = true;\n", " } else {\n", " var s2b = [0,0];\n", " s2b[0] = time*(um + Math.sqrt(g*hm));\n", " s2b[1] = time*(ur + Math.sqrt(g*hr));\n", " s2b.sort(d3.ascending);\n", " }\n", " \n", " for (var ii=0; ii <2*iter; ii++){\n", " var xx = q1.data[ii][0];\n", " // Calculate plot solution for 1 charactersitic (shock or rarefaction)\n", " if (shock1) {\n", " if (xx <= s1){ qsol1[ii] = hl; qsol2[ii] = hul;}\n", " else { qsol1[ii] = hm; qsol2[ii] = hum;}\n", " }\n", " else {\n", " if (xx <= s1b[0]) { qsol1[ii] = hl; qsol2[ii] = hul;}\n", " else if (xx <= s1b[1]) {\n", " var m1 = (hm - hl)/(s1b[1] - s1b[0]);\n", " var m2 = (hum - hul)/(s1b[1] - s1b[0]);\n", " qsol1[ii] = m1*(xx - s1b[0]) + hl;\n", " qsol2[ii] = m2*(xx - s1b[0]) + hul;}\n", " else { qsol1[ii] = hm; qsol2[ii] = hum;}\n", " }\n", " // Calculate plot solution for 2 charactersitic (shock or rarefaction)\n", " if (shock2) {\n", " if (xx > s2) {qsol1[ii] = hr; qsol2[ii] = hur;}\n", " }\n", " else{\n", " if (xx > s2b[0] && xx < s2b[1]) {\n", " var m1 = (hr - hm)/(s2b[1] - s2b[0]);\n", " var m2 = (hur - hum)/(s2b[1] - s2b[0]);\n", " qsol1[ii] = m1*(xx - s2b[0]) + hm;\n", " qsol2[ii] = m2*(xx - s2b[0]) + hum;}\n", " else if (xx > s2b[1]) { qsol1[ii] = hr; qsol2[ii] = hur;}\n", " }\n", " }\n", " var solution = [qsol1,qsol2];\n", " return solution;\n", " }\n", "\n", " // Function to update middle state \n", " function update_midstate(){\n", " var hstar = 0.05*(hl + hr);\n", " var solution = newton(hstar,hl,hr,hul,hur);\n", " hmfinal = solution[0];\n", " humfinal = solution[1]*hmfinal;\n", " var xx = obj.ax.x(hmfinal);\n", " var yy = obj.ax.y(humfinal);\n", " // Update middle state point and marker position\n", " midpoint.elements().transition().duration(5)\n", " .attr(\"transform\", \"translate(\" + [xx, yy] + \")\");\n", " // Move marker\n", " qmm.elements().transition().duration(1)\n", " .attr(\"transform\", \"translate(\" + [xx + 0.7*offx, yy + 0.7*offy] + \")\");\n", " }\n", " \n", " // Functon to update xt-plane\n", " function update_xtplane() {\n", " // Calculate shcok speeds from R-H conditions\n", " var lam1 = (humfinal - hul)/(hmfinal - hl); \n", " var lam2 = (hur - humfinal)/(hr - hmfinal); \n", " var lam1m = humfinal/hmfinal - Math.sqrt(g*hmfinal);\n", " var lam2m = humfinal/hmfinal + Math.sqrt(g*hmfinal);\n", " var lam1l = hul/hl - Math.sqrt(g*hl);\n", " var lam2r = hur/hr + Math.sqrt(g*hr);\n", " var color1 = \"red\";\n", " var color2 = \"red\";\n", " for (var ii=0; ii= hl) {\n", " shock1.data[ii][1] = shock1.data[ii][0]/lam1;\n", " rar1.data[ii][1] = rar1.data[ii][0]/lam1;\n", " color1 = \"red\";\n", " } else {\n", " shock1.data[ii][1] = shock1.data[ii][0]/lam1l;\n", " rar1.data[ii][1] = rar1.data[ii][0]/lam1m;\n", " color1 = \"blue\";\n", " }\n", " if (hmfinal >= hr) {\n", " shock2.data[ii][2] = shock2.data[ii][0]/lam2;\n", " rar2.data[ii][2] = rar2.data[ii][0]/lam2;\n", " color2 = \"red\";\n", " } else {\n", " shock2.data[ii][2] = shock2.data[ii][0]/lam2r;\n", " rar2.data[ii][2] = rar2.data[ii][0]/lam2m;\n", " color2 = \"blue\";\n", " } \n", " }\n", " // Do transitions\n", " shock1.elements().transition().duration(5)\n", " .attr(\"d\", shock1.datafunc(shock1.data))\n", " .style(\"stroke\", color1);\n", " shock2.elements().transition().duration(5)\n", " .attr(\"d\", shock2.datafunc(shock2.data))\n", " .style(\"stroke\", color2);\n", " rar1.elements().transition().duration(5)\n", " .attr(\"d\", rar1.datafunc(rar1.data))\n", " .style(\"stroke\", color1);\n", " rar2.elements().transition().duration(5)\n", " .attr(\"d\", rar2.datafunc(rar2.data))\n", " .style(\"stroke\", color2);\n", " }\n", " \n", " // Function to update solution plots\n", " function update_solplots() {\n", " var qsol = solplot(hl,hmfinal,hr,hul,humfinal,hur);\n", " for (var ii=0; ii<2*iter; ii++){\n", " q1.data[ii][1] = qsol[0][ii]; \n", " q2.data[ii][2] = qsol[1][ii]; \n", " }\n", " // Do transitions\n", " q1.elements().transition().duration(5)\n", " .attr(\"d\", q1.datafunc(q1.data));\n", " q2.elements().transition().duration(5)\n", " .attr(\"d\", q2.datafunc(q2.data));\n", " }\n", " \n", " // Initialize solution with given initial states before interacting\n", " update_midstate();\n", " update_xtplane();\n", " update_solplots();\n", " \n", " // Begin drag function\n", " function dragstarted(d) {\n", " d3.event.sourceEvent.stopPropagation();\n", " d3.select(this).classed(\"dragging\", true);\n", " }\n", " \n", " // The drag function called while dragging is happening (meat of code here)\n", " function dragged(d,i) {\n", " if (i == 0 || i ==1) {\n", " // Convert mouse coordinates in drag event (d3.event) to python coordinates d\n", " d[0] = obj.ax.x.invert(d3.event.x);\n", " d[1] = obj.ax.y.invert(d3.event.y); \n", " // Move ql and qr stored in obj (they have been selected in drag)\n", " d3.select(this)\n", " .attr(\"transform\", \"translate(\" + [d3.event.x,d3.event.y] + \")\");\n", " // If obj corresponds to ql, move all the other left elements \n", " }\n", " if (i==0){\n", " // Move marker\n", " qlm.elements().transition().duration(1)\n", " .attr(\"transform\", \"translate(\" + [d3.event.x + offx, d3.event.y + offy] + \")\"); \n", " // Re-calculate inital left variables when dragging \n", " hl = obj.offsets[0][0];\n", " hul = obj.offsets[0][1];\n", " \n", " // Draw Hugoniot loci through left state \n", " for (var ii=0; ii=0) {\n", " timeline.data[0][1] = ty;\n", " timeline.data[1][1] = ty;\n", " timeline.data[2][1] = ty;\n", " time = ty;\n", " }\n", " // Do transitions\n", " timeline.elements().transition().duration(5)\n", " .attr(\"d\", timeline.datafunc(timeline.data));\n", " }\n", " \n", " // Calculate middle state \n", " update_midstate();\n", " \n", " // Calculate solution plots of h and hu\n", " update_solplots();\n", " \n", " // Update characteristic in x-t plane\n", " update_xtplane();\n", " \n", " }\n", " // End dragging\n", " function dragended(d) {\n", " d3.select(this).classed(\"dragging\", false);\n", " }\n", " }\n", " mpld3.register_plugin(\"drag\", PPlaneNLPlugin); \n", " \"\"\"\n", "\n", " def __init__(self, points, midpoint, g, iters, iter_charac, time,\n", " hugol, hugor, intcl, intcr,\n", " hugol2, hugor2, intcl2, intcr2, \n", " qlmarker, qmmarker ,qrmarker, \n", " shock1, shock2, rar1, rar2, \n", " timedot, timeline,\n", " q1, q2):\n", " if isinstance(points, mpl.lines.Line2D):\n", " suffix = \"pts\"\n", " else:\n", " suffix = None\n", "\n", " self.dict_ = {\"type\": \"drag\",\n", " \"id\": utils.get_id(points, suffix),\n", " \"idmpoint\": utils.get_id(midpoint,suffix),\n", " \"idg\" : g,\n", " \"iditer\" : iters,\n", " \"iditer_charac\" : iter_charac,\n", " \"idtime\" : time,\n", " \"idhugol\" : utils.get_id(hugol),\n", " \"idhugor\" : utils.get_id(hugor),\n", " \"idintcl\" : utils.get_id(intcl),\n", " \"idintcr\" : utils.get_id(intcr),\n", " \"idhugol2\" : utils.get_id(hugol2),\n", " \"idhugor2\" : utils.get_id(hugor2),\n", " \"idintcl2\" : utils.get_id(intcl2),\n", " \"idintcr2\" : utils.get_id(intcr2),\n", " \"idqlm\": utils.get_id(qlmarker,suffix),\n", " \"idqmm\": utils.get_id(qmmarker,suffix),\n", " \"idqrm\": utils.get_id(qrmarker,suffix),\n", " \"idshock1\" : utils.get_id(shock1),\n", " \"idshock2\" : utils.get_id(shock2),\n", " \"idrar1\" : utils.get_id(rar1),\n", " \"idrar2\" : utils.get_id(rar2),\n", " \"idtimedot\" : utils.get_id(timedot,suffix),\n", " \"idtimeline\" : utils.get_id(timeline),\n", " \"idq1\" : utils.get_id(q1),\n", " \"idq2\" : utils.get_id(q2)\n", " }\n", "\n", " \n", "def hugoloci(g,hm,hside,huside,sign):\n", " term1 = hm*(hm - hside)\n", " hm[hm==0] = hm[hm==0] + 0.00000001\n", " term2 = np.sqrt(0.5*g*(1.0/hm + 1.0/hside))\n", " hloci = hm*huside/hside + sign*term1*term2\n", " return hloci\n", "\n", "def intcurve(g,hm,hside,huside,sign):\n", " term1 = 2.0*hm*np.sqrt(g*hside);\n", " term2 = 2.0*hm*np.sqrt(g*hm);\n", " intcurve = hm*huside/hside - sign*(term1 - term2);\n", " return intcurve\n", " \n", "\n", "# Create figure\n", "# Create a figure\n", "fig, axfull = plt.subplots(2,2, figsize=(13, 12))\n", "fig.subplots_adjust(left=0.05, right=0.9, bottom=0.1, top=0.9,\n", " hspace=0.3, wspace=0.15)\n", "\n", "# Have to do this to avoid issue with mpld3\n", "ax = axfull[0] # First row of plots\n", "axsol = axfull[1] # Second row of plost \n", "\n", "# Calculate ql and qr randomly\n", "qL = np.array([2.0, 5.0])\n", "qR = np.array([2.0, -5.0])\n", "dq = np.array([qR[0]-qL[0], qR[1]-qL[1]])\n", "g = 1.0\n", "iters = 100\n", "iter_charac = 2\n", "time = 2\n", "# eps required for bug in mpld3 (qL[0] cannot be the same than qR[0])\n", "eps = 0.00000001\n", "\n", "# PLOT PHASE PLANE\n", "xxL = np.linspace(0,qL[0],iters)\n", "xxL2 = np.linspace(qL[0],10,iters)\n", "xxR = np.linspace(0,qR[0]+eps,iters)\n", "xxR2 = np.linspace(qR[0]+eps,10,iters)\n", "\n", "#Plot midpoint\n", "qm = -1 +0.0*qL\n", "midpoint = ax[0].plot(qm[0],qm[1],'ob', alpha=0.9, markersize=8, markeredgewidth=1)\n", "\n", "# Plot hugoloci initial state\n", "yyL = hugoloci(g,xxL,qL[0],qL[1],-1)\n", "yyL2 = hugoloci(g,xxL2,qL[0],qL[1],-1)\n", "yyR = hugoloci(g,xxR,qR[0],qR[1],1)\n", "yyR2 = hugoloci(g,xxR2,qR[0],qR[1],1)\n", "hugol = ax[0].plot(xxL,yyL, '--r', linewidth=1.5)\n", "hugol2 = ax[0].plot(xxL2,yyL2, '-r', linewidth=2, label = 'Hugoniot Loci')\n", "hugor = ax[0].plot(xxR,yyR, '--r', linewidth=1.5, label = 'Hugoniot Loci (unphysical)')\n", "hugor2 = ax[0].plot(xxR2,yyR2, '-r', linewidth=2)\n", "\n", "# Plot integral curve initial state\n", "yyL = intcurve(g,xxL,qL[0],qL[1],-1)\n", "yyL2 = intcurve(g,xxL2,qL[0],qL[1],-1)\n", "yyR = intcurve(g,xxR,qR[0],qR[1],1)\n", "yyR2 = intcurve(g,xxR2,qR[0],qR[1],1)\n", "intcl = ax[0].plot(xxL,yyL, '-b', linewidth=2, label = 'Integral Curves')\n", "intcl2 = ax[0].plot(xxL2,yyL2, '--b', linewidth=1.5, label = 'Integral Curves (unphysical)')\n", "intcr = ax[0].plot(xxR,yyR, '-b', linewidth=2)\n", "intcr2 = ax[0].plot(xxR2,yyR2, '--b', linewidth=1.5)\n", "\n", "# Plot ql and qr\n", "points = ax[0].plot([qL[0],qR[0]], [qL[1], qR[1]], 'or', alpha=0.7, markersize=15, markeredgewidth=1)\n", "#data = [\"q_l\", \"q_r\"] \n", "\n", "# Plot markers\n", "offsetx = 0.3\n", "offsety = -3*offsetx\n", "qlmarker = ax[0].plot(qL[0]+offsetx, qL[1]+offsety, 'ok', marker=(r\"$ q_l $\"), markersize=15)\n", "qmmarker = ax[0].plot(qm[0]+offsetx, qm[1]+offsety, 'ok', marker=(r\"$ q_m $\"), markersize=15)\n", "qrmarker = ax[0].plot(qR[0]+offsetx, qR[1]+offsety, 'ok', marker=(r\"$ q_r $\"), markersize=15)\n", "\n", "# Set axis 1 properties\n", "ax[0].set_title(\"Phase plane\", fontsize=18)\n", "ax[0].axis([0,10,-15,15])\n", "ax[0].set_xlabel('h', fontsize=17)\n", "ax[0].set_ylabel('hu', fontsize=17)\n", "#ax[0].set_aspect('equal')\n", "ax[0].grid(alpha=0.1,color='k', linestyle='--')\n", "legend = ax[0].legend(loc='upper left', shadow=True)\n", "\n", "\n", "# PLOT x-t PLANE\n", "x_xtp = np.linspace(-10,10,iter_charac)\n", "x_xtp2 = np.linspace(-11,11,iter_charac)\n", "# Shock speeds lam1 and lam2\n", "lam2 = qL[1]/qL[0] - np.sqrt(g*qL[0])\n", "lam1 = qR[1]/qR[0] + np.sqrt(g*qR[0])\n", "char1 = x_xtp/lam1\n", "char2 = x_xtp/lam2\n", "char3 = x_xtp2/lam1\n", "char4 = x_xtp2/lam2\n", "shock1 = ax[1].plot(x_xtp, char1, '-k', linewidth=4, label=\"1-wave\")\n", "shock2 = ax[1].plot(x_xtp, char2, '-k', linewidth=2, label=\"2-wave\")\n", "rar1 = ax[1].plot(x_xtp2, char3, '-k', linewidth=4)\n", "rar2 = ax[1].plot(x_xtp2, char4, '-k', linewidth=2)\n", "timedot = ax[1].plot([100000,1000000,9], [-10,-10,time], 'ok' , markersize=15)\n", "timeline = ax[1].plot([-12,0,12], [time, time, time], '--k', linewidth = 3, label=\"time\")\n", "\n", "# Set axis 2 properties\n", "ax[1].set_title(\"x-t plane\", fontsize=18)\n", "ax[1].set_xlabel('x', fontsize=17)\n", "ax[1].set_ylabel('t', fontsize=17)\n", "ax[1].axis([-10,10,-0,5])\n", "ax[1].grid(alpha=0.1,color='k', linestyle='--')\n", "legend = ax[1].legend(loc='upper center', shadow=True)\n", "\n", "# PLOT SOLUTIONS\n", "xsol = np.linspace(-10,10,2*iters)\n", "hsol = 0*xsol + qL[0]\n", "husol = 0*xsol + qL[1]\n", "q1 = axsol[0].plot(xsol,hsol, '-b', linewidth = 4, alpha = 0.5)\n", "q2 = axsol[1].plot(xsol,husol, '-b', linewidth = 4, alpha = 0.5)\n", "\n", "def solplot(xsol,qL,qR,qM,g):\n", " hl = qL[0]; hm = qM[0]; hr = qR[0]\n", " ul = qL[1]/hl; um = qM[1]/hm; ur = qR[1]/hr \n", " lam = np.empty(4, dtype=float)\n", "\n", "# Set axis 3 properties\n", "axsol[0].set_title(\"Height h at time = %s\" %time, fontsize=18)\n", "axsol[0].set_xlabel('x', fontsize=17)\n", "axsol[0].set_ylabel('h', fontsize=17)\n", "axsol[0].axis([-10,10,0,10])\n", "axsol[0].grid(alpha=0.1,color='k', linestyle='--')\n", "\n", "# Set axis 4 properties\n", "axsol[1].set_title(\"Momentum hu at time = %s\" %time, fontsize=18)\n", "axsol[1].set_xlabel('x', fontsize=17)\n", "axsol[1].set_ylabel('hu', fontsize=17)\n", "axsol[1].axis([-10,10,-15,15])\n", "axsol[1].grid(alpha=0.1,color='k', linestyle='--')\n", "\n", "\n", "# Remove defult mpld3 plugins\n", "plugins.clear(fig)\n", "\n", "# Call mpld3 custom PPLane plugin to interact with plot\n", "plugins.connect(fig, PPlaneNLPlugin(points[0],midpoint[0], g,iters,iter_charac,time,\n", " hugol[0],hugor[0],intcl[0],intcr[0],\n", " hugol2[0],hugor2[0],intcl2[0],intcr2[0],\n", " qlmarker[0],qmmarker[0],qrmarker[0],\n", " shock1[0],shock2[0],rar1[0],rar2[0],\n", " timedot[0],timeline[0],\n", " q1[0],q2[0]))\n", "\n", "mpld3.display()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 0 }