(***********************************************************************

                    Mathematica-Compatible Notebook

This notebook can be used on any computer system with Mathematica 3.0,
MathReader 3.0, or any compatible application. The data for the notebook 
starts with the line of stars above.

To get the notebook into a Mathematica-compatible application, do one of 
the following:

* Save the data starting with the line of stars above into a file
  with a name ending in .nb, then open the file inside the application;

* Copy the data starting with the line of stars above to the
  clipboard, then use the Paste menu command inside the application.

Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode.  Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).

NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing the 
word CacheID, otherwise Mathematica-compatible applications may try to 
use invalid cache data.

For more information on notebooks and Mathematica-compatible 
applications, contact Wolfram Research:
  web: http://www.wolfram.com
  email: info@wolfram.com
  phone: +1-217-398-0700 (U.S.)

Notebook reader applications are available free of charge from 
Wolfram Research.
***********************************************************************)

(*CacheID: 232*)


(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[     49937,       1077]*)
(*NotebookOutlinePosition[     50776,       1104]*)
(*  CellTagsIndexPosition[     50732,       1100]*)
(*WindowFrame->Normal*)



Notebook[{
Cell[BoxData[
    RowBox[{
      RowBox[{"(*", " ", 
        RowBox[{\(lagrunconstrained1  .0\), "  ", "\n", "\t", "\n", 
          RowBox[{
            StyleBox["Mathematica",
              FontSlant->"Italic"], " ", "3.05"}], "\n", 
          \(D.\ B.\ Siano\ \ \(email : \ \ dimona@home.com\), \ \ Jan.\ 2, \ 
          1999\)}], "     ", "*)"}], "\n", "\t", "\n", 
      \( (*\tGets\ the\ Lagrange\ equations\ for\ the\ unconstrained\ movement
            \ of\ a\ trebuchet\ for\ input\ to\ a\ Runge - 
          Kutta\ \((or\ other\ method)\)\ of\ solving\ the\ differential\ 
            equations\ numerically\ in\ some\ other, \ 
        speedier\ language.\ First\ sets\ up\ the\ Lagrangian\  = \ 
          \(L\  = \ 
            \(T - V\  = \ KE - PE\ \ for\ the\ general\ trebuchet\)\), \ 
        which\ includes\ a\ beam, \ counterweight, \ and\ sling, \ 
        then\ uses\ Lagranges\ equations\ to\ get\ the\ equations\ of\ motion
          \ which\ are\ solved\ numerically\ for\ the\ times\ of\ 
          \(interest.\)\ *) \), " ", "\n", "\t", "\n", "\t", 
      \( (*\ Three\ joined\ rigid\ rods\ with\ two\ masses\ at\ either\ end\n
        have\ a\ pivot\ in\ the\ middle\ one; \ 
        they\ are\ released\ with\ an\ arbitrary\ initial\ configuration\ of\ 
          position\ and\ velocity.\n\nThe\ short\ arm\ of\ the\ pivoted\ 
          central\ rod\ \((the\ beam)\)\ is\ l1\ long, \ the\n
        long\ arm\ has\ a\ length\ l2; \ \ 
        The\ counterweight\ is\ at\ the\ end\ of\ an\ arm\ l4\ long, \ the\ \n
        
        projectile\ is\ at\ the\ end\ of\ the\ arm\ \((sling)\)\ l3\ 
          long.\ The\ beam\ and\ arms\ are\ assumed\ to\ be\ rigid\n
        \tm1\ is\ the\ mass\ of\ the\ counterweight, 
        m2\ the\ mass\ of\ the\ projectile.\ \ These\ are\ taken\ to\ be\ 
          point\ masses.\ The\ mass\ of\ the\ beam\ is\ \(mb.\)*) \), "\n", 
      "\n", \( (*There\ are\ three\ degrees\ of\ freedom\ to\ the\ problem, \ 
        which\ are\ the\ three\ angles\ required\ to\ specify\ an\ arbitrary\ 
          \(configuration : \n\t\t\n
            th\ is\ the\ angle\ the\ middle\ arm\ makes\ with\ the\ vertical
              \); \n
        phi\ is\ the\ angle\ the\ counterweight\ arm\ makes\ with\ the\n
        pivoted\ arm; \n\ psi\ is\ the\ angle\ the\ projectile\ arm\n
        makes\ with\ the\ pivoted\ arm.\ \nIn\ the\ usual\ initial\ 
          configuration, \ th\ is\ an\ angle\ greater\ than\ 90\ degrees, \ 
        while\ phi\ and\ psi\ are\ acute.\ \n\nThe\ origin\ of\ the\ 
          coordinate\ system\ is\ the\ axle\ about\ which\ the\ beam\ rotates
          \ in\ \ a\ plane.\ \ The\ positive\ directions\ are\ up\ and\ to\ 
          the\ right, \ 
        as\ usual.\n\t\t\nThe\ mass\ of\ the\ beam\ is\ distributed\ along\ 
          the\ beam, \ 
        with\ a\ center\ of\ mass\ measured\ from\ the\ axle\ 
          \((along\ the\ beam)\), \ rc, \ 
        which\ is\ positive\ when\ it\ is\ on\ the\ short\ end\ of\ the\ 
          beam.\ \ The\ radius\ of\ gyrartion\ of\ the\ beam, \ rg, \ 
        is\ the\ position\ at\ which\ all\ of\ its\ mass\ would\ be\ in\ order
          \ for\ it\ to\ have\ an\ equivalent\ rotational\ 
          energy.\ \ \n\tFor\ a\ uniform\ beam, 
        we\ would\ have\ \ rg^2 = 
          \(\((1/3)\)\ \((l2^2 - l1\ l2\  + l1^2)\)\ and\ rc\  = \ 
            \((l1 - l2)\)/2\)\n*) \), "\n", "\n", \( (*constants*) \), "\n", 
      \(cn = {l1, l2, l3, l4, m1, m2, mb, rg, rc}; \nl1 =. ; l2 =. ; l3 =. ; 
      l4 =. ; m1 =. ; m2 =. ; m3 =. ; mb =. ; g =. ; rg =. ; \nth =. ; 
      phi =. ; psi =. ; \n\n
       (*coords\ of\ the\ short\ end\ of\ the\ \(beam : \)*) \n
      x1[th_] := \ l1\ Sin[th]; \ny1[th_] := \(-l1\)\ Cos[th]; \n
       (*coords\ of\ long\ end\ of\ \ \(beam : \)*) \n
      x2[th_] := \(-l2\)\ Sin[th]; \ny2[th_] := l2\ Cos[th]; \n
       (*coords\ of\ center\ of\ mass\ of\ the\ counterweight*) \n
      x4[th_, phi_] := l1\ Sin[th] - l4\ Sin[phi + th]; \n
      y4[th_, phi_] := \(-l1\)\ Cos[th] + l4\ Cos[phi + th]; \n
       (*\ \ \ \ coords\ of\ the\ projectile\ end\ of\ the\ sling\ *) \n
      x3[th_, psi_] := \(-\((l3\ Sin[psi - th]\  + l2\ Sin[th])\)\); \n
      y3[th_, psi_] := \(-\((l3\ Cos[psi - th]\  - l2\ Cos[th])\)\); \)}]], 
  "Input"],

Cell[BoxData[
    \(\( (*\ potential\ energy\ of\ the\ three\ \(parts : \)*) \n
    vt[th_, phi_, psi_]\  := 
      m1\ g\ y4[th, phi] + m2\ g\ y3[th, psi] - mb\ g\ rc\ Cos[th]; \n\n
     (*\ kinetic\ energy\ of\ the\ three\ \(parts : \)*) \n
    ket[th_, phi_, psi_] := 
      m1/2\ \((\((Dt[x4[th, phi], t, Constants -> cn])\)^2 + 
              \((Dt[y4[th, phi], t, Constants -> cn])\)^2\n)\) + 
        m2/2\ \((
            \((Dt[x3[th, psi], t, Constants -> cn])\)^2 + 
              \((Dt[y3[th, psi, Constants -> cn], t])\)^2)\)\  + \ \n
        \((mb/2)\)\ \((rg^2)\)\ Dt[th, t, Constants -> cn]^2; \n\n
     (*\ the\ lagrangian*) \n
    lagrt[th_, phi_, psi_] := ket[th, phi, psi] - vt[th, phi, psi]; \)\)], 
  "Input"],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*\ define\ what\ is\ constant, \ 
      and\ the\ derivatives\ of\ the\ angles, \ 
      denoted\ by\ the\ \*"\""d \*"\"\< suffix\>:"*) \n
    ltrr = lagrt[th, phi, psi] /. {
          Dt[th, t, Constants \[Rule] {l1, l2, l3, l4, m1, m2, mb, rg, rc}]
               -> thd, \n\t
          Dt[phi, t, Constants \[Rule] {l1, l2, l3, l4, m1, m2, mb, rg, rc}]
               -> phid, \n\t\t
          Dt[psi, t, Constants \[Rule] {l1, l2, l3, l4, m1, m2, mb, rg, rc}]
               -> psid, Dt[th, t] -> thd, Dt[phi, t] -> phid, 
          Dt[psi, t] -> psid}\)\)], "Input"],

Cell[BoxData[
    \(1\/2\ mb\ rg\^2\ thd\^2 + g\ mb\ rc\ Cos[th] - 
      g\ m2\ \((\(-l3\)\ Cos[psi - th] + l2\ Cos[th])\) - 
      g\ m1\ \((\(-l1\)\ Cos[th] + l4\ Cos[phi + th])\) + 
      1\/2\ m2\ \((
        \((\(-l3\)\ \((psid - thd)\)\ Cos[psi - th] - l2\ thd\ Cos[th])\)\^2
             + \((l3\ psid\ Sin[psi - th] + 
                thd\ \((\(-l3\)\ Sin[psi - th] - l2\ Sin[th])\))\)\^2)\) + 
      1\/2\ m1\ \((
        \((l1\ thd\ Cos[th] - l4\ \((phid + thd)\)\ Cos[phi + th])\)\^2 + 
          \((l1\ thd\ Sin[th] - l4\ \((phid + thd)\)\ Sin[phi + th])\)\^2)
        \)\)], "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*lagrange' s\ eqn\ for\ th\ *) \n
    eqthsim = 
      Simplify[\ 
        Dt[D[ltrr, thd], t] - D[ltrr, th] == 0 /. {Dt[l1, t] -> 0, 
            Dt[l2, t] -> 0, Dt[l3, t] -> 0, Dt[l4, t] -> 0, Dt[mb, t] -> 0, 
            Dt[m1, t] -> 0, Dt[m2, t] -> 0, Dt[g, t] -> 0, Dt[rg, t] -> 0, 
            Dt[rc, t] -> 0, \n\t\tDt[th, t] -> thd, Dt[phi, t] -> phid, 
            Dt[psi, t] -> psid, \t\tDt[thd, t] -> thdd, Dt[phid, t] -> phidd, 
            Dt[psid, t] -> psidd}]\)\)], "Input"],

