1function [v_CE, dA, tau_m, F_m, F_CE] = muscle(m, phi, l_CE, A, S)
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26N = 1.5;
27K = 5;
28eps_ref = 0.04;
29w = 0.56*m.l_opt;
30c = 0.05;
31t_ecc = 0.01;
32
33
34
35r = [m.r0(1); ...
36 m.r0(2) .* cos(phi(2) - m.phi_max(2)); ...
37 m.r0(3) .* cos(phi(3) - m.phi_max(3))];
38dl_MTU = m.dir .* (...
39 m.rho .* m.r0(1) .* (phi(1)-m.phi_ref(1)) + ...
40 m.rho .* m.r0(2) .* (sin(phi(2) - m.phi_max(2)) - sin(m.phi_ref(2) - m.phi_max(2))) + ...
41 m.rho .* m.r0(3) .* (sin(phi(3) - m.phi_max(3)) - sin(m.phi_ref(3) - m.phi_max(3))));
42
43l_MTU = m.l_opt + m.l_slack + dl_MTU;
44l_SE = l_MTU - l_CE;
45
46
47dA = (S-A)/t_ecc;
48
49
50eps_SE = (l_SE-m.l_slack)/m.l_slack;
51
52if eps_SE > 0
53 f_SE = (eps_SE/eps_ref).^2;
54else
55 f_SE = 0;
56end
57F_SE = f_SE * m.F_max;
58
59
60l_min = m.l_opt - w;
61eps_BE = w/2;
62F_BE = m.F_max*(min(0,(l_CE-l_min))/(eps_BE))^2;
63
64
65eps_PE = w;
66F_PEi = m.F_max*((max(0,l_CE-m.l_opt))/(eps_PE))^2;
67
68
69
70f_l = exp(log(c).*abs((l_CE-m.l_opt)./w).^3);
71
72
73f_v = (F_SE+F_BE)/(A*m.F_max*f_l+F_PEi);
74
75
76
77
78
79
80
81
82
83if f_v <= 1
84 v_CE = m.v_max.*(1-f_v)./(K*f_v+1);
85elseif (f_v > 1) && (f_v <= N)
86 f_v1 = (f_v-N)/(N-1);
87 v_CE = m.v_max.*(1+f_v1)./(7.56*K*f_v1-1);
88else
89
90 v_CE = (f_v-N)*0.01-m.v_max;
91end
92
93F_CE = A * f_v * f_l * m.F_max;
94
95
96F_m = F_SE;
97tau_m = r .* (F_m .* m.dir);
98
99end
100