Force-length relation

[HeikoStark]

Hi!

In the last time I have worked on a second active contraction model. And I think its ready for posting.

It use the definition from Gordon (1966) and it should be able to use it for different animal/muscle/fibertype tissues.

I have used the framework described in the thread below:
(HOWTO) Insertion of a new fiber/contraction model

Feel free to discuss.

	else if (m_model == 2) // Gordon 1966
	{
		// strain energy derivative
		double Wp = 0;
		if (lamd > 1)
		{
			double lamdi = 1.0/lamd;
			double Wl;
			if (lamd < m_lam1)
			{
				Wl = lamdi*m_c3*(exp(m_c4*(lamd - 1)) - 1);
			}
			else
			{
				double c6 = m_c3*(exp(m_c4*(m_lam1-1))-1) - m_c5*m_lam1;
				Wl = lamdi*(m_c5*lamd + c6);
			}
			Wp  = 0.5*lamdi*Wl;
		}
		else 
		{
			Wp = 0;
		}

		// --- active contraction contribution ---		
		double Wa = 0;
		if (m_lcna >= 0)
		{
	                double Wl;	    
			if (lamd < m_ax) Wl = 0;
			else if (lamd < m_bx) Wl = m_ay + (lamd-m_ax)*((m_by-m_ay)/(m_bx-m_ax));
			else if (lamd < m_cx) Wl = m_by + (lamd-m_bx)*((m_cy-m_by)/(m_cx-m_bx));
			else if (lamd < m_dx) Wl = m_cy + (lamd-m_cx)*((m_dy-m_cy)/(m_dx-m_cx));
			else if (lamd < m_ex) Wl = m_dy + (lamd-m_dx)*((m_ey-m_dy)/(m_ex-m_dx));
			else Wl = 0;
			
			double alpha = m_plc->Value();		

			Wa = alpha * m_smax * Wl;
			
		}
		
		double W4 = Wp + Wa;

		// calculate dyad of a: AxA = (a x a)
		mat3ds AxA = dyad(a);

		// calculate FdWf/dCFt = I4*W4*(a x a)
		mat3ds T = (W4*I4) * AxA;
	
		// calculate stress: 
		s = T.dev()*twoJi;

		//-----------------------------------------	  
	}

[maas]

Hi Heiko,

Thanks for sharing your code. If you have any simple example problems, can you please upload those as well. I think that would be very useful for users who would like to try out your code.

Thanks,

Steve.