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

#define N 3

#define Ca  f[0]
#define IP3 f[1]
#define z   f[2]

double f[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.5;
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 d_Ca = 0.001;
//double d_IP3 = 0.12;

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 UllahJung(int i, double f[N])
{
    switch (i)
    {
    case 0: // Ca
        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) ) - ( v_3 * pow(Ca, 2) / (pow(k_3, 2) + pow(Ca, 2)) ) + ( c_1 * v_2 * ( c_0 / c_1 - (1 + 1 / c_1) * Ca ) ) + ( v_5 + v_6 * pow(IP3, 2) / (pow(k_2, 2) + pow(IP3, 2)) ) - k_1 * Ca;

    case 1: // IP3
        return ( IP3_star - IP3 ) / tau_IP3 + v_4 * (Ca + (1 - alpha) * k_4) / (Ca + k_4);

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

void RungeKutta(double dt, double f[N], double f_next[N])
{
    double k[N][4];

    // k1
    for (int i = 0; i < N; i++)
        k[i][0] = UllahJung(i, f) * dt;

    double phi_k1[N];
    for (int i = 0; i < N; i++)
        phi_k1[i] = f[i] + k[i][0] / 2;

    // k2
    for (int i = 0; i < N; i++)
        k[i][1] = UllahJung(i, phi_k1) * dt;

    double phi_k2[N];
    for (int i = 0; i < N; i++)
        phi_k2[i] = f[i] + k[i][1] / 2;

    // k3
    for (int i = 0; i < N; i++)
        k[i][2] = UllahJung(i, phi_k2) * dt;

    double phi_k3[N];
    for (int i = 0; i < N; i++)
        phi_k3[i] = f[i] + k[i][2] / 2;

    // k4
    for (int i = 0; i < N; i++)
        k[i][3] = UllahJung(i, phi_k3) * dt;

    for (int i = 0; i < N; i++)
        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;

    //Initial values at t = 0
    f[0] = Ca_0;
    f[1] = IP3_0;
    f[2] = z_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   = 180; // 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_Ca_IP3_z.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);
}