diff --git a/practical_work/rc_circuit/Makefile b/practical_work/rc_circuit/Makefile index 7580a51a5dd9711480606667893c38b22f6234d7..75b9932f33432aec4cd25d2d691ad7a86dee23c0 100644 --- a/practical_work/rc_circuit/Makefile +++ b/practical_work/rc_circuit/Makefile @@ -1,5 +1,5 @@ -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 diff --git a/practical_work/rc_circuit/main.c b/practical_work/rc_circuit/main.c index 6486b81f0016cb7bb9aa8da3578677fcd10d1c7c..9ff69495ee2e15c4d9f228648035183a4eb5ed0e 100644 --- a/practical_work/rc_circuit/main.c +++ b/practical_work/rc_circuit/main.c @@ -1,59 +1,229 @@ -#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();