C
C There are a total of 6 entries in the algebraic variable array.
C There are a total of 4 entries in each of the rate and state variable arrays.
C There are a total of 20 entries in the constant variable array.
C
C
C VOI is time in component environment (ms).
C ALGBRC(4) is Ca_i in component intracellular_ion_concentrations (uM).
C ALGBRC(1) is mtime in component intracellular_ion_concentrations (dimensionless).
C STATES(1) is xb in component crossbridges (dimensionless).
C STATES(2) is TRPN in component troponin (dimensionless).
C CONSTS(1) is k_xb in component crossbridges (per_ms).
C CONSTS(2) is nperm in component crossbridges (dimensionless).
C CONSTS(3) is perm50 in component crossbridges (dimensionless).
C ALGBRC(2) is permtot in component crossbridges (dimensionless).
C CONSTS(4) is Ca_50ref in component troponin (uM).
C CONSTS(5) is beta_1 in component troponin (dimensionless).
C CONSTS(6) is k_off in component troponin (per_ms).
C CONSTS(7) is n_TRPN in component troponin (dimensionless).
C CONSTS(17) is lambda_m in component filament_overlap (dimensionless).
C CONSTS(18) is Ca_50 in component troponin (uM).
C CONSTS(15) is lambda in component Myofilaments (dimensionless).
C CONSTS(16) is dlambdadt in component Myofilaments (per_ms).
C CONSTS(20) is overlap in component filament_overlap (dimensionless).
C CONSTS(8) is beta_0 in component filament_overlap (dimensionless).
C CONSTS(19) is lambda_s in component filament_overlap (dimensionless).
C CONSTS(9) is T_ref in component isometric_tension (kPa).
C ALGBRC(3) is T_0 in component isometric_tension (kPa).
C ALGBRC(5) is Q in component dynamic_stiffness (dimensionless).
C CONSTS(10) is a in component dynamic_stiffness (dimensionless).
C STATES(3) is Q_1 in component dynamic_stiffness (dimensionless).
C STATES(4) is Q_2 in component dynamic_stiffness (dimensionless).
C CONSTS(11) is A_1 in component dynamic_stiffness (dimensionless).
C CONSTS(12) is A_2 in component dynamic_stiffness (dimensionless).
C CONSTS(13) is alpha_1 in component dynamic_stiffness (per_ms).
C CONSTS(14) is alpha_2 in component dynamic_stiffness (per_ms).
C ALGBRC(6) is Tension in component dynamic_stiffness (kPa).
C RATES(1) is d/dt xb in component crossbridges (dimensionless).
C RATES(2) is d/dt TRPN in component troponin (dimensionless).
C RATES(3) is d/dt Q_1 in component dynamic_stiffness (dimensionless).
C RATES(4) is d/dt Q_2 in component dynamic_stiffness (dimensionless).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      STATES(1) = 0.00046
      STATES(2) = 0.0752
      CONSTS(1) = 0.1
      CONSTS(2) = 5
      CONSTS(3) = 0.35
      CONSTS(4) = 0.8
      CONSTS(5) = -1.5
      CONSTS(6) = 0.1
      CONSTS(7) = 2
      CONSTS(8) = 1.65
      CONSTS(9) = 120
      CONSTS(10) = 0.35
      STATES(3) = 0
      STATES(4) = 0
      CONSTS(11) = -29
      CONSTS(12) = 116
      CONSTS(13) = 0.1
      CONSTS(14) = 0.5
      CONSTS(15) = 1.00000
      CONSTS(16) = 0.00000
      CONSTS(17) = TERNRY(CONSTS(15).GT.1.20000, 1.20000, CONSTS(15))
      CONSTS(18) =  CONSTS(4)*(1.00000+ CONSTS(5)*(CONSTS(17) - 1.00000))
      CONSTS(19) = TERNRY(CONSTS(17).GE.0.870000, 0.870000, CONSTS(17))
      CONSTS(20) = 1.00000+ CONSTS(8)*((CONSTS(17)+CONSTS(19)) - 1.87000)
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      RATES(3) =  CONSTS(11)*CONSTS(16) -  CONSTS(13)*STATES(3)
      RATES(4) =  CONSTS(12)*CONSTS(16) -  CONSTS(14)*STATES(4)
      ALGBRC(2) =  STATES(2)/CONSTS(3) ** CONSTS(2) ** (1.0 / 2)
      RATES(1) =  CONSTS(1)*( ALGBRC(2)*(1.00000 - STATES(1)) -  (1.00000/ALGBRC(2))*STATES(1))
      ALGBRC(1) = (VOI -  167.000*INT(VOI/167.000))/1.00000
      ALGBRC(4) = TERNRY(ALGBRC(1).GE.1.17000.AND.ALGBRC(1).LT.30.8400,  1.00000*1.85358e-05*ALGBRC(1) ** 3.00000+ - 0.00159034*ALGBRC(1) ** 2.00000+ 0.0436459*ALGBRC(1) ** 1.00000+0.167079, TERNRY(ALGBRC(1).GE.30.8400, (( 1.00000*- 5.74585e-08*ALGBRC(1) ** 3.00000+ 3.11222e-05*ALGBRC(1) ** 2.00000) -  0.00661849*ALGBRC(1) ** 1.00000)+0.720442, 0.216000)
      RATES(2) =  CONSTS(6)*( ALGBRC(4)/CONSTS(18) ** CONSTS(7)*(1.00000 - STATES(2)) - STATES(2))
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(2) =  STATES(2)/CONSTS(3) ** CONSTS(2) ** (1.0 / 2)
      ALGBRC(1) = (VOI -  167.000*INT(VOI/167.000))/1.00000
      ALGBRC(4) = TERNRY(ALGBRC(1).GE.1.17000.AND.ALGBRC(1).LT.30.8400,  1.00000*1.85358e-05*ALGBRC(1) ** 3.00000+ - 0.00159034*ALGBRC(1) ** 2.00000+ 0.0436459*ALGBRC(1) ** 1.00000+0.167079, TERNRY(ALGBRC(1).GE.30.8400, (( 1.00000*- 5.74585e-08*ALGBRC(1) ** 3.00000+ 3.11222e-05*ALGBRC(1) ** 2.00000) -  0.00661849*ALGBRC(1) ** 1.00000)+0.720442, 0.216000)
      ALGBRC(3) =  CONSTS(9)*STATES(1)*CONSTS(20)
      ALGBRC(5) = STATES(3)+STATES(4)
      ALGBRC(6) = TERNRY(ALGBRC(5).LT.0.00000, ( ALGBRC(3)*( CONSTS(10)*ALGBRC(5)+1.00000))/(1.00000 - ALGBRC(5)), ( ALGBRC(3)*(1.00000+ (CONSTS(10)+2.00000)*ALGBRC(5)))/(1.00000+ALGBRC(5)))
      RETURN
      END
      REAL FUNCTION TERNRY(TEST, VALA, VALB)
      LOGICAL TEST
      REAL VALA, VALB
      IF (TEST) THEN
        TERNRY = VALA
      ELSE
        TERNRY = VALB
      ENDIF
      RETURN
      END