(*A mathematica notebook to calculate length dependent 
muscle properties (fiber length, sarcomere length and 
fiber angles) from a reference configuration and current 
muscle length.  Copyright 1996, Thomas Burkholder
You are licensed to use this package for non-commercial
purposes.  I'd like to hear from you if you do use this, or
need tech support at tburkhol@ucsd.edu*)
(*Assume rigid aponeurosis and tendon; Constant muscle
projection area; Uniform pennation, fiber angle etc
Notation:
	Lf1		Fiber Length in muscle refrence state
	La		Aponeurosis Length in refrence state
	alpha		aponeurosis angle
	beta		Fiber angle
	ML		Muscle length in current state
	MV		Muscle velocity
	theta		Joint angle or current pennation
	ZeroLength  Muscle length at refrence

depends on a list data structure of refrence:
	muscle length,	aponeurosis length,
	fiber length,	PCSA,
	aponeurosis angle,fiber angle,
	pennation angle,	moment arm function(unused)
	kilosarcomere #,	MA,
	fast fiber fraction,
	tendon length,	tendon "a",
	tendon "b"*)

BeginPackage["Muscle`"]

SLengthTension::usage =
		"SLengthTension[length] gives the normalized force
		produced by a human sarcomere at length in µm"

SVelocityTension::usage =
		"SVelocityTension[vin, Vmax] gives the normalized
		force of a sarcomere extending at vin, subject to
		specified Vmax."

ForceLength::usage =
		"ForceLength[fl, csa,n] gives the force of a muscle
		with PCSA csa, fiber length fl, and n sarcomeres in
		series."
								
ForceVelocity::usage =
		"ForceVelocity[svel, fast%] gives scaling factor
		for a muscle with fast% fast fibers lengthening at
		a velocity of fvel when n sarcomeres are in series.
		Requires definitions of global variables Fast and Slow,
		the Vmax for fast and slow fibers respecitvely."
			 
NewArch::usage =
		"NewArch[Muscle,ML] returns {aponeurosis angle, fiber
		angle, fiber length} for Muscle at new length ML"

Vf::usage =
		"Vf[MV,ML,Muscle] returns fiber velocity and fiber length
		of Muscle at length ML contracting at MV"

Begin["`Private`"]

SLengthTension[length_]:=Which[
					length < 1.05, 0,
					length < 1.65, length-1.05,
					length < 2.47, -0.2052 + 0.488 length,
					length <=2.81 , 1,
					length <4.33, 0.658 (4.33-length),
					True, 0]
	(*Human theoretical Cutts ('88) J Anat 160:79-88*)
	
SVelocityTension[vin_, Vmax_]:=If[vin <= 0,
					If[-vin <= Vmax, 
						(Vmax+vin)/(Vmax-4vin),
						0],
					   (1.8-0.8(Vmax-vin)/
							(Vmax+7.56vin))]
	(*Shortening curve from Close 1965 (Nature 206:718)
	  Lengthening from Katz 1939 (J Phys Lond 94:45)*)
	  
ForceLength[fl_, csa_,n_]:= 2.5 csa*
					SLengthTension[fl/n]
	(*based on 2.5 kg/cm^2, so units are kg*)
								
ForceVelocity[svel_, fast_]:=
			fast SVelocityTension[svel, Global`Fast] +
			 (1-fast) SVelocityTension[svel, Global`Slow]
	(*forces of fiber subpopulations add linearly*)
	
MuscleArea[Muscle_List]:= Muscle[[1]] Sin[Muscle[[5]]]
	(*conservation of projection area is used to find
		alpha.  this is a reduced form assuming const.
		aponeurosis length*)

NewArch[Muscle_List,ML_]:= Block[{a$ = ArcSin[
			MuscleArea[Muscle]/ML],
			b$ = ArcTan[
				Muscle[[2]] Sin[a$]/(ML -
				Muscle[[2]] Cos[a$])]},
			Return[{a$,b$,Muscle[[2]]/Sin[b$]*
							Sin[a$]}]]
	(*The fundamental model function:
			conservation of area gives alpha
			trig gives beta (arctan[height/length])
			law of sines gives fiber length*)
							
Vf[MV_,ML_,Muscle_List]:= Block[{Angles$2 = NewArch[
							Muscle, ML]},
		Return[{MV Cos[Angles$2[[1]]+Angles$2[[2]]]/
						Cos[Angles$2[[1]]],
				Angles$2[[3]]}]]
	(*from zuurbier and Huijing, with aponeurosis
		velocity = 0*)

End[]
EndPackage[]