Cell[BoxData[
    \(l4\^2\ m1\ phidd - l3\^2\ m2\ psidd + l1\^2\ m1\ thdd + 
        l4\^2\ m1\ thdd + l2\^2\ m2\ thdd + l3\^2\ m2\ thdd + 
        mb\ rg\^2\ thdd - l1\ l4\ m1\ \((phidd + 2\ thdd)\)\ Cos[phi] + 
        l2\ l3\ m2\ \((psidd - 2\ thdd)\)\ Cos[psi] + 
        l1\ l4\ m1\ phid\^2\ Sin[phi] + 2\ l1\ l4\ m1\ phid\ thd\ Sin[phi] - 
        l2\ l3\ m2\ psid\^2\ Sin[psi] + 2\ l2\ l3\ m2\ psid\ thd\ Sin[psi] - 
        g\ l3\ m2\ Sin[psi - th] + g\ l1\ m1\ Sin[th] - g\ l2\ m2\ Sin[th] + 
        g\ mb\ rc\ Sin[th] - g\ l4\ m1\ Sin[phi + th] == 0\)], "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[{
    \( (*similarly\ for\ phi\ and\ psi*) \n
    eqphisim = \ 
      Simplify[Dt[D[ltrr, phid], t] - D[ltrr, phi] == 0 /. {Dt[l1, t] -> 0, 
            Dt[l2, t] -> 0, Dt[l3, t] -> 0, Dt[l4, t] -> 0, Dt[mb, t] -> 0, 
            Dt[m1, t] -> 0, Dt[m2, t] -> 0, Dt[g, t] -> 0, Dt[rg, t] -> 0, 
            Dt[rc, t] -> 0, \n\t\tDt[th, t] -> thd, Dt[phi, t] -> phid, 
            Dt[psi, t] -> psid, \t\tDt[thd, t] -> thdd, Dt[phid, t] -> phidd, 
            Dt[psid, t] -> psidd}]\n\), 
    \(eqpsisim = \ 
      Simplify[Dt[D[ltrr, psid], t] - D[ltrr, psi] == 0 /. {Dt[l1, t] -> 0, 
            Dt[l2, t] -> 0, Dt[l3, t] -> 0, Dt[l4, t] -> 0, Dt[mb, t] -> 0, 
            Dt[m1, t] -> 0, Dt[m2, t] -> 0, Dt[g, t] -> 0, Dt[rg, t] -> 0, 
            Dt[rc, t] -> 0, \n\t\tDt[th, t] -> thd, Dt[phi, t] -> phid, 
            Dt[psi, t] -> psid, \t\tDt[thd, t] -> thdd, Dt[phid, t] -> phidd, 
            Dt[psid, t] -> psidd}]\)}], "Input"],

Cell[BoxData[
    \(l4\ m1\ \((
        l4\ phidd + l4\ thdd - l1\ thdd\ Cos[phi] - l1\ thd\^2\ Sin[phi] - 
          g\ Sin[phi + th])\) == 0\)], "Output"],

Cell[BoxData[
    \(l3\ m2\ \((
        l3\ psidd - l3\ thdd + l2\ thdd\ Cos[psi] - l2\ thd\^2\ Sin[psi] + 
          g\ Sin[psi - th])\) == 0\)], "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*\ 
      now\ solve\ the\ three\ equations\ for\ the\ three\ second\ 
        derivatives*) \n
    sol = Solve[{eqthsim, eqphisim, eqpsisim}, {thdd, phidd, psidd}]\)\)], 
  "Input"],

Cell[BoxData[
    \({{phidd \[Rule] 
          \(-\(\(\(-l1\)\ thd\^2\ Sin[phi] - g\ Sin[phi + th]\)\/l4\)\) + 
            \((\((l4 - l1\ Cos[phi])\)\ 
                  \((\(-l1\)\ l4\ m1\ phid\^2\ Sin[phi] - 
                      2\ l1\ l4\ m1\ phid\ thd\ Sin[phi] - 
                      l1\ l4\ m1\ thd\^2\ Sin[phi] + 
                      l1\^2\ m1\ thd\^2\ Cos[phi]\ Sin[phi] + 
                      l2\ l3\ m2\ psid\^2\ Sin[psi] - 
                      2\ l2\ l3\ m2\ psid\ thd\ Sin[psi] + 
                      l2\ l3\ m2\ thd\^2\ Sin[psi] - 
                      l2\^2\ m2\ thd\^2\ Cos[psi]\ Sin[psi] + 
                      g\ l2\ m2\ Cos[psi]\ Sin[psi - th] - 
                      g\ l1\ m1\ Sin[th] + g\ l2\ m2\ Sin[th] - 
                      g\ mb\ rc\ Sin[th] + g\ l1\ m1\ Cos[phi]\ Sin[phi + th])
                    \))\)/\((
                l4\ \((\(-l1\^2\)\ m1 - l2\^2\ m2 - mb\ rg\^2 + 
                    l1\^2\ m1\ Cos[phi]\^2 + l2\^2\ m2\ Cos[psi]\^2)\))\), 
        thdd \[Rule] 
          \(-\(\((\(-l1\)\ l4\ m1\ phid\^2\ Sin[phi] - 
                  2\ l1\ l4\ m1\ phid\ thd\ Sin[phi] - 
                  l1\ l4\ m1\ thd\^2\ Sin[phi] + 
                  l1\^2\ m1\ thd\^2\ Cos[phi]\ Sin[phi] + 
                  l2\ l3\ m2\ psid\^2\ Sin[psi] - 
                  2\ l2\ l3\ m2\ psid\ thd\ Sin[psi] + 
                  l2\ l3\ m2\ thd\^2\ Sin[psi] - 
                  l2\^2\ m2\ thd\^2\ Cos[psi]\ Sin[psi] + 
                  g\ l2\ m2\ Cos[psi]\ Sin[psi - th] - g\ l1\ m1\ Sin[th] + 
                  g\ l2\ m2\ Sin[th] - g\ mb\ rc\ Sin[th] + 
                  g\ l1\ m1\ Cos[phi]\ Sin[phi + th])\)/
              \((\(-l1\^2\)\ m1 - l2\^2\ m2 - mb\ rg\^2 + 
                  l1\^2\ m1\ Cos[phi]\^2 + l2\^2\ m2\ Cos[psi]\^2)\)\)\), 
        psidd \[Rule] 
          \(-\(\((\(-l3\)\ m2\ \((
                    \(-l4\)\ m1\ \((l4 - l1\ Cos[phi])\)\ 
                        \((l4\^2\ m1 - l1\ l4\ m1\ Cos[phi])\) + 
                      l4\^2\ m1\ \((
                        l1\^2\ m1 + l4\^2\ m1 + l2\^2\ m2 + l3\^2\ m2 + 
                          mb\ rg\^2 - 2\ l1\ l4\ m1\ Cos[phi] - 
                          2\ l2\ l3\ m2\ Cos[psi])\))\)\ 
                    \((\(-l2\)\ thd\^2\ Sin[psi] + g\ Sin[psi - th])\) + 
                  l3\ m2\ \((\(-l3\) + l2\ Cos[psi])\)\ 
                    \((\(-l4\)\ m1\ \((l4\^2\ m1 - l1\ l4\ m1\ Cos[phi])\)\ 
                          \((\(-l1\)\ thd\^2\ Sin[phi] - g\ Sin[phi + th])
                            \) + l4\^2\ m1\ \((
                          l1\ l4\ m1\ phid\^2\ Sin[phi] + 
                            2\ l1\ l4\ m1\ phid\ thd\ Sin[phi] - 
                            l2\ l3\ m2\ psid\^2\ Sin[psi] + 
                            2\ l2\ l3\ m2\ psid\ thd\ Sin[psi] - 
                            g\ l3\ m2\ Sin[psi - th] + g\ l1\ m1\ Sin[th] - 
                            g\ l2\ m2\ Sin[th] + g\ mb\ rc\ Sin[th] - 
                            g\ l4\ m1\ Sin[phi + th])\))\))\)/
              \((\(-l1\^2\)\ l3\^2\ l4\^2\ m1\^2\ m2 - 
                  l2\^2\ l3\^2\ l4\^2\ m1\ m2\^2 - 
                  l3\^2\ l4\^2\ m1\ m2\ mb\ rg\^2 + 
                  l1\^2\ l3\^2\ l4\^2\ m1\^2\ m2\ Cos[phi]\^2 + 
                  l2\^2\ l3\^2\ l4\^2\ m1\ m2\^2\ Cos[psi]\^2)\)\)\)}}\)], 
  "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*\ get\ the\ fortran\ form\ for\ each\ of\ these, \ 
      solving\ for\ the\ numerator\ and\ denominator\ separately\ \ *) \n
    s2 = \(Flatten[sol]\)[\([2]\)]; \nthddsol = thdd /. s2; \n
    numthdd = FortranForm[Numerator[thddsol]]\n\)\)], "Input"],

Cell["\<\
        l1*l4*m1*phid**2*Sin(phi) + 2*l1*l4*m1*phid*thd*Sin(phi) + 
     -  l1*l4*m1*thd**2*Sin(phi) - l1**2*m1*thd**2*Cos(phi)*Sin(phi) - 
     -  l2*l3*m2*psid**2*Sin(psi) + 2*l2*l3*m2*psid*thd*Sin(psi) - 
     -  l2*l3*m2*thd**2*Sin(psi) + l2**2*m2*thd**2*Cos(psi)*Sin(psi) - 
     -  g*l2*m2*Cos(psi)*Sin(psi - th) + g*l1*m1*Sin(th) - g*l2*m2*Sin(th) + 
     -  g*mb*rc*Sin(th) - g*l1*m1*Cos(phi)*Sin(phi + th)\
\>", "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\(\ndenomthdd = FortranForm[Denominator[thddsol]]\)\)], "Input"],

