C
C There are a total of 16 entries in the algebraic variable array.
C There are a total of 8 entries in each of the rate and state variable arrays.
C There are a total of 37 entries in the constant variable array.
C
C
C VOI is Time in component Environmental (second).
C CONSTS(1) is S in component Parameters (hertz).
C CONSTS(2) is Zeta in component Parameters (dimensionless).
C CONSTS(3) is R0 in component Parameters (mm).
C CONSTS(4) is Cwall in component Parameters (mm_per_mmHg).
C CONSTS(5) is Kne in component Parameters (mmHg_per_mm).
C CONSTS(6) is K1 in component Parameters (mmHg_per_mm).
C CONSTS(7) is K2 in component Parameters (mmHg_per_mm).
C CONSTS(8) is K3 in component Parameters (mmHg_per_mm).
C CONSTS(9) is Bwall in component Parameters (mmHg_s_per_sq_mm).
C CONSTS(10) is B1 in component Parameters (mmHg_s_per_mm).
C CONSTS(11) is B2 in component Parameters (mmHg_s_per_mm).
C CONSTS(12) is B3 in component Parameters (mmHg_s_per_mm).
C CONSTS(13) is Tsmax in component Parameters (AU).
C CONSTS(14) is Tpmax in component Parameters (AU).
C CONSTS(15) is Tsmin in component Parameters (AU).
C CONSTS(16) is Tpmin in component Parameters (AU).
C CONSTS(17) is Gcns in component Parameters (dimensionless).
C CONSTS(18) is Gs in component Parameters (per_Hertz).
C CONSTS(19) is Gp in component Parameters (per_Hertz).
C CONSTS(20) is tau_nor in component Parameters (second).
C CONSTS(21) is tau_ach in component Parameters (second).
C CONSTS(22) is tau_HR_nor in component Parameters (second).
C CONSTS(23) is tau_HR_ach in component Parameters (second).
C CONSTS(24) is HRo in component Parameters (Beats_per_min).
C CONSTS(25) is HRmax in component Parameters (Beats_per_min).
C CONSTS(26) is HRmin in component Parameters (Beats_per_min).
C CONSTS(27) is Beta in component Parameters (dimensionless).
C CONSTS(28) is delta_th in component Parameters (dimensionless).
C CONSTS(29) is q_nor in component Parameters (per_s).
C CONSTS(30) is q_ach in component Parameters (per_s).
C CONSTS(31) is K_nor in component Parameters (AU).
C CONSTS(32) is K_ach in component Parameters (AU).
C CONSTS(33) is Gamma in component Parameters (dimensionless).
C ALGBRC(11) is alpha_cns in component Nervous_System (hertz).
C ALGBRC(9) is n in component Nervous_System (hertz).
C ALGBRC(7) is Delta in component Coupling_Dynamics (dimensionless).
C CONSTS(34) is alpha_s0 in component Nervous_System (hertz).
C CONSTS(35) is alpha_p0 in component Nervous_System (hertz).
C STATES(1) is A in component Aortic_Wall (mm_sq).
C ALGBRC(1) is P in component Aortic_Wall (mmHg).
C ALGBRC(2) is R in component Aortic_Wall (mm).
C STATES(2) is Eps_1 in component Coupling_Dynamics (dimensionless).
C STATES(3) is Eps_2 in component Coupling_Dynamics (dimensionless).
C STATES(4) is Eps_3 in component Coupling_Dynamics (dimensionless).
C ALGBRC(5) is Eps_wall in component Coupling_Dynamics (dimensionless).
C ALGBRC(13) is Ts in component PNS_tones (AU).
C ALGBRC(14) is Tp in component PNS_tones (AU).
C STATES(5) is c_nor in component Norepinephrine (AU).
C STATES(6) is C_ach in component Acetylcholine (AU).
C ALGBRC(3) is delta_HR_ss in component Heart_Response_Nor (Beats_per_min).
C CONSTS(36) is delta_HR_smax in component Heart_Response_Nor (Beats_per_min).
C STATES(7) is delta_HR_s in component Heart_Response_Nor (Beats_per_min).
C ALGBRC(4) is delta_HR_ps in component HR_ach (Beats_per_min).
C CONSTS(37) is delta_HR_pmax in component HR_ach (Beats_per_min).
C ALGBRC(6) is delta_HR_pfast in component HR_ach (Beats_per_min).
C STATES(8) is delta_HR_pslow in component HR_ach (Beats_per_min).
C ALGBRC(8) is delta_HR_p in component HR_ach (Beats_per_min).
C ALGBRC(15) is HR in component HR_Combined (Beats_per_min).
C ALGBRC(12) is HR_p in component HR_Combined (Beats_per_min).
C ALGBRC(10) is HR_s in component HR_Combined (Beats_per_min).
C ALGBRC(16) is Period in component HR_Combined (Sec_per_Beat).
C RATES(1) is d/dt A in component Aortic_Wall (mm_sq).
C RATES(2) is d/dt Eps_1 in component Coupling_Dynamics (dimensionless).
C RATES(3) is d/dt Eps_2 in component Coupling_Dynamics (dimensionless).
C RATES(4) is d/dt Eps_3 in component Coupling_Dynamics (dimensionless).
C RATES(5) is d/dt c_nor in component Norepinephrine (AU).
C RATES(6) is d/dt C_ach in component Acetylcholine (AU).
C RATES(7) is d/dt delta_HR_s in component Heart_Response_Nor (Beats_per_min).
C RATES(8) is d/dt delta_HR_pslow in component HR_ach (Beats_per_min).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      CONSTS(1) = 480
      CONSTS(2) = 1
      CONSTS(3) = 1.6
      CONSTS(4) = 0.006
      CONSTS(5) = 1
      CONSTS(6) = 1.5
      CONSTS(7) = 3.75
      CONSTS(8) = 1.05
      CONSTS(9) = 1
      CONSTS(10) = 1
      CONSTS(11) = 10
      CONSTS(12) = 206.973
      CONSTS(13) = 4.12
      CONSTS(14) = 4.994
      CONSTS(15) = 0.5
      CONSTS(16) = 1.6
      CONSTS(17) = 1
      CONSTS(18) = 0.178
      CONSTS(19) = 0.492
      CONSTS(20) = 9.1
      CONSTS(21) = 0.2
      CONSTS(22) = 2.1
      CONSTS(23) = 2.5
      CONSTS(24) = 282.648
      CONSTS(25) = 483.218
      CONSTS(26) = 226.238
      CONSTS(27) = 0.175
      CONSTS(28) = 0
      CONSTS(29) = 0.1099
      CONSTS(30) = 5
      CONSTS(31) = 1.12
      CONSTS(32) = 0.65
      CONSTS(33) = 0.75
      CONSTS(34) = 58.6
      CONSTS(35) = 76.019
      STATES(1) = 15.20531
      STATES(2) = 0.2042
      STATES(3) = 0.183
      STATES(4) = 0.161
      STATES(5) = 1.441
      STATES(6) = 1.0
      STATES(7) = 0
      STATES(8) = 0
      CONSTS(36) = CONSTS(25) - CONSTS(24)
      CONSTS(37) = CONSTS(24) - CONSTS(26)
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(1) = 100.000+ 10.0000* sin( 5.00000*VOI)
      RATES(1) = (- ( (STATES(1)/ 3.14159265358979) ** (1.0 / 2) - CONSTS(3))/CONSTS(4)+ALGBRC(1))/CONSTS(9)
      ALGBRC(3) = ( CONSTS(36)*STATES(5) ** 2.00000)/(CONSTS(31) ** 2.00000+STATES(5) ** 2.00000)
      RATES(7) = (- STATES(7)+ALGBRC(3))/CONSTS(22)
      ALGBRC(4) = ( CONSTS(37)*STATES(6) ** 2.00000)/(CONSTS(32) ** 2.00000+STATES(6) ** 2.00000)
      RATES(8) = (- STATES(8)+ (1.00000 - CONSTS(33))*ALGBRC(4))/CONSTS(23)
      ALGBRC(2) = STATES(1)/ 3.14159265358979 ** 0.500000
      ALGBRC(5) = (ALGBRC(2) - CONSTS(3))/CONSTS(3)
