{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Formulas for a 4-bar Mechanism"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This worksheet computes the formulas needed for a Graphical User Interface,\n",
    "to model a 4-bar mechanism.  The mechanism is driven by a crank.\n",
    "The equations express the conditions that the 4-bar mechanism remains rigid as we turn the crank."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Defining a System of Polynomial Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 4-bar mechanism has four parameters:\n",
    "\n",
    "1. The first bar starts at (0,0) and ends at $(a,0)$, for some $a > 0$.\n",
    "\n",
    "2. The length of the crank is $L$, the crank is the second bar.\n",
    "\n",
    "3. The third bar connects the crank to the point $(x, y)$ and has length $r$.\n",
    "\n",
    "4. The fourth bar connects the point $(x, y)$ to $(a,0)$ and has length $R$.\n",
    "\n",
    "The polynomials in the system are then derived as follows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. The crank $C$ of length $L$ and angle $t$ \n",
    "   has coordinates $(L \\cos(t), L \\sin(t))$.\n",
    "   \n",
    "2. The point $D(x,y)$ is at distance $r$ from $B(a,0)$,\n",
    "   expressed by the equation $f$,\n",
    "   $f = (x-a)^2 + y^2 - r^2 = 0$.\n",
    "   \n",
    "3. The point $D(x,y)$ is at distance $R$ from \n",
    "   $C(L \\cos(t),L \\sin(t))$, expressed by\n",
    "   the equation $g$, \n",
    "   $g = (x - L \\cos(t))^2 + (y - L \\sin(t))^2 - R^2 = 0$.\n",
    "   \n",
    "4. Replacing $\\cos(t)$ and $\\sin(t)$ by $c$ and $s$ respectively,\n",
    "   adding $h = c^2 + s^2 - 1 = 0$ we obtain a system of algebraic equations\n",
    "   in $(x,y)$ and with parameters $(a,r,R,L)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(a - x)^2 - r^2 + y^2, (L*c - x)^2 + (L*s - y)^2 - R^2, c^2 + s^2 - 1]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x,y,a,r,t,L,R = var('x,y,a,r,t,L,R')\n",
    "f = (x - a)^2 + y^2 - r^2\n",
    "g = (x - L*cos(t))^2 + (y - L*sin(t))^2 - R^2\n",
    "c,s = var('c,s')\n",
    "h = c^2 + s^2 - 1\n",
    "gcs = g.subs(cos(t) == c).subs(sin(t) == s)\n",
    "[f,gcs,h]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'sage.symbolic.expression.Expression'>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Computing a Groebner basis with Singular"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now go to singular, declaring a polynomial ring."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Multivariate Polynomial Ring in x, y, a, r, L, R, c, s over Rational Field"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = PolynomialRing(QQ,8,names='xyarLRcs', order='lex')\n",
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We cast the original polynomials into the ring:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[x^2 - 2*x*a + y^2 + a^2 - r^2, x^2 - 2*x*L*c + y^2 - 2*y*L*s + L^2*c^2 + L^2*s^2 - R^2, c^2 + s^2 - 1]\n"
     ]
    }
   ],
   "source": [
    "F = P(f); G = P(gcs); H = P(h); print([F,G,H])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute a lexicographic Groebner basis with Singular"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[x^2 - 2*x*L*c + y^2 - 2*y*L*s + L^2 - R^2, x*y*L*s - 1/2*x*r^2 - 1/2*x*L^2*s^2 + 1/2*x*R^2 + y^2*a - y^2*L*c - 3/2*y*a*L*s + 1/2*y*L^2*c*s + 1/4*a^3 - 3/4*a^2*L*c - 1/4*a*r^2 + 3/4*a*L^2 - 3/4*a*R^2 + 3/4*r^2*L*c - 1/4*L^3*c + 1/4*L*R^2*c, x*a - x*L*c - y*L*s - 1/2*a^2 + 1/2*r^2 + 1/2*L^2 - 1/2*R^2, y^2*a^2 - 2*y^2*a*L*c + y^2*L^2 - y*a^2*L*s + 2*y*a*L^2*c*s - y*r^2*L*s - y*L^3*s + y*L*R^2*s + 1/4*a^4 - a^3*L*c - 1/2*a^2*r^2 - a^2*L^2*s^2 + 3/2*a^2*L^2 - 1/2*a^2*R^2 + a*r^2*L*c - a*L^3*c + a*L*R^2*c + 1/4*r^4 + r^2*L^2*s^2 - 1/2*r^2*L^2 - 1/2*r^2*R^2 + 1/4*L^4 - 1/2*L^2*R^2 + 1/4*R^4, c^2 + s^2 - 1]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I = P.ideal(F,G,H)\n",
    "gb = I.groebner_basis('libsingular:slimgb')\n",
    "gb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Processing the Groebner basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second element in the lexicographic Groebner basis is\n",
    "linear in x and we can express the x coordinate in terms of y\n",
    "and the other parameters of the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gb[2] = x*a - x*L*c - y*L*s - 1/2*a^2 + 1/2*r^2 + 1/2*L^2 - 1/2*R^2\n",
      "<class 'sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular'>\n"
     ]
    }
   ],
   "source": [
    "print('gb[2] =', gb[2])\n",
    "print(type(gb[2]))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "To solve the expression for x, we have to convert the singular\n",
    "polynomial back to a symbolic Sage expression.\n",
    "This happens with the application of the `SR()`\n",
    "where `SR` abbreviates `Symbolic Ring`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-L*c*x - L*s*y + 1/2*L^2 - 1/2*R^2 - 1/2*a^2 + 1/2*r^2 + a*x <class 'sage.symbolic.expression.Expression'>\n"
     ]
    }
   ],
   "source": [
    "px = SR(gb[2])\n",
    "print(px, type(px))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[\n",
      "x == -1/2*(2*L*s*y - L^2 + R^2 + a^2 - r^2)/(L*c - a)\n",
      "]\n",
      "x = -1/2*(2*L*s*y - L^2 + R^2 + a^2 - r^2)/(L*c - a)\n"
     ]
    }
   ],
   "source": [
    "xsol = solve([px==0],x); print(xsol)\n",
    "print('x =', xsol[0].rhs())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The third element in the lexicographic Groebner basis depends only\n",
    "on y and the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "y^2*a^2 - 2*y^2*a*L*c + y^2*L^2 - y*a^2*L*s + 2*y*a*L^2*c*s - y*r^2*L*s - y*L^3*s + y*L*R^2*s + 1/4*a^4 - a^3*L*c - 1/2*a^2*r^2 - a^2*L^2*s^2 + 3/2*a^2*L^2 - 1/2*a^2*R^2 + a*r^2*L*c - a*L^3*c + a*L*R^2*c + 1/4*r^4 + r^2*L^2*s^2 - 1/2*r^2*L^2 - 1/2*r^2*R^2 + 1/4*L^4 - 1/2*L^2*R^2 + 1/4*R^4"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gb[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We convert again to a symbolic Sage expression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y = 1/2*((2*L^2*a*c - L^3 + L*R^2 - L*a^2 - L*r^2)*s + sqrt(-L^6 + 2*L^4*R^2 - L^2*R^4 - a^6 - (7*L^2 - 2*R^2)*a^4 + (2*L*a*c - L^2 - a^2)*r^4 - (7*L^4 - 4*L^2*R^2 + R^4)*a^2 - 8*(L^2*a^4 + (L^4 - L^2*R^2)*a^2)*c^2 + 2*(4*L^2*a^2*c^2 + L^4 + L^2*R^2 + a^4 + (2*L^2 + R^2)*a^2 - 2*(2*L*a^3 + (2*L^3 + L*R^2)*a)*c)*r^2 + (4*L^4*a^2*c^2 + L^6 - 2*L^4*R^2 + L^2*R^4 + 5*L^2*a^4 + L^2*r^4 + 2*(3*L^4 - L^2*R^2)*a^2 + 2*(2*L^3*a*c - L^4 - L^2*R^2 - L^2*a^2)*r^2 - 4*(3*L^3*a^3 + (L^5 - L^3*R^2)*a)*c)*s^2 + 2*(3*L*a^5 + 2*(5*L^3 - 2*L*R^2)*a^3 + (3*L^5 - 4*L^3*R^2 + L*R^4)*a)*c))/(2*L*a*c - L^2 - a^2)\n"
     ]
    }
   ],
   "source": [
    "py = SR(gb[3])\n",
    "ysol = solve([py==0],y); print('y =', ysol[0].rhs())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We extract the discriminant of the solution for y:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1/2*((2*L^2*a*c - L^3 + L*R^2 - L*a^2 - L*r^2)*s + sqrt(-L^6 + 2*L^4*R^2 - L^2*R^4 - a^6 - (7*L^2 - 2*R^2)*a^4 + (2*L*a*c - L^2 - a^2)*r^4 - (7*L^4 - 4*L^2*R^2 + R^4)*a^2 - 8*(L^2*a^4 + (L^4 - L^2*R^2)*a^2)*c^2 + 2*(4*L^2*a^2*c^2 + L^4 + L^2*R^2 + a^4 + (2*L^2 + R^2)*a^2 - 2*(2*L*a^3 + (2*L^3 + L*R^2)*a)*c)*r^2 + (4*L^4*a^2*c^2 + L^6 - 2*L^4*R^2 + L^2*R^4 + 5*L^2*a^4 + L^2*r^4 + 2*(3*L^4 - L^2*R^2)*a^2 + 2*(2*L^3*a*c - L^4 - L^2*R^2 - L^2*a^2)*r^2 - 4*(3*L^3*a^3 + (L^5 - L^3*R^2)*a)*c)*s^2 + 2*(3*L*a^5 + 2*(5*L^3 - 2*L*R^2)*a^3 + (3*L^5 - 4*L^3*R^2 + L*R^4)*a)*c))/(2*L*a*c - L^2 - a^2)\n",
      "<function mul_vararg at 0x7f77c327ba30>\n",
      "[1/(2*L*a*c - L^2 - a^2), (2*L^2*a*c - L^3 + L*R^2 - L*a^2 - L*r^2)*s + sqrt(-L^6 + 2*L^4*R^2 - L^2*R^4 - a^6 - (7*L^2 - 2*R^2)*a^4 + (2*L*a*c - L^2 - a^2)*r^4 - (7*L^4 - 4*L^2*R^2 + R^4)*a^2 - 8*(L^2*a^4 + (L^4 - L^2*R^2)*a^2)*c^2 + 2*(4*L^2*a^2*c^2 + L^4 + L^2*R^2 + a^4 + (2*L^2 + R^2)*a^2 - 2*(2*L*a^3 + (2*L^3 + L*R^2)*a)*c)*r^2 + (4*L^4*a^2*c^2 + L^6 - 2*L^4*R^2 + L^2*R^4 + 5*L^2*a^4 + L^2*r^4 + 2*(3*L^4 - L^2*R^2)*a^2 + 2*(2*L^3*a*c - L^4 - L^2*R^2 - L^2*a^2)*r^2 - 4*(3*L^3*a^3 + (L^5 - L^3*R^2)*a)*c)*s^2 + 2*(3*L*a^5 + 2*(5*L^3 - 2*L*R^2)*a^3 + (3*L^5 - 4*L^3*R^2 + L*R^4)*a)*c), 1/2]\n"
     ]
    }
   ],
   "source": [
    "y0 = ysol[0].rhs(); print(y0)\n",
    "print(y0.operator()); print(y0.operands())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "s = (2*L^2*a*c - L^3 + L*R^2 - L*a^2 - L*r^2)*s + sqrt(-L^6 + 2*L^4*R^2 - L^2*R^4 - a^6 - (7*L^2 - 2*R^2)*a^4 + (2*L*a*c - L^2 - a^2)*r^4 - (7*L^4 - 4*L^2*R^2 + R^4)*a^2 - 8*(L^2*a^4 + (L^4 - L^2*R^2)*a^2)*c^2 + 2*(4*L^2*a^2*c^2 + L^4 + L^2*R^2 + a^4 + (2*L^2 + R^2)*a^2 - 2*(2*L*a^3 + (2*L^3 + L*R^2)*a)*c)*r^2 + (4*L^4*a^2*c^2 + L^6 - 2*L^4*R^2 + L^2*R^4 + 5*L^2*a^4 + L^2*r^4 + 2*(3*L^4 - L^2*R^2)*a^2 + 2*(2*L^3*a*c - L^4 - L^2*R^2 - L^2*a^2)*r^2 - 4*(3*L^3*a^3 + (L^5 - L^3*R^2)*a)*c)*s^2 + 2*(3*L*a^5 + 2*(5*L^3 - 2*L*R^2)*a^3 + (3*L^5 - 4*L^3*R^2 + L*R^4)*a)*c)\n"
     ]
    }
   ],
   "source": [
    "s = str(y0.operands()[1]); print('s =', s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L = ['(2*L^2*a*c - L^3 + L*R^2 - L*a^2 - L*r^2)*s + ', '(-L^6 + 2*L^4*R^2 - L^2*R^4 - a^6 - (7*L^2 - 2*R^2)*a^4 + (2*L*a*c - L^2 - a^2)*r^4 - (7*L^4 - 4*L^2*R^2 + R^4)*a^2 - 8*(L^2*a^4 + (L^4 - L^2*R^2)*a^2)*c^2 + 2*(4*L^2*a^2*c^2 + L^4 + L^2*R^2 + a^4 + (2*L^2 + R^2)*a^2 - 2*(2*L*a^3 + (2*L^3 + L*R^2)*a)*c)*r^2 + (4*L^4*a^2*c^2 + L^6 - 2*L^4*R^2 + L^2*R^4 + 5*L^2*a^4 + L^2*r^4 + 2*(3*L^4 - L^2*R^2)*a^2 + 2*(2*L^3*a*c - L^4 - L^2*R^2 - L^2*a^2)*r^2 - 4*(3*L^3*a^3 + (L^5 - L^3*R^2)*a)*c)*s^2 + 2*(3*L*a^5 + 2*(5*L^3 - 2*L*R^2)*a^3 + (3*L^5 - 4*L^3*R^2 + L*R^4)*a)*c)']\n"
     ]
    }
   ],
   "source": [
    "L = s.split('sqrt'); print('L =', L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "discriminant (as string) = (-L^6 + 2*L^4*R^2 - L^2*R^4 - a^6 - (7*L^2 - 2*R^2)*a^4 + (2*L*a*c - L^2 - a^2)*r^4 - (7*L^4 - 4*L^2*R^2 + R^4)*a^2 - 8*(L^2*a^4 + (L^4 - L^2*R^2)*a^2)*c^2 + 2*(4*L^2*a^2*c^2 + L^4 + L^2*R^2 + a^4 + (2*L^2 + R^2)*a^2 - 2*(2*L*a^3 + (2*L^3 + L*R^2)*a)*c)*r^2 + (4*L^4*a^2*c^2 + L^6 - 2*L^4*R^2 + L^2*R^4 + 5*L^2*a^4 + L^2*r^4 + 2*(3*L^4 - L^2*R^2)*a^2 + 2*(2*L^3*a*c - L^4 - L^2*R^2 - L^2*a^2)*r^2 - 4*(3*L^3*a^3 + (L^5 - L^3*R^2)*a)*c)*s^2 + 2*(3*L*a^5 + 2*(5*L^3 - 2*L*R^2)*a^3 + (3*L^5 - 4*L^3*R^2 + L*R^4)*a)*c)\n"
     ]
    }
   ],
   "source": [
    "d = L[1]; print('discriminant (as string) =', d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.5",
   "language": "sage",
   "name": "sagemath"
  },
  "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}