Select Git revision
main.c 6.68 KiB
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <plplot.h>
#define NSIZE 5000
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;
}
// TODO: try to do the high pass filter.
// double solve_vr(rc_circuit rc, double v0, double dv_in_dt, double dt)
// {
// return (dv_in_dt * (rc.c * rc.r) - 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.01 * rc.r * rc.c);
double omega2 = 1.0 / (2.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) + v_in * sin(2 * PI * omega2 * 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(1, 3);
// 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, 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");
compute_sin(NSIZE, t, vc, vr, v_input, rc, v0, v_in);
plot_data(NSIZE, t, vc, vr, v_input, "t", "v", "The discharge of a capacitor");
// Close PLplot library
plend();
return EXIT_SUCCESS;
}