Cell["\<\
-(l1**2*m1) - l2**2*m2 - mb*rg**2 + l1**2*m1*Cos(phi)**2 + \
l2**2*m2*Cos(psi)**2\
\>", "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*\ and\ similarly\ for\ the\ other\ two\ angles*) \n
    s1 = \(Flatten[sol]\)[\([1]\)]; \nphiddsol = Together[phidd /. s1]; \n\n
    numphidd = FortranForm[Numerator[phiddsol]]\n\)\)], "Input"],

Cell["\<\
        -(l1*l4**2*m1*phid**2*Sin(phi)) - \
2*l1*l4**2*m1*phid*thd*Sin(phi) - 
     -  l1**3*m1*thd**2*Sin(phi) - l1*l4**2*m1*thd**2*Sin(phi) - 
     -  l1*l2**2*m2*thd**2*Sin(phi) - l1*mb*rg**2*thd**2*Sin(phi) + 
     -  l1**2*l4*m1*phid**2*Cos(phi)*Sin(phi) + 
     -  2*l1**2*l4*m1*phid*thd*Cos(phi)*Sin(phi) + 
     -  2*l1**2*l4*m1*thd**2*Cos(phi)*Sin(phi) + 
     -  l1*l2**2*m2*thd**2*Cos(psi)**2*Sin(phi) + \
l2*l3*l4*m2*psid**2*Sin(psi) - 
     -  2*l2*l3*l4*m2*psid*thd*Sin(psi) + l2*l3*l4*m2*thd**2*Sin(psi) - 
     -  l1*l2*l3*m2*psid**2*Cos(phi)*Sin(psi) + 
     -  2*l1*l2*l3*m2*psid*thd*Cos(phi)*Sin(psi) - 
     -  l1*l2*l3*m2*thd**2*Cos(phi)*Sin(psi) - 
     -  l2**2*l4*m2*thd**2*Cos(psi)*Sin(psi) + 
     -  l1*l2**2*m2*thd**2*Cos(phi)*Cos(psi)*Sin(psi) + 
     -  g*l2*l4*m2*Cos(psi)*Sin(psi - th) - 
     -  g*l1*l2*m2*Cos(phi)*Cos(psi)*Sin(psi - th) - g*l1*l4*m1*Sin(th) + 
     -  g*l2*l4*m2*Sin(th) - g*l4*mb*rc*Sin(th) + g*l1**2*m1*Cos(phi)*Sin(th) \
- 
     -  g*l1*l2*m2*Cos(phi)*Sin(th) + g*l1*mb*rc*Cos(phi)*Sin(th) - 
     -  g*l1**2*m1*Sin(phi + th) - g*l2**2*m2*Sin(phi + th) - 
     -  g*mb*rg**2*Sin(phi + th) + g*l1*l4*m1*Cos(phi)*Sin(phi + th) + 
     -  g*l2**2*m2*Cos(psi)**2*Sin(phi + th)\
\>", "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(denomphidd = FortranForm[Denominator[phiddsol]]\)], "Input"],

Cell["\<\
        l4*(-(l1**2*m1) - l2**2*m2 - mb*rg**2 + \
l1**2*m1*Cos(phi)**2 + 
     -    l2**2*m2*Cos(psi)**2)\
\>", "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\(s3 = \(Flatten[sol]\)[\([3]\)]; \npsiddsol = Together[psidd /. s3]; \n
    numpsidd = FortranForm[Numerator[psiddsol]]\n\)\)], "Input"],

Cell["\<\
        l1*l3*l4*m1*phid**2*Sin(phi) + \
2*l1*l3*l4*m1*phid*thd*Sin(phi) + 
     -  l1*l3*l4*m1*thd**2*Sin(phi) - l1**2*l3*m1*thd**2*Cos(phi)*Sin(phi) - 
     -  l1*l2*l4*m1*phid**2*Cos(psi)*Sin(phi) - 
     -  2*l1*l2*l4*m1*phid*thd*Cos(psi)*Sin(phi) - 
     -  l1*l2*l4*m1*thd**2*Cos(psi)*Sin(phi) + 
     -  l1**2*l2*m1*thd**2*Cos(phi)*Cos(psi)*Sin(phi) - 
     -  l2*l3**2*m2*psid**2*Sin(psi) + 2*l2*l3**2*m2*psid*thd*Sin(psi) - 
     -  l1**2*l2*m1*thd**2*Sin(psi) - l2**3*m2*thd**2*Sin(psi) - 
     -  l2*l3**2*m2*thd**2*Sin(psi) - l2*mb*rg**2*thd**2*Sin(psi) + 
     -  l1**2*l2*m1*thd**2*Cos(phi)**2*Sin(psi) + 
     -  l2**2*l3*m2*psid**2*Cos(psi)*Sin(psi) - 
     -  2*l2**2*l3*m2*psid*thd*Cos(psi)*Sin(psi) + 
     -  2*l2**2*l3*m2*thd**2*Cos(psi)*Sin(psi) + g*l1**2*m1*Sin(psi - th) + 
     -  g*l2**2*m2*Sin(psi - th) + g*mb*rg**2*Sin(psi - th) - 
     -  g*l1**2*m1*Cos(phi)**2*Sin(psi - th) - 
     -  g*l2*l3*m2*Cos(psi)*Sin(psi - th) + g*l1*l3*m1*Sin(th) - 
     -  g*l2*l3*m2*Sin(th) + g*l3*mb*rc*Sin(th) - g*l1*l2*m1*Cos(psi)*Sin(th) \
+ 
     -  g*l2**2*m2*Cos(psi)*Sin(th) - g*l2*mb*rc*Cos(psi)*Sin(th) - 
     -  g*l1*l3*m1*Cos(phi)*Sin(phi + th) + 
     -  g*l1*l2*m1*Cos(phi)*Cos(psi)*Sin(phi + th)\
\>", "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(denompsidd = FortranForm[Denominator[psiddsol]]\)], "Input"],

Cell["\<\
        -(l1**2*l3*m1) - l2**2*l3*m2 - l3*mb*rg**2 + \
l1**2*l3*m1*Cos(phi)**2 + 
     -  l2**2*l3*m2*Cos(psi)**2\
\>", "Output"]
}, Open  ]],

Cell[BoxData[
    \( (*\ 
      these\ six\ results\ can\ be\ copied\ and\ pasted\ into\ a\ word\ 
        processor, \ 
      where\ substitutions\ and\ adjustments\ made\ for\ the\ translation\ 
        into\ other\ languages\ are\ made.\ \ To\ speed\ up\ the\ 
        calculations, \ 
      repeated\ terms\ are\ substituted\ for\ so\ that\ they\ are\ only\ 
        calculated\ \(once : \ \ e.g.\ Sin \((th)\)\ is\ replaced\ by\ sth\), 
      \ \(etc.\)\ \ *) \)], "Input"],

Cell[BoxData[
    \( (*\ 
      a\ similar\ set\ of\ equations\ are\ derived\ for\ the\ constrained\ 
        portion\ of\ the\ throw, \ 
      where\ the\ projectile\ slides\ along\ the\ trough\ at\ the\ 
        start.\ \ Now\ only\ two\ equations\ are\ required, \ 
      for\ the\ two\ variables\ th\ an\ phi.\ \ The\ solution\ for\ psi\ is\ 
        given\ in\ terms\ of\ the\ other\ \(two.\)*) \)], "Input"],

Cell[CellGroupData[{

Cell[BoxData[{
    \(\n\t (*\ the\ constraint\ function\ is*) \n
    fconstr[th_, phi_, psi_] := 
      psi - th + Pi/2 - ArcSin[\((l2/l3)\)\ *\ \((Sin[psis] + Cos[th])\)]; \n
     (*\ so\ the\ derivatives\ are*) \nath = D[fconstr[th, phi, psi], th]\), 
    \(aphi = D[fconstr[th, phi, psi], phi]\), 
    \(apsi = D[fconstr[th, phi, psi], psi]\n\)}], "Input"],

Cell[BoxData[
    \(\(-1\) + 
      \(l2\ Sin[th]\)\/\(l3\ 
          \@\(1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2\)\)\)], 
  "Output"],

Cell[BoxData[
    \(0\)], "Output"],

Cell[BoxData[
    \(1\)], "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*the\ lagrange\ multiplier\ is\ given\ by\ *) \n
    lam = \n\ \ \ \ \ \ 
      Dt[D[ltrr, psid], t] - D[ltrr, psi] /. {Dt[l1, t] -> 0, Dt[l2, t] -> 0, 
          Dt[l3, t] -> 0, Dt[l4, t] -> 0, Dt[mb, t] -> 0, Dt[m1, t] -> 0, 
          Dt[m2, t] -> 0, Dt[g, t] -> 0, \n\t\t\tDt[l1, {t, 2}] -> 0, \t
          Dt[l2, {t, 2}] -> 0, \tDt[l3, {t, 2}] -> 0, Dt[l4, {t, 2}] -> 0, \n
          \t\t\tDt[rc, t] -> 0, Dt[rg, t] -> 0, \n\t\tDt[th, t] -> thd, 
          Dt[phi, t] -> phid, Dt[psi, t] -> psid, \t\tDt[thd, t] -> thdd, 
          Dt[phid, t] -> phidd, Dt[psid, t] -> psidd}\)\)], "Input"],

