1 #ifndef OPENSIM_MUSCLE_H_
2 #define OPENSIM_MUSCLE_H_
3 /* -------------------------------------------------------------------------- *
4  *                             OpenSim:  Muscle.h                             *
5  * -------------------------------------------------------------------------- *
6  * The OpenSim API is a toolkit for musculoskeletal modeling and simulation.  *
7  * See http://opensim.stanford.edu and the NOTICE file for more information.  *
8  * OpenSim is developed at Stanford University and supported by the US        *
9  * National Institutes of Health (U54 GM072970, R24 HD065690) and by DARPA    *
10  * through the Warrior Web program.                                           *
11  *                                                                            *
12  * Copyright (c) 2005-2017 Stanford University and the Authors                *
13  * Author(s): Ajay Seth, Matthew Millard                                      *
14  *                                                                            *
15  * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
16  * not use this file except in compliance with the License. You may obtain a  *
17  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
18  *                                                                            *
19  * Unless required by applicable law or agreed to in writing, software        *
20  * distributed under the License is distributed on an "AS IS" BASIS,          *
21  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
22  * See the License for the specific language governing permissions and        *
23  * limitations under the License.                                             *
24  * -------------------------------------------------------------------------- */
25 
26 
27 // INCLUDE
28 #include "PathActuator.h"
29 
30 #ifdef SWIG
31     #ifdef OSIMSIMULATION_API
32         #undef OSIMSIMULATION_API
33         #define OSIMSIMULATION_API
34     #endif
35 #endif
36 
37 namespace OpenSim {
38 
39 //==============================================================================
40 //                              Muscle Exceptions
41 //==============================================================================
42 class MuscleCannotEquilibrate : public Exception {
43 public:
MuscleCannotEquilibrate(const std::string & file,size_t line,const std::string & func,const Object & obj,const std::string & detail)44     MuscleCannotEquilibrate(const std::string& file,
45                             size_t line,
46                             const std::string& func,
47                             const Object& obj,
48                             const std::string& detail) :
49         Exception(file, line, func, obj, detail) {
50         std::string msg = "Unable to compute equilibrium for this muscle.\n";
51         msg += "Please verify that the initial activation is valid and that ";
52         msg += "the length of the musculotendon actuator doesn't produce a ";
53         msg += "pennation angle of 90 degrees or a negative fiber length.";
54         addMessage(msg);
55     }
56 };
57 
58 //==============================================================================
59 //                                    Muscle
60 //==============================================================================
61 /**
62  * A base class for modeling a muscle-tendon actuator. It defines muscle parameters
63  * and methods to PathActuator, but does not implement all of the necessary methods,
64  * and remains an abstract class. The path information for a muscle is contained
65  * in PathActuator, and the force-generating behavior should be defined in
66  * the derived classes.
67  *
68  * This class defines a subset of muscle models that include an active fiber
69  * (contractile element) in series with a tendon. This class defines common
70  * data members and handles the geometry of a unipennate fiber in connection
71  * with a tendon. No states are assumed, but concrete classes are free to
72  * add whatever states are necessary to describe the specific behavior of a
73  * muscle.
74  *
75  * @author Ajay Seth, Matt Millard
76  *
77  * (Based on earlier work by Peter Loan and Frank C. Anderson.)
78  */
79 class OSIMSIMULATION_API Muscle : public PathActuator {
80 OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator);
81 public:
82 //=============================================================================
83 // PROPERTIES
84 //=============================================================================
85     OpenSim_DECLARE_PROPERTY(max_isometric_force, double,
86         "Maximum isometric force that the fibers can generate");
87     OpenSim_DECLARE_PROPERTY(optimal_fiber_length, double,
88         "Optimal length of the muscle fibers");
89     OpenSim_DECLARE_PROPERTY(tendon_slack_length, double,
90         "Resting length of the tendon");
91     OpenSim_DECLARE_PROPERTY(pennation_angle_at_optimal, double,
92         "Angle between tendon and fibers at optimal fiber length expressed in radians");
93     OpenSim_DECLARE_PROPERTY(max_contraction_velocity, double,
94         "Maximum contraction velocity of the fibers, in optimal fiberlengths/second");
95     OpenSim_DECLARE_PROPERTY(ignore_tendon_compliance, bool,
96         "Compute muscle dynamics ignoring tendon compliance. Tendon is assumed to be rigid.");
97     OpenSim_DECLARE_PROPERTY(ignore_activation_dynamics, bool,
98         "Compute muscle dynamics ignoring activation dynamics. Activation is equivalent to excitation.");
99 
100 //=============================================================================
101 // OUTPUTS
102 //=============================================================================
103     OpenSim_DECLARE_OUTPUT(excitation, double, getExcitation,
104             SimTK::Stage::Dynamics);
105     OpenSim_DECLARE_OUTPUT(activation, double, getActivation,
106             SimTK::Stage::Dynamics);
107     OpenSim_DECLARE_OUTPUT(fiber_length, double, getFiberLength,
108             SimTK::Stage::Position);
109     OpenSim_DECLARE_OUTPUT(pennation_angle, double, getPennationAngle,
110             SimTK::Stage::Position);
111     OpenSim_DECLARE_OUTPUT(cos_pennation_angle, double, getCosPennationAngle,
112             SimTK::Stage::Position);
113     OpenSim_DECLARE_OUTPUT(tendon_length, double, getTendonLength,
114             SimTK::Stage::Position);
115     OpenSim_DECLARE_OUTPUT(normalized_fiber_length, double,
116             getNormalizedFiberLength, SimTK::Stage::Position);
117     OpenSim_DECLARE_OUTPUT(fiber_length_along_tendon, double,
118             getFiberLengthAlongTendon, SimTK::Stage::Position);
119     OpenSim_DECLARE_OUTPUT(tendon_strain, double, getTendonStrain,
120             SimTK::Stage::Position);
121     OpenSim_DECLARE_OUTPUT(passive_force_multiplier, double,
122             getPassiveForceMultiplier, SimTK::Stage::Position);
123     OpenSim_DECLARE_OUTPUT(active_force_length_multiplier, double,
124             getActiveForceLengthMultiplier, SimTK::Stage::Position);
125     OpenSim_DECLARE_OUTPUT(fiber_velocity, double, getFiberVelocity,
126             SimTK::Stage::Velocity);
127     OpenSim_DECLARE_OUTPUT(normalized_fiber_velocity, double,
128             getNormalizedFiberVelocity, SimTK::Stage::Velocity);
129     OpenSim_DECLARE_OUTPUT(fiber_velocity_along_tendon, double,
130             getFiberVelocityAlongTendon, SimTK::Stage::Velocity);
131     OpenSim_DECLARE_OUTPUT(tendon_velocity, double, getTendonVelocity,
132             SimTK::Stage::Velocity);
133     OpenSim_DECLARE_OUTPUT(force_velocity_multiplier, double,
134             getForceVelocityMultiplier, SimTK::Stage::Velocity);
135     OpenSim_DECLARE_OUTPUT(pennation_angular_velocity, double,
136             getPennationAngularVelocity, SimTK::Stage::Velocity);
137     OpenSim_DECLARE_OUTPUT(fiber_force, double, getFiberForce,
138             SimTK::Stage::Dynamics);
139     OpenSim_DECLARE_OUTPUT(fiber_force_along_tendon, double,
140             getFiberForceAlongTendon, SimTK::Stage::Dynamics);
141     OpenSim_DECLARE_OUTPUT(active_fiber_force, double, getActiveFiberForce,
142             SimTK::Stage::Dynamics);
143     OpenSim_DECLARE_OUTPUT(passive_fiber_force, double, getPassiveFiberForce,
144             SimTK::Stage::Dynamics);
145     OpenSim_DECLARE_OUTPUT(active_fiber_force_along_tendon, double,
146             getActiveFiberForceAlongTendon, SimTK::Stage::Dynamics);
147     OpenSim_DECLARE_OUTPUT(passive_fiber_force_along_tendon, double,
148             getPassiveFiberForceAlongTendon, SimTK::Stage::Dynamics);
149     OpenSim_DECLARE_OUTPUT(tendon_force, double, getTendonForce,
150             SimTK::Stage::Dynamics);
151     OpenSim_DECLARE_OUTPUT(fiber_stiffness, double, getFiberStiffness,
152             SimTK::Stage::Dynamics);
153     OpenSim_DECLARE_OUTPUT(fiber_stiffness_along_tendon, double,
154             getFiberStiffnessAlongTendon, SimTK::Stage::Dynamics);
155     OpenSim_DECLARE_OUTPUT(tendon_stiffness, double, getTendonStiffness,
156             SimTK::Stage::Dynamics);
157     OpenSim_DECLARE_OUTPUT(muscle_stiffness, double, getMuscleStiffness,
158             SimTK::Stage::Dynamics);
159     OpenSim_DECLARE_OUTPUT(fiber_active_power, double, getFiberActivePower,
160             SimTK::Stage::Dynamics);
161     OpenSim_DECLARE_OUTPUT(fiber_passive_power, double, getFiberPassivePower,
162             SimTK::Stage::Dynamics);
163     OpenSim_DECLARE_OUTPUT(tendon_power, double, getTendonPower,
164             SimTK::Stage::Dynamics);
165     OpenSim_DECLARE_OUTPUT(muscle_power, double, getMusclePower,
166             SimTK::Stage::Dynamics);
167 
168 //=============================================================================
169 // PUBLIC METHODS
170 //=============================================================================
171     /** @name Constructors and Destructor
172      */
173     //@{
174     /** Default constructor. */
175     Muscle();
176 
177     // default destructor, copy constructor and copy assignment
178 
179     //--------------------------------------------------------------------------
180     // MUSCLE PARAMETER ACCESSORS
181     //--------------------------------------------------------------------------
182     /** @name Muscle Parameters Access Methods
183      */
184     //@{
185     /** get/set the maximum isometric force (in N) that the fibers can generate */
186     double getMaxIsometricForce() const;
187     void setMaxIsometricForce(double maxIsometricForce);
188 
189     /** get/set the optimal length (in m) of the muscle fibers (lumped as a single fiber) */
190     double getOptimalFiberLength() const;
191     void setOptimalFiberLength(double optimalFiberLength);
192 
193     /** get/set the resting (slack) length (in m) of the tendon that is in series with the muscle fiber */
194     double getTendonSlackLength() const;
195     void setTendonSlackLength(double tendonSlackLength);
196 
197     /** get/set the angle (in radians) between fibers at their optimal fiber length and the tendon */
198     double getPennationAngleAtOptimalFiberLength() const;
199     void setPennationAngleAtOptimalFiberLength(double pennationAngle);
200 
201     /** get/set the maximum contraction velocity of the fibers, in optimal fiber-lengths per second */
202     double getMaxContractionVelocity() const;
203     void setMaxContractionVelocity(double maxContractionVelocity);
204 
205     // End of Muscle Parameter Accessors.
206     //@}
207 
208     //--------------------------------------------------------------------------
209     // MUSCLE STATE DEPENDENT ACCESSORS
210     //--------------------------------------------------------------------------
211     /** @name Muscle State Dependent Access Methods
212      *  Get quantities of interest common to all muscles
213      */
214     //@{
215 
216     /** Get/set Modeling (runtime) option to ignore tendon compliance when
217     computing muscle dynamics. This does not directly modify the persistent
218     property value. **/
219     bool getIgnoreTendonCompliance(const SimTK::State& s) const;
220     void setIgnoreTendonCompliance(SimTK::State& s, bool ignore) const;
221 
222     /** Get/set Modeling (runtime) option to ignore activation dynamics when
223     computing muscle dynamics. This does not directly modify the persistent
224     property value. **/
225     bool getIgnoreActivationDynamics(const SimTK::State& s) const;
226     void setIgnoreActivationDynamics(SimTK::State& s, bool ignore) const;
227 
228     /** get the activation level of the muscle, which modulates the active force
229         of the muscle and has a normalized (0 to 1) value
230         Note: method remains virtual to permit override by deprecated muscles. */
231     virtual double getActivation(const SimTK::State& s) const;
232 
233     /** get the current working fiber length (m) for the muscle */
234     double getFiberLength(const SimTK::State& s) const;
235     /** get the current pennation angle (radians) between the fiber and tendon at the current fiber length  */
236     double getPennationAngle(const SimTK::State& s) const;
237     /** get the cosine of the current pennation angle (radians) between the fiber and tendon at the current fiber length  */
238     double getCosPennationAngle(const SimTK::State& s) const;
239     /** get the current tendon length (m)  given the current joint angles and fiber length */
240     double getTendonLength(const SimTK::State& s) const;
241     /** get the current normalized fiber length (fiber_length/optimal_fiber_length) */
242     double getNormalizedFiberLength(const SimTK::State& s) const;
243     /** get the current fiber length (m) projected (*cos(pennationAngle)) onto the tendon direction */
244     double getFiberLengthAlongTendon(const SimTK::State& s) const;
245     /** get the current tendon strain (delta_l/tendon_slack_length is dimensionless)  */
246     double getTendonStrain(const SimTK::State& s) const;
247 
248     /** the potential energy (J) stored in the fiber due to its parallel elastic element */
249     double getFiberPotentialEnergy(const SimTK::State& s) const;
250     /** the potential energy (J) stored in the tendon */
251     double getTendonPotentialEnergy(const SimTK::State& s) const;
252     /** the total potential energy (J) stored in the muscle */
253     double getMusclePotentialEnergy(const SimTK::State& s) const;
254 
255     /** get the passive fiber (parallel elastic element) force multiplier */
256     double getPassiveForceMultiplier(const SimTK::State& s) const;
257     /** get the active fiber (contractile element) force multiplier due to current fiber length */
258     double getActiveForceLengthMultiplier(const SimTK::State& s) const;
259 
260     /** get current fiber velocity (m/s) positive is lengthening */
261     double getFiberVelocity(const SimTK::State& s) const;
262     /** get normalized fiber velocity. This is the fiber velocity in m/s divided by
263     the maximum contraction velocity expressed in m/s; therefore, this quantity is
264     dimensionless and generally lies in the range [-1, 1]. */
265     double getNormalizedFiberVelocity(const SimTK::State& s) const;
266     /** get the current fiber velocity (m/s) projected onto the tendon direction */
267     double getFiberVelocityAlongTendon(const SimTK::State& s) const;
268     /** get pennation angular velocity (radians/s) */
269     double getPennationAngularVelocity(const SimTK::State& s) const;
270     /** get the tendon velocity (m/s) positive is lengthening */
271     double getTendonVelocity(const SimTK::State& s) const;
272     /** get the dimensionless multiplier resulting from the fiber's force-velocity curve */
273     double getForceVelocityMultiplier(const SimTK::State& s) const;
274 
275     /** get the current fiber force (N) applied to the tendon */
276     double getFiberForce(const SimTK::State& s) const;
277     /**get the force of the fiber (N/m) along the direction of the tendon*/
278     double getFiberForceAlongTendon(const SimTK::State& s) const;
279     /** get the current active fiber force (N) due to activation*force_length*force_velocity relationships */
280     double getActiveFiberForce(const SimTK::State& s) const;
281     /** get the total force applied by all passive elements in the fiber (N) */
282     double getPassiveFiberForce(const SimTK::State& s) const;
283     /** get the current active fiber force (N) projected onto the tendon direction */
284     double getActiveFiberForceAlongTendon(const SimTK::State& s) const;
285     /** get the total force applied by all passive elements in the fiber (N)
286         projected onto the tendon direction */
287     double getPassiveFiberForceAlongTendon(const SimTK::State& s) const;
288     /** get the current tendon force (N) applied to bones */
289     double getTendonForce(const SimTK::State& s) const;
290 
291     /** get the current fiber stiffness (N/m) defined as the partial derivative
292         of fiber force with respect to fiber length */
293     double getFiberStiffness(const SimTK::State& s) const;
294     /**get the stiffness of the fiber (N/m) along the direction of the tendon,
295     that is the partial derivative of the fiber force along the tendon with
296     respect to small changes in fiber length along the tendon*/
297     double getFiberStiffnessAlongTendon(const SimTK::State& s) const;
298     /** get the current tendon stiffness (N/m) defined as the partial derivative
299         of tendon force with respect to tendon length */
300     double getTendonStiffness(const SimTK::State& s) const;
301     /** get the current muscle stiffness (N/m) defined as the partial derivative
302         of muscle force with respect to muscle length */
303     double getMuscleStiffness(const SimTK::State& s) const;
304 
305     /** get the current active fiber power (W) */
306     double getFiberActivePower(const SimTK::State& s) const;
307     /** get the current passive fiber power (W) */
308     double getFiberPassivePower(const SimTK::State& s) const;
309 
310     /** get the current tendon power (W) */
311     double getTendonPower(const SimTK::State& s) const;
312     /** get the current muscle power (W) */
313     double getMusclePower(const SimTK::State& s) const;
314 
315     /** get the stress in the muscle (part of the Actuator interface as well) */
316     double getStress(const SimTK::State& s) const override;
317 
318     /** set the excitation (control) for this muscle. NOTE if controllers are connected to the
319         muscle and are adding in their controls, and setExcitation is called after the model's
320         computeControls(), then setExcitation will override the controller values. If called
321         before computeControls, then controller value(s) are added to the excitation set here. */
322     void setExcitation(SimTK::State& s, double excitation) const;
323     double getExcitation(const SimTK::State& s) const;
324 
325 
326     /** DEPRECATED: only for backward compatibility */
327     virtual void setActivation(SimTK::State& s, double activation) const = 0;
328 
329     // End of Muscle's State Dependent Accessors.
330     //@}
331 
332     /** Actuator interface for a muscle computes the tension in the muscle
333         and applied by the tendon to bones (i.e. not the fiber force) */
334     double computeActuation(const SimTK::State& s) const override = 0;
335 
336     /** @name Muscle initialization
337      */
338     //@{
339     /** Find and set the equilibrium state of the muscle (if any) */
computeEquilibrium(SimTK::State & s)340     void computeEquilibrium(SimTK::State& s) const override final {
341         return computeInitialFiberEquilibrium(s);
342     }
343     // End of Muscle's State Dependent Accessors.
344     //@}
345 
346     ///@cond
347     //--------------------------------------------------------------------------
348     // Estimate the muscle force for a given activation based on a rigid tendon
349     // assumption and neglecting passive fiber force. This provides a linear
350     // relationship between activation and force. This is used by CMC and
351     // StaticOptimization to solve the muscle force redundancy problem.
352     //--------------------------------------------------------------------------
353     virtual double calcInextensibleTendonActiveFiberForce(SimTK::State& s,
354                                                   double aActivation) const;
355     ///@endcond
356 //=============================================================================
357 // PROTECTED METHODS
358 //=============================================================================
359 protected:
360     struct MuscleLengthInfo;
361     struct FiberVelocityInfo;
362     struct MuscleDynamicsInfo;
363     struct MusclePotentialEnergyInfo;
364 
365     /** Developer Access to intermediate values calculate by the muscle model */
366     const MuscleLengthInfo& getMuscleLengthInfo(const SimTK::State& s) const;
367     MuscleLengthInfo& updMuscleLengthInfo(const SimTK::State& s) const;
368 
369     const FiberVelocityInfo& getFiberVelocityInfo(const SimTK::State& s) const;
370     FiberVelocityInfo& updFiberVelocityInfo(const SimTK::State& s) const;
371 
372     const MuscleDynamicsInfo& getMuscleDynamicsInfo(const SimTK::State& s) const;
373     MuscleDynamicsInfo& updMuscleDynamicsInfo(const SimTK::State& s) const;
374 
375     const MusclePotentialEnergyInfo& getMusclePotentialEnergyInfo(const SimTK::State& s) const;
376     MusclePotentialEnergyInfo& updMusclePotentialEnergyInfo(const SimTK::State& s) const;
377 
378     //--------------------------------------------------------------------------
379     // CALCULATIONS
380     //--------------------------------------------------------------------------
381     /** @name Muscle State Dependent Calculations
382      *  Developers must override these methods to implement the desired behavior
383      *  of their muscle models. Unless you are augmenting the behavior
384      *  of an existing muscle class or writing a new derived class, you do not
385      *  have access to these methods.
386      */
387     //@{
388     /** calculate muscle's position related values such fiber and tendon lengths,
389         normalized lengths, pennation angle, etc... */
390     virtual void calcMuscleLengthInfo(const SimTK::State& s,
391         MuscleLengthInfo& mli) const;
392 
393     /** calculate muscle's fiber velocity and pennation angular velocity, etc... */
394     virtual void calcFiberVelocityInfo(const SimTK::State& s,
395         FiberVelocityInfo& fvi) const;
396 
397     /** calculate muscle's active and passive force-length, force-velocity,
398         tendon force, relationships and their related values */
399     virtual void  calcMuscleDynamicsInfo(const SimTK::State& s,
400         MuscleDynamicsInfo& mdi) const;
401 
402     /** calculate muscle's fiber and tendon potential energy */
403     virtual void calcMusclePotentialEnergyInfo(const SimTK::State& s,
404         MusclePotentialEnergyInfo& mpei) const;
405 
406     /** This function modifies the fiber length in the supplied state such that
407     the fiber and tendon are developing the same force, taking activation and
408     velocity into account. This routine can assume that the state contains a
409     meaningful estimate of muscle activation, joint positions, and joint
410     velocities. For example, this can produce fiber lengths suited to
411     beginning a forward dynamics simulation.
412     computeFiberEquilibriumAtZeroVelocity(). */
413     virtual void computeInitialFiberEquilibrium(SimTK::State& s) const = 0;
414 
415     // End of Muscle's State Related Calculations.
416     //@}
417 
418     //--------------------------------------------------------------------------
419     // PARENT INTERFACES
420     //--------------------------------------------------------------------------
421     /** @name Interfaces imposed by parent classes
422      */
423     //@{
424 
425     /** Force interface applies tension to bodies, and Muscle also checks that
426         applied muscle tension is not negative. */
427     void computeForce(const SimTK::State& state,
428                       SimTK::Vector_<SimTK::SpatialVec>& bodyForces,
429                       SimTK::Vector& generalizedForce) const override;
430 
431     /** Potential energy stored by the muscle */
432     double computePotentialEnergy(const SimTK::State& state) const override;
433 
434     /** Override PathActuator virtual to calculate a preferred color for the
435     muscle path based on activation. **/
436     SimTK::Vec3 computePathColor(const SimTK::State& state) const override;
437 
438     /** Model Component creation interface */
439     void extendConnectToModel(Model& aModel) override;
440     void extendAddToSystem(SimTK::MultibodySystem& system) const override;
441     void extendSetPropertiesFromState(const SimTK::State &s) override;
442     void extendInitStateFromProperties(SimTK::State& state) const override;
443 
444     // Update the display geometry attached to the muscle
445     virtual void updateGeometry(const SimTK::State& s);
446     // End of Interfaces imposed by parent classes.
447     //@}
448 
449 
450 private:
451     void setNull();
452     void constructProperties();
453     void copyData(const Muscle &aMuscle);
454 
455     //--------------------------------------------------------------------------
456     // Implement Object interface.
457     //--------------------------------------------------------------------------
458     /** Override of the default implementation to account for versioning. */
459     void updateFromXMLNode(SimTK::Xml::Element& aNode, int versionNumber=-1) override;
460 
461 
462 //=============================================================================
463 // DATA
464 //=============================================================================
465 protected:
466 
467     /** The assumed fixed muscle-width from which the fiber pennation angle can
468         be calculated. */
469     double _muscleWidth;
470 
471  /**
472     The MuscleLengthInfo struct contains information about the muscle that is
473     strictly a function of the length of the fiber and the tendon, and the
474     orientation of the muscle fiber. w.r.t. a fixed orientation of the tendon.
475 
476     The function that populates this struct, calcMuscleLengthInfo, is
477     called at a point when only the position and orientation of the system are
478     known. This function is the first one that is called of the functions
479     calcMuscleLengthInfo, calcFiberVelocityInfo and calcMuscleDynamicInfo.
480     The velocity and acceleration of the muscle's path will not be known when
481     this function is called.
482 
483             NAME                    DIMENSION         UNITS
484              fiberLength              length            m
485              fiberLengthAlongTendon   length            m           [1]
486              normFiberLength          length/length     m/m         [2]
487 
488              tendonLength             length            m
489              normTendonLength         length/length     m/m         [3]
490              tendonStrain             length/length     m/m         [4]
491 
492              pennationAngle           angle             rad         [5]
493              cosPennationAngle        NA                NA
494              sinPennationAngle        NA                NA
495 
496              fiberPassiveForceLengthMultiplier   force/force     N/N      [6]
497              fiberActiveForceLengthMultiplier    force/force     N/N      [7]
498 
499             userDefinedLengthExtras     NA              NA            [8]
500 
501     [1] fiberLengthAlongTendon is the length of the muscle fiber as projected
502         on the tendon.
503 
504     [2] normFiberLength is the fiberLength normalized with respect to the
505         optimalFiberLength,
506 
507         normFiberLength = fiberLength / optimalFiberLength
508 
509         N.B. It is assumed that the optimalFiberLength of a muscle is also
510         its resting length.
511 
512     [3] normTendonLength is the tendonLength normalized with respect to the
513         tendonSlackLength
514 
515         normTendonLength = tendonLength / tendonSlackLength
516 
517     [4] Tendon strain is defined using the elongation of the material divided by
518         its resting length. This is identical to the engineering definition of
519         strain. Thus a tendonStrain of 0.01 means that the tendon is currently
520         1% longer than its resting length.
521 
522         tendonStrain = (tendonLength-tendonSlackLength)/tendonSlackLength
523 
524 
525     [5] The orientation of the muscle fiber is defined w.r.t. a fixed
526         orientation of the tendon. A pennation angle of 0 means that the fiber
527         is collinear to the tendon. It is normal for the pennation angle
528         to range between 0 and Pi/2 radians.
529 
530               Fiber                 Tendon
531         |===================||-----------------|   Pennation = 0
532 
533              ||-------------------|                Pennation = SimTK::Pi/3
534            //                                                   (60 degrees)
535           //
536          //
537         //
538 
539     [6] The fiberPassiveForceLengthMultiplier represents the elastic force the fiber
540         generates normalized w.r.t. the maximum isometric force of the fiber.
541         Is typically specified by a passiveForceLengthCurve.
542 
543 
544     [7] The fiberActiveForceLengthMultiplier is the scaling of the maximum force a fiber
545         can generate as a function of its length. This term usually follows a
546         curve that is zero at a normalized fiber length of 0.5, is 1 at a
547         normalized fiber length of 1, and then zero again at a normalized fiber
548         length of 1.5. This curve is generally an interpolation of experimental
549         data.
550 
551     [8] This vector is left for the muscle modeler to populate with any
552         computationally expensive quantities that are computed in
553         calcMuscleLengthInfo, and required for use in the user defined functions
554         calcFiberVelocityInfo and calcMuscleDynamicsInfo. None of the parent
555         classes make any assumptions about what is or isn't in this field
556         - use as necessary.
557 
558     */
559     struct MuscleLengthInfo{             //DIMENSION         Units
560         double fiberLength;              //length            m
561         double fiberLengthAlongTendon;   //length            m
562         double normFiberLength;          //length/length     m/m
563 
564         double tendonLength;             //length            m
565         double normTendonLength;         //length/length     m/m
566         double tendonStrain;             //length/length     m/m
567                                          //
568         double pennationAngle;           //angle             1/s (rads)
569         double cosPennationAngle;        //NA                NA
570         double sinPennationAngle;        //NA                NA
571 
572         double fiberPassiveForceLengthMultiplier;   //NA             NA
573         double fiberActiveForceLengthMultiplier;  //NA             NA
574 
575         SimTK::Vector userDefinedLengthExtras;//NA        NA
576 
MuscleLengthInfoMuscleLengthInfo577         MuscleLengthInfo():
578             fiberLength(SimTK::NaN),
579             fiberLengthAlongTendon(SimTK::NaN),
580             normFiberLength(SimTK::NaN),
581             tendonLength(SimTK::NaN),
582             normTendonLength(SimTK::NaN),
583             tendonStrain(SimTK::NaN),
584             pennationAngle(SimTK::NaN),
585             cosPennationAngle(SimTK::NaN),
586             sinPennationAngle(SimTK::NaN),
587             fiberPassiveForceLengthMultiplier(SimTK::NaN),
588             fiberActiveForceLengthMultiplier(SimTK::NaN),
589             userDefinedLengthExtras(0, SimTK::NaN){}
590         friend std::ostream& operator<<(std::ostream& o,
591             const MuscleLengthInfo& mli) {
592             o << "Muscle::MuscleLengthInfo should not be serialized!"
593               << std::endl;
594             return o;
595         }
596     };
597 
598     /**
599         FiberVelocityInfo contains velocity quantities related to the velocity
600         of the muscle (fiber + tendon) complex.
601 
602         The function that populates this struct, calcFiberVelocityInfo, is called
603         when position and velocity information is known. This function is the
604         second function that is called of these related
605         functions:calcMuscleLengthInfo,calcFiberVelocityInfo and
606         calcMuscleDynamicInfo. When calcFiberVelocity is called the acceleration
607         of the muscle path, and any forces the muscle experiences will not be
608         known.
609 
610             NAME                     DIMENSION             UNITS
611              fiberVelocity             length/time           m/s
612              fiberVelocityAlongTendon  length/time           m/s       [1]
613              normFiberVelocity         (length/time)/Vmax    NA        [2]
614 
615              pennationAngularVelocity  angle/time            rad/s     [3]
616 
617              tendonVelocity            length/time           m/s
618              normTendonVelocity        (length/time)/length  (m/s)/m   [4]
619 
620              fiberForceVelocityMultiplier force/force          NA        [5]
621 
622              userDefinedVelocityExtras    NA                   NA      [6]
623 
624         [1] fiberVelocityAlongTendon is the first derivative of the symbolic
625             equation that defines the fiberLengthAlongTendon.
626 
627         [2] normFiberVelocity is the fiberVelocity (in m/s) divided by
628             the optimal length of the fiber (in m) and by the maximum fiber
629             velocity (in optimal-fiber-lengths/s). normFiberVelocity has
630             units of 1/optimal-fiber-length.
631 
632         [3] The sign of the angular velocity is defined using the right
633             hand rule.
634 
635         [4] normTendonVelocity is the tendonVelocity (the lengthening velocity
636             of the tendon) divided by its resting length
637 
638         [5] The fiberForceVelocityMultiplier is the scaling factor that represents
639             how a muscle fiber's force generating capacity is modulated by the
640             contraction (concentric or eccentric) velocity of the fiber.
641             Generally this curve has a value of 1 at a fiber velocity of 0,
642             has a value of between 1.4-1.8 at the maximum eccentric contraction
643             velocity and a value of 0 at the maximum concentric contraction
644             velocity. The force velocity curve, which computes this term,
645             is usually an interpolation of an experimental curve.
646 
647         [6] This vector is left for the muscle modeler to populate with any
648             computationally expensive quantities that are computed in
649             calcFiberVelocityInfo, and required for use in the user defined
650             function calcMuscleDynamicsInfo. None of the parent classes make
651             any assumptions about what is or isn't in this field
652             - use as necessary.
653 
654     */
655     struct FiberVelocityInfo {              //DIMENSION             UNITS
656         double fiberVelocity;               //length/time           m/s
657         double fiberVelocityAlongTendon;    //length/time           m/s
658         double normFiberVelocity;           //(length/time)/Vmax    NA
659                                             //
660         double pennationAngularVelocity;    //angle/time            rad/s
661                                             //
662         double tendonVelocity;              //length/time           m/s
663         double normTendonVelocity;          //(length/time)/length  (m/s)/m
664 
665         double fiberForceVelocityMultiplier;     //force/force           NA
666 
667         SimTK::Vector userDefinedVelocityExtras;//NA                  NA
668 
FiberVelocityInfoFiberVelocityInfo669         FiberVelocityInfo():
670             fiberVelocity(SimTK::NaN),
671             fiberVelocityAlongTendon(SimTK::NaN),
672             normFiberVelocity(SimTK::NaN),
673             pennationAngularVelocity(SimTK::NaN),
674             tendonVelocity(SimTK::NaN),
675             normTendonVelocity(SimTK::NaN),
676             fiberForceVelocityMultiplier(SimTK::NaN),
677             userDefinedVelocityExtras(0,SimTK::NaN){};
678         friend std::ostream& operator<<(std::ostream& o,
679             const FiberVelocityInfo& fvi) {
680             o << "Muscle::FiberVelocityInfo should not be serialized!"
681               << std::endl;
682             return o;
683         }
684     };
685 
686     /**
687         MuscleDynamicsInfo contains quantities that are related to the forces
688         that the muscle generates.
689 
690         The function that populates this struct, calcMuscleDynamicsInfo, is
691         called when position and velocity information is known. This function
692         is the last function that is called of these related functions:
693         calcMuscleLengthInfo, calcFiberVelocityInfo and calcMuscleDynamicInfo.
694 
695 
696            NAME                         DIMENSION           UNITS
697             activation                  NA                   NA     [1]
698 
699             fiberForce                  force                N
700             fiberForceAlongTendon       force                N      [2]
701             normFiberForce              force/force          N/N    [3]
702             activeFiberForce            force                N      [4]
703             passiveFiberForce           force                N      [5]
704 
705             tendonForce                 force                N
706             normTendonForce             force/force          N/N    [6]
707 
708             fiberStiffness              force/length         N/m    [7]
709             fiberStiffnessAlongTendon   force/length         N/m    [8]
710             tendonStiffness             force/length         N/m    [9]
711             muscleStiffness             force/length         N/m    [10]
712 
713             fiberActivePower            force*velocity       W (N*m/s)
714             fiberPassivePower           force*velocity       W (N*m/s)
715             tendonPower                 force*velocity       W (N*m/s)
716             musclePower                 force*velocity       W (N*m/s)
717 
718             userDefinedDynamicsData     NA                   NA     [11]
719 
720         [1] This is a quantity that ranges between 0 and 1 that dictates how
721             on or activated a muscle is. This term may or may not have its own
722             time dependent behavior depending on the muscle model.
723 
724         [2] fiberForceAlongTendon is the fraction of the force that is developed
725             by the fiber that is transmitted to the tendon. This fraction
726             depends on the pennation model that is used for the muscle model
727 
728         [3] This is the force developed by the fiber scaled by the maximum
729             isometric contraction force. Note that the maximum isometric force
730             is defined as the maximum isometric force a muscle fiber develops
731             at its optimal pennation angle, and along the line of the fiber.
732 
733         [4] This is the portion of the fiber force that is created as a direct
734             consequence of the value of 'activation'.
735 
736         [5] This is the portion of the fiber force that is created by the
737             parallel elastic element within the fiber.
738 
739         [6] This is the tendonForce normalized by the maximum isometric
740             contraction force
741 
742         [7] fiberStiffness is defined as the partial derivative of fiber force
743             with respect to fiber length
744 
745         [8] fiberStiffnessAlongTendon is defined as the partial derivative of
746             fiber force along the tendon with respect to small changes in
747             the fiber length along the tendon. This quantity is normally
748             computed using the equations for fiberStiffness, and then using an
749             application of the chain rule to yield fiberStiffnessAlongTendon.
750 
751         [9] tendonStiffness is defined as the partial derivative of tendon
752             force with respect to tendon length
753 
754         [10] muscleStiffness is defined as the partial derivative of muscle force
755             with respect to changes in muscle length. This quantity can usually
756             be computed by noting that the tendon and the fiber are in series,
757             with the fiber at a pennation angle. Thus
758 
759             Kmuscle =   (Kfiber_along_tendon * Ktendon)
760                        /(Kfiber_along_tendon + Ktendon)
761 
762         [11] This vector is left for the muscle modeler to populate with any
763              computationally expensive quantities that might be of interest
764              after dynamics calculations are completed but maybe of use
765              in computing muscle derivatives or reporting values of interest.
766 
767     */
768     struct MuscleDynamicsInfo {     //DIMENSION             UNITS
769         double activation;              // NA                   NA
770                                         //
771         double fiberForce;              // force                N
772         double fiberForceAlongTendon;   // force                N
773         double normFiberForce;          // force/force          N/N
774         double activeFiberForce;        // force                N
775         double passiveFiberForce;       // force                N
776                                         //
777         double tendonForce;             // force                N
778         double normTendonForce;         // force/force          N/N
779                                         //
780         double fiberStiffness;          // force/length         N/m
781         double fiberStiffnessAlongTendon;//force/length         N/m
782         double tendonStiffness;         // force/length         N/m
783         double muscleStiffness;         // force/length         N/m
784                                         //
785         double fiberActivePower;        // force*velocity       W
786         double fiberPassivePower;       // force*velocity       W
787         double tendonPower;             // force*velocity       W
788         double musclePower;             // force*velocity       W
789 
790         SimTK::Vector userDefinedDynamicsExtras; //NA          NA
791 
MuscleDynamicsInfoMuscleDynamicsInfo792         MuscleDynamicsInfo():
793             activation(SimTK::NaN),
794             fiberForce(SimTK::NaN),
795             fiberForceAlongTendon(SimTK::NaN),
796             normFiberForce(SimTK::NaN),
797             activeFiberForce(SimTK::NaN),
798             passiveFiberForce(SimTK::NaN),
799             tendonForce(SimTK::NaN),
800             normTendonForce(SimTK::NaN),
801             fiberStiffness(SimTK::NaN),
802             fiberStiffnessAlongTendon(SimTK::NaN),
803             tendonStiffness(SimTK::NaN),
804             muscleStiffness(SimTK::NaN),
805             fiberActivePower(SimTK::NaN),
806             fiberPassivePower(SimTK::NaN),
807             tendonPower(SimTK::NaN),
808             musclePower(SimTK::NaN),
809             userDefinedDynamicsExtras(0, SimTK::NaN){};
810         friend std::ostream& operator<<(std::ostream& o,
811             const MuscleDynamicsInfo& mdi) {
812             o << "Muscle::MuscleDynamicsInfo should not be serialized!"
813               << std::endl;
814             return o;
815         }
816     };
817 
818     /**
819         MusclePotentialEnergyInfo contains quantities related to the potential
820         energy of the muscle (fiber + tendon) complex.
821 
822         The function that populates this struct, calcMusclePotentialEnergyInfo, can
823         be called when position information is known. This function is
824         dependent on calcMuscleLengthInfo.
825 
826         NAME                     DIMENSION              UNITS
827         fiberPotentialEnergy      force*distance         J (Nm)   [1]
828         tendonPotentialEnergy     force*distance         J (Nm)   [2]
829         musclePotentialEnergy     force*distance         J (Nm)   [3]
830 
831         userDefinedPotentialEnergyExtras                         [4]
832 
833         [4] This vector is left for the muscle modeler to populate with any
834             computationally expensive quantities that are computed in
835             calcMusclePotentialEnergyInfo, that might be useful for others to
836             access.
837 
838     */
839     struct MusclePotentialEnergyInfo {              //DIMENSION             UNITS
840         double fiberPotentialEnergy;     //force*distance    J (Nm)
841         double tendonPotentialEnergy;    //force*distance    J (Nm)
842         double musclePotentialEnergy;    //force*distance    J (Nm)
843 
844         SimTK::Vector userDefinedPotentialEnergyExtras;//NA                  NA
845 
MusclePotentialEnergyInfoMusclePotentialEnergyInfo846         MusclePotentialEnergyInfo():
847             fiberPotentialEnergy(SimTK::NaN),
848             tendonPotentialEnergy(SimTK::NaN),
849             musclePotentialEnergy(SimTK::NaN),
850             userDefinedPotentialEnergyExtras(0,SimTK::NaN){};
851         friend std::ostream& operator<<(std::ostream& o,
852             const MusclePotentialEnergyInfo& fvi) {
853             o << "Muscle::MusclePotentialEnergyInfo should not be serialized!"
854               << std::endl;
855             return o;
856         }
857     };
858 
859 
860     /** to support deprecated muscles */
861     double _maxIsometricForce;
862     double _optimalFiberLength;
863     double _pennationAngleAtOptimal;
864     double _tendonSlackLength;
865 
866 //=============================================================================
867 };  // END of class Muscle
868 //=============================================================================
869 //=============================================================================
870 
871 } // end of namespace OpenSim
872 
873 #endif // OPENSIM_MUSCLE_H_
874