From 6965d1dd06e7a6b8728926967c78c0cfd5724960 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 13:44:53 +0200
Subject: [PATCH 01/23] Add empty FiberForceLengthActiveDeGroote2016 class

---
 sympy/physics/_biomechanics/characteristic.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/sympy/physics/_biomechanics/characteristic.py b/sympy/physics/_biomechanics/characteristic.py
index ba335161ed75..29e05784c90e 100644
--- a/sympy/physics/_biomechanics/characteristic.py
+++ b/sympy/physics/_biomechanics/characteristic.py
@@ -7,6 +7,7 @@
 
 
 __all__ = [
+    'FiberForceLengthActiveDeGroote2016',
     'FiberForceLengthPassiveDeGroote2016',
     'FiberForceLengthPassiveInverseDeGroote2016',
     'TendonForceLengthDeGroote2016',
@@ -797,3 +798,7 @@ def _latex(self, printer):
         fl_M_pas = self.args[0]
         _fl_M_pas = printer._print(fl_M_pas)
         return r'\left( \operatorname{fl}^M_{pas} \right)^{-1} \left( %s \right)' % _fl_M_pas
+
+
+class FiberForceLengthActiveDeGroote2016(CharacteristicCurveFunction):
+    pass

From 2c7d48e21ebfb4b591fac8337db975d404ab9956 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 13:45:31 +0200
Subject: [PATCH 02/23] Add test class for FiberForceLengthActiveDeGroote2016
 with fixture

---
 .../tests/test_characteristic.py              | 24 +++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index cd6d577c54fc..93f758a9d10d 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -10,6 +10,7 @@
 from sympy.functions.elementary.exponential import exp, log
 from sympy.physics._biomechanics.characteristic import (
     CharacteristicCurveFunction,
+    FiberForceLengthActiveDeGroote2016,
     FiberForceLengthPassiveDeGroote2016,
     FiberForceLengthPassiveInverseDeGroote2016,
     TendonForceLengthDeGroote2016,
@@ -788,3 +789,26 @@ def test_lambdify_jax(self):
             1.2774998934,
         ])
         numpy.testing.assert_allclose(fl_M_pas_inv_callable(fl_M_pas), expected)
+
+
+class TestFiberForceLengthActiveDeGroote2016:
+
+    @pytest.fixture(autouse=True)
+    def _fiber_force_length_active_arguments_fixture(self):
+        self.l_M_tilde = Symbol('l_M_tilde')
+        self.c0 = Symbol('c_0')
+        self.c1 = Symbol('c_1')
+        self.c2 = Symbol('c_2')
+        self.c3 = Symbol('c_3')
+        self.c4 = Symbol('c_4')
+        self.c5 = Symbol('c_5')
+        self.c6 = Symbol('c_6')
+        self.c7 = Symbol('c_7')
+        self.c8 = Symbol('c_8')
+        self.c9 = Symbol('c_9')
+        self.c10 = Symbol('c_10')
+        self.c11 = Symbol('c_11')
+        self.constants = (
+            self.c0, self.c1, self.c2, self.c3, self.c4, self.c5,
+            self.c6, self.c7, self.c8, self.c9, self.c10, self.c11,
+        )

From 556b716d825b854212998d0e9259b5c748f00b5e Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 13:45:51 +0200
Subject: [PATCH 03/23] Add test case for inheritance

---
 sympy/physics/_biomechanics/tests/test_characteristic.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index 93f758a9d10d..eafdb57b9498 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -812,3 +812,9 @@ def _fiber_force_length_active_arguments_fixture(self):
             self.c0, self.c1, self.c2, self.c3, self.c4, self.c5,
             self.c6, self.c7, self.c8, self.c9, self.c10, self.c11,
         )
+
+    @staticmethod
+    def test_class():
+        assert issubclass(FiberForceLengthActiveDeGroote2016, Function)
+        assert issubclass(FiberForceLengthActiveDeGroote2016, CharacteristicCurveFunction)
+        assert FiberForceLengthActiveDeGroote2016.__name__ == 'FiberForceLengthActiveDeGroote2016'

From 999774d7aab62537d1444c8514e06adee0344bb3 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 13:46:08 +0200
Subject: [PATCH 04/23] Add test case for instantiation

---
 sympy/physics/_biomechanics/tests/test_characteristic.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index eafdb57b9498..7bf6414adbcc 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -818,3 +818,11 @@ def test_class():
         assert issubclass(FiberForceLengthActiveDeGroote2016, Function)
         assert issubclass(FiberForceLengthActiveDeGroote2016, CharacteristicCurveFunction)
         assert FiberForceLengthActiveDeGroote2016.__name__ == 'FiberForceLengthActiveDeGroote2016'
+
+    def test_instance(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        assert isinstance(fl_M_act, FiberForceLengthActiveDeGroote2016)
+        assert str(fl_M_act) == (
+            'FiberForceLengthActiveDeGroote2016(l_M_tilde, c_0, c_1, c_2, c_3, '
+            'c_4, c_5, c_6, c_7, c_8, c_9, c_10, c_11)'
+        )

From 15d60258befcd158fefa9143b5187cf55b25e2e3 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 13:46:38 +0200
Subject: [PATCH 05/23] Add test cases for doit method

---
 .../_biomechanics/tests/test_characteristic.py   | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index 7bf6414adbcc..6c29b92575e6 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -826,3 +826,19 @@ def test_instance(self):
             'FiberForceLengthActiveDeGroote2016(l_M_tilde, c_0, c_1, c_2, c_3, '
             'c_4, c_5, c_6, c_7, c_8, c_9, c_10, c_11)'
         )
+
+    def test_doit(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants).doit()
+        assert fl_M_act == (
+            self.c0*exp(-0.5*((self.l_M_tilde - self.c1)/(self.c2 + self.c3*self.l_M_tilde))**2)
+            + self.c4*exp(-0.5*((self.l_M_tilde - self.c5)/(self.c6 + self.c7*self.l_M_tilde))**2)
+            + self.c8*exp(-0.5*((self.l_M_tilde - self.c9)/(self.c10 + self.c11*self.l_M_tilde))**2)
+        )
+
+    def test_doit_evaluate_false(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants).doit(evaluate=False)
+        assert fl_M_act == (
+            self.c0*exp(-0.5*(UnevaluatedExpr(self.l_M_tilde - self.c1)/(self.c2 + self.c3*self.l_M_tilde))**2)
+            + self.c4*exp(-0.5*(UnevaluatedExpr(self.l_M_tilde - self.c5)/(self.c6 + self.c7*self.l_M_tilde))**2)
+            + self.c8*exp(-0.5*(UnevaluatedExpr(self.l_M_tilde - self.c9)/(self.c10 + self.c11*self.l_M_tilde))**2)
+        )

From 1ee2ed9670f25d0df0b3b59496cb804fedfa2c40 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 13:47:10 +0200
Subject: [PATCH 06/23] Add test case for with_default_constants alternate
 constructor

---
 .../tests/test_characteristic.py              | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index 6c29b92575e6..6a7dcd7a8646 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -842,3 +842,22 @@ def test_doit_evaluate_false(self):
             + self.c4*exp(-0.5*(UnevaluatedExpr(self.l_M_tilde - self.c5)/(self.c6 + self.c7*self.l_M_tilde))**2)
             + self.c8*exp(-0.5*(UnevaluatedExpr(self.l_M_tilde - self.c9)/(self.c10 + self.c11*self.l_M_tilde))**2)
         )
+
+    def test_with_default_constants(self):
+        constants = (
+            Float('0.814'),
+            Float('1.06'),
+            Float('0.162'),
+            Float('0.0633'),
+            Float('0.433'),
+            Float('0.717'),
+            Float('-0.0299'),
+            Rational(1, 5),
+            Rational(1, 10),
+            Integer(1),
+            Float('0.354'),
+            Integer(0),
+        )
+        fl_M_act_manual = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *constants)
+        fl_M_act_constants = FiberForceLengthActiveDeGroote2016.with_default_constants(self.l_M_tilde)
+        assert fl_M_act_manual == fl_M_act_constants

From 839532cbb5954ccdb3a5bd4ec47c6b98eeefd7ec Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 14:52:30 +0200
Subject: [PATCH 07/23] Add test cases for differentiation

---
 .../tests/test_characteristic.py              | 105 ++++++++++++++++++
 1 file changed, 105 insertions(+)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index 6a7dcd7a8646..cef220f2b6da 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -861,3 +861,108 @@ def test_with_default_constants(self):
         fl_M_act_manual = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *constants)
         fl_M_act_constants = FiberForceLengthActiveDeGroote2016.with_default_constants(self.l_M_tilde)
         assert fl_M_act_manual == fl_M_act_constants