Cell[BoxData[
    \(g\ l3\ m2\ Sin[psi - th] - 
      1\/2\ m2\ \((
        2\ l3\ \((psid - thd)\)\ 
            \((\(-l3\)\ \((psid - thd)\)\ Cos[psi - th] - l2\ thd\ Cos[th])
              \)\ Sin[psi - th] + 
          2\ \((l3\ psid\ Cos[psi - th] - l3\ thd\ Cos[psi - th])\)\ 
            \((l3\ psid\ Sin[psi - th] + 
                thd\ \((\(-l3\)\ Sin[psi - th] - l2\ Sin[th])\))\))\) + 
      1\/2\ m2\ \((
        2\ l3\ \((psid - thd)\)\ 
            \((\(-l3\)\ \((psid - thd)\)\ Cos[psi - th] - l2\ thd\ Cos[th])
              \)\ Sin[psi - th] - 
          2\ l3\ Cos[psi - th]\ 
            \((\(-l3\)\ \((psidd - thdd)\)\ Cos[psi - th] - 
                l2\ thdd\ Cos[th] + l3\ \((psid - thd)\)\^2\ Sin[psi - th] + 
                l2\ thd\^2\ Sin[th])\) + 
          2\ l3\ \((psid - thd)\)\ Cos[psi - th]\ 
            \((l3\ psid\ Sin[psi - th] + 
                thd\ \((\(-l3\)\ Sin[psi - th] - l2\ Sin[th])\))\) + 
          2\ l3\ Sin[psi - th]\ 
            \((l3\ psid\ \((psid - thd)\)\ Cos[psi - th] + 
                thd\ \((
                  \(-l3\)\ \((psid - thd)\)\ Cos[psi - th] - l2\ thd\ Cos[th])
                  \) + l3\ psidd\ Sin[psi - th] + 
                thdd\ \((\(-l3\)\ Sin[psi - th] - l2\ Sin[th])\))\))\)\)], 
  "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(lamath2 = 
      lam*ath /. 
        \((\(l2\ Sin[th]\)\/\(l3\ 
                \@\(1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2\)\))
            \) -> radc; \n\n
    eqthl = Simplify[
        \((Dt[D[ltrr, thd], t] - D[ltrr, th])\) /. {Dt[l1, t] -> 0, 
            Dt[l2, t] -> 0, Dt[l3, t] -> 0, Dt[l4, t] -> 0, Dt[mb, t] -> 0, 
            Dt[m1, t] -> 0, Dt[m2, t] -> 0, Dt[g, t] -> 0, \n\t\t\t
            Dt[l1, {t, 2}] -> 0, \tDt[l2, {t, 2}] -> 0, \t
            Dt[l3, {t, 2}] -> 0, \tDt[l4, {t, 2}] -> 0, \n\t\t\t\t\t\t\t
            Dt[rc, t] -> 0, Dt[rg, t] -> 0, \n\t\tDt[th, t] -> thd, 
            Dt[phi, t] -> phid, Dt[psi, t] -> psid, \t\tDt[thd, t] -> thdd, 
            Dt[phid, t] -> phidd, Dt[psid, t] -> psidd}]; \n\n
    eqth2 = ExpandAll[eqthl - lamath2 == 0]\)], "Input"],

Cell[BoxData[
    \(General::"spell1" \( : \ \) 
      "Possible spelling error: new symbol name \"\!\(eqthl\)\" is similar to \
existing symbol \"\!\(eqth\)\"."\)], "Message"],

Cell[BoxData[
    \(l4\^2\ m1\ phidd - l3\^2\ m2\ psidd + l1\^2\ m1\ thdd + 
        l4\^2\ m1\ thdd + l2\^2\ m2\ thdd + l3\^2\ m2\ thdd + 
        mb\ rg\^2\ thdd - l1\ l4\ m1\ phidd\ Cos[phi] - 
        2\ l1\ l4\ m1\ thdd\ Cos[phi] + l2\ l3\ m2\ psidd\ Cos[psi] - 
        2\ l2\ l3\ m2\ thdd\ Cos[psi] + l3\^2\ m2\ psidd\ Cos[psi - th]\^2 - 
        l3\^2\ m2\ psidd\ radc\ Cos[psi - th]\^2 - 
        l3\^2\ m2\ thdd\ Cos[psi - th]\^2 + 
        l3\^2\ m2\ radc\ thdd\ Cos[psi - th]\^2 + 
        l2\ l3\ m2\ thdd\ Cos[psi - th]\ Cos[th] - 
        l2\ l3\ m2\ radc\ thdd\ Cos[psi - th]\ Cos[th] + 
        l1\ l4\ m1\ phid\^2\ Sin[phi] + 2\ l1\ l4\ m1\ phid\ thd\ Sin[phi] - 
        l2\ l3\ m2\ psid\^2\ Sin[psi] + 2\ l2\ l3\ m2\ psid\ thd\ Sin[psi] - 
        g\ l3\ m2\ radc\ Sin[psi - th] - 
        l2\ l3\ m2\ thd\^2\ Cos[th]\ Sin[psi - th] + 
        l2\ l3\ m2\ radc\ thd\^2\ Cos[th]\ Sin[psi - th] + 
        l3\^2\ m2\ psidd\ Sin[psi - th]\^2 - 
        l3\^2\ m2\ psidd\ radc\ Sin[psi - th]\^2 - 
        l3\^2\ m2\ thdd\ Sin[psi - th]\^2 + 
        l3\^2\ m2\ radc\ thdd\ Sin[psi - th]\^2 + g\ l1\ m1\ Sin[th] - 
        g\ l2\ m2\ Sin[th] + g\ mb\ rc\ Sin[th] - 
        l2\ l3\ m2\ thd\^2\ Cos[psi - th]\ Sin[th] + 
        l2\ l3\ m2\ radc\ thd\^2\ Cos[psi - th]\ Sin[th] - 
        l2\ l3\ m2\ thdd\ Sin[psi - th]\ Sin[th] + 
        l2\ l3\ m2\ radc\ thdd\ Sin[psi - th]\ Sin[th] - 
        g\ l4\ m1\ Sin[phi + th] == 0\)], "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*\ and\ similarly\ for\ \(phi : \)\ \ *) \n
    eqphil = Simplify[
        \((Dt[D[ltrr, phid], t] - D[ltrr, phi])\) /. {Dt[l1, t] -> 0, 
            Dt[l2, t] -> 0, Dt[l3, t] -> 0, Dt[l4, t] -> 0, Dt[mb, t] -> 0, 
            Dt[m1, t] -> 0, Dt[m2, t] -> 0, Dt[g, t] -> 0, \n\t\t\t
            Dt[l1, {t, 2}] -> 0, \tDt[l2, {t, 2}] -> 0, \t
            Dt[l3, {t, 2}] -> 0, \tDt[l4, {t, 2}] -> 0, \n\t\t\t\t\t\t\t
            Dt[rc, t] -> 0, Dt[rg, t] -> 0, \n\t\tDt[th, t] -> thd, 
            Dt[phi, t] -> phid, Dt[psi, t] -> psid, \t\tDt[thd, t] -> thdd, 
            Dt[phid, t] -> phidd, Dt[psid, t] -> psidd}]; \n\n
    eqphi2 = eqphil - lam*aphi == 0\)\)], "Input"],

Cell[BoxData[
    \(l4\ m1\ \((
        l4\ phidd + l4\ thdd - l1\ thdd\ Cos[phi] - l1\ thd\^2\ Sin[phi] - 
          g\ Sin[phi + th])\) == 0\)], "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[{
    \( (*\ the\ derivatgive\ for\ psi\ \(are : \)*) \n
    psi = th - Pi/2 + ArcSin[\((l2/l3)\)*\((Sin[psis] + Cos[th])\)]\), 
    \(\tpsid2 = \t\t
      Dt[psi, t] /. {Dt[l2, t] -> 0, Dt[l3, t] -> 0, \ Dt[psis, t] -> 0, 
          Dt[th, t] -> thd}\), 
    \(\tpsidd2 = 
      Dt[psid2, t] /. {Dt[l2, t] -> 0, Dt[l3, t] -> 0, \ Dt[psis, t] -> 0, 
          Dt[th, t] -> thd, Dt[thd, t] -> thdd}\)}], "Input"],

Cell[BoxData[
    \(\(-\(\[Pi]\/2\)\) + th + 
      ArcSin[\(l2\ \((Cos[th] + Sin[psis])\)\)\/l3]\)], "Output"],

Cell[BoxData[
    \(thd - 
      \(l2\ thd\ Sin[th]\)\/\(l3\ 
          \@\(1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2\)\)\)], 
  "Output"],

Cell[BoxData[
    \(thdd - 
      \(l2\ thd\^2\ Cos[th]\)\/\(l3\ 
          \@\(1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2\)\) - 
      \(l2\ thdd\ Sin[th]\)\/\(l3\ 
          \@\(1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2\)\) + 
      \(l2\^3\ thd\^2\ \((Cos[th] + Sin[psis])\)\ Sin[th]\^2\)\/\(l3\^3\ 
          \((1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2)
            \)\^\(3/2\)\)\)], "Output"]
}, Open  ]],

Cell[BoxData[
    \(\( (*
      define\ these\ intermediate\ variables\ to\ simplify\ the\ 
        expressions*) \n
    psid3 = psid2
         /. {1/\@\(1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2\) -> 
            1/rad}; \n
    psidd3 = psidd2
         /. {1/\@\(1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2\) -> 
            1/rad, 1/
              \((1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2)\)^
                \((3/2)\) -> 1/rad^3}; \)\)], "Input"],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*\ 
      the\ two\ equations\ to\ be\ solved\ for\ thdd\ and\ phidd\ are\ 
        \(then : \)*) \n
    eq4 = {\({eqth2, eqphi2} /. {psidd -> psidd3, psid -> psid3, 
              Cos[th + ArcSin[\(l2\ \((Cos[th] + Sin[psis])\)\)\/l3]] -> 
                Expand[Cos[
                    th + ArcSin[\(l2\ \((Cos[th] + Sin[psis])\)\)\/l3]], 
                  Trig -> True], \n\t
              Sin[th + ArcSin[\(l2\ \((Cos[th] + Sin[psis])\)\)\/l3]] -> 
                Expand[Sin[
                    th + ArcSin[\(l2\ \((Cos[th] + Sin[psis])\)\)\/l3]], 
                  Trig -> True], 
              \@\(1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2\) -> rad}
              \) /. \@\(1 - 
                \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2\) -> rad}\)\)],
   "Input"],

