#define _USE_MATH_DEFINES
#include "math.h"
#include <stdlib.h>
#include <stdio.h>
#include <locale.h>
#include <time.h>
#include <stdbool.h>

#define N 8

#define Ca  f[0]
#define IP3 f[1]
#define z   f[2]
#define V   f[3]
#define m   f[4]
#define n   f[5]
#define h   f[6]
#define G   f[7]

double f[N];

double k[N][4];
double phi_k1[N];
double phi_k2[N];
double phi_k3[N];

double c_0 = 2;
double c_1 = 0.185;
double v_1 = 6;
double v_2 = 0.11;
double v_3 = 2.2;
double v_4 = 0.0;
double v_5 = 0.025;
double v_6 = 0.2;
double k_1 = 0.5;
double k_2 = 1;
double k_3 = 0.1;
double k_4 = 1.1;
double a_2 = 0.14;
double d_1 = 0.13;
double d_2 = 1.049;
double d_3 = 0.9434;
double d_5 = 0.082;
double alpha_G = 25;
double beta_G = 500;
double alpha = 0.8;
double tau_IP3 = 7.143;
double IP3_star = 0.16;
double alpha_Glu = 2;
//double d_Ca = 0.001;
//double d_IP3 = 0.12;

double C = 1; // muF/cm^2
double g_K  =  35; // mS/cm^2
double g_Na =  40; // mS/cm^2
double g_L  = 0.3; // mS/cm^2
double E_K  = -77; // mV
double E_Na =  55; // mV
double E_L  = -65; // mV
double I_app = 1.2;

int RandomI(int min, int max)
{
    return ((double)rand() / (RAND_MAX - 1)) * (max - min) + min;
}

double RandomD(double min, double max)
{
    return ((double)rand() / RAND_MAX) * (max - min) + min;
}

double J_channel(double f[N])
{
    return c_1 * v_1 * pow(IP3, 3) * pow(Ca, 3) * pow(z, 3) * (c_0 / c_1 - (1 + 1 / c_1) * Ca) / pow((IP3 + d_1) * (Ca + d_5), 3);
}

double J_pump(double f[N])
{
    return v_3 * pow(Ca, 2) / (pow(k_3, 2) + pow(Ca, 2));
}

double J_leak(double f[N])
{
    return c_1 * v_2 * (c_0 / c_1 - (1 + 1 / c_1) * Ca);
}

double J_in(double f[N])
{
    return v_5 + v_6 * pow(IP3, 2) / (pow(k_2, 2) + pow(IP3, 2));
}

double J_out(double f[N])
{
    return k_1 * Ca;
}

double J_PLC(double f[N])
{
    return v_4 * (Ca + (1 - alpha) * k_4) / (Ca + k_4);
}

double J_Glu(double f[N])
{
    return alpha_Glu / (1 + exp(-(G - 0.25) / 0.01));
}

double alpha_n(double f[N])
{
    return 0.02 * (V - 25) / (1 - exp(-(V - 25) / 9));
}

double beta_n(double f[N])
{
    return -0.002 * (V - 25) / (1 - exp((V - 25) / 9));
}

double alpha_m(double f[N])
{
    return 0.182 * (V + 35) / (1 - exp(-(V + 35) / 9));
}

double beta_m(double f[N])
{
    return -0.124 * (V + 35) / (1 - exp((V + 35) / 9));
}

double alpha_h(double f[N])
{
    return 0.25 * exp(-(V + 90) / 12);
}

double beta_h(double f[N])
{
    return 0.25 * exp((V + 62) / 6) / exp((V + 90) / 12);
}

double UllahJung_HodgkinHuxley(int i, double f[N])
{
    switch (i)
    {
    case 0: // Ca
        return J_channel(f) - J_pump(f) + J_leak(f) + J_in(f) - J_out(f);

    case 1: // IP3
        return (IP3_star - IP3) / tau_IP3 + J_PLC(f) + J_Glu(f);

    case 2: // z
        return a_2 * ( d_2 * ((IP3 + d_1) / (IP3 + d_3)) * (1 - z) - Ca * z );

    case 3: // V
        return 1000 * ((g_Na * pow(m, 3) * h * (E_Na - V) + g_K * n * (E_K - V) + g_L * (E_L - V) + I_app) / C);

    case 4: // m
        return 1000 * (alpha_m(f) * (1 - m) - beta_m(f) * m);

    case 5: // n
        return 1000 * (alpha_n(f) * (1 - n) - beta_n(f) * n);

    case 6: // h
        return 1000 * (alpha_h(f) * (1 - h) - beta_h(f) * h);

    case 7: // G
        return -alpha_G * G + beta_G * (1 / (1 + exp(-V / 0.5)));
    }

    return 0;
}