CALL minimize(minfunc_0, CONSTS, VARIABLES, <VAR>)
      ALGBRC(7) = ALGBRC(5) - STATES(2)
      ALGBRC(9) =  CONSTS(1)*(ALGBRC(7) -  CONSTS(2)*CONSTS(28))
      ALGBRC(11) =  CONSTS(17)*ALGBRC(9)
      ALGBRC(13) = CONSTS(15)+(CONSTS(13) - CONSTS(15))/(EXP( CONSTS(18)*(ALGBRC(11) - CONSTS(34)))+1.00000)
      RATES(5) = - (STATES(5)/CONSTS(20))+ CONSTS(29)*ALGBRC(13)
      ALGBRC(14) = CONSTS(16)+(CONSTS(14) - CONSTS(16))/(EXP( - CONSTS(19)*(ALGBRC(11) - CONSTS(35)))+1.00000)
      RATES(6) = - (STATES(6)/CONSTS(21))+ CONSTS(30)*ALGBRC(14)
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(1) = 100.000+ 10.0000* sin( 5.00000*VOI)
      ALGBRC(3) = ( CONSTS(36)*STATES(5) ** 2.00000)/(CONSTS(31) ** 2.00000+STATES(5) ** 2.00000)
      ALGBRC(4) = ( CONSTS(37)*STATES(6) ** 2.00000)/(CONSTS(32) ** 2.00000+STATES(6) ** 2.00000)
      ALGBRC(2) = STATES(1)/ 3.14159265358979 ** 0.500000
      ALGBRC(5) = (ALGBRC(2) - CONSTS(3))/CONSTS(3)
      ALGBRC(7) = ALGBRC(5) - STATES(2)
      ALGBRC(9) =  CONSTS(1)*(ALGBRC(7) -  CONSTS(2)*CONSTS(28))
      ALGBRC(11) =  CONSTS(17)*ALGBRC(9)
      ALGBRC(13) = CONSTS(15)+(CONSTS(13) - CONSTS(15))/(EXP( CONSTS(18)*(ALGBRC(11) - CONSTS(34)))+1.00000)
      ALGBRC(14) = CONSTS(16)+(CONSTS(14) - CONSTS(16))/(EXP( - CONSTS(19)*(ALGBRC(11) - CONSTS(35)))+1.00000)
      ALGBRC(6) =  CONSTS(33)*ALGBRC(4)
      ALGBRC(8) = ALGBRC(6)+STATES(8)
      ALGBRC(10) = CONSTS(24)+STATES(7)
      ALGBRC(12) = CONSTS(24) - ALGBRC(8)
      ALGBRC(15) = ALGBRC(12)+( (ALGBRC(10) - CONSTS(24))*(ALGBRC(12) -  CONSTS(27)*CONSTS(26)))/(CONSTS(24) -  CONSTS(27)*CONSTS(26))
      ALGBRC(16) = 60.0000/ALGBRC(15)
      RETURN
      END
REAL FUNCTION minfunc_0(CONSTS, VARIABLES)
            REAL CONSTS(*), VARIABLES(*) 
      minfunc_0[1] = abs(RATES(2) - (( CONSTS(5)*(ALGBRC(5) - STATES(2)) -  CONSTS(6)*(STATES(2) - STATES(3)))/CONSTS(10)+RATES(3)))
      minfunc_0[2] = abs(RATES(3) - (( CONSTS(6)*(STATES(2) - STATES(3)) -  CONSTS(7)*(STATES(3) - STATES(4)))+ CONSTS(10)*RATES(2)+ CONSTS(11)*RATES(4))/(CONSTS(10)+CONSTS(11)))
      minfunc_0[3] = abs(RATES(4) - (( CONSTS(7)*(STATES(3) - STATES(4)) -  CONSTS(8)*STATES(4))+ CONSTS(11)*RATES(3))/(CONSTS(11)+CONSTS(12)))
      RETURN
      END