+
+    def test_differentiate_wrt_l_M_tilde(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            self.c0*(
+                self.c3*(self.l_M_tilde - self.c1)**2/(self.c2 + self.c3*self.l_M_tilde)**3
+                + (self.l_M_tilde - self.c1)/((self.c2 + self.c3*self.l_M_tilde)**2)
+            )*exp((self.c1 - self.l_M_tilde)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
+            + self.c4*(
+                self.c7*(self.l_M_tilde - self.c5)**2/(self.c6 + self.c7*self.l_M_tilde)**3
+                + (self.l_M_tilde - self.c5)/((self.c6 + self.c7*self.l_M_tilde)**2)
+            )*exp((self.c5 - self.l_M_tilde)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
+            + self.c8*(
+                self.c11*(self.l_M_tilde - self.c9)**2/(self.c10 + self.c11*self.l_M_tilde)**3
+                + (self.l_M_tilde - self.c9)/((self.c10 + self.c11*self.l_M_tilde)**2)
+            )*exp((self.c9 - self.l_M_tilde)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
+        )
+        assert fl_M_act.diff(self.l_M_tilde) == expected
+
+    def test_differentiate_wrt_c0(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = exp((self.c1 - self.l_M_tilde)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
+        assert fl_M_act.diff(self.c0) == expected
+
+    def test_differentiate_wrt_c1(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            self.c0*(self.c1 - self.l_M_tilde)/(self.c2 + self.c3*self.l_M_tilde)**2
+            *exp((self.c1 - self.l_M_tilde)**2 /(2*(self.c2 + self.c3*self.l_M_tilde)**2))
+        )
+        assert fl_M_act.diff(self.c1) == expected
+
+    def test_differentiate_wrt_c2(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            self.c0*(self.l_M_tilde - self.c1)**2/(self.c2 + self.c3*self.l_M_tilde)**3
+            *exp((self.c1 - self.l_M_tilde)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
+        )
+        assert fl_M_act.diff(self.c2) == expected
+
+    def test_differentiate_wrt_c3(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            self.c0*self.l_M_tilde*(self.l_M_tilde - self.c1)**2/(self.c2 + self.c3*self.l_M_tilde)**3
+            *exp((self.c1 - self.l_M_tilde)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
+        )
+        assert fl_M_act.diff(self.c3) == expected
+
+    def test_differentiate_wrt_c4(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = exp((self.c5 - self.l_M_tilde)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
+        assert fl_M_act.diff(self.c4) == expected
+
+    def test_differentiate_wrt_c5(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            self.c4*(self.c5 - self.l_M_tilde)/(self.c6 + self.c7*self.l_M_tilde)**2
+            *exp((self.c5 - self.l_M_tilde)**2 /(2*(self.c6 + self.c7*self.l_M_tilde)**2))
+        )
+        assert fl_M_act.diff(self.c5) == expected
+
+    def test_differentiate_wrt_c6(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            self.c4*(self.l_M_tilde - self.c5)**2/(self.c6 + self.c7*self.l_M_tilde)**3
+            *exp((self.c5 - self.l_M_tilde)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
+        )
+        assert fl_M_act.diff(self.c6) == expected
+
+    def test_differentiate_wrt_c7(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            self.c4*self.l_M_tilde*(self.l_M_tilde - self.c5)**2/(self.c6 + self.c7*self.l_M_tilde)**3
+            *exp((self.c5 - self.l_M_tilde)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
+        )
+        assert fl_M_act.diff(self.c7) == expected
+
+    def test_differentiate_wrt_c8(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = exp((self.c9 - self.l_M_tilde)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
+        assert fl_M_act.diff(self.c8) == expected
+
+    def test_differentiate_wrt_c9(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            self.c8*(self.c9 - self.l_M_tilde)/(self.c10 + self.c11*self.l_M_tilde)**2
+            *exp((self.c9 - self.l_M_tilde)**2 /(2*(self.c10 + self.c11*self.l_M_tilde)**2))
+        )
+        assert fl_M_act.diff(self.c9) == expected
+
+    def test_differentiate_wrt_c10(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            self.c8*(self.l_M_tilde - self.c9)**2/(self.c10 + self.c11*self.l_M_tilde)**3
+            *exp((self.c9 - self.l_M_tilde)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
+        )
+        assert fl_M_act.diff(self.c10) == expected
+
+    def test_differentiate_wrt_c11(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            self.c8*self.l_M_tilde*(self.l_M_tilde - self.c9)**2/(self.c10 + self.c11*self.l_M_tilde)**3
+            *exp((self.c9 - self.l_M_tilde)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
+        )
+        assert fl_M_act.diff(self.c11) == expected

From f4a0e915d62b44e9afdfa20ecc0545951c42155a Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 14:52:57 +0200
Subject: [PATCH 08/23] Simplify expected expression returned by doit method

---
 .../_biomechanics/tests/test_characteristic.py       | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index cef220f2b6da..8bf91b9dfe61 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -830,17 +830,17 @@ def test_instance(self):
     def test_doit(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants).doit()
         assert fl_M_act == (
-            self.c0*exp(-0.5*((self.l_M_tilde - self.c1)/(self.c2 + self.c3*self.l_M_tilde))**2)
-            + self.c4*exp(-0.5*((self.l_M_tilde - self.c5)/(self.c6 + self.c7*self.l_M_tilde))**2)
-            + self.c8*exp(-0.5*((self.l_M_tilde - self.c9)/(self.c10 + self.c11*self.l_M_tilde))**2)
+            self.c0*exp((((self.c1 - self.l_M_tilde)/(self.c2 + self.c3*self.l_M_tilde))**2)/2)
+            + self.c4*exp((((self.c5 - self.l_M_tilde)/(self.c6 + self.c7*self.l_M_tilde))**2)/2)
+            + self.c8*exp((((self.c9 - self.l_M_tilde)/(self.c10 + self.c11*self.l_M_tilde))**2)/2)
         )
 
     def test_doit_evaluate_false(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants).doit(evaluate=False)
         assert fl_M_act == (
-            self.c0*exp(-0.5*(UnevaluatedExpr(self.l_M_tilde - self.c1)/(self.c2 + self.c3*self.l_M_tilde))**2)
-            + self.c4*exp(-0.5*(UnevaluatedExpr(self.l_M_tilde - self.c5)/(self.c6 + self.c7*self.l_M_tilde))**2)
-            + self.c8*exp(-0.5*(UnevaluatedExpr(self.l_M_tilde - self.c9)/(self.c10 + self.c11*self.l_M_tilde))**2)
+            self.c0*exp(((UnevaluatedExpr(self.c1 - self.l_M_tilde)/(self.c2 + self.c3*self.l_M_tilde))**2)/2)
+            + self.c4*exp(((UnevaluatedExpr(self.c5 - self.l_M_tilde)/(self.c6 + self.c7*self.l_M_tilde))**2)/2)
+            + self.c8*exp(((UnevaluatedExpr(self.c9 - self.l_M_tilde)/(self.c10 + self.c11*self.l_M_tilde))**2)/2)
         )
 
     def test_with_default_constants(self):

From 5bda35b70307bfe148ea8129bbfa0f8fefa7082f Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 15:11:56 +0200
Subject: [PATCH 09/23] Add test cases for LaTeX printing

---
 .../_biomechanics/tests/test_characteristic.py     | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index 8bf91b9dfe61..d17cb51d5f79 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -966,3 +966,17 @@ def test_differentiate_wrt_c11(self):
             *exp((self.c9 - self.l_M_tilde)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.c11) == expected
+
+    def test_function_print_latex(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = r'\operatorname{fl}^M_{act} \left( l_{M tilde} \right)'
+        assert LatexPrinter().doprint(fl_M_act) == expected
+
+    def test_expression_print_latex(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
+        expected = (
+            r'c_{0} e^{\frac{\left(c_{1} - l_{M tilde}\right)^{2}}{2 \left(c_{2} + c_{3} l_{M tilde}\right)^{2}}} '
+            r'+ c_{4} e^{\frac{\left(c_{5} - l_{M tilde}\right)^{2}}{2 \left(c_{6} + c_{7} l_{M tilde}\right)^{2}}} '
+            r'+ c_{8} e^{\frac{\left(c_{9} - l_{M tilde}\right)^{2}}{2 \left(c_{10} + c_{11} l_{M tilde}\right)^{2}}}'
+        )
+        assert LatexPrinter().doprint(fl_M_act.doit()) == expected

From 05e21426569119c03ad18e16bb122e4e4d5cab2e Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 15:27:06 +0200
Subject: [PATCH 10/23] Add test cases for code printing

---
 .../tests/test_characteristic.py              | 156 ++++++++++++++++++
 1 file changed, 156 insertions(+)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index d17cb51d5f79..681cc5d710c0 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -980,3 +980,159 @@ def test_expression_print_latex(self):
             r'+ c_{8} e^{\frac{\left(c_{9} - l_{M tilde}\right)^{2}}{2 \left(c_{10} + c_{11} l_{M tilde}\right)^{2}}}'
         )
         assert LatexPrinter().doprint(fl_M_act.doit()) == expected
+
+    @pytest.mark.parametrize(
+        'code_printer, expected',
+        [
+            (
+                C89CodePrinter,
+                (
+                    '0.81399999999999995*exp(19.051973784484073*pow(1.0600000000000001 - l_M_tilde, 2)/pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*exp((25.0/2.0)*pow(0.71699999999999997 - l_M_tilde, 2)/pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*exp(3.9899134986753491*pow(1 - l_M_tilde, 2))'
+                ),
+            ),
+            (
+                C99CodePrinter,
+                (
+                    '0.81399999999999995*exp(19.051973784484073*pow(1.0600000000000001 - l_M_tilde, 2)/pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*exp((25.0/2.0)*pow(0.71699999999999997 - l_M_tilde, 2)/pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*exp(3.9899134986753491*pow(1 - l_M_tilde, 2))'
+                ),
+            ),
+            (
+                C11CodePrinter,
+                (
+                    '0.81399999999999995*exp(19.051973784484073*pow(1.0600000000000001 - l_M_tilde, 2)/pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*exp((25.0/2.0)*pow(0.71699999999999997 - l_M_tilde, 2)/pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*exp(3.9899134986753491*pow(1 - l_M_tilde, 2))'
+                ),
+            ),
+            (
+                CXX98CodePrinter,
+                (
+                    '0.81399999999999995*exp(19.051973784484073*std::pow(1.0600000000000001 - l_M_tilde, 2)/std::pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*exp((25.0/2.0)*std::pow(0.71699999999999997 - l_M_tilde, 2)/std::pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*exp(3.9899134986753491*std::pow(1 - l_M_tilde, 2))'
+                ),
+            ),
+            (
+                CXX11CodePrinter,
+                (
+                    '0.81399999999999995*std::exp(19.051973784484073*std::pow(1.0600000000000001 - l_M_tilde, 2)/std::pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*std::exp((25.0/2.0)*std::pow(0.71699999999999997 - l_M_tilde, 2)/std::pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*std::exp(3.9899134986753491*std::pow(1 - l_M_tilde, 2))'
+                ),
+            ),
+            (
+                CXX17CodePrinter,
+                (
+                    '0.81399999999999995*std::exp(19.051973784484073*std::pow(1.0600000000000001 - l_M_tilde, 2)/std::pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*std::exp((25.0/2.0)*std::pow(0.71699999999999997 - l_M_tilde, 2)/std::pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*std::exp(3.9899134986753491*std::pow(1 - l_M_tilde, 2))'
+                ),
+            ),
+            (
+                FCodePrinter,
+                (
+                    '      0.814d0*exp(19.051973784484073d0*(1.06d0 - l_M_tilde)**2/(\n'
+                    '      @ 0.39074074074074072d0*l_M_tilde + 1.0d0)**2) + 0.433d0*exp(\n'
+                    '     @ 12.5d0\n'
+                    '      @ *(0.717d0 - l_M_tilde)**2/(l_M_tilde - 0.14949999999999999d0)**2\n'
+                    '      @ ) + (1.0d0/10.0d0)*exp(3.9899134986753491d0*(1 - l_M_tilde)**2)'
+                ),
+            ),
+            (
+                OctaveCodePrinter,
+                (
+                    '0.814*exp(19.0519737844841*(1.06 - l_M_tilde).^2./(0.390740740740741*l_M_tilde + 1).^2) '
+                    '+ 0.433*exp(25*(0.717 - l_M_tilde).^2./(2*(l_M_tilde - 0.1495).^2)) '
+                    '+ exp(3.98991349867535*(1 - l_M_tilde).^2)/10'
+                ),
+            ),
+            (
+                PythonCodePrinter,
+                (
+                    '0.814*math.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*math.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*math.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                ),
+            ),
+            (
+                NumPyPrinter,
+                (
+                    '0.814*numpy.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*numpy.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*numpy.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                ),
+            ),
+            (
+                SciPyPrinter,
+                (
+                    '0.814*numpy.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*numpy.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*numpy.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                ),
+            ),
+            (
+                CuPyPrinter,
+                (
+                    '0.814*cupy.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*cupy.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*cupy.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                ),
+            ),
+            (
+                JaxPrinter,
+                (
+                    '0.814*jax.numpy.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*jax.numpy.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*jax.numpy.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                ),
+            ),
+            (
+                MpmathPrinter,
+                (
+                    'mpmath.mpf((0, 7331860193359167, -53, 53))'
+                    '*mpmath.exp(mpmath.mpf((0, 5362653877279683, -48, 53))'
+                    '*(mpmath.mpf((0, 2386907802506363, -51, 52)) - l_M_tilde)**2'
+                    '/(mpmath.mpf((0, 3519479708796943, -53, 52))*l_M_tilde + 1)**2) '
+                    '+ mpmath.mpf((0, 7800234554605699, -54, 53))'
+                    '*mpmath.exp((mpmath.mpf(25)/mpmath.mpf(2))'
+                    '*(mpmath.mpf((0, 6458161865649291, -53, 53)) - l_M_tilde)**2'
+                    '/(l_M_tilde + mpmath.mpf((1, 5386305154335113, -55, 53)))**2) '
+                    '+ (mpmath.mpf(1)/mpmath.mpf(10))'
+                    '*mpmath.exp(mpmath.mpf((0, 8984486472937407, -51, 53))'
+                    '*(1 - l_M_tilde)**2)'
+                ),
+            ),
+            (
+                LambdaPrinter,
+                (
+                    '0.814*math.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*math.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*math.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                ),
+            ),
+        ]
+    )
+    def test_print_code(self, code_printer, expected):
+        fl_M_act = FiberForceLengthActiveDeGroote2016.with_default_constants(self.l_M_tilde)
+        assert code_printer().doprint(fl_M_act) == expected
+
+    def test_derivative_print_code(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016.with_default_constants(self.l_M_tilde)
+        fl_M_act_dl_M_tilde = fl_M_act.diff(self.l_M_tilde)
+        expected = (
+            '(0.79798269973507*l_M_tilde - 0.79798269973507)'
+            '*math.exp(3.98991349867535*(1 - l_M_tilde)**2) '
+            '+ (31.0166133211401*(l_M_tilde - 1.06)/(0.390740740740741*l_M_tilde + 1)**2 '
+            '+ 13.6174190361677*(0.943396226415094*l_M_tilde - 1)**2'
+            '/(0.390740740740741*l_M_tilde + 1)**3)'
+            '*math.exp(21.4067977442463*(1 - 0.943396226415094*l_M_tilde)**2'
+            '/(0.390740740740741*l_M_tilde + 1)**2) '
+            '+ (10.825*(l_M_tilde - 0.717)**2/(l_M_tilde - 0.1495)**3 '
+            '+ 10.825*(l_M_tilde - 0.717)/(l_M_tilde - 0.1495)**2)'
+            '*math.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2)'
+        )
+        assert PythonCodePrinter().doprint(fl_M_act_dl_M_tilde) == expected

From e68b05ef813cd882dcf9aa2241114c663a4869f2 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 15:47:53 +0200
Subject: [PATCH 11/23] Fix mistake in doit method

---
 .../_biomechanics/tests/test_characteristic.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index 681cc5d710c0..fa438d83cdb0 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -830,17 +830,17 @@ def test_instance(self):
     def test_doit(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants).doit()
         assert fl_M_act == (
-            self.c0*exp((((self.c1 - self.l_M_tilde)/(self.c2 + self.c3*self.l_M_tilde))**2)/2)
-            + self.c4*exp((((self.c5 - self.l_M_tilde)/(self.c6 + self.c7*self.l_M_tilde))**2)/2)
-            + self.c8*exp((((self.c9 - self.l_M_tilde)/(self.c10 + self.c11*self.l_M_tilde))**2)/2)
+            self.c0*exp(-(((self.l_M_tilde - self.c1)/(self.c2 + self.c3*self.l_M_tilde))**2)/2)
+            + self.c4*exp(-(((self.l_M_tilde - self.c5)/(self.c6 + self.c7*self.l_M_tilde))**2)/2)
+            + self.c8*exp(-(((self.l_M_tilde - self.c9)/(self.c10 + self.c11*self.l_M_tilde))**2)/2)
         )
 
     def test_doit_evaluate_false(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants).doit(evaluate=False)
         assert fl_M_act == (
-            self.c0*exp(((UnevaluatedExpr(self.c1 - self.l_M_tilde)/(self.c2 + self.c3*self.l_M_tilde))**2)/2)
-            + self.c4*exp(((UnevaluatedExpr(self.c5 - self.l_M_tilde)/(self.c6 + self.c7*self.l_M_tilde))**2)/2)
-            + self.c8*exp(((UnevaluatedExpr(self.c9 - self.l_M_tilde)/(self.c10 + self.c11*self.l_M_tilde))**2)/2)
+            self.c0*exp(-((UnevaluatedExpr(self.l_M_tilde - self.c1)/(self.c2 + self.c3*self.l_M_tilde))**2)/2)
+            + self.c4*exp(-((UnevaluatedExpr(self.l_M_tilde - self.c5)/(self.c6 + self.c7*self.l_M_tilde))**2)/2)
+            + self.c8*exp(-((UnevaluatedExpr(self.l_M_tilde - self.c9)/(self.c10 + self.c11*self.l_M_tilde))**2)/2)
         )
 
     def test_with_default_constants(self):
@@ -975,9 +975,9 @@ def test_function_print_latex(self):
     def test_expression_print_latex(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
         expected = (
-            r'c_{0} e^{\frac{\left(c_{1} - l_{M tilde}\right)^{2}}{2 \left(c_{2} + c_{3} l_{M tilde}\right)^{2}}} '
-            r'+ c_{4} e^{\frac{\left(c_{5} - l_{M tilde}\right)^{2}}{2 \left(c_{6} + c_{7} l_{M tilde}\right)^{2}}} '
-            r'+ c_{8} e^{\frac{\left(c_{9} - l_{M tilde}\right)^{2}}{2 \left(c_{10} + c_{11} l_{M tilde}\right)^{2}}}'
+            r'c_{0} e^{- \frac{\left(- c_{1} + l_{M tilde}\right)^{2}}{2 \left(c_{2} + c_{3} l_{M tilde}\right)^{2}}} '
+            r'+ c_{4} e^{- \frac{\left(- c_{5} + l_{M tilde}\right)^{2}}{2 \left(c_{6} + c_{7} l_{M tilde}\right)^{2}}} '
+            r'+ c_{8} e^{- \frac{\left(- c_{9} + l_{M tilde}\right)^{2}}{2 \left(c_{10} + c_{11} l_{M tilde}\right)^{2}}}'
         )
         assert LatexPrinter().doprint(fl_M_act.doit()) == expected
 

From a501db88a9d24ebe60b2ae33a1067654c3e00b12 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 15:49:33 +0200
Subject: [PATCH 12/23] Add test cases for lambdify

---
 .../tests/test_characteristic.py              | 33 +++++++++++++++++++
 1 file changed, 33 insertions(+)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index fa438d83cdb0..8cc41091e017 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -1136,3 +1136,36 @@ def test_derivative_print_code(self):
             '*math.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2)'
         )
         assert PythonCodePrinter().doprint(fl_M_act_dl_M_tilde) == expected
+
+    def test_lambdify(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016.with_default_constants(self.l_M_tilde)
+        fl_M_act_callable = lambdify(self.l_M_tilde, fl_M_act)
+        assert fl_M_act_callable(1.0) == pytest.approx(0.9941398866)
+
+    @pytest.mark.skipif(numpy is None, reason='NumPy not installed')
+    def test_lambdify_numpy(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016.with_default_constants(self.l_M_tilde)
+        fl_M_act_callable = lambdify(self.l_M_tilde, fl_M_act, 'numpy')
+        l_M_tilde = numpy.array([0.0, 0.5, 1.0, 1.5, 2.0])
+        expected = numpy.array([
+            0.0018501319,
+            0.0529122812,
+            0.9941398866,
+            0.2312431531,
+            0.0069595432,
+        ])
+        numpy.testing.assert_allclose(fl_M_act_callable(l_M_tilde), expected)
+
+    @pytest.mark.skipif(jax is None, reason='JAX not installed')
+    def test_lambdify_jax(self):
+        fl_M_act = FiberForceLengthActiveDeGroote2016.with_default_constants(self.l_M_tilde)
+        fl_M_act_callable = jax.jit(lambdify(self.l_M_tilde, fl_M_act, 'jax'))
+        l_M_tilde = jax.numpy.array([0.0, 0.5, 1.0, 1.5, 2.0])
+        expected = jax.numpy.array([
+            0.0018501319,
+            0.0529122812,
+            0.9941398866,
+            0.2312431531,
+            0.0069595432,
+        ])
+        numpy.testing.assert_allclose(fl_M_act_callable(l_M_tilde), expected)

From 1ad5323eed58c861020cb82efabd699c16e4367d Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:20:23 +0200
Subject: [PATCH 13/23] Fix mistakes in expected results from fdiff method

---
 .../tests/test_characteristic.py              | 44 +++++++++----------
 1 file changed, 22 insertions(+), 22 deletions(-)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index 8cc41091e017..938d4fbf6c1e 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -867,29 +867,29 @@ def test_differentiate_wrt_l_M_tilde(self):
         expected = (
             self.c0*(
                 self.c3*(self.l_M_tilde - self.c1)**2/(self.c2 + self.c3*self.l_M_tilde)**3
-                + (self.l_M_tilde - self.c1)/((self.c2 + self.c3*self.l_M_tilde)**2)
-            )*exp((self.c1 - self.l_M_tilde)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
+                + (self.c1 - self.l_M_tilde)/((self.c2 + self.c3*self.l_M_tilde)**2)
+            )*exp(-(self.l_M_tilde - self.c1)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
             + self.c4*(
                 self.c7*(self.l_M_tilde - self.c5)**2/(self.c6 + self.c7*self.l_M_tilde)**3
-                + (self.l_M_tilde - self.c5)/((self.c6 + self.c7*self.l_M_tilde)**2)
-            )*exp((self.c5 - self.l_M_tilde)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
+                + (self.c5 - self.l_M_tilde)/((self.c6 + self.c7*self.l_M_tilde)**2)
+            )*exp(-(self.l_M_tilde - self.c5)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
             + self.c8*(
                 self.c11*(self.l_M_tilde - self.c9)**2/(self.c10 + self.c11*self.l_M_tilde)**3
-                + (self.l_M_tilde - self.c9)/((self.c10 + self.c11*self.l_M_tilde)**2)
-            )*exp((self.c9 - self.l_M_tilde)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
+                + (self.c9 - self.l_M_tilde)/((self.c10 + self.c11*self.l_M_tilde)**2)
+            )*exp(-(self.l_M_tilde - self.c9)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.l_M_tilde) == expected
 
     def test_differentiate_wrt_c0(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
-        expected = exp((self.c1 - self.l_M_tilde)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
-        assert fl_M_act.diff(self.c0) == expected
+        expected = exp(-(self.l_M_tilde - self.c1)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
+        assert fl_M_act.doit().diff(self.c0) == expected
 
     def test_differentiate_wrt_c1(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
         expected = (
-            self.c0*(self.c1 - self.l_M_tilde)/(self.c2 + self.c3*self.l_M_tilde)**2
-            *exp((self.c1 - self.l_M_tilde)**2 /(2*(self.c2 + self.c3*self.l_M_tilde)**2))
+            self.c0*(self.l_M_tilde - self.c1)/(self.c2 + self.c3*self.l_M_tilde)**2
+            *exp(-(self.l_M_tilde - self.c1)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.c1) == expected
 
@@ -897,7 +897,7 @@ def test_differentiate_wrt_c2(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
         expected = (
             self.c0*(self.l_M_tilde - self.c1)**2/(self.c2 + self.c3*self.l_M_tilde)**3
-            *exp((self.c1 - self.l_M_tilde)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
+            *exp(-(self.l_M_tilde - self.c1)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.c2) == expected
 
@@ -905,20 +905,20 @@ def test_differentiate_wrt_c3(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
         expected = (
             self.c0*self.l_M_tilde*(self.l_M_tilde - self.c1)**2/(self.c2 + self.c3*self.l_M_tilde)**3
-            *exp((self.c1 - self.l_M_tilde)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
+            *exp(-(self.l_M_tilde - self.c1)**2/(2*(self.c2 + self.c3*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.c3) == expected
 
     def test_differentiate_wrt_c4(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
-        expected = exp((self.c5 - self.l_M_tilde)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
+        expected = exp(-(self.l_M_tilde - self.c5)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
         assert fl_M_act.diff(self.c4) == expected
 
     def test_differentiate_wrt_c5(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
         expected = (
-            self.c4*(self.c5 - self.l_M_tilde)/(self.c6 + self.c7*self.l_M_tilde)**2
-            *exp((self.c5 - self.l_M_tilde)**2 /(2*(self.c6 + self.c7*self.l_M_tilde)**2))
+            self.c4*(self.l_M_tilde - self.c5)/(self.c6 + self.c7*self.l_M_tilde)**2
+            *exp(-(self.l_M_tilde - self.c5)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.c5) == expected
 
@@ -926,7 +926,7 @@ def test_differentiate_wrt_c6(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
         expected = (
             self.c4*(self.l_M_tilde - self.c5)**2/(self.c6 + self.c7*self.l_M_tilde)**3
-            *exp((self.c5 - self.l_M_tilde)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
+            *exp(-(self.l_M_tilde - self.c5)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.c6) == expected
 
@@ -934,20 +934,20 @@ def test_differentiate_wrt_c7(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
         expected = (
             self.c4*self.l_M_tilde*(self.l_M_tilde - self.c5)**2/(self.c6 + self.c7*self.l_M_tilde)**3
-            *exp((self.c5 - self.l_M_tilde)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
+            *exp(-(self.l_M_tilde - self.c5)**2/(2*(self.c6 + self.c7*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.c7) == expected
 
     def test_differentiate_wrt_c8(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
-        expected = exp((self.c9 - self.l_M_tilde)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
+        expected = exp(-(self.l_M_tilde - self.c9)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
         assert fl_M_act.diff(self.c8) == expected
 
     def test_differentiate_wrt_c9(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
         expected = (
-            self.c8*(self.c9 - self.l_M_tilde)/(self.c10 + self.c11*self.l_M_tilde)**2
-            *exp((self.c9 - self.l_M_tilde)**2 /(2*(self.c10 + self.c11*self.l_M_tilde)**2))
+            self.c8*(self.l_M_tilde - self.c9)/(self.c10 + self.c11*self.l_M_tilde)**2
+            *exp(-(self.l_M_tilde - self.c9)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.c9) == expected
 
@@ -955,7 +955,7 @@ def test_differentiate_wrt_c10(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
         expected = (
             self.c8*(self.l_M_tilde - self.c9)**2/(self.c10 + self.c11*self.l_M_tilde)**3
-            *exp((self.c9 - self.l_M_tilde)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
+            *exp(-(self.l_M_tilde - self.c9)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.c10) == expected
 
@@ -963,7 +963,7 @@ def test_differentiate_wrt_c11(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016(self.l_M_tilde, *self.constants)
         expected = (
             self.c8*self.l_M_tilde*(self.l_M_tilde - self.c9)**2/(self.c10 + self.c11*self.l_M_tilde)**3
-            *exp((self.c9 - self.l_M_tilde)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
+            *exp(-(self.l_M_tilde - self.c9)**2/(2*(self.c10 + self.c11*self.l_M_tilde)**2))
         )
         assert fl_M_act.diff(self.c11) == expected
 

From 058abc96ca7513d6b0906d3b554dd95335adfd78 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:27:53 +0200
Subject: [PATCH 14/23] Fix mistakes in expected code printing

---
 .../tests/test_characteristic.py              | 158 +++++++++++-------
 1 file changed, 101 insertions(+), 57 deletions(-)

diff --git a/sympy/physics/_biomechanics/tests/test_characteristic.py b/sympy/physics/_biomechanics/tests/test_characteristic.py
index 938d4fbf6c1e..9138daa77bd6 100644
--- a/sympy/physics/_biomechanics/tests/test_characteristic.py
+++ b/sympy/physics/_biomechanics/tests/test_characteristic.py
@@ -987,131 +987,175 @@ def test_expression_print_latex(self):
             (
                 C89CodePrinter,
                 (
-                    '0.81399999999999995*exp(19.051973784484073*pow(1.0600000000000001 - l_M_tilde, 2)/pow(0.39074074074074072*l_M_tilde + 1, 2)) '
-                    '+ 0.433*exp((25.0/2.0)*pow(0.71699999999999997 - l_M_tilde, 2)/pow(l_M_tilde - 0.14949999999999999, 2)) '
-                    '+ (1.0/10.0)*exp(3.9899134986753491*pow(1 - l_M_tilde, 2))'
+                    '0.81399999999999995*exp(-19.051973784484073'
+                    '*pow(l_M_tilde - 1.0600000000000001, 2)'
+                    '/pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*exp(-25.0/2.0'
+                    '*pow(l_M_tilde - 0.71699999999999997, 2)'
+                    '/pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*exp(-3.9899134986753491'
+                    '*pow(l_M_tilde - 1, 2))'
                 ),
             ),
             (
                 C99CodePrinter,
                 (
-                    '0.81399999999999995*exp(19.051973784484073*pow(1.0600000000000001 - l_M_tilde, 2)/pow(0.39074074074074072*l_M_tilde + 1, 2)) '
-                    '+ 0.433*exp((25.0/2.0)*pow(0.71699999999999997 - l_M_tilde, 2)/pow(l_M_tilde - 0.14949999999999999, 2)) '
-                    '+ (1.0/10.0)*exp(3.9899134986753491*pow(1 - l_M_tilde, 2))'
+                    '0.81399999999999995*exp(-19.051973784484073'
+                    '*pow(l_M_tilde - 1.0600000000000001, 2)'
+                    '/pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*exp(-25.0/2.0'
+                    '*pow(l_M_tilde - 0.71699999999999997, 2)'
+                    '/pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*exp(-3.9899134986753491'
+                    '*pow(l_M_tilde - 1, 2))'
                 ),
             ),
             (
                 C11CodePrinter,
                 (
-                    '0.81399999999999995*exp(19.051973784484073*pow(1.0600000000000001 - l_M_tilde, 2)/pow(0.39074074074074072*l_M_tilde + 1, 2)) '
-                    '+ 0.433*exp((25.0/2.0)*pow(0.71699999999999997 - l_M_tilde, 2)/pow(l_M_tilde - 0.14949999999999999, 2)) '
-                    '+ (1.0/10.0)*exp(3.9899134986753491*pow(1 - l_M_tilde, 2))'
+                    '0.81399999999999995*exp(-19.051973784484073'
+                    '*pow(l_M_tilde - 1.0600000000000001, 2)'
+                    '/pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*exp(-25.0/2.0'
+                    '*pow(l_M_tilde - 0.71699999999999997, 2)'
+                    '/pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*exp(-3.9899134986753491'
+                    '*pow(l_M_tilde - 1, 2))'
                 ),
             ),
             (
                 CXX98CodePrinter,
                 (
-                    '0.81399999999999995*exp(19.051973784484073*std::pow(1.0600000000000001 - l_M_tilde, 2)/std::pow(0.39074074074074072*l_M_tilde + 1, 2)) '
-                    '+ 0.433*exp((25.0/2.0)*std::pow(0.71699999999999997 - l_M_tilde, 2)/std::pow(l_M_tilde - 0.14949999999999999, 2)) '
-                    '+ (1.0/10.0)*exp(3.9899134986753491*std::pow(1 - l_M_tilde, 2))'
+                    '0.81399999999999995*exp(-19.051973784484073'
+                    '*std::pow(l_M_tilde - 1.0600000000000001, 2)'
+                    '/std::pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*exp(-25.0/2.0'
+                    '*std::pow(l_M_tilde - 0.71699999999999997, 2)'
+                    '/std::pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*exp(-3.9899134986753491'
+                    '*std::pow(l_M_tilde - 1, 2))'
                 ),
             ),
             (
                 CXX11CodePrinter,
                 (
-                    '0.81399999999999995*std::exp(19.051973784484073*std::pow(1.0600000000000001 - l_M_tilde, 2)/std::pow(0.39074074074074072*l_M_tilde + 1, 2)) '
-                    '+ 0.433*std::exp((25.0/2.0)*std::pow(0.71699999999999997 - l_M_tilde, 2)/std::pow(l_M_tilde - 0.14949999999999999, 2)) '
-                    '+ (1.0/10.0)*std::exp(3.9899134986753491*std::pow(1 - l_M_tilde, 2))'
+                    '0.81399999999999995*std::exp(-19.051973784484073'
+                    '*std::pow(l_M_tilde - 1.0600000000000001, 2)'
+                    '/std::pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*std::exp(-25.0/2.0'
+                    '*std::pow(l_M_tilde - 0.71699999999999997, 2)'
+                    '/std::pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*std::exp(-3.9899134986753491'
+                    '*std::pow(l_M_tilde - 1, 2))'
                 ),
             ),
             (
                 CXX17CodePrinter,
                 (
-                    '0.81399999999999995*std::exp(19.051973784484073*std::pow(1.0600000000000001 - l_M_tilde, 2)/std::pow(0.39074074074074072*l_M_tilde + 1, 2)) '
-                    '+ 0.433*std::exp((25.0/2.0)*std::pow(0.71699999999999997 - l_M_tilde, 2)/std::pow(l_M_tilde - 0.14949999999999999, 2)) '
-                    '+ (1.0/10.0)*std::exp(3.9899134986753491*std::pow(1 - l_M_tilde, 2))'
+                    '0.81399999999999995*std::exp(-19.051973784484073'
+                    '*std::pow(l_M_tilde - 1.0600000000000001, 2)'
+                    '/std::pow(0.39074074074074072*l_M_tilde + 1, 2)) '
+                    '+ 0.433*std::exp(-25.0/2.0'
+                    '*std::pow(l_M_tilde - 0.71699999999999997, 2)'
+                    '/std::pow(l_M_tilde - 0.14949999999999999, 2)) '
+                    '+ (1.0/10.0)*std::exp(-3.9899134986753491'
+                    '*std::pow(l_M_tilde - 1, 2))'
                 ),
             ),
             (
                 FCodePrinter,
                 (
-                    '      0.814d0*exp(19.051973784484073d0*(1.06d0 - l_M_tilde)**2/(\n'
+                    '      0.814d0*exp(-19.051973784484073d0*(l_M_tilde - 1.06d0)**2/(\n'
                     '      @ 0.39074074074074072d0*l_M_tilde + 1.0d0)**2) + 0.433d0*exp(\n'
-                    '     @ 12.5d0\n'
-                    '      @ *(0.717d0 - l_M_tilde)**2/(l_M_tilde - 0.14949999999999999d0)**2\n'
-                    '      @ ) + (1.0d0/10.0d0)*exp(3.9899134986753491d0*(1 - l_M_tilde)**2)'
+                    '      @ -12.5d0*(l_M_tilde - 0.717d0)**2/(l_M_tilde -\n'
+                    '      @ 0.14949999999999999d0)**2) + (1.0d0/10.0d0)*exp(\n'
+                    '      @ -3.9899134986753491d0*(l_M_tilde - 1)**2)'
                 ),
             ),
             (
                 OctaveCodePrinter,
                 (
-                    '0.814*exp(19.0519737844841*(1.06 - l_M_tilde).^2./(0.390740740740741*l_M_tilde + 1).^2) '
-                    '+ 0.433*exp(25*(0.717 - l_M_tilde).^2./(2*(l_M_tilde - 0.1495).^2)) '
-                    '+ exp(3.98991349867535*(1 - l_M_tilde).^2)/10'
+                    '0.814*exp(-19.0519737844841*(l_M_tilde - 1.06).^2'
+                    './(0.390740740740741*l_M_tilde + 1).^2) '
+                    '+ 0.433*exp(-25*(l_M_tilde - 0.717).^2'
+                    './(2*(l_M_tilde - 0.1495).^2)) '
+                    '+ exp(-3.98991349867535*(l_M_tilde - 1).^2)/10'
                 ),
             ),
             (
                 PythonCodePrinter,
                 (
-                    '0.814*math.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
-                    '+ 0.433*math.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
-                    '+ (1/10)*math.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                    '0.814*math.exp(-19.0519737844841*(l_M_tilde - 1.06)**2'
+                    '/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*math.exp(-25/2*(l_M_tilde - 0.717)**2'
+                    '/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*math.exp(-3.98991349867535*(l_M_tilde - 1)**2)'
                 ),
             ),
             (
                 NumPyPrinter,
                 (
-                    '0.814*numpy.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
-                    '+ 0.433*numpy.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
-                    '+ (1/10)*numpy.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                    '0.814*numpy.exp(-19.0519737844841*(l_M_tilde - 1.06)**2'
+                    '/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*numpy.exp(-25/2*(l_M_tilde - 0.717)**2'
+                    '/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*numpy.exp(-3.98991349867535*(l_M_tilde - 1)**2)'
                 ),
             ),
             (
                 SciPyPrinter,
                 (
-                    '0.814*numpy.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
-                    '+ 0.433*numpy.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
-                    '+ (1/10)*numpy.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                    '0.814*numpy.exp(-19.0519737844841*(l_M_tilde - 1.06)**2'
+                    '/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*numpy.exp(-25/2*(l_M_tilde - 0.717)**2'
+                    '/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*numpy.exp(-3.98991349867535*(l_M_tilde - 1)**2)'
                 ),
             ),
             (
                 CuPyPrinter,
                 (
-                    '0.814*cupy.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
-                    '+ 0.433*cupy.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
-                    '+ (1/10)*cupy.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                    '0.814*cupy.exp(-19.0519737844841*(l_M_tilde - 1.06)**2'
+                    '/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*cupy.exp(-25/2*(l_M_tilde - 0.717)**2'
+                    '/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*cupy.exp(-3.98991349867535*(l_M_tilde - 1)**2)'
                 ),
             ),
             (
                 JaxPrinter,
                 (
-                    '0.814*jax.numpy.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
-                    '+ 0.433*jax.numpy.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
-                    '+ (1/10)*jax.numpy.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                    '0.814*jax.numpy.exp(-19.0519737844841*(l_M_tilde - 1.06)**2'
+                    '/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*jax.numpy.exp(-25/2*(l_M_tilde - 0.717)**2'
+                    '/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*jax.numpy.exp(-3.98991349867535*(l_M_tilde - 1)**2)'
                 ),
             ),
             (
                 MpmathPrinter,
                 (
                     'mpmath.mpf((0, 7331860193359167, -53, 53))'
-                    '*mpmath.exp(mpmath.mpf((0, 5362653877279683, -48, 53))'
-                    '*(mpmath.mpf((0, 2386907802506363, -51, 52)) - l_M_tilde)**2'
+                    '*mpmath.exp(-mpmath.mpf((0, 5362653877279683, -48, 53))'
+                    '*(l_M_tilde + mpmath.mpf((1, 2386907802506363, -51, 52)))**2'
                     '/(mpmath.mpf((0, 3519479708796943, -53, 52))*l_M_tilde + 1)**2) '
                     '+ mpmath.mpf((0, 7800234554605699, -54, 53))'
-                    '*mpmath.exp((mpmath.mpf(25)/mpmath.mpf(2))'
-                    '*(mpmath.mpf((0, 6458161865649291, -53, 53)) - l_M_tilde)**2'
+                    '*mpmath.exp(-mpmath.mpf(25)/mpmath.mpf(2)'
+                    '*(l_M_tilde + mpmath.mpf((1, 6458161865649291, -53, 53)))**2'
                     '/(l_M_tilde + mpmath.mpf((1, 5386305154335113, -55, 53)))**2) '
                     '+ (mpmath.mpf(1)/mpmath.mpf(10))'
-                    '*mpmath.exp(mpmath.mpf((0, 8984486472937407, -51, 53))'
-                    '*(1 - l_M_tilde)**2)'
+                    '*mpmath.exp(-mpmath.mpf((0, 8984486472937407, -51, 53))'
+                    '*(l_M_tilde - 1)**2)'
                 ),
             ),
             (
                 LambdaPrinter,
                 (
-                    '0.814*math.exp(19.0519737844841*(1.06 - l_M_tilde)**2/(0.390740740740741*l_M_tilde + 1)**2) '
-                    '+ 0.433*math.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2) '
-                    '+ (1/10)*math.exp(3.98991349867535*(1 - l_M_tilde)**2)'
+                    '0.814*math.exp(-19.0519737844841*(l_M_tilde - 1.06)**2'
+                    '/(0.390740740740741*l_M_tilde + 1)**2) '
+                    '+ 0.433*math.exp(-25/2*(l_M_tilde - 0.717)**2'
+                    '/(l_M_tilde - 0.1495)**2) '
+                    '+ (1/10)*math.exp(-3.98991349867535*(l_M_tilde - 1)**2)'
                 ),
             ),
         ]
@@ -1124,16 +1168,16 @@ def test_derivative_print_code(self):
         fl_M_act = FiberForceLengthActiveDeGroote2016.with_default_constants(self.l_M_tilde)
         fl_M_act_dl_M_tilde = fl_M_act.diff(self.l_M_tilde)
         expected = (
-            '(0.79798269973507*l_M_tilde - 0.79798269973507)'
-            '*math.exp(3.98991349867535*(1 - l_M_tilde)**2) '
-            '+ (31.0166133211401*(l_M_tilde - 1.06)/(0.390740740740741*l_M_tilde + 1)**2 '
+            '(0.79798269973507 - 0.79798269973507*l_M_tilde)'
+            '*math.exp(-3.98991349867535*(l_M_tilde - 1)**2) '
+            '+ (10.825*(0.717 - l_M_tilde)/(l_M_tilde - 0.1495)**2 '
+            '+ 10.825*(l_M_tilde - 0.717)**2/(l_M_tilde - 0.1495)**3)'
+            '*math.exp(-25/2*(l_M_tilde - 0.717)**2/(l_M_tilde - 0.1495)**2) '
+            '+ (31.0166133211401*(1.06 - l_M_tilde)/(0.390740740740741*l_M_tilde + 1)**2 '
             '+ 13.6174190361677*(0.943396226415094*l_M_tilde - 1)**2'
             '/(0.390740740740741*l_M_tilde + 1)**3)'
-            '*math.exp(21.4067977442463*(1 - 0.943396226415094*l_M_tilde)**2'
-            '/(0.390740740740741*l_M_tilde + 1)**2) '
-            '+ (10.825*(l_M_tilde - 0.717)**2/(l_M_tilde - 0.1495)**3 '
-            '+ 10.825*(l_M_tilde - 0.717)/(l_M_tilde - 0.1495)**2)'
-            '*math.exp((25/2)*(0.717 - l_M_tilde)**2/(l_M_tilde - 0.1495)**2)'
+            '*math.exp(-21.4067977442463*(0.943396226415094*l_M_tilde - 1)**2'
+            '/(0.390740740740741*l_M_tilde + 1)**2)'
         )
         assert PythonCodePrinter().doprint(fl_M_act_dl_M_tilde) == expected
 

From 46f6bfb9cef7c7da31e8cd27f7d7f315433b58a5 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:38:57 +0200
Subject: [PATCH 15/23] Add class-level docstring

---
 sympy/physics/_biomechanics/characteristic.py | 32 ++++++++++++++++++-
 1 file changed, 31 insertions(+), 1 deletion(-)

diff --git a/sympy/physics/_biomechanics/characteristic.py b/sympy/physics/_biomechanics/characteristic.py
index 29e05784c90e..1647bfbc5ae8 100644
--- a/sympy/physics/_biomechanics/characteristic.py
+++ b/sympy/physics/_biomechanics/characteristic.py
@@ -801,4 +801,34 @@ def _latex(self, printer):
 
 
 class FiberForceLengthActiveDeGroote2016(CharacteristicCurveFunction):
-    pass
+    r"""Active muscle fiber force-length curve based on De Groote et al., 2016
+    [1].
+
+    Explanation
+    ===========
+
+    The function is defined by the equation:
+
+    $fl^M_{act} = c_0 \exp{-\frac{1}{2}\left(\frac{\tilde{l}^M - c_1}{\left(c_2 + c_3 \tilde{l}^M\right)}\right)^2}
+        + c_4 \exp{-\frac{1}{2}\left(\frac{\tilde{l}^M - c_5}{\left(c_6 + c_7 \tilde{l}^M\right)}\right)^2}
+        + c_8 \exp{-\frac{1}{2}\left(\frac{\tilde{l}^M - c_9}{\left(c_{10} + c_{11} \tilde{l}^M\right)}\right)^2}$
+
+    with constant values of $c0 = 0.814$, $c1 = 1.06$, $c2 = 0.162$,
+    $c3 = 0.0633$, $c4 = 0.433$, $c5 = 0.717$, $c6 = -0.0299$, $c7 = 0.2$,
+    $c8 = 0.1$, $c9 = 1.0$, $c10 = 0.354$, and $c11 = 0.0$.
+
+    While it is possible to change the constant values, these were carefully
+    selected in the original publication to give the characteristic curve
+    specific and required properties. For example, the function produces a
+    active fiber force of 1 at a normalized fiber length of 1, and an active
+    fiber force of 0 at normalized fiber lengths of 0 and 2.
+
+    References
+    ==========
+
+    .. [1] De Groote, F., Kinney, A. L., Rao, A. V., & Fregly, B. J., Evaluation
+           of direct collocation optimal control problem formulations for
+           solving the muscle redundancy problem, Annals of biomedical
+           engineering, 44(10), (2016) pp. 2922-2936
+
+    """

From 6257e10f35e4c1ace771e629ef84be617af379d8 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:39:23 +0200
Subject: [PATCH 16/23] Implement with_default_constants alternate constructor

---
 sympy/physics/_biomechanics/characteristic.py | 48 +++++++++++++++++++
 1 file changed, 48 insertions(+)

diff --git a/sympy/physics/_biomechanics/characteristic.py b/sympy/physics/_biomechanics/characteristic.py
index 1647bfbc5ae8..fd9bab2e4200 100644
--- a/sympy/physics/_biomechanics/characteristic.py
+++ b/sympy/physics/_biomechanics/characteristic.py
@@ -832,3 +832,51 @@ class FiberForceLengthActiveDeGroote2016(CharacteristicCurveFunction):
            engineering, 44(10), (2016) pp. 2922-2936
 
     """
+
+    @classmethod
+    def with_default_constants(cls, l_M_tilde):
+        r"""Recommended constructor that will use the published constants.
+
+        Explanation
+        ===========
+
+        Returns a new instance of the inverse muscle fiber act force-length
+        function using the four constant values specified in the original
+        publication.
+
+        These have the values:
+
+        $c0 = 0.814$
+        $c1 = 1.06$
+        $c2 = 0.162$
+        $c3 = 0.0633$
+        $c4 = 0.433$
+        $c5 = 0.717$
+        $c6 = -0.0299$
+        $c7 = 0.2$
+        $c8 = 0.1$
+        $c9 = 1.0$
+        $c10 = 0.354$
+        $c11 = 0.0$
+
+        Parameters
+        ==========
+
+        fl_M_act : Any (sympifiable)
+            Normalized passive muscle fiber force as a function of muscle fiber
+            length.
+
+        """
+        c0 = Float('0.814')
+        c1 = Float('1.06')
+        c2 = Float('0.162')
+        c3 = Float('0.0633')
+        c4 = Float('0.433')
+        c5 = Float('0.717')
+        c6 = Float('-0.0299')
+        c7 = Rational(1, 5)
+        c8 = Rational(1, 10)
+        c9 = Integer(1)
+        c10 = Float('0.354')
+        c11 = Integer(0)
+        return cls(l_M_tilde, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11)

From b1bf54095c53dd3b909af259711dca291fb1fec8 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:39:44 +0200
Subject: [PATCH 17/23] Implement eval and _eval_evalf evaluation methods

---
 sympy/physics/_biomechanics/characteristic.py | 53 +++++++++++++++++++
 1 file changed, 53 insertions(+)

diff --git a/sympy/physics/_biomechanics/characteristic.py b/sympy/physics/_biomechanics/characteristic.py
index fd9bab2e4200..399e637469eb 100644
--- a/sympy/physics/_biomechanics/characteristic.py
+++ b/sympy/physics/_biomechanics/characteristic.py
@@ -880,3 +880,56 @@ def with_default_constants(cls, l_M_tilde):
         c10 = Float('0.354')
         c11 = Integer(0)
         return cls(l_M_tilde, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11)
+
+    @classmethod
+    def eval(cls, l_M_tilde, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11):
+        """Evaluation of basic inputs.
+
+        Parameters
+        ==========
+
+        l_M_tilde : Any (sympifiable)
+            Normalized muscle fiber length.
+        c0 : Any (sympifiable)
+            The first constant in the characteristic equation. The published
+            value is ``0.814``.
+        c1 : Any (sympifiable)
+            The second constant in the characteristic equation. The published
+            value is ``1.06``.
+        c2 : Any (sympifiable)
+            The third constant in the characteristic equation. The published
+            value is ``0.162``.
+        c3 : Any (sympifiable)
+            The fourth constant in the characteristic equation. The published
+            value is ``0.0633``.
+        c4 : Any (sympifiable)
+            The fifth constant in the characteristic equation. The published
+            value is ``0.433``.
+        c5 : Any (sympifiable)
+            The sixth constant in the characteristic equation. The published
+            value is ``0.717``.
+        c6 : Any (sympifiable)
+            The seventh constant in the characteristic equation. The published
+            value is ``-0.0299``.
+        c7 : Any (sympifiable)
+            The eighth constant in the characteristic equation. The published
+            value is ``0.2``.
+        c8 : Any (sympifiable)
+            The ninth constant in the characteristic equation. The published
+            value is ``0.1``.
+        c9 : Any (sympifiable)
+            The tenth constant in the characteristic equation. The published
+            value is ``1.0``.
+        c10 : Any (sympifiable)
+            The eleventh constant in the characteristic equation. The published
+            value is ``0.354``.
+        c11 : Any (sympifiable)
+            The tweflth constant in the characteristic equation. The published
+            value is ``0.0``.
+
+        """
+        pass
+
+    def _eval_evalf(self, prec):
+        """Evaluate the expression numerically using ``evalf``."""
+        return self.doit(deep=False, evaluate=False)._eval_evalf(prec)

From 4ec0afb8926f9e1154f5f2b58cda744bcb52f240 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:40:05 +0200
Subject: [PATCH 18/23] Implement doit method

---
 sympy/physics/_biomechanics/characteristic.py | 39 +++++++++++++++++++
 1 file changed, 39 insertions(+)

diff --git a/sympy/physics/_biomechanics/characteristic.py b/sympy/physics/_biomechanics/characteristic.py
index 399e637469eb..3042ed634a54 100644
--- a/sympy/physics/_biomechanics/characteristic.py
+++ b/sympy/physics/_biomechanics/characteristic.py
@@ -933,3 +933,42 @@ def eval(cls, l_M_tilde, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11):
     def _eval_evalf(self, prec):
         """Evaluate the expression numerically using ``evalf``."""
         return self.doit(deep=False, evaluate=False)._eval_evalf(prec)
+
+    def doit(self, deep=True, evaluate=True, **hints):
+        """Evaluate the expression defining the function.
+
+        Parameters
+        ==========
+
+        deep : bool
+            Whether ``doit`` should be recursively called. Default is ``True``.
+        evaluate : bool.
+            Whether the SymPy expression should be evaluated as it is
+            constructed. If ``False``, then no constant folding will be
+            conducted which will leave the expression in a more numerically-
+            stable for values of ``l_M_tilde`` that correspond to a sensible
+            operating range for a musculotendon. Default is ``True``.
+        **kwargs : dict[str, Any]
+            Additional keyword argument pairs to be recursively passed to
+            ``doit``.
+
+        """
+        l_M_tilde, *constants = self.args
+        if deep:
+            hints['evaluate'] = evaluate
+            l_M_tilde = l_M_tilde.doit(deep=deep, **hints)
+            constants = [c.doit(deep=deep, **hints) for c in constants]
+        c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11 = constants
+
+        if evaluate:
+            return (
+                c0*exp(-(((l_M_tilde - c1)/(c2 + c3*l_M_tilde))**2)/2)
+                + c4*exp(-(((l_M_tilde - c5)/(c6 + c7*l_M_tilde))**2)/2)
+                + c8*exp(-(((l_M_tilde - c9)/(c10 + c11*l_M_tilde))**2)/2)
+            )
+
+        return (
+            c0*exp(-((UnevaluatedExpr(l_M_tilde - c1)/(c2 + c3*l_M_tilde))**2)/2)
+            + c4*exp(-((UnevaluatedExpr(l_M_tilde - c5)/(c6 + c7*l_M_tilde))**2)/2)
+            + c8*exp(-((UnevaluatedExpr(l_M_tilde - c9)/(c10 + c11*l_M_tilde))**2)/2)
+        )

From 969273cf6066159e1c065c30848a7fff10aad73c Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:40:24 +0200
Subject: [PATCH 19/23] Implement fdiff method for differentiation

---
 sympy/physics/_biomechanics/characteristic.py | 82 +++++++++++++++++++
 1 file changed, 82 insertions(+)

diff --git a/sympy/physics/_biomechanics/characteristic.py b/sympy/physics/_biomechanics/characteristic.py
index 3042ed634a54..86a709467518 100644
--- a/sympy/physics/_biomechanics/characteristic.py
+++ b/sympy/physics/_biomechanics/characteristic.py
@@ -972,3 +972,85 @@ def doit(self, deep=True, evaluate=True, **hints):
             + c4*exp(-((UnevaluatedExpr(l_M_tilde - c5)/(c6 + c7*l_M_tilde))**2)/2)
             + c8*exp(-((UnevaluatedExpr(l_M_tilde - c9)/(c10 + c11*l_M_tilde))**2)/2)
         )
+
+    def fdiff(self, argindex=1):
+        """Derivative of the function with respect to a single argument.
+
+        Parameters
+        ==========
+
+        argindex : int
+            The index of the function's arguments with respect to which the
+            derivative should be taken. Argument indexes start at ``1``.
+            Default is ``1``.
+
+        """
+        l_M_tilde, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11 = self.args
+        if argindex == 1:
+            return (
+                c0*(
+                    c3*(l_M_tilde - c1)**2/(c2 + c3*l_M_tilde)**3
+                    + (c1 - l_M_tilde)/((c2 + c3*l_M_tilde)**2)
+                )*exp(-(l_M_tilde - c1)**2/(2*(c2 + c3*l_M_tilde)**2))
+                + c4*(
+                    c7*(l_M_tilde - c5)**2/(c6 + c7*l_M_tilde)**3
+                    + (c5 - l_M_tilde)/((c6 + c7*l_M_tilde)**2)
+                )*exp(-(l_M_tilde - c5)**2/(2*(c6 + c7*l_M_tilde)**2))
+                + c8*(
+                    c11*(l_M_tilde - c9)**2/(c10 + c11*l_M_tilde)**3
+                    + (c9 - l_M_tilde)/((c10 + c11*l_M_tilde)**2)
+                )*exp(-(l_M_tilde - c9)**2/(2*(c10 + c11*l_M_tilde)**2))
+            )
+        elif argindex == 2:
+            return exp(-(l_M_tilde - c1)**2/(2*(c2 + c3*l_M_tilde)**2))
+        elif argindex == 3:
+            return (
+                c0*(l_M_tilde - c1)/(c2 + c3*l_M_tilde)**2
+                *exp(-(l_M_tilde - c1)**2 /(2*(c2 + c3*l_M_tilde)**2))
+            )
+        elif argindex == 4:
+            return (
+                c0*(l_M_tilde - c1)**2/(c2 + c3*l_M_tilde)**3
+                *exp(-(l_M_tilde - c1)**2/(2*(c2 + c3*l_M_tilde)**2))
+            )
+        elif argindex == 5:
+            return (
+                c0*l_M_tilde*(l_M_tilde - c1)**2/(c2 + c3*l_M_tilde)**3
+                *exp(-(l_M_tilde - c1)**2/(2*(c2 + c3*l_M_tilde)**2))
+            )
+        elif argindex == 6:
+            return exp(-(l_M_tilde - c5)**2/(2*(c6 + c7*l_M_tilde)**2))
+        elif argindex == 7:
+            return (
+                c4*(l_M_tilde - c5)/(c6 + c7*l_M_tilde)**2
+                *exp(-(l_M_tilde - c5)**2 /(2*(c6 + c7*l_M_tilde)**2))
+            )
+        elif argindex == 8:
+            return (
+                c4*(l_M_tilde - c5)**2/(c6 + c7*l_M_tilde)**3
+                *exp(-(l_M_tilde - c5)**2/(2*(c6 + c7*l_M_tilde)**2))
+            )
+        elif argindex == 9:
+            return (
+                c4*l_M_tilde*(l_M_tilde - c5)**2/(c6 + c7*l_M_tilde)**3
+                *exp(-(l_M_tilde - c5)**2/(2*(c6 + c7*l_M_tilde)**2))
+            )
+        elif argindex == 10:
+            return exp(-(l_M_tilde - c9)**2/(2*(c10 + c11*l_M_tilde)**2))
+        elif argindex == 11:
+            return (
+                c8*(l_M_tilde - c9)/(c10 + c11*l_M_tilde)**2
+                *exp(-(l_M_tilde - c9)**2 /(2*(c10 + c11*l_M_tilde)**2))
+            )
+        elif argindex == 12:
+            return (
+                c8*(l_M_tilde - c9)**2/(c10 + c11*l_M_tilde)**3
+                *exp(-(l_M_tilde - c9)**2/(2*(c10 + c11*l_M_tilde)**2))
+            )
+        elif argindex == 13:
+            return (
+                c8*l_M_tilde*(l_M_tilde - c9)**2/(c10 + c11*l_M_tilde)**3
+                *exp(-(l_M_tilde - c9)**2/(2*(c10 + c11*l_M_tilde)**2))
+            )
+
+        raise ArgumentIndexError(self, argindex)

From a105a349cfa5edc2c99e8278f9577f14500ba880 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:40:39 +0200
Subject: [PATCH 20/23] Implement _latex method for LaTeX printing

---
 sympy/physics/_biomechanics/characteristic.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/sympy/physics/_biomechanics/characteristic.py b/sympy/physics/_biomechanics/characteristic.py
index 86a709467518..aa309bac4dc3 100644
--- a/sympy/physics/_biomechanics/characteristic.py
+++ b/sympy/physics/_biomechanics/characteristic.py
@@ -1054,3 +1054,17 @@ def fdiff(self, argindex=1):
             )
 
         raise ArgumentIndexError(self, argindex)
+
+    def _latex(self, printer):
+        """Print a LaTeX representation of the function defining the curve.
+
+        Parameters
+        ==========
+
+        printer : Printer
+            The printer to be used to print the LaTeX string representation.
+
+        """
+        l_M_tilde = self.args[0]
+        _l_M_tilde = printer._print(l_M_tilde)
+        return r'\operatorname{fl}^M_{act} \left( %s \right)' % _l_M_tilde

From 179775dbdf39fb6f38b67f8d20cd0e18c7525ea7 Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:41:11 +0200
Subject: [PATCH 21/23] Add FiberForceLengthActiveDeGroote2016 class to module
 namespace

---
 sympy/physics/_biomechanics/__init__.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/sympy/physics/_biomechanics/__init__.py b/sympy/physics/_biomechanics/__init__.py
index a81a06b0dda7..ddacbee2c30b 100644
--- a/sympy/physics/_biomechanics/__init__.py
+++ b/sympy/physics/_biomechanics/__init__.py
@@ -7,6 +7,7 @@
 """
 
 from .characteristic import (
+   FiberForceLengthActiveDeGroote2016,
    FiberForceLengthPassiveDeGroote2016,
    FiberForceLengthPassiveInverseDeGroote2016,
    TendonForceLengthDeGroote2016,
@@ -16,6 +17,7 @@
 
 __all__ = [
    # Musculotendon characteristic curve functions
+   'FiberForceLengthActiveDeGroote2016',
    'FiberForceLengthPassiveDeGroote2016',
    'FiberForceLengthPassiveInverseDeGroote2016',
    'TendonForceLengthDeGroote2016',

From b79c0175b52bc5b734b72adb877a7f398e19afaa Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:41:33 +0200
Subject: [PATCH 22/23] Add test case for FiberForceLengthActiveDeGroote2016
 class

---
 sympy/core/tests/test_args.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index 5725199e41da..7c917e4fe884 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -3390,6 +3390,12 @@ def test_sympy__physics___biomechanics__characteristic__FiberForceLengthPassiveI
     assert _test_args(FiberForceLengthPassiveInverseDeGroote2016(fl_M_pas, c0, c1))
 
 
+def test_sympy__physics___biomechanics__characteristic__FiberForceLengthActiveDeGroote2016():
+    from sympy.physics._biomechanics import FiberForceLengthActiveDeGroote2016
+    l_M_tilde, c0, c1 = symbols('l_M_tilde, c0, c1')
+    assert _test_args(FiberForceLengthActiveDeGroote2016(l_M_tilde, c0, c1))
+
+
 def test_sympy__physics__paulialgebra__Pauli():
     from sympy.physics.paulialgebra import Pauli
     assert _test_args(Pauli(1))

From 5293584172169b15fe1548d63b368211407cd46d Mon Sep 17 00:00:00 2001
From: Sam Brockie <sambrockie@icloud.com>
Date: Wed, 16 Aug 2023 16:57:34 +0200
Subject: [PATCH 23/23] Fix error in FiberForceLengthActiveDeGroote2016 test

---
 sympy/core/tests/test_args.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index 7c917e4fe884..6397a4004853 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -3392,8 +3392,8 @@ def test_sympy__physics___biomechanics__characteristic__FiberForceLengthPassiveI
 
 def test_sympy__physics___biomechanics__characteristic__FiberForceLengthActiveDeGroote2016():
     from sympy.physics._biomechanics import FiberForceLengthActiveDeGroote2016
-    l_M_tilde, c0, c1 = symbols('l_M_tilde, c0, c1')
-    assert _test_args(FiberForceLengthActiveDeGroote2016(l_M_tilde, c0, c1))
+    l_M_tilde, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11 = symbols('l_M_tilde, c0:12')
+    assert _test_args(FiberForceLengthActiveDeGroote2016(l_M_tilde, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11))
 
 
 def test_sympy__physics__paulialgebra__Pauli():