/*
 * The Pole-Zero Constitutive Material Law
 * 
 * Model Status
 * 
 * This is the original unchecked version of the model imported
 * from the previous CellML model repository, 24-Jan-2006.
 * 
 * Model Structure
 * 
 * Myocardial tissue consists of layers of interconnected sheets
 * of tissue, separated by cleavage planes. The muscle fibres lie
 * in the plane of a sheet, and adjacent fibres are coupled more
 * strongly in the plane of the sheet than traverse to it. This
 * results in three microstructural axes:
 * 
 * one along the fibre direction called the fibre axis;
 * 
 * one orthogonal to the fibre axis but also in the plane of the
 * sheet called the sheet axis; and
 * 
 * one which is orthogonal to these two, directed across the cleavage
 * planes called the sheet normal.
 * 
 * The most obvious difference in the material behaviour of the
 * myocardial tisue along these three axes is the limiting strain
 * for an elastic response. Another feature of the biaxial tests
 * is that the stress-strain curve along one axis is almost independent
 * of the degree of lateral stretch. To account for these microstructural
 * properties, a strain energy function called the pole-zero law
 * has been proposed (see below). By using a separate pole for
 * each microstructurally defined axis, the different strain limiting
 * behaviour along each axis is accounted for without the need
 * for numerically unstable, large exponents which are needed in
 * exponential or power law expressions.
 * 
 * The original paper reference is cited below:
 * 
 * Computational mechanics of the heart: from tissue structure
 * to ventricular function, M P Nash and P J Hunter, Journal of
 * Elasticity 61(1/3): 113-141, 2000.
 */

import nsrunit;
unit conversion on;
unit strain=1  dimensionless;
unit stress=1  dimensionless;
unit pole=1  dimensionless;
unit curvature=1  dimensionless;
unit scale=1  dimensionless;

math main {
	//Warning:  the following variables were set 'extern' or given
	//  an initial value of '0' because the model would otherwise be
	//  underdetermined:  E11, E22, E33, E12, E13, E23
	extern real E11 strain;
	extern real E22 strain;
	extern real E33 strain;
	extern real E12 strain;
	extern real E13 strain;
	extern real E23 strain;
	real k11 scale;
	k11=0;
	real a11 pole;
	a11=0;
	real b11 curvature;
	b11=0;
	real k22 scale;
	k22=0;
	real a22 pole;
	a22=0;
	real b22 curvature;
	b22=0;
	real k33 scale;
	k33=0;
	real a33 pole;
	a33=0;
	real b33 curvature;
	b33=0;
	real k12 scale;
	k12=0;
	real a12 pole;
	a12=0;
	real b12 curvature;
	b12=0;
	real k13 scale;
	k13=0;
	real a13 pole;
	a13=0;
	real b13 curvature;
	b13=0;
	real k23 scale;
	k23=0;
	real a23 pole;
	a23=0;
	real b23 curvature;
	b23=0;
	real Tdev11 stress;
	real Tdev22 stress;
	real Tdev33 stress;
	real Tdev12 stress;
	real Tdev13 stress;
	real Tdev23 stress;
	real Elim11 strain;
	real Elim22 strain;
	real Elim33 strain;
	real Elim12 strain;
	real Elim13 strain;
	real Elim23 strain;
	real Tol dimensionless;
	Tol=0.9;
	real Sign12 dimensionless;
	real Sign13 dimensionless;
	real Sign23 dimensionless;
	real Eabs12 strain;
	real Eabs13 strain;
	real Eabs23 strain;

	// <component name="interface">

	// <component name="equations">
	Elim11=(Tol*a11);
	Elim22=(Tol*a22);
	Elim33=(Tol*a33);
	Elim12=(Tol*a12);
	Elim13=(Tol*a13);
	Elim23=(Tol*a23);
	Sign12=(if (E12>=0) 1 else -1);
	Sign13=(if (E13>=0) 1 else -1);
	Sign23=(if (E23>=0) 1 else -1);
	Eabs12=(if (E12>=0) E12 else (-1)*E12);
	Eabs13=(if (E13>=0) E13 else (-1)*E13);
	Eabs23=(if (E23>=0) E23 else (-1)*E23);
	Tdev11=(if (E11<=0) 2*k11*E11/a11^b11 else if (E11>Elim11) k11*Elim11/(a11-Elim11)^b11*(2+b11*Elim11/(a11-Elim11))+k11/(a11-Elim11)^b11*(2+4*b11*Elim11/(a11-Elim11)+b11*Elim11^2*(b11+1)/(a11-Elim11)^2)*(E11-Elim11) else k11*E11/(a11-E11)^b11*(2+b11*E11/(a11-E11)));
	Tdev22=(if (E22<=0) 2*k22*E22/a22^b22 else if (E22>Elim22) k22*Elim22/(a22-Elim22)^b22*(2+b22*Elim22/(a22-Elim22))+k22/(a22-Elim22)^b22*(2+4*b22*Elim22/(a22-Elim22)+b22*Elim22^2*(b22+1)/(a22-Elim22)^2)*(E22-Elim22) else k22*E22/(a22-E22)^b22*(2+b22*E22/(a22-E22)));
	Tdev33=(if (E33<=0) 2*k33*E33/a33^b33 else if (E33>Elim33) k33*Elim33/(a33-Elim33)^b33*(2+b33*Elim33/(a33-Elim33))+k33/(a33-Elim33)^b33*(2+4*b33*Elim33/(a33-Elim33)+b33*Elim33^2*(b33+1)/(a33-Elim33)^2)*(E33-Elim33) else k33*E33/(a33-E33)^b33*(2+b33*E33/(a33-E33)));
	Tdev12=(Sign12*(if (Eabs12>Elim12) k12*Elim12/(a12-Elim12)^b12*(2+b12*Elim12/(a12-Elim12))+k12/(a12-Elim12)^b12*(2+4*b12*Elim12/(a12-Elim12)+b12*Elim12^2*(b12+1)/(a12-Elim12)^2)*(Eabs12-Elim12) else k12*Eabs12/(a12-Eabs12)^b12*(2+b12*Eabs12/(a12-Eabs12))));
	Tdev13=(Sign13*(if (Eabs13>Elim13) k13*Elim13/(a13-Elim13)^b13*(2+b13*Elim13/(a13-Elim13))+k13/(a13-Elim13)^b13*(2+4*b13*Elim13/(a13-Elim13)+b13*Elim13^2*(b13+1)/(a13-Elim13)^2)*(Eabs13-Elim13) else k13*Eabs13/(a13-Eabs13)^b13*(2+b13*Eabs13/(a13-Eabs13))));
	Tdev23=(Sign23*(if (Eabs23>Elim23) k23*Elim23/(a23-Elim23)^b23*(2+b23*Elim23/(a23-Elim23))+k23/(a23-Elim23)^b23*(2+4*b23*Elim23/(a23-Elim23)+b23*Elim23^2*(b23+1)/(a23-Elim23)^2)*(Eabs23-Elim23) else k23*Eabs23/(a23-Eabs23)^b23*(2+b23*Eabs23/(a23-Eabs23))));
}