## 2 Materials and Methods
### Wrapping Muscles
and the associated surfaces \(\mathbf{x}(u,v)\); then it calculates the geodesics starting from the two intersections \(P\) and \(Q\) on the surface, \(\mathbf{c}_{P}\) and \(\mathbf{c}_{Q}\), varying the attachment points location until the original muscle \(AB\) is split in two straight-line muscle units \(AP\) and \(QB\) tangent to the surface while the related geodesics are collinear to each other in correspondence with the geodesics' closest points \(P^{*}\) and \(Q^{*}\)[34]. From the differential geometry, the tangent, normal and binormal unit vectors, \(\hat{\mathbf{t}}\), \(\hat{\mathbf{n}}\) and \(\hat{\mathbf{b}}\), respectively, are given in (24). \[ \begin{array}{l}\hat{\mathbf{t}}=\mathbf{c}^{\prime}=\left[\begin{array}{cc}\mathbf{ x}_{u}&\mathbf{x}_{v}\\ \hline\mathbf{x}^{\prime}\end{array}\right]\left[\begin{array}{l}u^{\prime}\\ v^{\prime}\end{array}\right]\\ \hat{\mathbf{n}}=\frac{\
abla f}{\left|\
abla f\right|}\\ \hat{\mathbf{b}}=\hat{\mathbf{t}}\times\hat{\mathbf{n}}\end{array}  \tag{24}\]  
The objective is to find the unknown geodesics initial positions, the attachment points locations described by the surface parameters \(\left(u_{0P},v_{0P}\right)\) and \(\left(u_{0Q},v_{0Q}\right)\), that set to zero the function \(f\) in (25), in which the first two components impose the surface's tangency of the split muscles, while the last  
Figure 5: Geodesics on a wrapping surface. two components guarantee the geodesics' collinearity. The non-linear closed problem can be solved using the [PERSON] iterative method, numerically calculating numerically the Jacobian function \(\mathbf{J}_{f}\). \[ \begin{split}\mathbf{\varepsilon}=\begin{bmatrix}u_{0P}\\ u_{0Q}\\ v_{0Q}\\ \end{bmatrix}\to\begin{cases}r_{P}=\mathbf{x}(u_{0P},v_{0P})\\ r_{Q}=\mathbf{x}(u_{0Q},v_{0Q})\end{cases}\to\begin{cases}\hat{\mathbf{t}}_{P}=\frac{r _{P}-r_{A}}{|r_{Q}-r_{B}|}\\ \hat{\mathbf{t}}_{Q}=\frac{|r_{Q}-r_{B}|}{|r_{Q}-r_{B}|}\end{cases}\\ \mathbf{f}(\mathbf{\varepsilon})=\begin{bmatrix}\hat{\mathbf{t}}_{P}^{T}\hat{\mathbf{n}}_{P}\\ \hat{\mathbf{t}}_{P}^{T}\hat{\mathbf{n}}_{Q}\end{cases}\mathbf{\Omega}_{Q}\\ \mathbf{\Omega}_{P}^{T}\left(r_{P}\cdot\mathbf{t}_{Q}\right)\end{split} \mathbf{0}\to \mathbf{\varepsilon}^{(k+1)}=\mathbf{\varepsilon}^{(k)}-\mathbf{J}_{f}^{-1}\Big{(}\mathbf{ \varepsilon}^{(k)}\Big{)}\mathbf{f}(\mathbf{\varepsilon}^{(k)})\end{split} \tag{25}\]  
Once the Newton cycle has reached convergence, the two geodesics \(\mathbf{c}_{P}\) and \(\mathbf{c}_{Q}\) are joined in correspondence of the collinearity point, resulting in the definitive muscle path \(c_{PQ}\). The surface types analysed in this work are the cylinder \(x_{c}\), with radius \(r\) and height \(h\), and the ellipsoid \(\mathbf{x}_{e}\), with semiaxes \(s_{x}\), \(s_{y}\) and \(s_{z}\), given in (26). \[ \begin{split}\mathbf{x}_{c}(u,v)=\begin{bmatrix}r\cos u\\ r\sin u\\ v\end{bmatrix}&0\leq u\leq 2\pi\\ -\frac{h}{2}\leq v\leq\frac{h}{2}&\mathbf{x}_{e}(u,v)=\begin{bmatrix}s_{x}\cos v \cos u\\ s_{y}\cos v\sin u\\ s_{z}\sin v\end{bmatrix}&0\leq u\leq 2\pi\\ -\frac{\pi}{2}\leq v\leq\frac{\pi}{2}\end{split}  \tag{26}\] \[f_{c}(\mathbf{x})=x^{2}+y^{2}-r^{2}=0\quad\quad f_{e}(\mathbf{x})= \left(\frac{x}{s_{x}}\right)^{2}+\left(\frac{y}{s_{y}}\right)^{2}+\left(\frac {z}{s_{z}}\right)^{2}-1=0\]