Cell[BoxData[
    \({{l4\^2\ m1\ phidd + g\ l3\ m2\ rad\ radc + l1\^2\ m1\ thdd + 
            l4\^2\ m1\ thdd + l2\^2\ m2\ thdd + l3\^2\ m2\ thdd + 
            mb\ rg\^2\ thdd - l1\ l4\ m1\ phidd\ Cos[phi] - 
            2\ l1\ l4\ m1\ thdd\ Cos[phi] + 
            l2\ l3\ m2\ rad\ thd\^2\ Cos[th] - 
            l2\ l3\ m2\ rad\ radc\ thd\^2\ Cos[th] + 
            l1\ l4\ m1\ phid\^2\ Sin[phi] + 
            2\ l1\ l4\ m1\ phid\ thd\ Sin[phi] + 
            l2\^2\ m2\ thdd\ Cos[th]\ \((Cos[th] + Sin[psis])\) - 
            l2\^2\ m2\ radc\ thdd\ Cos[th]\ \((Cos[th] + Sin[psis])\) - 
            l2\^2\ m2\ thdd\ \((Cos[th] + Sin[psis])\)\^2 + 
            l2\^2\ m2\ radc\ thdd\ \((Cos[th] + Sin[psis])\)\^2 - 
            l3\^2\ m2\ thdd\ \((
              1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2)\) + 
            l3\^2\ m2\ radc\ thdd\ \((
              1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2)\) + 
            g\ l1\ m1\ Sin[th] - g\ l2\ m2\ Sin[th] + g\ mb\ rc\ Sin[th] + 
            l2\ l3\ m2\ rad\ thdd\ Sin[th] - 
            l2\ l3\ m2\ rad\ radc\ thdd\ Sin[th] - 
            l2\^2\ m2\ thd\^2\ \((Cos[th] + Sin[psis])\)\ Sin[th] + 
            l2\^2\ m2\ radc\ thd\^2\ \((Cos[th] + Sin[psis])\)\ Sin[th] - 
            2\ l2\ l3\ m2\ thdd\ \((
              \(l2\ Cos[th]\^2\)\/l3 + \(l2\ Cos[th]\ Sin[psis]\)\/l3 + 
                rad\ Sin[th])\) - 
            2\ l2\ l3\ m2\ thd\ \((thd - \(l2\ thd\ Sin[th]\)\/\(l3\ rad\))
              \)\ \((rad\ Cos[th] - \(l2\ Cos[th]\ Sin[th]\)\/l3 - 
                  \(l2\ Sin[psis]\ Sin[th]\)\/l3)\) + 
            l2\ l3\ m2\ \((thd - \(l2\ thd\ Sin[th]\)\/\(l3\ rad\))\)\^2\ 
              \((rad\ Cos[th] - \(l2\ Cos[th]\ Sin[th]\)\/l3 - 
                  \(l2\ Sin[psis]\ Sin[th]\)\/l3)\) - 
            l3\^2\ m2\ \((
              thdd - \(l2\ thd\^2\ Cos[th]\)\/\(l3\ rad\) - 
                \(l2\ thdd\ Sin[th]\)\/\(l3\ rad\) + 
                \(l2\^3\ thd\^2\ \((Cos[th] + Sin[psis])\)\ 
                    Sin[th]\^2\)\/\(l3\^3\ rad\^3\))\) + 
            l2\^2\ m2\ \((Cos[th] + Sin[psis])\)\^2\ 
              \((thdd - \(l2\ thd\^2\ Cos[th]\)\/\(l3\ rad\) - 
                  \(l2\ thdd\ Sin[th]\)\/\(l3\ rad\) + 
                  \(l2\^3\ thd\^2\ \((Cos[th] + Sin[psis])\)\ 
                      Sin[th]\^2\)\/\(l3\^3\ rad\^3\))\) - 
            l2\^2\ m2\ radc\ \((Cos[th] + Sin[psis])\)\^2\ 
              \((thdd - \(l2\ thd\^2\ Cos[th]\)\/\(l3\ rad\) - 
                  \(l2\ thdd\ Sin[th]\)\/\(l3\ rad\) + 
                  \(l2\^3\ thd\^2\ \((Cos[th] + Sin[psis])\)\ 
                      Sin[th]\^2\)\/\(l3\^3\ rad\^3\))\) + 
            l3\^2\ m2\ \((1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2)
              \)\ \((thdd - \(l2\ thd\^2\ Cos[th]\)\/\(l3\ rad\) - 
                  \(l2\ thdd\ Sin[th]\)\/\(l3\ rad\) + 
                  \(l2\^3\ thd\^2\ \((Cos[th] + Sin[psis])\)\ 
                      Sin[th]\^2\)\/\(l3\^3\ rad\^3\))\) - 
            l3\^2\ m2\ radc\ \((
              1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2)\)\ 
              \((thdd - \(l2\ thd\^2\ Cos[th]\)\/\(l3\ rad\) - 
                  \(l2\ thdd\ Sin[th]\)\/\(l3\ rad\) + 
                  \(l2\^3\ thd\^2\ \((Cos[th] + Sin[psis])\)\ 
                      Sin[th]\^2\)\/\(l3\^3\ rad\^3\))\) + 
            l2\ l3\ m2\ \((
              \(l2\ Cos[th]\^2\)\/l3 + \(l2\ Cos[th]\ Sin[psis]\)\/l3 + 
                rad\ Sin[th])\)\ 
              \((thdd - \(l2\ thd\^2\ Cos[th]\)\/\(l3\ rad\) - 
                  \(l2\ thdd\ Sin[th]\)\/\(l3\ rad\) + 
                  \(l2\^3\ thd\^2\ \((Cos[th] + Sin[psis])\)\ 
                      Sin[th]\^2\)\/\(l3\^3\ rad\^3\))\) - 
            g\ l4\ m1\ Sin[phi + th] == 0, 
        l4\ m1\ \((
            l4\ phidd + l4\ thdd - l1\ thdd\ Cos[phi] - 
              l1\ thd\^2\ Sin[phi] - g\ Sin[phi + th])\) == 0}}\)], "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*\ this\ is\ a\ mess, \ but\ can\ be\ simplified\ to\ *) \n
    eq5 = ExpandAll[eq4] /. radc -> l2*Sin[th]/\((l3*rad)\)\)\)], "Input"],

Cell[BoxData[
    \({{l4\^2\ m1\ phidd + l1\^2\ m1\ thdd + l4\^2\ m1\ thdd + 
            l2\^2\ m2\ thdd + mb\ rg\^2\ thdd - l1\ l4\ m1\ phidd\ Cos[phi] - 
            2\ l1\ l4\ m1\ thdd\ Cos[phi] - 
            \(l2\^3\ m2\ thd\^2\ Cos[th]\^3\)\/\(l3\ rad\) + 
            l1\ l4\ m1\ phid\^2\ Sin[phi] + 
            2\ l1\ l4\ m1\ phid\ thd\ Sin[phi] - 
            \(l2\^3\ m2\ thd\^2\ Cos[th]\^2\ Sin[psis]\)\/\(l3\ rad\) + 
            g\ l1\ m1\ Sin[th] + g\ mb\ rc\ Sin[th] - 
            2\ l2\^2\ m2\ thd\^2\ Cos[th]\ Sin[th] + 
            \(l2\^2\ m2\ thd\^2\ Cos[th]\ Sin[th]\)\/rad\^2 - 
            \(2\ l2\^3\ m2\ thdd\ Cos[th]\^2\ Sin[th]\)\/\(l3\ rad\) - 
            \(2\ l2\^3\ m2\ thdd\ Cos[th]\ Sin[psis]\ Sin[th]\)\/\(l3\ rad
                \) - 2\ l2\^2\ m2\ thdd\ Sin[th]\^2 + 
            \(l2\^2\ m2\ thdd\ Sin[th]\^2\)\/rad\^2 + 
            \(2\ l2\^3\ m2\ thd\^2\ Cos[th]\ Sin[th]\^2\)\/\(l3\ rad\) + 
            \(l2\^5\ m2\ thd\^2\ Cos[th]\^3\ Sin[th]\^2\)\/\(l3\^3\ 
                rad\^3\) + 
            \(l2\^3\ m2\ thd\^2\ Sin[psis]\ Sin[th]\^2\)\/\(l3\ rad\) + 
            \(2\ l2\^5\ m2\ thd\^2\ Cos[th]\^2\ Sin[psis]\ 
                Sin[th]\^2\)\/\(l3\^3\ rad\^3\) + 
            \(l2\^5\ m2\ thd\^2\ Cos[th]\ Sin[psis]\^2\ 
                Sin[th]\^2\)\/\(l3\^3\ rad\^3\) - 
            \(l2\^4\ m2\ thd\^2\ Cos[th]\ Sin[th]\^3\)\/\(l3\^2\ rad\^4\) - 
            \(l2\^4\ m2\ thd\^2\ Sin[psis]\ Sin[th]\^3\)\/\(l3\^2\ rad\^4\) - 
            g\ l4\ m1\ Sin[phi + th] == 0, 
        l4\^2\ m1\ phidd + l4\^2\ m1\ thdd - l1\ l4\ m1\ thdd\ Cos[phi] - 
            l1\ l4\ m1\ thd\^2\ Sin[phi] - g\ l4\ m1\ Sin[phi + th] == 0}}\)],
   "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*\ and\ finally\ solved\ for\ thdd\ and\ \(phidd : \)*) \n
    sol5 = Solve[eq5, {thdd, phidd}]\)\)], "Input"],

