Add curvature calculation for sphere, sinc, ellipsoid and scaled sinc

bf82c992 · Niklas Bönisch · e4e491d2 · bf82c992 · bf82c992 · bf82c992
Commit bf82c992 authored 4 years ago by Niklas Bönisch
--- a/src/implicitSurfaces/curvature.ipynb
+++ b/src/implicitSurfaces/curvature.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from sympy import *\n",
+    "from sympy.printing import print_ccode\n",
+    "from sympy.parsing.sympy_parser import parse_expr\n",
+    "from sympy.plotting import plot\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "x0, x1, x2 = symbols(\"x0 x1 x2\")\n",
+    "c0,c1,c2 = symbols(\"O0:3\")\n",
+    "a0, a1, a2 = symbols(\"a0:3\")\n",
+    "R = symbols(\"R\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Curvature of an ellipsoid"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "-1 + (-O2 + x2)**2/a2**2 + (-O1 + x1)**2/a1**2 + (-O0 + x0)**2/a0**2\n",
+      "PYTHON:\n",
+      "1.0/(a2**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O2 + 2*x2)*(-O2 + x2)/(a2**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2)) + 1.0/(a1**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O1 + 2*x1)*(-O1 + x1)/(a1**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2)) + 1.0/(a0**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O0 + 2*x0)*(-O0 + x0)/(a0**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2))\n",
+      "C CODE:\n",
+      "1.0/(pow(a2, 2)*sqrt(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4))) - 2.0*(-2*O2 + 2*x2)*(-O2 + x2)/(pow(a2, 6)*pow(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4), 3.0/2.0)) + 1.0/(pow(a1, 2)*sqrt(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4))) - 2.0*(-2*O1 + 2*x1)*(-O1 + x1)/(pow(a1, 6)*pow(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4), 3.0/2.0)) + 1.0/(pow(a0, 2)*sqrt(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4))) - 2.0*(-2*O0 + 2*x0)*(-O0 + x0)/(pow(a0, 6)*pow(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4), 3.0/2.0))\n"
+     ]
+    },
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\frac{1.0}{a_{2}^{2} \\sqrt{\\frac{4 \\left(- O_{2} + x_{2}\\right)^{2}}{a_{2}^{4}} + \\frac{4 \\left(- O_{1} + x_{1}\\right)^{2}}{a_{1}^{4}} + \\frac{4 \\left(- O_{0} + x_{0}\\right)^{2}}{a_{0}^{4}}}} - \\frac{2.0 \\left(- 2 O_{2} + 2 x_{2}\\right) \\left(- O_{2} + x_{2}\\right)}{a_{2}^{6} \\left(\\frac{4 \\left(- O_{2} + x_{2}\\right)^{2}}{a_{2}^{4}} + \\frac{4 \\left(- O_{1} + x_{1}\\right)^{2}}{a_{1}^{4}} + \\frac{4 \\left(- O_{0} + x_{0}\\right)^{2}}{a_{0}^{4}}\\right)^{\\frac{3}{2}}} + \\frac{1.0}{a_{1}^{2} \\sqrt{\\frac{4 \\left(- O_{2} + x_{2}\\right)^{2}}{a_{2}^{4}} + \\frac{4 \\left(- O_{1} + x_{1}\\right)^{2}}{a_{1}^{4}} + \\frac{4 \\left(- O_{0} + x_{0}\\right)^{2}}{a_{0}^{4}}}} - \\frac{2.0 \\left(- 2 O_{1} + 2 x_{1}\\right) \\left(- O_{1} + x_{1}\\right)}{a_{1}^{6} \\left(\\frac{4 \\left(- O_{2} + x_{2}\\right)^{2}}{a_{2}^{4}} + \\frac{4 \\left(- O_{1} + x_{1}\\right)^{2}}{a_{1}^{4}} + \\frac{4 \\left(- O_{0} + x_{0}\\right)^{2}}{a_{0}^{4}}\\right)^{\\frac{3}{2}}} + \\frac{1.0}{a_{0}^{2} \\sqrt{\\frac{4 \\left(- O_{2} + x_{2}\\right)^{2}}{a_{2}^{4}} + \\frac{4 \\left(- O_{1} + x_{1}\\right)^{2}}{a_{1}^{4}} + \\frac{4 \\left(- O_{0} + x_{0}\\right)^{2}}{a_{0}^{4}}}} - \\frac{2.0 \\left(- 2 O_{0} + 2 x_{0}\\right) \\left(- O_{0} + x_{0}\\right)}{a_{0}^{6} \\left(\\frac{4 \\left(- O_{2} + x_{2}\\right)^{2}}{a_{2}^{4}} + \\frac{4 \\left(- O_{1} + x_{1}\\right)^{2}}{a_{1}^{4}} + \\frac{4 \\left(- O_{0} + x_{0}\\right)^{2}}{a_{0}^{4}}\\right)^{\\frac{3}{2}}}$"
+      ],
+      "text/plain": [
+       "1.0/(a2**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O2 + 2*x2)*(-O2 + x2)/(a2**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2)) + 1.0/(a1**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O1 + 2*x1)*(-O1 + x1)/(a1**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2)) + 1.0/(a0**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O0 + 2*x0)*(-O0 + x0)/(a0**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2))"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "a, b, c = symbols(\"a, b, c\")\n",
+    "ellipsoid = (x0-c0)**2 / a0**2 + (x1-c1)**2 / a1**2 + (x2-c2)**2 / a2**2 - 1\n",
+    "print(ellipsoid)\n",
+    "\n",
+    "grad = [simplify(diff(ellipsoid, x0)),\n",
+    "        simplify(diff(ellipsoid, x1)),\n",
+    "        simplify(diff(ellipsoid, x2))]\n",
+    "\n",
+    "# print(\"Gradient:    {}\".format(grad))\n",
+    "\n",
+    "\n",
+    "norm = sqrt(grad[0]**2+grad[1]**2+grad[2]**2)\n",
+    "grad_norm = [None, None, None]\n",
+    "grad_norm[0] = grad[0]/norm\n",
+    "grad_norm[1] = grad[1]/norm\n",
+    "grad_norm[2] = grad[2]/norm\n",
+    "# print(\"Norm:        {}\".format(norm))\n",
+    "# print(\"GradientNorm:{}\".format(grad_norm))\n",
+    "div0 = diff(grad_norm[0], x0)\n",
+    "# print('div0')\n",
+    "div1 = diff(grad_norm[1], x1)\n",
+    "# print('div1')\n",
+    "div2 = diff(grad_norm[2], x2)\n",
+    "# print('div2')\n",
+    "div = div0 + div1 + div2\n",
+    "# print('div')\n",
+    "# div = diff(grad_norm[0], x0) + diff(grad_norm[1], x1) + diff(grad_norm[2],x2)\n",
+    "# div = simplify(div)\n",
+    "div = 0.5*div\n",
+    "# print(\"DIV:         {}\".format(div))\n",
+    "k = div\n",
+    "print(\"PYTHON:\")\n",
+    "print(k)\n",
+    "print(\"C CODE:\")\n",
+    "print_ccode(k)\n",
+    "k"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Curvature of the sinc"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "PYTHON:\n",
+      "0.5*A*(O0 - x0)*(5*O0 - 5*x0)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(7/2)*sqrt(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)) + 0.5*A*(O0 - x0)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-A**2*(-2*O0 + 2*x0)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) - A**2*(O0 - x0)**2*(10*O0 - 10*x0)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**6) - A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-2*w**2*(-O0 + x0)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-3*O0 + 3*x0)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-O0 + x0)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + 2*(2*O0 - 2*x0)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) - A**2*(10*O0 - 10*x0)*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**6) - A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-2*w**2*(-O0 + x0)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-3*O0 + 3*x0)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-O0 + x0)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + 2*(2*O0 - 2*x0)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(5/2)*(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)**(3/2)) + 0.5*A*(O0 - x0)*(-w**2*(-O0 + x0)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + w*(-3*O0 + 3*x0)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + w*(-O0 + x0)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + (2*O0 - 2*x0)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(5/2)*sqrt(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)) + 0.5*A*(O1 - x1)*(5*O1 - 5*x1)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(7/2)*sqrt(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)) + 0.5*A*(O1 - x1)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-A**2*(O0 - x0)**2*(10*O1 - 10*x1)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**6) - A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-2*w**2*(-O1 + x1)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-3*O1 + 3*x1)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-O1 + x1)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + 2*(2*O1 - 2*x1)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) - A**2*(-2*O1 + 2*x1)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) - A**2*(O1 - x1)**2*(10*O1 - 10*x1)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**6) - A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-2*w**2*(-O1 + x1)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-3*O1 + 3*x1)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-O1 + x1)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + 2*(2*O1 - 2*x1)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(5/2)*(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)**(3/2)) + 0.5*A*(O1 - x1)*(-w**2*(-O1 + x1)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + w*(-3*O1 + 3*x1)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + w*(-O1 + x1)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + (2*O1 - 2*x1)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(5/2)*sqrt(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)) - 1.0*A*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(5/2)*sqrt(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1))\n",
+      "C CODE:\n",
+      "0.5*A*(O0 - x0)*(5*O0 - 5*x0)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 7.0/2.0)*sqrt(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1)) + 0.5*A*(O0 - x0)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-1.0/2.0*pow(A, 2)*(-2*O0 + 2*x0)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) - 1.0/2.0*pow(A, 2)*pow(O0 - x0, 2)*(10*O0 - 10*x0)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 6)) - 1.0/2.0*pow(A, 2)*pow(O0 - x0, 2)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-2*pow(w, 2)*(-O0 + x0)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-3*O0 + 3*x0)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-O0 + x0)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + 2*(2*O0 - 2*x0)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) - 1.0/2.0*pow(A, 2)*(10*O0 - 10*x0)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 6)) - 1.0/2.0*pow(A, 2)*pow(O1 - x1, 2)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-2*pow(w, 2)*(-O0 + x0)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-3*O0 + 3*x0)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-O0 + x0)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + 2*(2*O0 - 2*x0)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5.0/2.0)*pow(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1, 3.0/2.0)) + 0.5*A*(O0 - x0)*(-pow(w, 2)*(-O0 + x0)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + w*(-3*O0 + 3*x0)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + w*(-O0 + x0)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + (2*O0 - 2*x0)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5.0/2.0)*sqrt(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1)) + 0.5*A*(O1 - x1)*(5*O1 - 5*x1)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 7.0/2.0)*sqrt(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1)) + 0.5*A*(O1 - x1)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-1.0/2.0*pow(A, 2)*pow(O0 - x0, 2)*(10*O1 - 10*x1)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 6)) - 1.0/2.0*pow(A, 2)*pow(O0 - x0, 2)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-2*pow(w, 2)*(-O1 + x1)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-3*O1 + 3*x1)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-O1 + x1)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + 2*(2*O1 - 2*x1)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) - 1.0/2.0*pow(A, 2)*(-2*O1 + 2*x1)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) - 1.0/2.0*pow(A, 2)*pow(O1 - x1, 2)*(10*O1 - 10*x1)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 6)) - 1.0/2.0*pow(A, 2)*pow(O1 - x1, 2)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-2*pow(w, 2)*(-O1 + x1)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-3*O1 + 3*x1)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-O1 + x1)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + 2*(2*O1 - 2*x1)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5.0/2.0)*pow(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1, 3.0/2.0)) + 0.5*A*(O1 - x1)*(-pow(w, 2)*(-O1 + x1)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + w*(-3*O1 + 3*x1)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + w*(-O1 + x1)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + (2*O1 - 2*x1)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5.0/2.0)*sqrt(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1)) - 1.0*A*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5.0/2.0)*sqrt(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1))\n"
+     ]
+    }
+   ],
+   "source": [
+    "omega, amp = symbols('w A')\n",
+    "sinc = x2-c2-amp*sin(omega*sqrt((x0-c0)**2+(x1-c1)**2))/(omega*sqrt((x0-c0)**2+(x1-c1)**2))\n",
+    "grad = [simplify(diff(sinc, x0)),\n",
+    "        simplify(diff(sinc, x1)),\n",
+    "        simplify(diff(sinc, x2))]\n",
+    "\n",
+    "# print(\"Gradient:    {}\".format(grad))\n",
+    "\n",
+    "\n",
+    "norm = sqrt(grad[0]**2+grad[1]**2+grad[2]**2)\n",
+    "grad_norm = [None, None, None]\n",
+    "grad_norm[0] = grad[0]/norm\n",
+    "grad_norm[1] = grad[1]/norm\n",
+    "grad_norm[2] = grad[2]/norm\n",
+    "# print(\"Norm:        {}\".format(norm))\n",
+    "# print(\"GradientNorm:{}\".format(grad_norm))\n",
+    "div0 = diff(grad_norm[0], x0)\n",
+    "# print('div0')\n",
+    "div1 = diff(grad_norm[1], x1)\n",
+    "# print('div1')\n",
+    "div2 = diff(grad_norm[2], x2)\n",
+    "# print('div2')\n",
+    "div = div0 + div1 + div2\n",
+    "# print('div')\n",
+    "div = 0.5*div\n",
+    "# print(\"DIV:         {}\".format(div))\n",
+    "k = div\n",
+    "print(\"PYTHON:\")\n",
+    "print(k)\n",
+    "print(\"C CODE:\")\n",
+    "print_ccode(k)\n",
+    "# k"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Curvature of the scaled sinc"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.8.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
+%% Cell type:code id: tags:
+``` python
+from sympy import *
+from sympy.printing import print_ccode
+from sympy.parsing.sympy_parser import parse_expr
+from sympy.plotting import plot
+import matplotlib.pyplot as plt
+import numpy as np
+x0, x1, x2 = symbols("x0 x1 x2")
+c0,c1,c2 = symbols("O0:3")
+a0, a1, a2 = symbols("a0:3")
+R = symbols("R")
+```
+%% Cell type:markdown id: tags:
+# Curvature of an ellipsoid
+%% Cell type:code id: tags:
+``` python
+a, b, c = symbols("a, b, c")
+ellipsoid = (x0-c0)**2 / a0**2 + (x1-c1)**2 / a1**2 + (x2-c2)**2 / a2**2 - 1
+print(ellipsoid)
+grad = [simplify(diff(ellipsoid, x0)),
+        simplify(diff(ellipsoid, x1)),
+        simplify(diff(ellipsoid, x2))]
+# print("Gradient:    {}".format(grad))
+norm = sqrt(grad[0]**2+grad[1]**2+grad[2]**2)
+grad_norm = [None, None, None]
+grad_norm[0] = grad[0]/norm
+grad_norm[1] = grad[1]/norm
+grad_norm[2] = grad[2]/norm
+# print("Norm:        {}".format(norm))
+# print("GradientNorm:{}".format(grad_norm))
+div0 = diff(grad_norm[0], x0)
+# print('div0')
+div1 = diff(grad_norm[1], x1)
+# print('div1')
+div2 = diff(grad_norm[2], x2)
+# print('div2')
+div = div0 + div1 + div2
+# print('div')
+# div = diff(grad_norm[0], x0) + diff(grad_norm[1], x1) + diff(grad_norm[2],x2)
+# div = simplify(div)
+div = 0.5*div
+# print("DIV:         {}".format(div))
+k = div
+print("PYTHON:")
+print(k)
+print("C CODE:")
+print_ccode(k)
+k
+```
+%% Output
+    -1 + (-O2 + x2)**2/a2**2 + (-O1 + x1)**2/a1**2 + (-O0 + x0)**2/a0**2
+    PYTHON:
+    1.0/(a2**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O2 + 2*x2)*(-O2 + x2)/(a2**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2)) + 1.0/(a1**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O1 + 2*x1)*(-O1 + x1)/(a1**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2)) + 1.0/(a0**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O0 + 2*x0)*(-O0 + x0)/(a0**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2))
+    C CODE:
+    1.0/(pow(a2, 2)*sqrt(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4))) - 2.0*(-2*O2 + 2*x2)*(-O2 + x2)/(pow(a2, 6)*pow(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4), 3.0/2.0)) + 1.0/(pow(a1, 2)*sqrt(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4))) - 2.0*(-2*O1 + 2*x1)*(-O1 + x1)/(pow(a1, 6)*pow(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4), 3.0/2.0)) + 1.0/(pow(a0, 2)*sqrt(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4))) - 2.0*(-2*O0 + 2*x0)*(-O0 + x0)/(pow(a0, 6)*pow(4*pow(-O2 + x2, 2)/pow(a2, 4) + 4*pow(-O1 + x1, 2)/pow(a1, 4) + 4*pow(-O0 + x0, 2)/pow(a0, 4), 3.0/2.0))
+    $\displaystyle \frac{1.0}{a_{2}^{2} \sqrt{\frac{4 \left(- O_{2} + x_{2}\right)^{2}}{a_{2}^{4}} + \frac{4 \left(- O_{1} + x_{1}\right)^{2}}{a_{1}^{4}} + \frac{4 \left(- O_{0} + x_{0}\right)^{2}}{a_{0}^{4}}}} - \frac{2.0 \left(- 2 O_{2} + 2 x_{2}\right) \left(- O_{2} + x_{2}\right)}{a_{2}^{6} \left(\frac{4 \left(- O_{2} + x_{2}\right)^{2}}{a_{2}^{4}} + \frac{4 \left(- O_{1} + x_{1}\right)^{2}}{a_{1}^{4}} + \frac{4 \left(- O_{0} + x_{0}\right)^{2}}{a_{0}^{4}}\right)^{\frac{3}{2}}} + \frac{1.0}{a_{1}^{2} \sqrt{\frac{4 \left(- O_{2} + x_{2}\right)^{2}}{a_{2}^{4}} + \frac{4 \left(- O_{1} + x_{1}\right)^{2}}{a_{1}^{4}} + \frac{4 \left(- O_{0} + x_{0}\right)^{2}}{a_{0}^{4}}}} - \frac{2.0 \left(- 2 O_{1} + 2 x_{1}\right) \left(- O_{1} + x_{1}\right)}{a_{1}^{6} \left(\frac{4 \left(- O_{2} + x_{2}\right)^{2}}{a_{2}^{4}} + \frac{4 \left(- O_{1} + x_{1}\right)^{2}}{a_{1}^{4}} + \frac{4 \left(- O_{0} + x_{0}\right)^{2}}{a_{0}^{4}}\right)^{\frac{3}{2}}} + \frac{1.0}{a_{0}^{2} \sqrt{\frac{4 \left(- O_{2} + x_{2}\right)^{2}}{a_{2}^{4}} + \frac{4 \left(- O_{1} + x_{1}\right)^{2}}{a_{1}^{4}} + \frac{4 \left(- O_{0} + x_{0}\right)^{2}}{a_{0}^{4}}}} - \frac{2.0 \left(- 2 O_{0} + 2 x_{0}\right) \left(- O_{0} + x_{0}\right)}{a_{0}^{6} \left(\frac{4 \left(- O_{2} + x_{2}\right)^{2}}{a_{2}^{4}} + \frac{4 \left(- O_{1} + x_{1}\right)^{2}}{a_{1}^{4}} + \frac{4 \left(- O_{0} + x_{0}\right)^{2}}{a_{0}^{4}}\right)^{\frac{3}{2}}}$
+    1.0/(a2**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O2 + 2*x2)*(-O2 + x2)/(a2**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2)) + 1.0/(a1**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O1 + 2*x1)*(-O1 + x1)/(a1**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2)) + 1.0/(a0**2*sqrt(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)) - 2.0*(-2*O0 + 2*x0)*(-O0 + x0)/(a0**6*(4*(-O2 + x2)**2/a2**4 + 4*(-O1 + x1)**2/a1**4 + 4*(-O0 + x0)**2/a0**4)**(3/2))
+%% Cell type:markdown id: tags:
+# Curvature of the sinc
+%% Cell type:code id: tags:
+``` python
+omega, amp = symbols('w A')
+sinc = x2-c2-amp*sin(omega*sqrt((x0-c0)**2+(x1-c1)**2))/(omega*sqrt((x0-c0)**2+(x1-c1)**2))
+grad = [simplify(diff(sinc, x0)),
+        simplify(diff(sinc, x1)),
+        simplify(diff(sinc, x2))]
+# print("Gradient:    {}".format(grad))
+norm = sqrt(grad[0]**2+grad[1]**2+grad[2]**2)
+grad_norm = [None, None, None]
+grad_norm[0] = grad[0]/norm
+grad_norm[1] = grad[1]/norm
+grad_norm[2] = grad[2]/norm
+# print("Norm:        {}".format(norm))
+# print("GradientNorm:{}".format(grad_norm))
+div0 = diff(grad_norm[0], x0)
+# print('div0')
+div1 = diff(grad_norm[1], x1)
+# print('div1')
+div2 = diff(grad_norm[2], x2)
+# print('div2')
+div = div0 + div1 + div2
+# print('div')
+div = 0.5*div
+# print("DIV:         {}".format(div))
+k = div
+print("PYTHON:")
+print(k)
+print("C CODE:")
+print_ccode(k)
+# k
+```
+%% Output
+    PYTHON:
+    0.5*A*(O0 - x0)*(5*O0 - 5*x0)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(7/2)*sqrt(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)) + 0.5*A*(O0 - x0)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-A**2*(-2*O0 + 2*x0)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) - A**2*(O0 - x0)**2*(10*O0 - 10*x0)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**6) - A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-2*w**2*(-O0 + x0)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-3*O0 + 3*x0)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-O0 + x0)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + 2*(2*O0 - 2*x0)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) - A**2*(10*O0 - 10*x0)*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**6) - A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-2*w**2*(-O0 + x0)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-3*O0 + 3*x0)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-O0 + x0)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + 2*(2*O0 - 2*x0)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(5/2)*(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)**(3/2)) + 0.5*A*(O0 - x0)*(-w**2*(-O0 + x0)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + w*(-3*O0 + 3*x0)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + w*(-O0 + x0)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + (2*O0 - 2*x0)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(5/2)*sqrt(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)) + 0.5*A*(O1 - x1)*(5*O1 - 5*x1)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(7/2)*sqrt(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)) + 0.5*A*(O1 - x1)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-A**2*(O0 - x0)**2*(10*O1 - 10*x1)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**6) - A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-2*w**2*(-O1 + x1)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-3*O1 + 3*x1)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-O1 + x1)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + 2*(2*O1 - 2*x1)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) - A**2*(-2*O1 + 2*x1)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) - A**2*(O1 - x1)**2*(10*O1 - 10*x1)*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**6) - A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))*(-2*w**2*(-O1 + x1)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-3*O1 + 3*x1)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + 2*w*(-O1 + x1)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + 2*(2*O1 - 2*x1)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(2*w**2*((O0 - x0)**2 + (O1 - x1)**2)**5))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(5/2)*(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)**(3/2)) + 0.5*A*(O1 - x1)*(-w**2*(-O1 + x1)*((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + w*(-3*O1 + 3*x1)*sqrt((O0 - x0)**2 + (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) + w*(-O1 + x1)*(-(O0 - x0)**2 - (O1 - x1)**2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2))/sqrt((O0 - x0)**2 + (O1 - x1)**2) + (2*O1 - 2*x1)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(5/2)*sqrt(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1)) - 1.0*A*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))/(w*((O0 - x0)**2 + (O1 - x1)**2)**(5/2)*sqrt(A**2*(O0 - x0)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + A**2*(O1 - x1)**2*(w*((O0 - x0)**2 + (O1 - x1)**2)**(3/2)*cos(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)) - ((O0 - x0)**2 + (O1 - x1)**2)*sin(w*sqrt((O0 - x0)**2 + (O1 - x1)**2)))**2/(w**2*((O0 - x0)**2 + (O1 - x1)**2)**5) + 1))
+    C CODE:
+    0.5*A*(O0 - x0)*(5*O0 - 5*x0)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 7.0/2.0)*sqrt(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1)) + 0.5*A*(O0 - x0)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-1.0/2.0*pow(A, 2)*(-2*O0 + 2*x0)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) - 1.0/2.0*pow(A, 2)*pow(O0 - x0, 2)*(10*O0 - 10*x0)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 6)) - 1.0/2.0*pow(A, 2)*pow(O0 - x0, 2)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-2*pow(w, 2)*(-O0 + x0)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-3*O0 + 3*x0)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-O0 + x0)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + 2*(2*O0 - 2*x0)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) - 1.0/2.0*pow(A, 2)*(10*O0 - 10*x0)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 6)) - 1.0/2.0*pow(A, 2)*pow(O1 - x1, 2)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-2*pow(w, 2)*(-O0 + x0)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-3*O0 + 3*x0)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-O0 + x0)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + 2*(2*O0 - 2*x0)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5.0/2.0)*pow(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1, 3.0/2.0)) + 0.5*A*(O0 - x0)*(-pow(w, 2)*(-O0 + x0)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + w*(-3*O0 + 3*x0)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + w*(-O0 + x0)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + (2*O0 - 2*x0)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5.0/2.0)*sqrt(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1)) + 0.5*A*(O1 - x1)*(5*O1 - 5*x1)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 7.0/2.0)*sqrt(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1)) + 0.5*A*(O1 - x1)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-1.0/2.0*pow(A, 2)*pow(O0 - x0, 2)*(10*O1 - 10*x1)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 6)) - 1.0/2.0*pow(A, 2)*pow(O0 - x0, 2)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-2*pow(w, 2)*(-O1 + x1)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-3*O1 + 3*x1)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-O1 + x1)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + 2*(2*O1 - 2*x1)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) - 1.0/2.0*pow(A, 2)*(-2*O1 + 2*x1)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) - 1.0/2.0*pow(A, 2)*pow(O1 - x1, 2)*(10*O1 - 10*x1)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 6)) - 1.0/2.0*pow(A, 2)*pow(O1 - x1, 2)*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))*(-2*pow(w, 2)*(-O1 + x1)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-3*O1 + 3*x1)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + 2*w*(-O1 + x1)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + 2*(2*O1 - 2*x1)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5.0/2.0)*pow(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1, 3.0/2.0)) + 0.5*A*(O1 - x1)*(-pow(w, 2)*(-O1 + x1)*(pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + w*(-3*O1 + 3*x1)*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) + w*(-O1 + x1)*(-pow(O0 - x0, 2) - pow(O1 - x1, 2))*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)))/sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2)) + (2*O1 - 2*x1)*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5.0/2.0)*sqrt(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1)) - 1.0*A*(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))))/(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5.0/2.0)*sqrt(pow(A, 2)*pow(O0 - x0, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + pow(A, 2)*pow(O1 - x1, 2)*pow(w*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 3.0/2.0)*cos(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))) - (pow(O0 - x0, 2) + pow(O1 - x1, 2))*sin(w*sqrt(pow(O0 - x0, 2) + pow(O1 - x1, 2))), 2)/(pow(w, 2)*pow(pow(O0 - x0, 2) + pow(O1 - x1, 2), 5)) + 1))
+%% Cell type:markdown id: tags:
+# Curvature of the scaled sinc
+%% Cell type:code id: tags:
+``` python
+```
--- a/src/implicitSurfaces/implicitSurfaces.C
+++ b/src/implicitSurfaces/implicitSurfaces.C
--- a/src/implicitSurfaces/implicitSurfaces.H
+++ b/src/implicitSurfaces/implicitSurfaces.H
@@ -71,6 +71,8 @@ namespace Foam {
            virtual scalar value(const vector&) const = 0;
            virtual scalar operator()(const vector&) const = 0;
            virtual vector grad(const vector&) const = 0;
+            virtual scalar curvature(const vector&) const = 0;
    };
    class plane : public implicitSurface
@@ -96,6 +98,8 @@ namespace Foam {
            virtual vector grad(const vector& x) const;
+            virtual scalar curvature(const vector& x) const;
            vector position() const;
            vector normal() const;
@@ -124,6 +128,8 @@ namespace Foam {
            virtual vector grad(const vector& x) const;
+            virtual scalar curvature(const vector& x) const;
            vector center() const;
            scalar radius() const;
@@ -153,6 +159,8 @@ namespace Foam {
            virtual vector grad(const vector& x) const;
+            virtual scalar curvature(const vector& x) const;
            vector center() const;
            vector axes() const;
@@ -180,6 +188,8 @@ namespace Foam {
            virtual vector grad(const vector& x) const;
+            virtual scalar curvature(const vector& x) const;
            vector origin() const;
            scalar amplitude() const;
@@ -209,6 +219,8 @@ namespace Foam {
            virtual vector grad(const vector& x) const;
+            virtual scalar curvature(const vector& x) const;
            vector origin() const;
            scalar amplitude() const;