n62.c

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

void initializeConstants(double *states, double *variables)
{
    states[0] = 0.01;
    states[1] = 0.8;
    states[2] = 0.01;
    states[3] = -87.0;
    variables[0] = -60.0;
    variables[1] = 0.075;
    variables[2] = 12.0;
    variables[3] = 40.0;
    variables[4] = 400.0;
}

void computeComputedConstants(double *variables)
{
}

void computeRates(double voi, double *states, double *rates, double *variables)
{
    variables[8] = 0.1*(-states[3]-48.0)/(exp((-states[3]-48.0)/15.0)-1.0);
    variables[9] = 0.12*(states[3]+8.0)/(exp((states[3]+8.0)/5.0)-1.0);
    rates[0] = variables[8]*(1.0-states[0])-variables[9]*states[0];
    variables[10] = 0.17*exp((-states[3]-90.0)/20.0);
    variables[11] = 1.0/(1.0+exp((-states[3]-42.0)/10.0));
    rates[1] = variables[10]*(1.0-states[1])-variables[11]*states[1];
    variables[14] = 0.0001*(-states[3]-50.0)/(exp((-states[3]-50.0)/10.0)-1.0);
    variables[15] = 0.002*exp((-states[3]-90.0)/80.0);
    rates[2] = variables[14]*(1.0-states[2])-variables[15]*states[2];
    variables[6] = pow(states[0], 3.0)*states[1]*variables[4];
    variables[7] = (variables[6]+0.14)*(states[3]-variables[3]);
    variables[5] = variables[1]*(states[3]-variables[0]);
    variables[12] = 1.2*exp((-states[3]-90.0)/50.0)+0.015*exp((states[3]+90.0)/60.0);
    variables[13] = 1.2*pow(states[2], 4.0);
    variables[16] = (variables[12]+variables[13])*(states[3]+100.0);
    rates[3] = -(variables[7]+variables[16]+variables[5])/variables[2];
}

void computeVariables(double voi, double *states, double *rates, double *variables)
{
    variables[5] = variables[1]*(states[3]-variables[0]);
    variables[6] = pow(states[0], 3.0)*states[1]*variables[4];
    variables[7] = (variables[6]+0.14)*(states[3]-variables[3]);
    variables[8] = 0.1*(-states[3]-48.0)/(exp((-states[3]-48.0)/15.0)-1.0);
    variables[9] = 0.12*(states[3]+8.0)/(exp((states[3]+8.0)/5.0)-1.0);
    variables[10] = 0.17*exp((-states[3]-90.0)/20.0);
    variables[11] = 1.0/(1.0+exp((-states[3]-42.0)/10.0));
    variables[12] = 1.2*exp((-states[3]-90.0)/50.0)+0.015*exp((states[3]+90.0)/60.0);
    variables[13] = 1.2*pow(states[2], 4.0);
    variables[14] = 0.0001*(-states[3]-50.0)/(exp((-states[3]-50.0)/10.0)-1.0);
    variables[15] = 0.002*exp((-states[3]-90.0)/80.0);
    variables[16] = (variables[12]+variables[13])*(states[3]+100.0);
}

#define NB_OF_STATES 4
#define NB_OF_VARIABLES 17
#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 < 100000; ++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) % 100 == 0) {
            computeVariables(voi, states, rates, variables);

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

    return 0;
}