Cell[BoxData[
    \({{thdd \[Rule] 
          \(-\(\((l2\^3\ l3\^2\ m2\ rad\^3\ thd\^2\ Cos[th]\^3 - 
                  l1\ l3\^3\ l4\ m1\ phid\^2\ rad\^4\ Sin[phi] - 
                  2\ l1\ l3\^3\ l4\ m1\ phid\ rad\^4\ thd\ Sin[phi] - 
                  l1\ l3\^3\ l4\ m1\ rad\^4\ thd\^2\ Sin[phi] + 
                  l1\^2\ l3\^3\ m1\ rad\^4\ thd\^2\ Cos[phi]\ Sin[phi] + 
                  l2\^3\ l3\^2\ m2\ rad\^3\ thd\^2\ Cos[th]\^2\ Sin[psis] - 
                  g\ l1\ l3\^3\ m1\ rad\^4\ Sin[th] - 
                  g\ l3\^3\ mb\ rad\^4\ rc\ Sin[th] - 
                  l2\^2\ l3\^3\ m2\ rad\^2\ thd\^2\ Cos[th]\ Sin[th] + 
                  2\ l2\^2\ l3\^3\ m2\ rad\^4\ thd\^2\ Cos[th]\ Sin[th] - 
                  2\ l2\^3\ l3\^2\ m2\ rad\^3\ thd\^2\ Cos[th]\ Sin[th]\^2 - 
                  l2\^5\ m2\ rad\ thd\^2\ Cos[th]\^3\ Sin[th]\^2 - 
                  l2\^3\ l3\^2\ m2\ rad\^3\ thd\^2\ Sin[psis]\ Sin[th]\^2 - 
                  2\ l2\^5\ m2\ rad\ thd\^2\ Cos[th]\^2\ Sin[psis]\ 
                    Sin[th]\^2 - 
                  l2\^5\ m2\ rad\ thd\^2\ Cos[th]\ Sin[psis]\^2\ Sin[th]\^2 + 
                  l2\^4\ l3\ m2\ thd\^2\ Cos[th]\ Sin[th]\^3 + 
                  l2\^4\ l3\ m2\ thd\^2\ Sin[psis]\ Sin[th]\^3 + 
                  g\ l1\ l3\^3\ m1\ rad\^4\ Cos[phi]\ Sin[phi + th])\)/
              \((l3\^2\ rad\^2\ \((
                  \(-l1\^2\)\ l3\ m1\ rad\^2 - l2\^2\ l3\ m2\ rad\^2 - 
                    l3\ mb\ rad\^2\ rg\^2 + 
                    l1\^2\ l3\ m1\ rad\^2\ Cos[phi]\^2 + 
                    2\ l2\^3\ m2\ rad\ Cos[th]\^2\ Sin[th] + 
                    2\ l2\^3\ m2\ rad\ Cos[th]\ Sin[psis]\ Sin[th] - 
                    l2\^2\ l3\ m2\ Sin[th]\^2 + 
                    2\ l2\^2\ l3\ m2\ rad\^2\ Sin[th]\^2)\))\)\)\), 
        phidd \[Rule] 
          \(-\(\((\(-\((l1\^2\ m1 + l4\^2\ m1 + l2\^2\ m2 + mb\ rg\^2 - 
                          2\ l1\ l4\ m1\ Cos[phi] - 
                          \(2\ l2\^3\ m2\ Cos[th]\^2\ Sin[th]\)\/\(l3\ rad
                              \) - \(2\ l2\^3\ m2\ Cos[th]\ Sin[psis]\ 
                              Sin[th]\)\/\(l3\ rad\) - 
                          2\ l2\^2\ m2\ Sin[th]\^2 + 
                          \(l2\^2\ m2\ Sin[th]\^2\)\/rad\^2)\)\)\ 
                    \((\(-l1\)\ l4\ m1\ thd\^2\ Sin[phi] - 
                        g\ l4\ m1\ Sin[phi + th])\) + 
                  \((l4\^2\ m1 - l1\ l4\ m1\ Cos[phi])\)\ 
                    \((\(-\(\(l2\^3\ m2\ thd\^2\ Cos[th]\^3\)\/\(l3\ rad
                                \)\)\) + l1\ l4\ m1\ phid\^2\ Sin[phi] + 
                        2\ l1\ l4\ m1\ phid\ thd\ Sin[phi] - 
                        \(l2\^3\ m2\ thd\^2\ Cos[th]\^2\ Sin[psis]\)\/\(l3\ 
                            rad\) + g\ l1\ m1\ Sin[th] + g\ mb\ rc\ Sin[th] - 
                        2\ l2\^2\ m2\ thd\^2\ Cos[th]\ Sin[th] + 
                        \(l2\^2\ m2\ thd\^2\ Cos[th]\ Sin[th]\)\/rad\^2 + 
                        \(2\ l2\^3\ m2\ thd\^2\ Cos[th]\ Sin[th]\^2\)\/\(l3\ 
                            rad\) + 
                        \(l2\^5\ m2\ thd\^2\ Cos[th]\^3\ 
                            Sin[th]\^2\)\/\(l3\^3\ rad\^3\) + 
                        \(l2\^3\ m2\ thd\^2\ Sin[psis]\ Sin[th]\^2\)\/\(l3\ 
                            rad\) + 
                        \(2\ l2\^5\ m2\ thd\^2\ Cos[th]\^2\ Sin[psis]\ 
                            Sin[th]\^2\)\/\(l3\^3\ rad\^3\) + 
                        \(l2\^5\ m2\ thd\^2\ Cos[th]\ Sin[psis]\^2\ 
                            Sin[th]\^2\)\/\(l3\^3\ rad\^3\) - 
                        \(l2\^4\ m2\ thd\^2\ Cos[th]\ Sin[th]\^3\)\/\(l3\^2\ 
                            rad\^4\) - 
                        \(l2\^4\ m2\ thd\^2\ Sin[psis]\ 
                            Sin[th]\^3\)\/\(l3\^2\ rad\^4\) - 
                        g\ l4\ m1\ Sin[phi + th])\))\)/
              \((\((l4\^2\ m1 - l1\ l4\ m1\ Cos[phi])\)\^2 - 
                  l4\^2\ m1\ \((
                    l1\^2\ m1 + l4\^2\ m1 + l2\^2\ m2 + mb\ rg\^2 - 
                      2\ l1\ l4\ m1\ Cos[phi] - 
                      \(2\ l2\^3\ m2\ Cos[th]\^2\ Sin[th]\)\/\(l3\ rad\) - 
                      \(2\ l2\^3\ m2\ Cos[th]\ Sin[psis]\ Sin[th]\)\/\(l3\ rad
                          \) - 2\ l2\^2\ m2\ Sin[th]\^2 + 
                      \(l2\^2\ m2\ Sin[th]\^2\)\/rad\^2)\))\)\)\)}}\)], 
  "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[{
    \(thddsolnum = Expand[Numerator[thdd /. \(Flatten[sol5]\)[\([1]\)]]]\), 
    \(thddsoldenom = 
      Expand[Denominator[thdd /. \(Flatten[sol5]\)[\([1]\)]]]\)}], "Input"],

Cell[BoxData[
    \(\(-l2\^3\)\ l3\^2\ m2\ rad\^3\ thd\^2\ Cos[th]\^3 + 
      l1\ l3\^3\ l4\ m1\ phid\^2\ rad\^4\ Sin[phi] + 
      2\ l1\ l3\^3\ l4\ m1\ phid\ rad\^4\ thd\ Sin[phi] + 
      l1\ l3\^3\ l4\ m1\ rad\^4\ thd\^2\ Sin[phi] - 
      l1\^2\ l3\^3\ m1\ rad\^4\ thd\^2\ Cos[phi]\ Sin[phi] - 
      l2\^3\ l3\^2\ m2\ rad\^3\ thd\^2\ Cos[th]\^2\ Sin[psis] + 
      g\ l1\ l3\^3\ m1\ rad\^4\ Sin[th] + g\ l3\^3\ mb\ rad\^4\ rc\ Sin[th] + 
      l2\^2\ l3\^3\ m2\ rad\^2\ thd\^2\ Cos[th]\ Sin[th] - 
      2\ l2\^2\ l3\^3\ m2\ rad\^4\ thd\^2\ Cos[th]\ Sin[th] + 
      2\ l2\^3\ l3\^2\ m2\ rad\^3\ thd\^2\ Cos[th]\ Sin[th]\^2 + 
      l2\^5\ m2\ rad\ thd\^2\ Cos[th]\^3\ Sin[th]\^2 + 
      l2\^3\ l3\^2\ m2\ rad\^3\ thd\^2\ Sin[psis]\ Sin[th]\^2 + 
      2\ l2\^5\ m2\ rad\ thd\^2\ Cos[th]\^2\ Sin[psis]\ Sin[th]\^2 + 
      l2\^5\ m2\ rad\ thd\^2\ Cos[th]\ Sin[psis]\^2\ Sin[th]\^2 - 
      l2\^4\ l3\ m2\ thd\^2\ Cos[th]\ Sin[th]\^3 - 
      l2\^4\ l3\ m2\ thd\^2\ Sin[psis]\ Sin[th]\^3 - 
      g\ l1\ l3\^3\ m1\ rad\^4\ Cos[phi]\ Sin[phi + th]\)], "Output"],

Cell[BoxData[
    \(\(-l1\^2\)\ l3\^3\ m1\ rad\^4 - l2\^2\ l3\^3\ m2\ rad\^4 - 
      l3\^3\ mb\ rad\^4\ rg\^2 + l1\^2\ l3\^3\ m1\ rad\^4\ Cos[phi]\^2 + 
      2\ l2\^3\ l3\^2\ m2\ rad\^3\ Cos[th]\^2\ Sin[th] + 
      2\ l2\^3\ l3\^2\ m2\ rad\^3\ Cos[th]\ Sin[psis]\ Sin[th] - 
      l2\^2\ l3\^3\ m2\ rad\^2\ Sin[th]\^2 + 
      2\ l2\^2\ l3\^3\ m2\ rad\^4\ Sin[th]\^2\)], "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[{
    \(\n (*\ similarly, \ *) \n\t
    phiddsolnum = Expand[Numerator[phidd /. \(Flatten[sol5]\)[\([2]\)]]]\), 
    \(phiddsoldenom = 
      Expand[Denominator[phidd /. \(Flatten[sol5]\)[\([2]\)]]]\)}], "Input"],

