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

#define NSIZE 500

double PI = 4.0 * atan(1.0);

typedef struct _rc_circuit
{
    double r, c;
} rc_circuit;

rc_circuit rc_circuit_create(double r, double c)
{
    rc_circuit tmp = {.r = r, .c = c};
    return tmp;
}

double rc_circuit_charge(double v_in, rc_circuit rc, double t)
{
    return v_in * (1.0 - exp(-t / (rc.c * rc.r)));
}

double rc_circuit_discharge(double v_in, rc_circuit rc, double t)
{
    return v_in * exp(-t / (rc.c * rc.r));
}

double solve_vc(rc_circuit rc, double v0, double v_in, double dt)
{
    return (v_in - v0) / (rc.c * rc.r) * dt + v0;
}

double solve_vr(double v_in, double vc)
{
    return v_in - vc;
}

PLFLT rc_min(int size, PLFLT x[size])
{
    PLFLT val = x[0];
    for (int i = 0; i < size; ++i)
    {
        val = x[i] < val ? x[i] : val;
    }
    return val;
}

PLFLT rc_max(int size, PLFLT x[size])
{
    PLFLT val = x[0];
    for (int i = 0; i < size; ++i)
    {
        val = x[i] > val ? x[i] : val;
    }
    return val;
}

void plot_data(int size, PLFLT x[size], PLFLT y[size], PLFLT z[size], PLFLT ref[size], char *xaxis, char *yaxis, char *title)
{
    PLFLT xmin = rc_min(size, x);
    PLFLT xmax = rc_max(size, x);
    PLFLT ymin = fmin(rc_min(size, ref), fmin(rc_min(size, z), rc_min(size, y)));
    PLFLT ymax = fmax(rc_max(size, ref), fmax(rc_max(size, z), rc_max(size, y)));

    PLINT nlegend = 3;
    PLCHAR_VECTOR text[3], symbols[3];
    PLINT opt_array[3];
    PLINT text_colors[3];
    PLINT line_colors[3];
    PLINT line_styles[3];
    PLFLT line_widths[3];
    PLINT symbol_numbers[3], symbol_colors[3];
    PLFLT symbol_scales[3];
    PLFLT legend_width, legend_height;

    // Draw a legend
    // First legend entry.
    opt_array[0] = PL_LEGEND_LINE;
    text_colors[0] = 3;
    text[0] = "V_c";
    line_colors[0] = 3;
    line_styles[0] = 1;
    line_widths[0] = 1.;
    // note from the above opt_array the first symbol (and box) indices
    // do not have to be specified, at least in C.
    symbols[0] = "";
    // Second legend entry.
    opt_array[1] = PL_LEGEND_LINE;
    text_colors[1] = 4;
    text[1] = "V_r";
    line_colors[1] = 4;
    line_styles[1] = 1;
    line_widths[1] = 1.;
    // note from the above opt_array the first symbol (and box) indices
    // do not have to be specified, at least in C.
    symbols[1] = "";
    // First legend entry.
    opt_array[2] = PL_LEGEND_LINE;
    text_colors[2] = 5;
    text[2] = "V_in";
    line_colors[2] = 5;
    line_styles[2] = 1;
    line_widths[2] = 1.;
    // note from the above opt_array the first symbol (and box) indices
    // do not have to be specified, at least in C.
    symbols[2] = "";

    // Create a labelled box to hold the plot.
    plcol0(2);
    plenv(xmin, xmax, ymin, ymax, 0, 0);
    pllab(xaxis, yaxis, title);

    // Plot the data that was prepared above.
    plcol0(3);
    plline(size, x, y);
    plcol0(4);
    plline(size, x, z);
    plcol0(5);
    plline(size, x, ref);

    pllegend(&legend_width, &legend_height,
             PL_LEGEND_BACKGROUND | PL_LEGEND_BOUNDING_BOX, 0,
             0.0, 0.0, 0.1, 15,
             1, 1, 0, 0,
             nlegend, opt_array,
             1.0, 1.0, 2.0,
             1., text_colors, (const char **)text,
             NULL, NULL, NULL, NULL,
             line_colors, line_styles, line_widths,
             symbol_colors, symbol_scales, symbol_numbers, (const char **)symbols);
    plflush();
}

