C
C There are a total of 10 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 24 entries in the constant variable array.
C
C
C VOI is time in component environment (millisecond).
C STATES(1) is V in component membrane (millivolt).
C CONSTS(1) is Cm in component membrane (femtofarad).
C ALGBRC(5) is Ica in component Ica (femtoampere).
C ALGBRC(7) is Is1 in component Is1 (femtoampere).
C ALGBRC(10) is Is2 in component Is2 (femtoampere).
C ALGBRC(9) is Il in component Il (femtoampere).
C ALGBRC(8) is Ik in component Ik (femtoampere).
C CONSTS(2) is gCa in component Ica (picosiemens).
C CONSTS(3) is VCa in component Ica (millivolt).
C ALGBRC(1) is m_infinity in component m (dimensionless).
C CONSTS(4) is vm in component m (millivolt).
C CONSTS(5) is sm in component m (millivolt).
C CONSTS(6) is gs1 in component Is1 (picosiemens).
C CONSTS(7) is VK in component Ik (millivolt).
C STATES(2) is s1 in component s1 (dimensionless).
C ALGBRC(2) is s1_infinity in component s1 (dimensionless).
C CONSTS(8) is autos1 in component s1 (dimensionless).
C CONSTS(9) is s1knot in component s1 (dimensionless).
C CONSTS(10) is tau_s1 in component s1 (millisecond).
C CONSTS(11) is vs1 in component s1 (millivolt).
C CONSTS(12) is ss1 in component s1 (millivolt).
C CONSTS(13) is gK in component Ik (picosiemens).
C STATES(3) is n in component n (dimensionless).
C ALGBRC(3) is n_infinity in component n (dimensionless).
C CONSTS(14) is tau_n_bar in component n (millisecond).
C ALGBRC(6) is tau_n in component n (millisecond).
C CONSTS(15) is vn in component n (millivolt).
C CONSTS(16) is sn in component n (millivolt).
C CONSTS(17) is gl in component Il (picosiemens).
C CONSTS(18) is Vl in component Il (millivolt).
C CONSTS(19) is gs2 in component Is2 (picosiemens).
C STATES(4) is s2 in component s2 (dimensionless).
C ALGBRC(4) is s2_infinity in component s2 (dimensionless).
C CONSTS(20) is autos2 in component s2 (dimensionless).
C CONSTS(21) is s2knot in component s2 (dimensionless).
C CONSTS(22) is tau_s2 in component s2 (millisecond).
C CONSTS(23) is vs2 in component s2 (millivolt).
C CONSTS(24) is ss2 in component s2 (millivolt).
C RATES(1) is d/dt V in component membrane (millivolt).
C RATES(2) is d/dt s1 in component s1 (dimensionless).
C RATES(3) is d/dt n in component n (dimensionless).
C RATES(4) is d/dt s2 in component s2 (dimensionless).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      STATES(1) = -40.0
      CONSTS(1) = 4525.0
      CONSTS(2) = 280.0
      CONSTS(3) = 100.0
      CONSTS(4) = -22.0
      CONSTS(5) = 7.5
      CONSTS(6) = 20.0
      CONSTS(7) = -80.0
      STATES(2) = 0.9
      CONSTS(8) = 1
      CONSTS(9) = 1
      CONSTS(10) = 1000.0
      CONSTS(11) = -50.0
      CONSTS(12) = 5
      CONSTS(13) = 1300.0
      STATES(3) = 0.0
      CONSTS(14) = 8.25
      CONSTS(15) = -9.0
      CONSTS(16) = 10.0
      CONSTS(17) = 25.0
      CONSTS(18) = -40.0
      CONSTS(19) = 16
      STATES(4) = 0.5
      CONSTS(20) = 1
      CONSTS(21) = 0.49
      CONSTS(22) = 30000.0
      CONSTS(23) = -40.0
      CONSTS(24) = 15
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(2) = 1.00000/(1.00000+EXP((CONSTS(11) - STATES(1))/CONSTS(12)))
      RATES(2) =  CONSTS(8)*((ALGBRC(2) - STATES(2))/CONSTS(10))+ (1.00000 - CONSTS(8))*(CONSTS(9) - STATES(2))
      ALGBRC(4) = 1.00000/(1.00000+EXP((CONSTS(23) - STATES(1))/CONSTS(24)))
      RATES(4) =  CONSTS(20)*((ALGBRC(4) - STATES(4))/CONSTS(22))+ (1.00000 - CONSTS(20))*(CONSTS(21) - STATES(4))
      ALGBRC(3) = 1.00000/(1.00000+EXP((CONSTS(15) - STATES(1))/CONSTS(16)))
      ALGBRC(6) = CONSTS(14)/(1.00000+EXP((STATES(1) - CONSTS(15))/CONSTS(16)))
      RATES(3) = (ALGBRC(3) - STATES(3))/ALGBRC(6)
      ALGBRC(1) = 1.00000/(1.00000+EXP((CONSTS(4) - STATES(1))/CONSTS(5)))
      ALGBRC(5) =  CONSTS(2)*ALGBRC(1)*(STATES(1) - CONSTS(3))
      ALGBRC(7) =  CONSTS(6)*STATES(2)*(STATES(1) - CONSTS(7))
      ALGBRC(10) =  CONSTS(19)*STATES(4)*(STATES(1) - CONSTS(7))
      ALGBRC(9) =  CONSTS(17)*(STATES(1) - CONSTS(18))
      ALGBRC(8) =  CONSTS(13)*STATES(3)*(STATES(1) - CONSTS(7))
      RATES(1) = - (ALGBRC(5)+ALGBRC(7)+ALGBRC(10)+ALGBRC(9)+ALGBRC(8))/CONSTS(1)
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(2) = 1.00000/(1.00000+EXP((CONSTS(11) - STATES(1))/CONSTS(12)))
      ALGBRC(4) = 1.00000/(1.00000+EXP((CONSTS(23) - STATES(1))/CONSTS(24)))
      ALGBRC(3) = 1.00000/(1.00000+EXP((CONSTS(15) - STATES(1))/CONSTS(16)))
      ALGBRC(6) = CONSTS(14)/(1.00000+EXP((STATES(1) - CONSTS(15))/CONSTS(16)))
      ALGBRC(1) = 1.00000/(1.00000+EXP((CONSTS(4) - STATES(1))/CONSTS(5)))
      ALGBRC(5) =  CONSTS(2)*ALGBRC(1)*(STATES(1) - CONSTS(3))
      ALGBRC(7) =  CONSTS(6)*STATES(2)*(STATES(1) - CONSTS(7))
      ALGBRC(10) =  CONSTS(19)*STATES(4)*(STATES(1) - CONSTS(7))
      ALGBRC(9) =  CONSTS(17)*(STATES(1) - CONSTS(18))
      ALGBRC(8) =  CONSTS(13)*STATES(3)*(STATES(1) - CONSTS(7))
      RETURN
      END