Cell[BoxData[
    \(\(l2\^3\ l4\^2\ m1\ m2\ thd\^2\ Cos[th]\^3\)\/\(l3\ rad\) - 
      \(l1\ l2\^3\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Cos[th]\^3\)\/\(l3\ rad\) - 
      l1\ l4\^3\ m1\^2\ phid\^2\ Sin[phi] - 
      2\ l1\ l4\^3\ m1\^2\ phid\ thd\ Sin[phi] - 
      l1\^3\ l4\ m1\^2\ thd\^2\ Sin[phi] - 
      l1\ l4\^3\ m1\^2\ thd\^2\ Sin[phi] - 
      l1\ l2\^2\ l4\ m1\ m2\ thd\^2\ Sin[phi] - 
      l1\ l4\ m1\ mb\ rg\^2\ thd\^2\ Sin[phi] + 
      l1\^2\ l4\^2\ m1\^2\ phid\^2\ Cos[phi]\ Sin[phi] + 
      2\ l1\^2\ l4\^2\ m1\^2\ phid\ thd\ Cos[phi]\ Sin[phi] + 
      2\ l1\^2\ l4\^2\ m1\^2\ thd\^2\ Cos[phi]\ Sin[phi] + 
      \(l2\^3\ l4\^2\ m1\ m2\ thd\^2\ Cos[th]\^2\ Sin[psis]\)\/\(l3\ rad\) - 
      \(l1\ l2\^3\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Cos[th]\^2\ Sin[psis]\)\/\(l3
          \ rad\) - g\ l1\ l4\^2\ m1\^2\ Sin[th] - 
      g\ l4\^2\ m1\ mb\ rc\ Sin[th] + 
      g\ l1\^2\ l4\ m1\^2\ Cos[phi]\ Sin[th] + 
      g\ l1\ l4\ m1\ mb\ rc\ Cos[phi]\ Sin[th] + 
      2\ l2\^2\ l4\^2\ m1\ m2\ thd\^2\ Cos[th]\ Sin[th] - 
      \(l2\^2\ l4\^2\ m1\ m2\ thd\^2\ Cos[th]\ Sin[th]\)\/rad\^2 - 
      2\ l1\ l2\^2\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Cos[th]\ Sin[th] + 
      \(l1\ l2\^2\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Cos[th]\ Sin[th]\)\/rad\^2 + 
      \(2\ l1\ l2\^3\ l4\ m1\ m2\ thd\^2\ Cos[th]\^2\ Sin[phi]\ 
          Sin[th]\)\/\(l3\ rad\) + 
      \(2\ l1\ l2\^3\ l4\ m1\ m2\ thd\^2\ Cos[th]\ Sin[phi]\ Sin[psis]\ 
          Sin[th]\)\/\(l3\ rad\) - 
      \(2\ l2\^3\ l4\^2\ m1\ m2\ thd\^2\ Cos[th]\ Sin[th]\^2\)\/\(l3\ rad\) + 
      \(2\ l1\ l2\^3\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Cos[th]\ 
          Sin[th]\^2\)\/\(l3\ rad\) - 
      \(l2\^5\ l4\^2\ m1\ m2\ thd\^2\ Cos[th]\^3\ Sin[th]\^2\)\/\(l3\^3\ 
          rad\^3\) + 
      \(l1\ l2\^5\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Cos[th]\^3\ 
          Sin[th]\^2\)\/\(l3\^3\ rad\^3\) + 
      2\ l1\ l2\^2\ l4\ m1\ m2\ thd\^2\ Sin[phi]\ Sin[th]\^2 - 
      \(l1\ l2\^2\ l4\ m1\ m2\ thd\^2\ Sin[phi]\ Sin[th]\^2\)\/rad\^2 - 
      \(l2\^3\ l4\^2\ m1\ m2\ thd\^2\ Sin[psis]\ Sin[th]\^2\)\/\(l3\ rad\) + 
      \(l1\ l2\^3\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Sin[psis]\ Sin[th]\^2\)\/\(l3
          \ rad\) - 
      \(2\ l2\^5\ l4\^2\ m1\ m2\ thd\^2\ Cos[th]\^2\ Sin[psis]\ 
          Sin[th]\^2\)\/\(l3\^3\ rad\^3\) + 
      \(2\ l1\ l2\^5\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Cos[th]\^2\ Sin[psis]\ 
          Sin[th]\^2\)\/\(l3\^3\ rad\^3\) - 
      \(l2\^5\ l4\^2\ m1\ m2\ thd\^2\ Cos[th]\ Sin[psis]\^2\ 
          Sin[th]\^2\)\/\(l3\^3\ rad\^3\) + 
      \(l1\ l2\^5\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Cos[th]\ Sin[psis]\^2\ 
          Sin[th]\^2\)\/\(l3\^3\ rad\^3\) + 
      \(l2\^4\ l4\^2\ m1\ m2\ thd\^2\ Cos[th]\ Sin[th]\^3\)\/\(l3\^2\ 
          rad\^4\) - 
      \(l1\ l2\^4\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Cos[th]\ 
          Sin[th]\^3\)\/\(l3\^2\ rad\^4\) + 
      \(l2\^4\ l4\^2\ m1\ m2\ thd\^2\ Sin[psis]\ Sin[th]\^3\)\/\(l3\^2\ 
          rad\^4\) - 
      \(l1\ l2\^4\ l4\ m1\ m2\ thd\^2\ Cos[phi]\ Sin[psis]\ 
          Sin[th]\^3\)\/\(l3\^2\ rad\^4\) - 
      g\ l1\^2\ l4\ m1\^2\ Sin[phi + th] - 
      g\ l2\^2\ l4\ m1\ m2\ Sin[phi + th] - 
      g\ l4\ m1\ mb\ rg\^2\ Sin[phi + th] + 
      g\ l1\ l4\^2\ m1\^2\ Cos[phi]\ Sin[phi + th] + 
      \(2\ g\ l2\^3\ l4\ m1\ m2\ Cos[th]\^2\ Sin[th]\ Sin[phi + th]\)\/\(l3\ 
          rad\) + \(2\ g\ l2\^3\ l4\ m1\ m2\ Cos[th]\ Sin[psis]\ Sin[th]\ 
          Sin[phi + th]\)\/\(l3\ rad\) + 
      2\ g\ l2\^2\ l4\ m1\ m2\ Sin[th]\^2\ Sin[phi + th] - 
      \(g\ l2\^2\ l4\ m1\ m2\ Sin[th]\^2\ Sin[phi + th]\)\/rad\^2\)], "Output"],

Cell[BoxData[
    \(\(-l1\^2\)\ l4\^2\ m1\^2 - l2\^2\ l4\^2\ m1\ m2 - 
      l4\^2\ m1\ mb\ rg\^2 + l1\^2\ l4\^2\ m1\^2\ Cos[phi]\^2 + 
      \(2\ l2\^3\ l4\^2\ m1\ m2\ Cos[th]\^2\ Sin[th]\)\/\(l3\ rad\) + 
      \(2\ l2\^3\ l4\^2\ m1\ m2\ Cos[th]\ Sin[psis]\ Sin[th]\)\/\(l3\ rad\) + 
      2\ l2\^2\ l4\^2\ m1\ m2\ Sin[th]\^2 - 
      \(l2\^2\ l4\^2\ m1\ m2\ Sin[th]\^2\)\/rad\^2\)], "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(FortranForm[thddsolnum]\)], "Input"],

Cell["\<\
        -(l2**3*l3**2*m2*rad**3*thd**2*Cos(th)**3) + 
     -  l1*l3**3*l4*m1*phid**2*rad**4*Sin(phi) + 
     -  2*l1*l3**3*l4*m1*phid*rad**4*thd*Sin(phi) + 
     -  l1*l3**3*l4*m1*rad**4*thd**2*Sin(phi) - 
     -  l1**2*l3**3*m1*rad**4*thd**2*Cos(phi)*Sin(phi) - 
     -  l2**3*l3**2*m2*rad**3*thd**2*Cos(th)**2*Sin(psis) + 
     -  g*l1*l3**3*m1*rad**4*Sin(th) + g*l3**3*mb*rad**4*rc*Sin(th) + 
     -  l2**2*l3**3*m2*rad**2*thd**2*Cos(th)*Sin(th) - 
     -  2*l2**2*l3**3*m2*rad**4*thd**2*Cos(th)*Sin(th) + 
     -  2*l2**3*l3**2*m2*rad**3*thd**2*Cos(th)*Sin(th)**2 + 
     -  l2**5*m2*rad*thd**2*Cos(th)**3*Sin(th)**2 + 
     -  l2**3*l3**2*m2*rad**3*thd**2*Sin(psis)*Sin(th)**2 + 
     -  2*l2**5*m2*rad*thd**2*Cos(th)**2*Sin(psis)*Sin(th)**2 + 
     -  l2**5*m2*rad*thd**2*Cos(th)*Sin(psis)**2*Sin(th)**2 - 
     -  l2**4*l3*m2*thd**2*Cos(th)*Sin(th)**3 - 
     -  l2**4*l3*m2*thd**2*Sin(psis)*Sin(th)**3 - 
     -  g*l1*l3**3*m1*rad**4*Cos(phi)*Sin(phi + th)\
\>", "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(FortranForm[thddsoldenom]\)], "Input"],

Cell["\<\
        -(l1**2*l3**3*m1*rad**4) - l2**2*l3**3*m2*rad**4 - 
     -  l3**3*mb*rad**4*rg**2 + l1**2*l3**3*m1*rad**4*Cos(phi)**2 + 
     -  2*l2**3*l3**2*m2*rad**3*Cos(th)**2*Sin(th) + 
     -  2*l2**3*l3**2*m2*rad**3*Cos(th)*Sin(psis)*Sin(th) - 
     -  l2**2*l3**3*m2*rad**2*Sin(th)**2 + 2*l2**2*l3**3*m2*rad**4*Sin(th)**2\
\
\>", "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\(FortranForm[phiddsolnum]\n\)\)], "Input"],

