C
C There are a total of 29 entries in the algebraic variable array.
C There are a total of 9 entries in each of the rate and state variable arrays.
C There are a total of 48 entries in the constant variable array.
C
C
C VOI is time in component environment (millisecond).
C CONSTS(1) is isotonic_mode in component parameters (dimensionless).
C CONSTS(2) is alpha_1 in component parameters (per_um).
C CONSTS(3) is beta_1 in component parameters (mN).
C CONSTS(4) is alpha_2 in component parameters (per_um).
C CONSTS(5) is beta_2 in component parameters (mN).
C CONSTS(6) is alpha_3 in component parameters (per_um).
C CONSTS(7) is beta_3 in component parameters (mN).
C CONSTS(8) is lambda in component parameters (mN).
C CONSTS(9) is A_half in component parameters (dimensionless).
C CONSTS(10) is mu in component parameters (dimensionless).
C CONSTS(11) is chi in component parameters (dimensionless).
C CONSTS(12) is chi_0 in component parameters (dimensionless).
C CONSTS(13) is m_0 in component parameters (dimensionless).
C CONSTS(14) is v_max in component parameters (um_per_msec).
C CONSTS(15) is a in component parameters (dimensionless).
C CONSTS(16) is d_h in component parameters (dimensionless).
C CONSTS(17) is alpha_P in component parameters (dimensionless).
C CONSTS(18) is alpha_PG in component parameters (dimensionless).
C CONSTS(19) is S_0 in component parameters_izakov_et_al_1991 (um).
C ALGBRC(4) is q_v in component parameters_izakov_et_al_1991 (per_millisecond).
C CONSTS(20) is q_1 in component parameters_izakov_et_al_1991 (per_millisecond).
C CONSTS(21) is q_2 in component parameters_izakov_et_al_1991 (per_millisecond).
C CONSTS(22) is q_3 in component parameters_izakov_et_al_1991 (per_millisecond).
C CONSTS(46) is v_1 in component parameters_izakov_et_al_1991 (um_per_msec).
C CONSTS(23) is alpha_P in component parameters_izakov_et_al_1991 (per_um).
C CONSTS(24) is alpha_S in component parameters_izakov_et_al_1991 (per_um).
C CONSTS(25) is alpha_G in component parameters_izakov_et_al_1991 (dimensionless).
C CONSTS(26) is a_on in component parameters_izakov_et_al_1991 (per_millisecond).
C CONSTS(27) is a_off in component parameters_izakov_et_al_1991 (per_millisecond).
C CONSTS(28) is k_A in component parameters_izakov_et_al_1991 (dimensionless).
C STATES(1) is v in component CE_velocity (um_per_msec).
C ALGBRC(27) is F_CE in component force (mN).
C ALGBRC(6) is F_muscle in component force (mN).
C ALGBRC(1) is F_XSE in component force (mN).
C ALGBRC(2) is F_SE in component force (mN).
C ALGBRC(3) is F_PE in component force (mN).
C STATES(2) is N in component crossbridge_kinetics (dimensionless).
C ALGBRC(14) is k_P_vis in component CE_velocity (mN).
C ALGBRC(16) is k_S_vis in component PE_velocity (mN).
C STATES(3) is w in component PE_velocity (um_per_msec).
C STATES(4) is l_1 in component length (um).
C STATES(5) is l_2 in component length (um).
C STATES(6) is l_3 in component length (um).
C ALGBRC(26) is p_v in component average_crossbridge_force (dimensionless).
C ALGBRC(23) is K_chi in component crossbridge_kinetics (per_millisecond).
C ALGBRC(7) is M_A in component crossbridge_kinetics (dimensionless).
C ALGBRC(8) is n_1 in component crossbridge_kinetics (dimensionless).
C ALGBRC(9) is L_oz in component crossbridge_kinetics (dimensionless).
C ALGBRC(21) is k_p_v in component crossbridge_kinetics (per_millisecond).
C ALGBRC(22) is k_m_v in component crossbridge_kinetics (per_millisecond).
C STATES(7) is A in component calcium_handling (dimensionless).
C ALGBRC(20) is G_star in component average_crossbridge_force (dimensionless).
C ALGBRC(10) is dl_1_dt in component length (um_per_msec).
C ALGBRC(12) is dl_2_dt in component length (um_per_msec).
C ALGBRC(5) is dl_3_dt in component length (um_per_msec).
C ALGBRC(28) is phi_chi in component CE_velocity (um_per_msec2).
C ALGBRC(29) is p_prime_v in component average_crossbridge_force (per_millisecond).
C CONSTS(29) is alpha_P_lengthening in component CE_velocity (per_um).
C CONSTS(30) is beta_P_lengthening in component CE_velocity (mN).
C CONSTS(31) is alpha_P_shortening in component CE_velocity (per_um).
C CONSTS(32) is beta_P_shortening in component CE_velocity (mN).
C CONSTS(33) is alpha_S_lengthening in component PE_velocity (per_um).
C CONSTS(34) is beta_S_lengthening in component PE_velocity (mN).
C CONSTS(35) is alpha_S_shortening in component PE_velocity (per_um).
C CONSTS(36) is beta_S_shortening in component PE_velocity (mN).
C ALGBRC(19) is P_star in component average_crossbridge_force (dimensionless).
C ALGBRC(18) is gamma in component average_crossbridge_force (dimensionless).
C CONSTS(47) is case_1 in component average_crossbridge_force (per_millisecond).
C ALGBRC(24) is case_2 in component average_crossbridge_force (per_millisecond).
C CONSTS(48) is case_3 in component average_crossbridge_force (per_millisecond).
C ALGBRC(25) is case_4 in component average_crossbridge_force (per_millisecond).
C ALGBRC(15) is dA_dt in component calcium_handling (per_millisecond).
C ALGBRC(11) is N_A in component calcium_handling (dimensionless).
C ALGBRC(13) is pi_N_A in component calcium_handling (dimensionless).
C STATES(8) is B in component calcium_handling (dimensionless).
C ALGBRC(17) is dB_dt in component calcium_handling (per_millisecond).
C STATES(9) is Ca_C in component calcium_handling (dimensionless).
C CONSTS(37) is A_tot in component calcium_handling (dimensionless).
C CONSTS(38) is B_tot in component calcium_handling (dimensionless).
C CONSTS(39) is b_on in component calcium_handling (per_millisecond).
C CONSTS(40) is b_off in component calcium_handling (per_millisecond).
C CONSTS(41) is a_c in component calcium_handling (per_millisecond2).
C CONSTS(42) is b_c in component calcium_handling (per_millisecond2).
C CONSTS(43) is r_Ca in component calcium_handling (per_millisecond).
C CONSTS(44) is q_Ca in component calcium_handling (dimensionless).
C CONSTS(45) is t_d in component calcium_handling (millisecond).
C RATES(2) is d/dt N in component crossbridge_kinetics (dimensionless).
C RATES(4) is d/dt l_1 in component length (um).
C RATES(5) is d/dt l_2 in component length (um).
C RATES(6) is d/dt l_3 in component length (um).
C RATES(1) is d/dt v in component CE_velocity (um_per_msec).
C RATES(3) is d/dt w in component PE_velocity (um_per_msec).
C RATES(7) is d/dt A in component calcium_handling (dimensionless).
C RATES(8) is d/dt B in component calcium_handling (dimensionless).
C RATES(9) is d/dt Ca_C in component calcium_handling (dimensionless).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      CONSTS(1) = 1
      CONSTS(2) = 19
      CONSTS(3) = 2.9
      CONSTS(4) = 14.6
      CONSTS(5) = 0.01
      CONSTS(6) = 48
      CONSTS(7) = 0.1
      CONSTS(8) = 960
      CONSTS(9) = 0.6
      CONSTS(10) = 3
      CONSTS(11) = 0.705
      CONSTS(12) = 3
      CONSTS(13) = 0.9
      CONSTS(14) = 0.0056
      CONSTS(15) = 0.25
      CONSTS(16) = 0.5
      CONSTS(17) = 4
      CONSTS(18) = 1
      CONSTS(19) = 0.77
      CONSTS(20) = 0.017
      CONSTS(21) = 0.26
      CONSTS(22) = 0.03
      CONSTS(23) = 4
      CONSTS(24) = 4
      CONSTS(25) = 4
      CONSTS(26) = 2.9e-2
      CONSTS(27) = 0.2
      CONSTS(28) = 2.8
      STATES(1) = 0
      STATES(2) = 1
      STATES(3) = 0
      STATES(4) = 1
      STATES(5) = 1
      STATES(6) = 1
      STATES(7) = 0
      CONSTS(29) = 16
      CONSTS(30) = 15
      CONSTS(31) = 16
      CONSTS(32) = 15
      CONSTS(33) = 39
      CONSTS(34) = 80
      CONSTS(35) = 46
      CONSTS(36) = 60
      STATES(8) = 0
      STATES(9) = 0
      CONSTS(37) = 1
      CONSTS(38) = 0.4
      CONSTS(39) = 2.6
      CONSTS(40) = 0.182
      CONSTS(41) = 0.005
      CONSTS(42) = 0.03
      CONSTS(43) = 0.65
      CONSTS(44) = 50
      CONSTS(45) = 33
      CONSTS(46) = CONSTS(14)/10.0000
      CONSTS(47) = ( CONSTS(15)*(0.400000+ 0.400000*CONSTS(15)))/( CONSTS(14)* (CONSTS(15)+1.00000)*0.400000 ** 2.00000)
      CONSTS(48) = ( 0.400000*CONSTS(15)+1.00000)/( CONSTS(15)*CONSTS(14))
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(5) = TERNRY(CONSTS(1).EQ.1.00000, 0.00000, - STATES(3))
      RATES(6) = ALGBRC(5)
      ALGBRC(10) = STATES(1)
      RATES(4) = ALGBRC(10)
      ALGBRC(12) = STATES(3)
      RATES(5) = ALGBRC(12)
      ALGBRC(9) = (STATES(4)+CONSTS(19))/(0.460000+CONSTS(19))
      ALGBRC(11) = STATES(2)/( ALGBRC(9)*STATES(7))
      ALGBRC(13) = TERNRY(ALGBRC(11).GE.1.00000, 1.00000, 0.0200000 ** ALGBRC(11))
      ALGBRC(15) =  CONSTS(26)*(CONSTS(37) - STATES(7))*STATES(9) -  CONSTS(27)*EXP( - CONSTS(28)*STATES(7))*ALGBRC(13)*STATES(7)
      RATES(7) = ALGBRC(15)
      ALGBRC(17) =  CONSTS(39)*(CONSTS(38) - STATES(8))*STATES(9) -  CONSTS(40)*STATES(8)
      RATES(8) = ALGBRC(17)
      RATES(9) = TERNRY(VOI.LT.CONSTS(45),  CONSTS(42)*VOI*(1.00000 - EXP( - CONSTS(41)*VOI ** 2.00000))*EXP( - CONSTS(41)*VOI ** 2.00000), (- ALGBRC(15) - ALGBRC(17)) -  CONSTS(43)*EXP( - CONSTS(44)*STATES(9))*STATES(9))
      ALGBRC(7) = STATES(7) ** CONSTS(10)/(STATES(7) ** CONSTS(10)+CONSTS(9) ** CONSTS(10))
      ALGBRC(8) =  0.600000*STATES(4)+0.500000
      ALGBRC(4) = TERNRY(STATES(1).LE.0.00000, CONSTS(20) - ( CONSTS(21)*STATES(1))/CONSTS(14), CONSTS(22))
      ALGBRC(18) = ( CONSTS(15)*CONSTS(16)*CONSTS(46)/CONSTS(14) ** 2.00000)/( 3.00000*CONSTS(15)*CONSTS(16) - ( (CONSTS(15)+1.00000)*STATES(1))/CONSTS(14))
      ALGBRC(19) = TERNRY(STATES(1).LE.0.00000, ( CONSTS(15)*(1.00000+STATES(1)/CONSTS(14)))/(CONSTS(15) - STATES(1)/CONSTS(14)), (1.00000+CONSTS(16)) - ( CONSTS(16) ** 2.00000*CONSTS(15))/( (( CONSTS(15)*CONSTS(16))/ALGBRC(18))*STATES(1)/CONSTS(14) ** 2.00000+( (CONSTS(15)+1.00000)*STATES(1))/CONSTS(14)+ CONSTS(15)*CONSTS(16)))
      ALGBRC(20) = TERNRY(- CONSTS(14).LE.STATES(1).AND.STATES(1).LE.0.00000, 1.00000+( 0.600000*STATES(1))/CONSTS(14), TERNRY(0.00000.LT.STATES(1).AND.STATES(1).LE.CONSTS(46), ALGBRC(19)/(( (( 0.400000*CONSTS(15)+1.00000)/CONSTS(15))*STATES(1))/CONSTS(14)+1.00000), ( ALGBRC(19)*EXP( - CONSTS(25)*(STATES(1) - CONSTS(46))/CONSTS(14) ** CONSTS(18)))/(( (( 0.400000*CONSTS(15)+1.00000)/CONSTS(15))*STATES(1))/CONSTS(14)+1.00000))
      ALGBRC(21) =  CONSTS(11)*CONSTS(12)*ALGBRC(4)*CONSTS(13)*ALGBRC(20)
      ALGBRC(22) =  CONSTS(12)*ALGBRC(4)*(1.00000 -  CONSTS(11)*CONSTS(13)*ALGBRC(20))
      ALGBRC(23) =  ALGBRC(21)*ALGBRC(7)*ALGBRC(8)*ALGBRC(9)*(1.00000 - STATES(2)) -  ALGBRC(22)*STATES(2)
      RATES(2) = ALGBRC(23)
      ALGBRC(14) = TERNRY(ALGBRC(10).LE.0.00000,  CONSTS(30)*EXP( CONSTS(29)*STATES(4)),  CONSTS(32)*EXP( CONSTS(31)*STATES(4)))
      ALGBRC(26) = ALGBRC(19)/ALGBRC(20)
      ALGBRC(24) = ( CONSTS(15)*1.00000*(1.00000+ 0.400000*CONSTS(15)+( 1.20000*STATES(1))/CONSTS(14)+ 0.600000*STATES(1)/CONSTS(14) ** 2.00000))/( CONSTS(14)* (CONSTS(15) - STATES(1)/CONSTS(14))*(1.00000+( 0.600000*STATES(1))/CONSTS(14)) ** 2.00000)
      ALGBRC(25) =  (1.00000/CONSTS(14))*EXP( - CONSTS(25)*STATES(1)/CONSTS(14) - CONSTS(46)/CONSTS(14) ** CONSTS(18))*(( 0.400000*CONSTS(15)+1.00000)/CONSTS(15)+ CONSTS(25)*CONSTS(17)*(1.00000+( ( 0.400000*CONSTS(15)+1.00000)*STATES(1))/( CONSTS(15)*CONSTS(14))))*STATES(1)/CONSTS(14) - CONSTS(46)/CONSTS(14) ** CONSTS(18) - 1.00000
CALL minimize(minfunc_0, CONSTS, VARIABLES, <VAR>)
      RATES(1) = ALGBRC(28)
      ALGBRC(16) = TERNRY(ALGBRC(12).LE.ALGBRC(10),  CONSTS(34)*EXP( CONSTS(33)*(STATES(5) - STATES(4))),  CONSTS(36)*EXP( CONSTS(35)*(STATES(5) - STATES(4))))
      RATES(3) = TERNRY(CONSTS(1).EQ.1.00000, (( ALGBRC(16)*(ALGBRC(28) -  CONSTS(24)*STATES(3) - STATES(1) ** 2.00000) -  CONSTS(2)*CONSTS(3)*EXP( CONSTS(2)*(STATES(5) - STATES(4)))*(STATES(3) - STATES(1))*1.00000) -  CONSTS(4)*CONSTS(5)*EXP( CONSTS(4)*STATES(5))*STATES(3)*1.00000)/ALGBRC(16), (( ALGBRC(16)*(ALGBRC(28) -  CONSTS(24)*STATES(3) - STATES(1) ** 2.00000) -  CONSTS(2)*CONSTS(3)*EXP( CONSTS(2)*(STATES(5) - STATES(4)))*(STATES(3) - STATES(1))*1.00000) -  ( CONSTS(4)*CONSTS(5)*EXP( CONSTS(4)*STATES(5))+ CONSTS(6)*CONSTS(7)*EXP( CONSTS(6)*STATES(6)))*STATES(3)*1.00000)/ALGBRC(16))
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(5) = TERNRY(CONSTS(1).EQ.1.00000, 0.00000, - STATES(3))
      ALGBRC(10) = STATES(1)
      ALGBRC(12) = STATES(3)
      ALGBRC(9) = (STATES(4)+CONSTS(19))/(0.460000+CONSTS(19))
      ALGBRC(11) = STATES(2)/( ALGBRC(9)*STATES(7))
      ALGBRC(13) = TERNRY(ALGBRC(11).GE.1.00000, 1.00000, 0.0200000 ** ALGBRC(11))
      ALGBRC(15) =  CONSTS(26)*(CONSTS(37) - STATES(7))*STATES(9) -  CONSTS(27)*EXP( - CONSTS(28)*STATES(7))*ALGBRC(13)*STATES(7)
      ALGBRC(17) =  CONSTS(39)*(CONSTS(38) - STATES(8))*STATES(9) -  CONSTS(40)*STATES(8)
      ALGBRC(7) = STATES(7) ** CONSTS(10)/(STATES(7) ** CONSTS(10)+CONSTS(9) ** CONSTS(10))
      ALGBRC(8) =  0.600000*STATES(4)+0.500000
      ALGBRC(4) = TERNRY(STATES(1).LE.0.00000, CONSTS(20) - ( CONSTS(21)*STATES(1))/CONSTS(14), CONSTS(22))
      ALGBRC(18) = ( CONSTS(15)*CONSTS(16)*CONSTS(46)/CONSTS(14) ** 2.00000)/( 3.00000*CONSTS(15)*CONSTS(16) - ( (CONSTS(15)+1.00000)*STATES(1))/CONSTS(14))
      ALGBRC(19) = TERNRY(STATES(1).LE.0.00000, ( CONSTS(15)*(1.00000+STATES(1)/CONSTS(14)))/(CONSTS(15) - STATES(1)/CONSTS(14)), (1.00000+CONSTS(16)) - ( CONSTS(16) ** 2.00000*CONSTS(15))/( (( CONSTS(15)*CONSTS(16))/ALGBRC(18))*STATES(1)/CONSTS(14) ** 2.00000+( (CONSTS(15)+1.00000)*STATES(1))/CONSTS(14)+ CONSTS(15)*CONSTS(16)))
      ALGBRC(20) = TERNRY(- CONSTS(14).LE.STATES(1).AND.STATES(1).LE.0.00000, 1.00000+( 0.600000*STATES(1))/CONSTS(14), TERNRY(0.00000.LT.STATES(1).AND.STATES(1).LE.CONSTS(46), ALGBRC(19)/(( (( 0.400000*CONSTS(15)+1.00000)/CONSTS(15))*STATES(1))/CONSTS(14)+1.00000), ( ALGBRC(19)*EXP( - CONSTS(25)*(STATES(1) - CONSTS(46))/CONSTS(14) ** CONSTS(18)))/(( (( 0.400000*CONSTS(15)+1.00000)/CONSTS(15))*STATES(1))/CONSTS(14)+1.00000))
      ALGBRC(21) =  CONSTS(11)*CONSTS(12)*ALGBRC(4)*CONSTS(13)*ALGBRC(20)
      ALGBRC(22) =  CONSTS(12)*ALGBRC(4)*(1.00000 -  CONSTS(11)*CONSTS(13)*ALGBRC(20))
      ALGBRC(23) =  ALGBRC(21)*ALGBRC(7)*ALGBRC(8)*ALGBRC(9)*(1.00000 - STATES(2)) -  ALGBRC(22)*STATES(2)
      ALGBRC(14) = TERNRY(ALGBRC(10).LE.0.00000,  CONSTS(30)*EXP( CONSTS(29)*STATES(4)),  CONSTS(32)*EXP( CONSTS(31)*STATES(4)))
      ALGBRC(26) = ALGBRC(19)/ALGBRC(20)
      ALGBRC(24) = ( CONSTS(15)*1.00000*(1.00000+ 0.400000*CONSTS(15)+( 1.20000*STATES(1))/CONSTS(14)+ 0.600000*STATES(1)/CONSTS(14) ** 2.00000))/( CONSTS(14)* (CONSTS(15) - STATES(1)/CONSTS(14))*(1.00000+( 0.600000*STATES(1))/CONSTS(14)) ** 2.00000)
      ALGBRC(25) =  (1.00000/CONSTS(14))*EXP( - CONSTS(25)*STATES(1)/CONSTS(14) - CONSTS(46)/CONSTS(14) ** CONSTS(18))*(( 0.400000*CONSTS(15)+1.00000)/CONSTS(15)+ CONSTS(25)*CONSTS(17)*(1.00000+( ( 0.400000*CONSTS(15)+1.00000)*STATES(1))/( CONSTS(15)*CONSTS(14))))*STATES(1)/CONSTS(14) - CONSTS(46)/CONSTS(14) ** CONSTS(18) - 1.00000
      ALGBRC(16) = TERNRY(ALGBRC(12).LE.ALGBRC(10),  CONSTS(34)*EXP( CONSTS(33)*(STATES(5) - STATES(4))),  CONSTS(36)*EXP( CONSTS(35)*(STATES(5) - STATES(4))))
      ALGBRC(1) =  CONSTS(7)*(EXP( CONSTS(6)*STATES(6)) - 1.00000)
      ALGBRC(2) =  CONSTS(3)*(EXP( CONSTS(2)*(STATES(5) - STATES(4))) - 1.00000)
      ALGBRC(3) =  CONSTS(5)*(EXP( CONSTS(4)*STATES(5)) - 1.00000)
      ALGBRC(6) = ALGBRC(1)
      ALGBRC(27) =  CONSTS(8)*ALGBRC(26)*STATES(2)
      RETURN
      END
REAL FUNCTION minfunc_0(CONSTS, VARIABLES)
            REAL CONSTS(*), VARIABLES(*) 
      minfunc_0[1] = abs(ALGBRC(28) - (TERNRY(CONSTS(1).EQ.1.00000, ( CONSTS(8)*ALGBRC(23)*ALGBRC(26)*1.00000+ CONSTS(23)*ALGBRC(14)*STATES(1) ** 2.00000+ CONSTS(4)*CONSTS(5)*EXP( CONSTS(4)*STATES(5))*STATES(3)*1.00000)/( CONSTS(8)*STATES(2)*ALGBRC(29)*1.00000+ALGBRC(14)), ( CONSTS(8)*ALGBRC(23)*ALGBRC(26)*1.00000+ CONSTS(23)*ALGBRC(14)*STATES(1) ** 2.00000+ ( CONSTS(4)*CONSTS(5)*EXP( CONSTS(4)*STATES(5))+ CONSTS(6)*CONSTS(7)*EXP( CONSTS(6)*STATES(6)))*STATES(3)*1.00000)/( CONSTS(8)*STATES(2)*ALGBRC(29)*1.00000+ALGBRC(14)))))
      minfunc_0[2] = abs(ALGBRC(29) - (TERNRY(STATES(1).LE.- CONSTS(14),  CONSTS(47)*ALGBRC(28), TERNRY(- CONSTS(14).LT.STATES(1).AND.STATES(1).LE.0.00000,  ALGBRC(24)*ALGBRC(28), TERNRY(0.00000.LT.STATES(1).AND.STATES(1).LE.CONSTS(46),  CONSTS(48)*ALGBRC(28),  ALGBRC(25)*ALGBRC(28))))
      RETURN
      END
      REAL FUNCTION TERNRY(TEST, VALA, VALB)
      LOGICAL TEST
      REAL VALA, VALB
      IF (TEST) THEN
        TERNRY = VALA
      ELSE
        TERNRY = VALB
      ENDIF
      RETURN
      END