57   using namespace NEWMAT;
 
   62                       const ColumnVector & qpp, 
const ColumnVector & dq,
 
   63                       ColumnVector & ltorque, ColumnVector & dtorque)
 
   73    if(q.Ncols() != 1 || q.Nrows() != dof) error(
"q has wrong dimension");
 
   74    if(qp.Ncols() != 1 || qp.Nrows() != dof) error(
"qp has wrong dimension");
 
   75    if(qpp.Ncols() != 1 || qpp.Nrows() != dof) error(
"qpp has wrong dimension");
 
   76    if(dq.Ncols() != 1 || qpp.Nrows() != dof) error(
"dq has wrong dimension");
 
   77    ltorque = ColumnVector(dof);
 
   78    dtorque = ColumnVector(dof);
 
   83    z0(1) = 0.0; z0(2) = 0.0;  z0(3) = 1.0;
 
   89    for(i = 1; i <= dof; i++)
 
   92       p[i] = ColumnVector(3);
 
   93       p[i](1) = links[i].get_a();
 
   94       p[i](2) = links[i].get_d() * Rt(2,3);
 
   95       p[i](3) = links[i].get_d() * Rt(3,3);
 
   96       if(links[i].get_joint_type() != 0)
 
   98          dp[i] = ColumnVector(3);
 
  103       if(links[i].get_joint_type() == 0)
 
  105          w[i] = Rt*(w[i-1] + z0*qp(i));
 
  106          dw[i] = Rt*(dw[i-1] - Q*w[i-1]*dq(i));
 
  107          wp[i] = Rt*(wp[i-1] + z0*qpp(i)
 
  108                      + CrossProduct(w[i-1],z0*qp(i)));
 
  109          dwp[i] = Rt*(dwp[i-1]
 
  110                       + CrossProduct(dw[i-1],z0*qp(i))
 
  111                       - Q*(wp[i-1]+z0*qpp(i)+CrossProduct(w[i-1],z0*qp(i)))
 
  113          vp[i] = CrossProduct(wp[i],p[i])
 
  114                  + CrossProduct(w[i],CrossProduct(w[i],p[i]))
 
  116          dvp[i] = CrossProduct(dwp[i],p[i])
 
  117                   + CrossProduct(dw[i],CrossProduct(w[i],p[i]))
 
  118                   + CrossProduct(w[i],CrossProduct(dw[i],p[i]))
 
  119                   + Rt*(dvp[i-1] - Q*vp[i-1]*dq(i));
 
  126          dwp[i] = Rt*dwp[i-1];
 
  127          vp[i] = Rt*(vp[i-1] + z0*qpp(i)
 
  128                      + 2.0*CrossProduct(w[i],Rt*z0*qp(i)))
 
  129                  + CrossProduct(wp[i],p[i])
 
  130                  + CrossProduct(w[i],CrossProduct(w[i],p[i]));
 
  131          dvp[i] = Rt*(dvp[i-1]
 
  132                       + 2.0*CrossProduct(dw[i-1],z0*qp(i)))
 
  133                   + CrossProduct(dwp[i],p[i])
 
  134                   + CrossProduct(dw[i],CrossProduct(w[i],p[i]))
 
  135                   + CrossProduct(w[i],CrossProduct(dw[i],p[i]))
 
  136                   + (CrossProduct(wp[i],dp[i])
 
  137                      + CrossProduct(w[i],CrossProduct(w[i],dp[i])))
 
  140       a[i] = CrossProduct(wp[i],links[i].r)
 
  141              + CrossProduct(w[i],CrossProduct(w[i],links[i].r))
 
  143       da[i] = CrossProduct(dwp[i],links[i].r)
 
  144               + CrossProduct(dw[i],CrossProduct(w[i],links[i].r))
 
  145               + CrossProduct(w[i],CrossProduct(dw[i],links[i].r))
 
  149    for(i = dof; i >= 1; i--)
 
  151       F[i] = a[i] * links[i].m;
 
  152       N[i] = links[i].I*wp[i] + CrossProduct(w[i],links[i].I*w[i]);
 
  153       dF[i] = da[i] * links[i].m;
 
  154       dN[i] = links[i].I*dwp[i] + CrossProduct(dw[i],links[i].I*w[i])
 
  155               + CrossProduct(w[i],links[i].I*dw[i]);
 
  160          n[i] = CrossProduct(p[i],f[i])
 
  161                 + CrossProduct(links[i].r,F[i]) + N[i];
 
  162          dn[i] = CrossProduct(p[i],df[i])
 
  163                  + CrossProduct(links[i].r,dF[i]) + dN[i];
 
  164          if(links[i].get_joint_type() != 0)
 
  165             dn[i] += CrossProduct(dp[i],f[i])*dq(i);
 
  169          f[i] = links[i+1].R*f[i+1] + F[i];
 
  170          df[i] = links[i+1].R*df[i+1] + dF[i];
 
  171          if(links[i+1].get_joint_type() == 0)
 
  172             df[i] += Q*links[i+1].R*f[i+1]*dq(i+1);
 
  173          n[i] = links[i+1].R*n[i+1] + CrossProduct(p[i],f[i])
 
  174                 + CrossProduct(links[i].r,F[i]) + N[i];
 
  176          dn[i] = links[i+1].R*dn[i+1] + CrossProduct(p[i],df[i])
 
  177                  + CrossProduct(links[i].r,dF[i]) + dN[i];
 
  178          if(links[i+1].get_joint_type() == 0)
 
  179             dn[i] += Q*links[i+1].R*n[i+1]*dq(i+1);
 
  181             dn[i] += CrossProduct(dp[i],f[i])*dq(i);
 
  183       if(links[i].get_joint_type() == 0)
 
  185          temp = ((z0.t()*links[i].R)*n[i]);
 
  186          ltorque(i) = temp(1,1);
 
  187          temp = ((z0.t()*links[i].R)*dn[i]);
 
  188          dtorque(i) = temp(1,1);
 
  192          temp = ((z0.t()*links[i].R)*f[i]);
 
  193          ltorque(i) = temp(1,1);
 
  194          temp = ((z0.t()*links[i].R)*df[i]);
 
  195          dtorque(i) = temp(1,1);
 
  202                        const ColumnVector & qpp, 
const ColumnVector & dq,
 
  203                        ColumnVector & ltorque, ColumnVector & dtorque)
 
  213    if(q.Ncols() != 1 || q.Nrows() != dof) error(
"q has wrong dimension");
 
  214    if(qp.Ncols() != 1 || qp.Nrows() != dof) error(
"qp has wrong dimension");
 
  215    if(qpp.Ncols() != 1 || qpp.Nrows() != dof) error(
"qpp has wrong dimension");
 
  216    if(dq.Ncols() != 1 || dq.Nrows() != dof) error(
"dq has wrong dimension");
 
  217    ltorque = ColumnVector(dof);
 
  218    dtorque = ColumnVector(dof);
 
  223    z0(1) = 0.0; z0(2) = 0.0; z0(3) = 1.0;
 
  228    for(i = 1; i <= dof; i++)
 
  232       if(links[i].get_joint_type() != 0)
 
  234          dp[i] = ColumnVector(3);
 
  239       if(links[i].get_joint_type() == 0)
 
  241          w[i] = Rt*w[i-1] + z0*qp(i);
 
  242          dw[i] = Rt*dw[i-1] - Q*Rt*w[i-1]*dq(i);
 
  243          wp[i] = Rt*wp[i-1] + CrossProduct(Rt*w[i-1],z0*qp(i))
 
  245          dwp[i] = Rt*dwp[i-1] + CrossProduct(Rt*dw[i-1],z0*qp(i))
 
  246                   - (Q*Rt*wp[i-1] + CrossProduct(Q*Rt*w[i-1],z0*qp(i)))*dq(i);
 
  247          vp[i] = Rt*(CrossProduct(wp[i-1],p[i])
 
  248                      + CrossProduct(w[i-1],CrossProduct(w[i-1],p[i]))
 
  250          dvp[i] = Rt*(CrossProduct(dwp[i-1],p[i])
 
  251                       + CrossProduct(dw[i-1],CrossProduct(w[i-1],p[i]))
 
  252                       + CrossProduct(w[i-1],CrossProduct(dw[i-1],p[i])) + dvp[i-1])
 
  253                   - Q*Rt*(CrossProduct(wp[i-1],p[i])
 
  254                           + CrossProduct(w[i-1],CrossProduct(w[i-1],p[i])) + vp[i-1])*dq(i);
 
  261          dwp[i] = Rt*dwp[i-1];
 
  262          vp[i] = Rt*(vp[i-1] + CrossProduct(wp[i-1],p[i])
 
  263                      + CrossProduct(w[i-1],CrossProduct(w[i-1],p[i])))
 
  264                  + z0*qpp(i)+ 2.0*CrossProduct(w[i],z0*qp(i));
 
  265          dvp[i] = Rt*(CrossProduct(dwp[i-1],p[i])
 
  266                       + CrossProduct(dw[i-1],CrossProduct(w[i-1],p[i]))
 
  267                       + CrossProduct(w[i-1],CrossProduct(dw[i-1],p[i])) + dvp[i-1]
 
  268                       + (CrossProduct(wp[i-1],dp[i])
 
  269                          + CrossProduct(w[i-1],CrossProduct(w[i-1],dp[i])))*dq(i))
 
  270                   + 2*CrossProduct(dw[i],z0*qp(i));
 
  272       a[i] = CrossProduct(wp[i],links[i].r)
 
  273              + CrossProduct(w[i],CrossProduct(w[i],links[i].r))
 
  275       da[i] = CrossProduct(dwp[i],links[i].r)
 
  276               + CrossProduct(dw[i],CrossProduct(w[i],links[i].r))
 
  277               + CrossProduct(w[i],CrossProduct(dw[i],links[i].r))
 
  281    for(i = dof; i >= 1; i--) {
 
  282       F[i] = a[i] * links[i].m;
 
  283       N[i] = links[i].I*wp[i] + CrossProduct(w[i],links[i].I*w[i]);
 
  284       dF[i] = da[i] * links[i].m;
 
  285       dN[i] = links[i].I*dwp[i] + CrossProduct(dw[i],links[i].I*w[i])
 
  286               + CrossProduct(w[i],links[i].I*dw[i]);
 
  292          n[i] = CrossProduct(links[i].r,F[i]) + N[i];
 
  293          dn[i] = dN[i] + CrossProduct(links[i].r,dF[i]);
 
  297          f[i] = links[i+1].R*f[i+1] + F[i];
 
  298          df[i] = links[i+1].R*df[i+1] + dF[i];
 
  299          if(links[i+1].get_joint_type() == 0)
 
  300             df[i] += links[i+1].R*Q*f[i+1]*dq(i+1);
 
  302          n[i] = links[i+1].R*n[i+1] + CrossProduct(p[i+1],links[i+1].R*f[i+1])
 
  303                 + CrossProduct(links[i].r,F[i]) + N[i];
 
  304          dn[i] = links[i+1].R*dn[i+1] + CrossProduct(p[i+1],links[i+1].R*df[i+1])
 
  305                  + CrossProduct(links[i].r,dF[i]) + dN[i];
 
  306          if(links[i+1].get_joint_type() == 0)
 
  307             dn[i] += (links[i+1].R*Q*n[i+1] +
 
  308                       CrossProduct(p[i+1],links[i+1].R*Q*f[i+1]))*dq(i+1);
 
  310             dn[i] += CrossProduct(dp[i+1],links[i+1].R*f[i+1])*dq(i+1);
 
  313       if(links[i].get_joint_type() == 0)
 
  316          ltorque(i) = temp(1,1);
 
  318          dtorque(i) = temp(1,1);
 
  323          ltorque(i) = temp(1,1);
 
  325          dtorque(i) = temp(1,1);
 
  332                                 const ColumnVector & qpp, 
const ColumnVector & dq,
 
  333                                 ColumnVector & ltorque, ColumnVector & dtorque)
 
  343    if(q.Ncols() != 1 || q.Nrows() != dof) error(
"q has wrong dimension");
 
  344    if(qp.Ncols() != 1 || qp.Nrows() != dof) error(
"qp has wrong dimension");
 
  345    if(qpp.Ncols() != 1 || qpp.Nrows() != dof) error(
"qpp has wrong dimension");
 
  346    if(dq.Ncols() != 1 || dq.Nrows() != dof) error(
"dq has wrong dimension");
 
  347    ltorque = ColumnVector(dof);
 
  348    dtorque = ColumnVector(dof);
 
  353    z0(1) = 0.0; z0(2) = 0.0; z0(3) = 1.0;
 
  358    for(i = 1; i <= dof; i++)
 
  362       if(links[i].get_joint_type() != 0)
 
  364          dp[i] = ColumnVector(3);
 
  369       if(links[i].get_joint_type() == 0)
 
  371          w[i] = Rt*w[i-1] + z0*qp(i);
 
  372          dw[i] = Rt*dw[i-1] - Q*Rt*w[i-1]*dq(i);
 
  373          wp[i] = Rt*wp[i-1] + CrossProduct(Rt*w[i-1],z0*qp(i))
 
  375          dwp[i] = Rt*dwp[i-1] + CrossProduct(Rt*dw[i-1],z0*qp(i))
 
  376                   - (Q*Rt*wp[i-1] + CrossProduct(Q*Rt*w[i-1],z0*qp(i)))*dq(i);
 
  377          vp[i] = Rt*(CrossProduct(wp[i-1],p[i])
 
  378                      + CrossProduct(w[i-1],CrossProduct(w[i-1],p[i]))
 
  380          dvp[i] = Rt*(CrossProduct(dwp[i-1],p[i])
 
  381                       + CrossProduct(dw[i-1],CrossProduct(w[i-1],p[i]))
 
  382                       + CrossProduct(w[i-1],CrossProduct(dw[i-1],p[i])) + dvp[i-1])
 
  383                   - Q*Rt*(CrossProduct(wp[i-1],p[i])
 
  384                           + CrossProduct(w[i-1],CrossProduct(w[i-1],p[i])) + vp[i-1])*dq(i);
 
  391          dwp[i] = Rt*dwp[i-1];
 
  392          vp[i] = Rt*(vp[i-1] + CrossProduct(wp[i-1],p[i])
 
  393                      + CrossProduct(w[i-1],CrossProduct(w[i-1],p[i])))
 
  394                  + z0*qpp(i)+ 2.0*CrossProduct(w[i],z0*qp(i));
 
  395          dvp[i] = Rt*(CrossProduct(dwp[i-1],p[i])
 
  396                       + CrossProduct(dw[i-1],CrossProduct(w[i-1],p[i]))
 
  397                       + CrossProduct(w[i-1],CrossProduct(dw[i-1],p[i])) + dvp[i-1]
 
  398                       + (CrossProduct(wp[i-1],dp[i])
 
  399                          + CrossProduct(w[i-1],CrossProduct(w[i-1],dp[i])))*dq(i));
 
  403    for(i = dof; i >= 1; i--) {
 
  404       F[i] = vp[i]*links[i].m + CrossProduct(wp[i], links[i].mc) +
 
  405              CrossProduct(w[i], CrossProduct(w[i], links[i].mc));
 
  406       dF[i] = dvp[i]*links[i].m + CrossProduct(dwp[i],links[i].mc)
 
  407               + CrossProduct(dw[i],CrossProduct(w[i],links[i].mc))
 
  408               + CrossProduct(w[i],CrossProduct(dw[i],links[i].mc));
 
  409       N[i] = links[i].I*wp[i] + CrossProduct(w[i],links[i].I*w[i])
 
  410              - CrossProduct(vp[i], links[i].mc);
 
  411       dN[i] = links[i].I*dwp[i] + CrossProduct(dw[i],links[i].I*w[i])
 
  412               + CrossProduct(w[i],links[i].I*dw[i])
 
  413               - CrossProduct(dvp[i],links[i].mc);
 
  424          f[i] = links[i+1].R*f[i+1] + F[i];
 
  425          df[i] = links[i+1].R*df[i+1] + dF[i];
 
  426          if(links[i+1].get_joint_type() == 0)
 
  427             df[i] += links[i+1].R*Q*f[i+1]*dq(i+1);
 
  429          n[i] = links[i+1].R*n[i+1] + CrossProduct(p[i+1],links[i+1].R*f[i+1])
 
  431          dn[i] = links[i+1].R*dn[i+1] + CrossProduct(p[i+1],links[i+1].R*df[i+1])
 
  433          if(links[i+1].get_joint_type() == 0)
 
  434             dn[i] += (links[i+1].R*Q*n[i+1] +
 
  435                       CrossProduct(p[i+1],links[i+1].R*Q*f[i+1]))*dq(i+1);
 
  437             dn[i] += CrossProduct(dp[i+1],links[i+1].R*f[i+1])*dq(i+1);
 
  440       if(links[i].get_joint_type() == 0)
 
  443          ltorque(i) = temp(1,1);
 
  445          dtorque(i) = temp(1,1);
 
  450          ltorque(i) = temp(1,1);
 
  452          dtorque(i) = temp(1,1);