adedd plplot stuff and etsting a bit

practical_work/rc_circuit/Makefile

-CC=clang
-OPTS=-g -O3 -Wall -Wextra -fsanitize=address -fsanitize=leak -std=c11
+CC=gcc
+OPTS=-g -Wall -Wextra -fsanitize=address -fsanitize=leak -std=c11 #-O3
 INCLUDE=-I/usr/include/plplot/
 LINK=-lm -fsanitize=address -fsanitize=leak  -lplplot
 
@@ -18,9 +18,5 @@ ode_o1.o: ode_o1.c ode_o1.h
 rc.o: rc.c rc.h
 	$(CC) $(OPTS) -c $^
 
-test:
-	make -C tests test
-
 clean:
 	rm -f *.o main
-	make -C tests clean

practical_work/rc_circuit/main.c

-#include "rc.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #include <plplot.h>
 
-#define NSIZE 5000
+#define NSIZE 500
 
-int main(int argc, char *argv[]) {
-    rc_state rc = rc_state_create(1.0, 1.0, 1.0);
+double PI = 4.0 * atan(1.0);
 
+typedef struct _rc_circuit
+{
+    double r, c;
+} rc_circuit;
 
-// ============================================================ //
-// RC with constant != 0 tension at input and 0 initial tension //
-// ============================================================ //
-    double v, v_ini;
-    v = v_ini = 0.0;
+rc_circuit rc_circuit_create(double r, double c)
+{
+    rc_circuit tmp = {.r = r, .c = c};
+    return tmp;
+}
 
-    double dt = 0.001;
-    double max_t = 5.0;
+double rc_circuit_charge(double v_in, rc_circuit rc, double t)
+{
+    return v_in * (1.0 - exp(-t / (rc.c * rc.r)));
+}
 
-    PLFLT x[NSIZE], y[NSIZE];
-    PLFLT xmin = 0., xmax = max_t, ymin = rc.eps-1, ymax = rc.eps+1;
+double rc_circuit_discharge(double v_in, rc_circuit rc, double t)
+{
+    return v_in * exp(-t / (rc.c * rc.r));
+}
 
-    for (double t = 0.0; t < max_t; t += dt) {
-        double ve = exact_solution(t, v_ini, &rc);
-        printf("t = %f, v = %f, v_e = %f, diff = %f\n", t, v, ve, ve - v);
-        v = rc_advance(v, dt, &rc);
+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);
 
-    v = v_ini = 1.0;
-    double omega_low = 1.0;
-    double omega = 100.0;
     int i = 0;
-    for (double t = 0.0; t < max_t; t += dt) {
-        rc.eps = v_ini * (1.0 + cos(omega * t) + cos(omega_low * t));
-        double ve = exact_solution(t, v_ini, &rc);
-        printf("t = %f, v = %f, v_e = %f, diff = %f\n", t, v, ve, ve - v);
-        x[i] = t;
-        y[i] = v;
-        v = rc_advance(v, dt, &rc);
+    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 );
+    plparseopts(&argc, argv, PL_PARSE_FULL);
+    plsetopt("dev", "wxwidgets");
+    plssub(2, 2);
 
     // Initialize plplot
     plinit();
 
-    // Create a labelled box to hold the plot.
-    plenv( xmin, xmax, ymin, ymax, 0, 0 );
-    pllab( "t", "v", "Simple PLplot demo of of the charge of a capacitor" );
+    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");
 
-    // Plot the data that was prepared above.
-    plline( NSIZE, x, y );
+    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();