void RungeKutta(double dt, double f[N], double f_next[N])
{
    #pragma omp parallel for
    for (int i = 0; i < N; i++)
    {
        k[i][0] = UllahJung_HodgkinHuxley(i, f) * dt;
        phi_k1[i] = f[i] + k[i][0] / 2;
        k[i][1] = UllahJung_HodgkinHuxley(i, phi_k1) * dt;
        phi_k2[i] = f[i] + k[i][1] / 2;
        k[i][2] = UllahJung_HodgkinHuxley(i, phi_k2) * dt;
        phi_k3[i] = f[i] + k[i][2] / 2;
        k[i][3] = UllahJung_HodgkinHuxley(i, phi_k3) * dt;

        f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
    }
}

void CopyArray(double source[N], double target[N])
{
    for (int i = 0; i < N; i++)
        target[i] = source[i];
}

bool Approximately(double a, double b)
{
    if (a < 0)
        a = -a;

    if (b < 0)
        b = -b;

    return a - b <= 0.000001;
}

int main()
//main(int argc, char *argv[])
{
    //sscanf(argv[1], "%lf", &v_4);

    FILE *fp0;
    srand(time(NULL));

    double Ca_0  = 0.07;
    double IP3_0 = 0.16;
    double z_0   = 0.67;

    double V_0 = -58.7085;
    double m_0 = 0.0953;
    double n_0 = 0.000913;
    double h_0 = 0.3662;

    double G_0 = 0;

    //Initial values at t = 0
    f[0] = Ca_0;
    f[1] = IP3_0;
    f[2] = z_0;

    f[3] = V_0;
    f[4] = m_0;
    f[5] = n_0;
    f[6] = h_0;

    f[7] = G_0;

    /*fp0 = fopen("last_values.txt", "r");
    for (int i = 0; i < N; i++)
    {
        fscanf(fp0, "%lf", &f[i]);
    }
    fclose(fp0);*/

    const double t_start = 0;
    const double t_max   = 120; // sec
    const double dt      = 0.00001;

    double t = t_start;

    //fp0 = fopen("v_4.txt", "w+");
    //fprintf(fp0, "%f\n", v_4);
    //fclose(fp0);

    fp0 = fopen("results.txt", "w+");
    //setlocale(LC_NUMERIC, "French_Canada.1252");

    clock_t start_rk4, end_rk4;
    start_rk4 = clock();
    int lastPercent = -1;

    while (t < t_max || Approximately(t, t_max))
    {
        fprintf(fp0, "%f\t", t);
        for (int i = 0; i < N; i++)
        {
            fprintf(fp0, i == N - 1 ? "%f" : "%f\t", f[i]);
        }
        fprintf(fp0, "\n");

        double f_next[N];

        RungeKutta(dt, f, f_next);
        CopyArray(f_next, f);

        t += dt;

        int percent = (int)(100 * (t - t_start) / (t_max - t_start));
        if (percent != lastPercent)
        {
            printf("Progress: %d%%\n", percent);
            lastPercent = percent;
        }
    }

    end_rk4 = clock();
    double extime_rk4 = (double)(end_rk4 - start_rk4) / CLOCKS_PER_SEC;
    int minutes = (int)extime_rk4 / 60;
    int seconds = (int)extime_rk4 % 60;
    printf("\nExecution time is: %d minutes %d seconds\n ", minutes, seconds);

    fclose(fp0);

    fp0 = fopen("time_exec.txt", "w+");
    fprintf(fp0, "%f\n", extime_rk4);
    fclose(fp0);

    /*fp0 = fopen("last_values.txt", "w+");
    for (int i = 0; i < N; i++)
    {
        fprintf(fp0, i == N - 1 ? "%f" : "%f\t", f[i]);
    }
    fclose(fp0);*/
}