(* See muscles.m for assumptions and notation.  Edited
 DynLength and DynLengthTen 10/7/96*)




BeginPackage["TendonForce`", "Muscle`"]

IsoM::usage =
	"IsoM[Muscle, ML] returns the force generated by
	isometric contraction of Muscle at length ML"

IsoK::usage = 
	"IsoK[Muscle, MV, ML] returns the force generated by
	Muscle contracting isokinetically at velocity MV and length ML";

NewLength::usage = 
	"NewLength[Muscle,OldLength, omega, dt] returns a
	list containing {Muscle length, Muscle velocity, Tendon stretch}
	given Muscle, OldLength = {Initial length, velocity, tendon
	stretch}, angular velocity of contraction, and angular step size";

StaticLength::usage =
	"StaticLength[Muscle, L] returns the final length of Muscle
	after tendon stretch from an isometric contraction"

DynLength::usage =
	"DynLength[Muscle, {OldL, OldV, Old dLt}, increase in Lm,
	and time step  and returns {Length, Velocity, Tendon Stretch}"

ForceHistory::usage =
	"ForceHistory[LHistory, Muscle, dTime] accepts an array of
	muscle length (without tendon) and returns a matrix of{true ML,
	Muscle Velocity, Tendon Stretch}"

Begin["`Private`"]

IsoM[Muscle_List, l_] := SLengthTension[
					NewArch[Muscle, l][[3]]/
							Muscle[[9]]] *
						2.5 Muscle[[4]]
		(*finds new fiber (sarcomere) length and
			multiplies by 2.5 PCSA - supercedes
			ForceLength*)

IsoK[Muscle_List,v_,l_]:= Block[{a$ = Vf[v, l, Muscle]},
			Return[SLengthTension[ a$[[2]]/
					Muscle[[9]]] *2.5 Muscle[[4]] *
				ForceVelocity[a$[[1]]/Muscle[[9]],
						Muscle[[11]]]]]
		(*finds fiber(sarcomere) length and velocity (in
			a$) looks up length-tension & velocity
			tension parameters*)


StaticLength[Muscle_List,l0_]:= Block[ {dl = 0,dl2 = 1, l=l0,
		a$ = 2.5 Muscle[[13]] Muscle[[4]],
		b$ = - Muscle[[12]]/Muscle[[14]],
		f},
			While[Abs[dl-dl2] > 10^-4 ,
	dl2 = dl;
	f = IsoM[Muscle, l];
	If[f ==0,
		dl = dl2,
		dl = b$ Log[1 + f / a$]];
	l = l0 + dl];Return[l]]
		(*iteratively finds muscle length appropriate for
			muscle with tendon in isometric contraction
			Based on tendon P/Po = a (Exp[b dl/l] - 1)
			or dl = l/b Log[1 + (P/Po)/a.
			units need to be strain, base E logs*) 


DynLength[Muscle_List,OldL_List, dL_, dT_]:= Block[
		{dl = dL, dl2 = OldL[[3]], l = OldL[[1]] + dL,
			v = dL/dT,
			Lo = OldL[[1]] + dL +OldL[[3]],
			Vo = dL/dT,
			a$ = 2.5 Muscle[[13]] Muscle[[4]],
			b$ = Muscle[[12]]/Muscle[[14]]},
	While[Abs[dl - dl2] > 10^-4,
		dl2 = dl;
		dl = b$ Log[1 + IsoK[
					Muscle,
					v,
					l] / a$];
		dl = (dl+dl2)/2;
	(*goes 1/2 way to point where muscle can stretch tendon-
			this makes it approach sol'n faster;
				greater stability going "down"*)
		l = Lo-dl;
		v = (dL - (dl-OldL[[3]])) / dT];
	Return[{l, v,  dl}]]
			(*Same iterative solution, but now incorp's
			velocity - including velocity changes
			associated with changes in tendon length*)

(*DynLengthTen is faster for long tendons, but has trouble
	with short.*)
DynLengthTen[Muscle_List,OldL_List, dL_, dT_]:= Block[
		{Q = OldL[[1]] + dL + OldL[[3]],
			R = (dL + OldL[[3]])/dT,
			t = OldL[[3]],
			a$ = Muscle[[13]] (2.5 Muscle[[4]])},
	t= t/.FindRoot[IsoK[Muscle, R - t/dT, Q - t] ==
			a$ (Exp[Muscle[[14]] t/
						Muscle[[12]]] - 1),
				{t, {OldL[[3]],OldL[[3]] - dL},
					0, Muscle[[12]]}];
	Return[{Q - t, R - t/dT, t}]]


NewLength[Muscle_List,OldLength_List, omega_, dt_]:=
	If[omega == 0.,
		Do[a$ = StaticLength[Muscle,
				OldLength[[1]] + OldLength[[3]] +
							Muscle[[10]] dt];
			Return[{a$, 0,
				OldLength[[1]]+Muscle[[10]] dt-a$+
				  OldLength[[3]]}],{1}],
		If[Muscle[[12]] < Muscle[[3]], 
			DynLength[Muscle,OldLength, Muscle[[10]] dt,
					dt/omega],
		DynLengthTen[Muscle,OldLength, Muscle[[10]] dt,
					dt/omega]]]
		(*Select which method to find muscle force
			unfortunately, developed while trying to
			keep muscle and moment arms in one block
			(hence the omega) just something to watch*)


ForceHistory[kinetic_List, Muscle_List,dt_] := Block[{
	il = N[StaticLength[Muscle, kinetic[[1]]]],
	size = Length[kinetic],
	result = Table[
			{0,0,0}, {i, 1, Length[kinetic]}]},
	
	result[[1]] = {il, (kinetic[[2]]-kinetic[[1]])/dt,
					kinetic[[1]] - il};
	For[ index = 2, index <= size, index++,
		If[IsoM[Muscle,
			  result[[index-1, 1]]] == 0,
			result[[index]] = {kinetic[[index]],
					(kinetic[[index]] -
						kinetic[[index-1]])/dt,
						result[[index-1,3]]},
			result[[index]] = N[NewLength[
						Muscle,
						result[[index-1]],
					(kinetic[[index]] -
						kinetic[[index-1]])/
							(dt Muscle[[8]]),
					(kinetic[[index]] -
						kinetic[[index-1]])/
							Muscle[[8]]]]]];
	Return[result]]
		(*takes a list of passive muscle lengths evenly
			spaced with time increment dt and returns
			appropriate forces by calls to NewLength*)
	
End[]
EndPackage[]