void compute_charge(int size, PLFLT time[size], PLFLT vc[size], PLFLT vr[size], PLFLT v_input[size], rc_circuit rc, double v0, double v_in)
{
    double dt = 5.0 * rc.r * rc.c / size;
    double max_t = size * dt;
    double v = v0;

    int i = 0;
    for (double t = 0.0; t < max_t && i < size; t += dt)
    {
        // double ve = rc_circuit_charge(v_in, rc, t);
        // printf("t = %f, v = %f, v_e = %f, diff = %f\n", t, v, ve, ve - v);
        time[i] = t;
        vc[i] = v;
        vr[i] = solve_vr(v_in, v);
        v_input[i] = v_in;
        v = solve_vc(rc, v, v_in, dt);
        i += 1;
    }
}

void compute_discharge(int size, PLFLT time[size], PLFLT vc[size], PLFLT vr[size], PLFLT v_input[size], rc_circuit rc, double v0, double v_in)

{
    double dt = 5.0 * rc.r * rc.c / size;
    double max_t = size * dt;
    double v = v0;

    int i = 0;
    for (double t = 0.0; t < max_t && i < size; t += dt)
    {
        // double ve = rc_circuit_discharge(v0, rc, t);
        // printf("i = %d, t = %f, v = %f, v_e = %f, diff = %f\n", i, t, v, ve, ve - v);
        time[i] = t;
        vc[i] = v;
        vr[i] = solve_vr(v_in, v);
        v_input[i] = v_in;
        v = solve_vc(rc, v, v_in, dt);
        i += 1;
    }
}

void compute_sin(int size, PLFLT time[size], PLFLT vc[size], PLFLT vr[size], PLFLT v_input[size], rc_circuit rc, double v0, double v_in)
{
    double dt = 5.0 * rc.r * rc.c / size;
    double max_t = size * dt;
    double v = v0;
    double v_in_time = v_in;
    double omega1 = 1.0 / (0.1 * rc.r * rc.c);
    double omega2 = 1.0 / (10.0 * rc.r * rc.c);

    int i = 0;
    for (double t = 0.0; t < max_t && i < size; t += dt)
    {
        // double ve = rc_circuit_charge(v_in, rc, t);
        // printf("t = %f, v = %f, v_e = %f, diff = %f\n", t, v, ve, ve - v);
        time[i] = t;
        vc[i] = v;
        vr[i] = solve_vr(v_in_time, v);
        v_input[i] = v_in_time;
        v = solve_vc(rc, v, v_in_time, dt);
        v_in_time = v_in * sin(2 * PI * omega1 * t);
        i += 1;
    }
}

int main(int argc, char *argv[])
{
    rc_circuit rc = rc_circuit_create(1.0, 1.0);

    // ============================================================ //
    // RC with constant != 0 tension at input and 0 initial tension //
    // ============================================================ //

    // Parse and process command line arguments
    plparseopts(&argc, argv, PL_PARSE_FULL);
    plsetopt("dev", "wxwidgets");
    plssub(2, 2);

    // Initialize plplot
    plinit();

    double v0 = 0.0;
    double v_in = 1.0;
    PLFLT t[NSIZE], vc[NSIZE], vr[NSIZE], v_input[NSIZE];
    // compute_charge(NSIZE, t, vc, vr, rc, v0, v_in);
    compute_sin(NSIZE, t, vc, vr, v_input, rc, v0, v_in);
    plot_data(NSIZE, t, vc, vr, v_input, "t", "v", "The charge of a capacitor");

    compute_discharge(NSIZE, t, vc, vr, v_input, rc, v_in, v0);
    plot_data(NSIZE, t, vc, vr, v_input, "t", "v", "The discharge of a capacitor");

    // Close PLplot library
    plend();

    return EXIT_SUCCESS;
}