hh52.c

#include <math.h>
#include <stdio.h>

void initializeConstants(double *states, double *variables)
{
    states[0] = 0.05;
    states[1] = 0.6;
    states[2] = 0.325;
    states[3] = 0.0;
    variables[0] = 0.3;
    variables[1] = 1.0;
    variables[2] = 0.0;
    variables[3] = 36.0;
    variables[4] = 120.0;
}

void computeComputedConstants(double *variables)
{
    variables[6] = variables[2]-10.613;
    variables[8] = variables[2]-115.0;
    variables[14] = variables[2]+12.0;
}

void computeRates(double voi, double *states, double *rates, double *variables)
{
    variables[10] = 0.1*(states[3]+25.0)/(exp((states[3]+25.0)/10.0)-1.0);
    variables[11] = 4.0*exp(states[3]/18.0);
    rates[0] = variables[10]*(1.0-states[0])-variables[11]*states[0];
    variables[12] = 0.07*exp(states[3]/20.0);
    variables[13] = 1.0/(exp((states[3]+30.0)/10.0)+1.0);
    rates[1] = variables[12]*(1.0-states[1])-variables[13]*states[1];
    variables[16] = 0.01*(states[3]+10.0)/(exp((states[3]+10.0)/10.0)-1.0);
    variables[17] = 0.125*exp(states[3]/80.0);
    rates[2] = variables[16]*(1.0-states[2])-variables[17]*states[2];
    variables[5] = ((voi >= 10.0) && (voi <= 10.5))?-20.0:0.0;
    variables[7] = variables[0]*(states[3]-variables[6]);
    variables[15] = variables[3]*pow(states[2], 4.0)*(states[3]-variables[14]);
    variables[9] = variables[4]*pow(states[0], 3.0)*states[1]*(states[3]-variables[8]);
    rates[3] = -(-variables[5]+variables[9]+variables[15]+variables[7])/variables[1];
}

void computeVariables(double voi, double *states, double *rates, double *variables)
{
    variables[7] = variables[0]*(states[3]-variables[6]);
    variables[9] = variables[4]*pow(states[0], 3.0)*states[1]*(states[3]-variables[8]);
    variables[10] = 0.1*(states[3]+25.0)/(exp((states[3]+25.0)/10.0)-1.0);
    variables[11] = 4.0*exp(states[3]/18.0);
    variables[12] = 0.07*exp(states[3]/20.0);
    variables[13] = 1.0/(exp((states[3]+30.0)/10.0)+1.0);
    variables[15] = variables[3]*pow(states[2], 4.0)*(states[3]-variables[14]);
    variables[16] = 0.01*(states[3]+10.0)/(exp((states[3]+10.0)/10.0)-1.0);
    variables[17] = 0.125*exp(states[3]/80.0);
}

#define NB_OF_STATES 4
#define NB_OF_VARIABLES 18
#define DELTA_T 0.01

void outputData(double voi, double *states, double *variables)
{
    printf("%f", voi);

    for (int i = 0; i < NB_OF_STATES; ++i) {
        printf (",%f", states[i]);
    }

    for (int i = 0; i < NB_OF_VARIABLES; ++i) {
        printf (",%f", variables[i]);
    }

    printf("\n");
}

int main()
{
    double states[NB_OF_STATES];
    double rates[NB_OF_STATES];
    double variables[NB_OF_VARIABLES];

    initializeConstants(states, variables);
    computeComputedConstants(variables);
    computeRates(0.0, states, rates, variables);
    computeVariables(0.0, states, rates, variables);

    outputData(0.0, states, variables);

    for (int i = 0; i < 5000; ++i) {
        double voi = i*DELTA_T;

        computeRates(voi, states, rates, variables);

        for (int j = 0; j < NB_OF_STATES; ++j) {
            states[j] += DELTA_T*rates[j];
        }

        if ((i+1) % 10 == 0) {
            computeVariables(voi, states, rates, variables);

            outputData((i+1)*DELTA_T, states, variables);
        }
    }

    return 0;
}