C
C There are a total of 4 entries in the algebraic variable array.
C There are a total of 3 entries in each of the rate and state variable arrays.
C There are a total of 19 entries in the constant variable array.
C
C
C VOI is time in component environment (hour).
C ALGBRC(1) is Qv in component v (per_second).
C STATES(1) is Vv in component v (mV).
C CONSTS(1) is tau_v in component v (second).
C CONSTS(2) is v_vm in component v (mV_second).
C CONSTS(3) is awake_start in component v (hour).
C CONSTS(4) is awake_end in component v (hour).
C ALGBRC(2) is Qm in component m (per_second).
C CONSTS(5) is Qmax in component model_parameters (per_second).
C ALGBRC(4) is D in component D (mV).
C CONSTS(6) is theta in component model_parameters (mV).
C CONSTS(7) is sigma in component model_parameters (mV).
C STATES(2) is Vm in component m (mV).
C CONSTS(8) is tau_m in component m (second).
C CONSTS(9) is v_mv in component m (mV_second).
C CONSTS(10) is A in component m (mV).
C CONSTS(11) is awake_start in component m (hour).
C CONSTS(12) is awake_end in component m (hour).
C STATES(3) is H in component H (nM).
C CONSTS(13) is chi in component H (hour).
C CONSTS(14) is mu in component H (nM_second).
C ALGBRC(3) is C in component D (dimensionless).
C CONSTS(15) is c0 in component D (dimensionless).
C CONSTS(19) is omega in component D (per_hour).
C CONSTS(16) is alpha in component D (hour).
C CONSTS(17) is v_vc in component D (mV).
C CONSTS(18) is v_vh in component D (mV_per_nM).
C RATES(1) is d/dt Vv in component v (mV).
C RATES(2) is d/dt Vm in component m (mV).
C RATES(3) is d/dt H in component H (nM).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      STATES(1) = -12.5
      CONSTS(1) = 10.0
      CONSTS(2) = -2.1
      CONSTS(3) = 167.04
      CONSTS(4) = 240.0
      CONSTS(5) = 100.0
      CONSTS(6) = 10.0
      CONSTS(7) = 3.0
      STATES(2) = 0.0
      CONSTS(8) = 10.0
      CONSTS(9) = -1.8
      CONSTS(10) = 1.3
      CONSTS(11) = 160.0
      CONSTS(12) = 240.0
      STATES(3) = 15.0
      CONSTS(13) = 45.0
      CONSTS(14) = 4.4
      CONSTS(15) = 4.5
      CONSTS(16) = 0
      CONSTS(17) = -2.9
      CONSTS(18) = 1.0
      CONSTS(19) = ( 2.00000* 3.14159265358979)/24.0000
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(1) = CONSTS(5)/(1.00000+EXP(- (STATES(1) - CONSTS(6))/CONSTS(7)))
      RATES(2) = TERNRY(VOI.GE.CONSTS(11).AND.VOI.LT.CONSTS(12), 0.00000, ((CONSTS(10)+ CONSTS(9)*ALGBRC(1)) - STATES(2))/(CONSTS(8)/3600.00))
      ALGBRC(2) = CONSTS(5)/(1.00000+EXP(- (STATES(2) - CONSTS(6))/CONSTS(7)))
      RATES(3) = ( CONSTS(14)*ALGBRC(2) - STATES(3))/CONSTS(13)
      ALGBRC(3) = CONSTS(15)+cos( CONSTS(19)*(VOI - CONSTS(16)))
      ALGBRC(4) =  CONSTS(17)*ALGBRC(3)+ CONSTS(18)*STATES(3)
      RATES(1) = TERNRY(VOI.GE.CONSTS(3).AND.VOI.LT.CONSTS(4), 0.00000, (( CONSTS(2)*ALGBRC(2)+ALGBRC(4)) - STATES(1))/(CONSTS(1)/3600.00))
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(1) = CONSTS(5)/(1.00000+EXP(- (STATES(1) - CONSTS(6))/CONSTS(7)))
      ALGBRC(2) = CONSTS(5)/(1.00000+EXP(- (STATES(2) - CONSTS(6))/CONSTS(7)))
      ALGBRC(3) = CONSTS(15)+cos( CONSTS(19)*(VOI - CONSTS(16)))
      ALGBRC(4) =  CONSTS(17)*ALGBRC(3)+ CONSTS(18)*STATES(3)
      RETURN
      END
      REAL FUNCTION TERNRY(TEST, VALA, VALB)
      LOGICAL TEST
      REAL VALA, VALB
      IF (TEST) THEN
        TERNRY = VALA
      ELSE
        TERNRY = VALB
      ENDIF
      RETURN
      END