Cell["\<\
        (l2**3*l4**2*m1*m2*thd**2*Cos(th)**3)/(l3*rad) - 
     -  (l1*l2**3*l4*m1*m2*thd**2*Cos(phi)*Cos(th)**3)/(l3*rad) - 
     -  l1*l4**3*m1**2*phid**2*Sin(phi) - 2*l1*l4**3*m1**2*phid*thd*Sin(phi) \
- 
     -  l1**3*l4*m1**2*thd**2*Sin(phi) - l1*l4**3*m1**2*thd**2*Sin(phi) - 
     -  l1*l2**2*l4*m1*m2*thd**2*Sin(phi) - l1*l4*m1*mb*rg**2*thd**2*Sin(phi) \
+ 
     -  l1**2*l4**2*m1**2*phid**2*Cos(phi)*Sin(phi) + 
     -  2*l1**2*l4**2*m1**2*phid*thd*Cos(phi)*Sin(phi) + 
     -  2*l1**2*l4**2*m1**2*thd**2*Cos(phi)*Sin(phi) + 
     -  (l2**3*l4**2*m1*m2*thd**2*Cos(th)**2*Sin(psis))/(l3*rad) - 
     -  (l1*l2**3*l4*m1*m2*thd**2*Cos(phi)*Cos(th)**2*Sin(psis))/(l3*rad) - 
     -  g*l1*l4**2*m1**2*Sin(th) - g*l4**2*m1*mb*rc*Sin(th) + 
     -  g*l1**2*l4*m1**2*Cos(phi)*Sin(th) + g*l1*l4*m1*mb*rc*Cos(phi)*Sin(th) \
+ 
     -  2*l2**2*l4**2*m1*m2*thd**2*Cos(th)*Sin(th) - 
     -  (l2**2*l4**2*m1*m2*thd**2*Cos(th)*Sin(th))/rad**2 - 
     -  2*l1*l2**2*l4*m1*m2*thd**2*Cos(phi)*Cos(th)*Sin(th) + 
     -  (l1*l2**2*l4*m1*m2*thd**2*Cos(phi)*Cos(th)*Sin(th))/rad**2 + 
     -  (2*l1*l2**3*l4*m1*m2*thd**2*Cos(th)**2*Sin(phi)*Sin(th))/(l3*rad) + 
     -  (2*l1*l2**3*l4*m1*m2*thd**2*Cos(th)*Sin(phi)*Sin(psis)*Sin(th))/
     -   (l3*rad) - (2*l2**3*l4**2*m1*m2*thd**2*Cos(th)*Sin(th)**2)/(l3*rad) \
+ 
     -  (2*l1*l2**3*l4*m1*m2*thd**2*Cos(phi)*Cos(th)*Sin(th)**2)/(l3*rad) - 
     -  (l2**5*l4**2*m1*m2*thd**2*Cos(th)**3*Sin(th)**2)/(l3**3*rad**3) + 
     -  (l1*l2**5*l4*m1*m2*thd**2*Cos(phi)*Cos(th)**3*Sin(th)**2)/
     -   (l3**3*rad**3) + 2*l1*l2**2*l4*m1*m2*thd**2*Sin(phi)*Sin(th)**2 - 
     -  (l1*l2**2*l4*m1*m2*thd**2*Sin(phi)*Sin(th)**2)/rad**2 - 
     -  (l2**3*l4**2*m1*m2*thd**2*Sin(psis)*Sin(th)**2)/(l3*rad) + 
     -  (l1*l2**3*l4*m1*m2*thd**2*Cos(phi)*Sin(psis)*Sin(th)**2)/(l3*rad) - 
     -  (2*l2**5*l4**2*m1*m2*thd**2*Cos(th)**2*Sin(psis)*Sin(th)**2)/
     -   (l3**3*rad**3) + (2*l1*l2**5*l4*m1*m2*thd**2*Cos(phi)*Cos(th)**2*
     -     Sin(psis)*Sin(th)**2)/(l3**3*rad**3) - 
     -  (l2**5*l4**2*m1*m2*thd**2*Cos(th)*Sin(psis)**2*Sin(th)**2)/
     -   (l3**3*rad**3) + \
(l1*l2**5*l4*m1*m2*thd**2*Cos(phi)*Cos(th)*Sin(psis)**2*
     -     Sin(th)**2)/(l3**3*rad**3) + 
     -  (l2**4*l4**2*m1*m2*thd**2*Cos(th)*Sin(th)**3)/(l3**2*rad**4) - 
     -  (l1*l2**4*l4*m1*m2*thd**2*Cos(phi)*Cos(th)*Sin(th)**3)/(l3**2*rad**4) \
+ 
     -  (l2**4*l4**2*m1*m2*thd**2*Sin(psis)*Sin(th)**3)/(l3**2*rad**4) - 
     -  (l1*l2**4*l4*m1*m2*thd**2*Cos(phi)*Sin(psis)*Sin(th)**3)/(l3**2*rad**\
4) - 
     -  g*l1**2*l4*m1**2*Sin(phi + th) - g*l2**2*l4*m1*m2*Sin(phi + th) - 
     -  g*l4*m1*mb*rg**2*Sin(phi + th) + 
     -  g*l1*l4**2*m1**2*Cos(phi)*Sin(phi + th) + 
     -  (2*g*l2**3*l4*m1*m2*Cos(th)**2*Sin(th)*Sin(phi + th))/(l3*rad) + 
     -  (2*g*l2**3*l4*m1*m2*Cos(th)*Sin(psis)*Sin(th)*Sin(phi + th))/(l3*rad) \
+ 
     -  2*g*l2**2*l4*m1*m2*Sin(th)**2*Sin(phi + th) - 
     -  (g*l2**2*l4*m1*m2*Sin(th)**2*Sin(phi + th))/rad**2\
\>", "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(FortranForm[phiddsoldenom]\)], "Input"],

Cell["\<\
        -(l1**2*l4**2*m1**2) - l2**2*l4**2*m1*m2 - \
l4**2*m1*mb*rg**2 + 
     -  l1**2*l4**2*m1**2*Cos(phi)**2 + 
     -  (2*l2**3*l4**2*m1*m2*Cos(th)**2*Sin(th))/(l3*rad) + 
     -  (2*l2**3*l4**2*m1*m2*Cos(th)*Sin(psis)*Sin(th))/(l3*rad) + 
     -  2*l2**2*l4**2*m1*m2*Sin(th)**2 - \
(l2**2*l4**2*m1*m2*Sin(th)**2)/rad**2\
\>", "Output"]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(\( (*\ where\ rad\ is\ *) \n\t\t
    FortranForm[\@\(1 - \(l2\^2\ \((Cos[th] + Sin[psis])\)\^2\)\/l3\^2\)]
      \)\)], "Input"],

Cell["Sqrt(1 - (l2**2*(Cos(th) + Sin(psis))**2)/l3**2)", "Output"]
}, Open  ]],

Cell[BoxData[
    \( (*\ 
      These\ last\ five\ expressions\ are\ further\ manipulated\ in\ a\ word\ 
        processor\ to\ conform\ to\ the\ target\ language.\ \n\tIf\ the\ 
        target\ languange\ is\ C\ or\ one\ of\ its\ relatives, \ 
      then\ the\ labor\ can\ be\ shortened\ somewhat\ by\ using\ CForm\ rather
        \ than\ FortranForm, \ of\ \(course.\)\ \ \n*) \)], "Input"],

Cell[BoxData[
    \( (*\n\t\n\t\tThus, \ 
      the\ equations\ for\ the\ second\ derivatives\ of\ the\ angles\ are\ all
          \ derived\ \((in\ each\ of\ the\ two\ movement\ regimes)\)\ for\ 
          input\ to\ a\ runge - kutta\ solving\ \(subroutine.\)\n\t\n\t*) \)],
   "Input"]
},
FrontEndVersion->"Macintosh 3.0",
ScreenRectangle->{{0, 832}, {0, 604}},
WindowSize->{667, 499},
WindowMargins->{{Automatic, 23}, {Automatic, 0}},
MacintoshSystemPageSetup->"\<\
00<0001804P000000]P2:?oQon82n@960dL5:0?l0080001804P000000]P2:001
0000I00000400`<300000BL?00400@0000000000000006P801T1T00000000000
00000000000000000000000000000000\>"
]


(***********************************************************************
Cached data follows.  If you edit this Notebook file directly, not using
Mathematica, you must remove the line containing CacheID at the top of 
the file.  The cache data will then be recreated when you save this file 
from within Mathematica.
***********************************************************************)

(*CellTagsOutline
CellTagsIndex->{}
*)

(*CellTagsIndex
CellTagsIndex->{}
*)

(*NotebookFileOutline
Notebook[{
Cell[1709, 49, 4301, 71, 1046, "Input"],
Cell[6013, 122, 730, 14, 191, "Input"],

Cell[CellGroupData[{
Cell[6768, 140, 595, 11, 116, "Input"],
Cell[7366, 153, 598, 11, 153, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[8001, 169, 513, 9, 101, "Input"],
Cell[8517, 180, 574, 8, 109, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[9128, 193, 950, 15, 191, "Input"],
Cell[10081, 210, 156, 3, 28, "Output"],
Cell[10240, 215, 156, 3, 28, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[10433, 223, 206, 5, 41, "Input"],
Cell[10642, 230, 3307, 57, 411, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[13986, 292, 276, 4, 101, "Input"],
Cell[14265, 298, 440, 7, 108, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[14742, 310, 84, 1, 41, "Input"],
Cell[14829, 313, 107, 3, 38, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[14973, 321, 218, 3, 101, "Input"],
Cell[15194, 326, 1253, 24, 304, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[16484, 355, 80, 1, 26, "Input"],
Cell[16567, 358, 131, 4, 40, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[16735, 367, 157, 2, 71, "Input"],
Cell[16895, 371, 1248, 23, 304, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[18180, 399, 80, 1, 26, "Input"],
Cell[18263, 402, 139, 4, 52, "Output"]
}, Open  ]],
Cell[18417, 409, 478, 9, 101, "Input"],
Cell[18898, 420, 413, 7, 86, "Input"],

Cell[CellGroupData[{
Cell[19336, 431, 359, 6, 146, "Input"],
Cell[19698, 439, 145, 4, 68, "Output"],
Cell[19846, 445, 35, 1, 26, "Output"],
Cell[19884, 448, 35, 1, 26, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[19956, 454, 619, 9, 161, "Input"],
Cell[20578, 465, 1275, 26, 227, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[21890, 496, 822, 15, 262, "Input"],
Cell[22715, 513, 176, 3, 37, "Message"],
Cell[22894, 518, 1457, 25, 243, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[24388, 548, 698, 11, 176, "Input"],
Cell[25089, 561, 156, 3, 28, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[25282, 569, 424, 8, 101, "Input"],
Cell[25709, 579, 111, 2, 41, "Output"],
Cell[25823, 583, 147, 4, 68, "Output"],
Cell[25973, 589, 424, 8, 125, "Output"]
}, Open  ]],
Cell[26412, 600, 478, 11, 149, "Input"],

Cell[CellGroupData[{
Cell[26915, 615, 813, 16, 266, "Input"],
Cell[27731, 633, 3899, 67, 821, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[31667, 705, 158, 2, 41, "Input"],
Cell[31828, 709, 1690, 28, 407, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[33555, 742, 134, 2, 41, "Input"],
Cell[33692, 746, 4355, 70, 701, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[38084, 821, 189, 3, 41, "Input"],
Cell[38276, 826, 1067, 17, 164, "Output"],
Cell[39346, 845, 383, 6, 62, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[39766, 856, 225, 4, 71, "Input"],
Cell[39994, 862, 3527, 62, 882, "Output"],
Cell[43524, 926, 395, 6, 98, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[43956, 937, 56, 1, 26, "Input"],
Cell[44015, 940, 990, 18, 262, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[45042, 963, 58, 1, 26, "Input"],
Cell[45103, 966, 349, 7, 94, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[45489, 978, 63, 1, 41, "Input"],
Cell[45555, 981, 2978, 52, 626, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[48570, 1038, 59, 1, 26, "Input"],
Cell[48632, 1041, 350, 8, 94, "Output"]
}, Open  ]],

Cell[CellGroupData[{
Cell[49019, 1054, 148, 3, 67, "Input"],
Cell[49170, 1059, 66, 0, 38, "Output"]
}, Open  ]],
Cell[49251, 1062, 392, 6, 101, "Input"],
Cell[49646, 1070, 287, 5, 116, "Input"]
}
]
*)




(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)