// Formas de onda de un circuito RLC serie o paralelo.
// Autor: Juan Claudio Regidor
// Versi—n: 1.0
// Fecha: 4/jun/02
//

title "Circuito RLC de 2do orden"

help
{@
Formas de onda de un circuito RLC serie o paralelo.
@}

variable R		// Resistencia
variable L		// Inductancia
variable C		// Capacitancia
variable i0		// Corriente inicial del inductor
variable v0		// Voltaje inicial del condensador
variable vl
variable vc
variable il
variable tmin
variable tmax
variable t
variable a		// tipo de circuito
variable esc	// escala variable?
variable kvar	//
variable w		// tama–o de la ventana
variable g

init (tmin, tmax, a, R, L, C, i0, v0, esc, kvar, g, w) = init

make (t, vl, vc, il) = solve(tmin, tmax, a, R, L, C, i0, v0)


menu "RLC serie" _checkmark(a == 1) a = 1
menu "RLC paralelo" _checkmark(a == 2) a = 2
separator
menu "Raices" _checkmark(g == 1) g = 1
menu "Potencia" _checkmark(g == 2) g = 2
menu "Energ’a almacenada" _checkmark(g == 3) g = 3
separator
menu "800x600" _checkmark(w == 1) w = 1
menu "1024x768" _checkmark(w == 2) w = 2
menu "1152x768" _checkmark(w == 3) w = 3


figure "Respuesta segundo orden"
	draw res(0, 10, a, R, L, C, i0, v0, esc, g, w)

figure "Par‡metros"
	draw select(a, R, L, C, i0, v0, esc, kvar)
	mousedrag (R, L, C, i0, v0, esc, kvar) = dragSlider(R, L, C, i0, v0, esc, kvar, _id, _nb, _x1)

figure "Raices"
	draw raices(R, L, C, esc)

figure "Potencia"
	draw potencia(t, tmin, tmax, a, esc, vl, vc, il)

figure "Energia almacenada"
	draw energia(t, tmin, tmax, a, esc, vc, il, L, C)



functions
{@
function (tmin, tmax, a, R, L, C, i0, v0, esc, kvar, g, w) = init
	tmin = 0;
	tmax = 10;
	a = 1;
	R = 1;
	L = 1;
	C = 1;
	i0 = 0.25;
	v0 = 0.8;
	esc = 0;
	kvar = 0;
	g = 1;
	w = 1;


function res(tmin, tmax, a, R, L, C, i0, v0, esc, g, w)
	switch g
	case 1
		subplots('Respuesta segundo orden\tRaices\nPar‡metros');
	case 2
		subplots('Respuesta segundo orden\tPotencia\nPar‡metros');
	case 3
		subplots('Respuesta segundo orden\tEnergia almacenada\nPar‡metros');
	end
	switch w
	case 1
		x = 800; y = 600;
	case 2
		x = 1024; y = 768;
	case 3
		x = 1152; y = 768;
	end
	x=x-26;y=y-80; ux=x/8; uy=y/8; uy3 = 19.5;
	subplotpos([0,ux/2,0,uy-uy3; ux/2,ux,0,uy-uy3; 0,ux,uy-uy3,uy;0,0,0,0]);
	//subplotpos([0,ux/2,0,uy/2+4; ux/2+1,ux,0,uy/2+4; 0,ux,uy/2+5,uy;0,0,0,0]);
	resp(tmin, tmax, a, R, L, C, i0, v0, esc);


function suma_exp(tmin, tmax, tau1, tau2, k1, k2, esc);
	u = 0;
	if esc == 0
		scale([tmin, tmax, -1, 1]);
		u = 1;
	end
	(t, x1) = sigExp(tmin, tmax, 0, tau1);
	(t, x2) = sigExp(tmin, tmax, 0, tau2);
	plot(t, k1*x1, 'g:');
	plot(t, k2*x2, 'm:');
	plot(t, k1*x1 + k2*x2, 'r');
	if u == 0
		u = max(max(k1*x1), max(k2*x2));
		u = max(u, max(k1*x1 + k2*x2));
	end
	text(tmax, 0.8*u, sprintf('%.2G*exp(-t/%.2G) %+.2G*exp(-t/%.2G)', k1, tau1, k2, tau2), 'r');
	line([1, 0], 0, 'c');
	line([0, 1], 0, 'c');


function t_exp(tmin, tmax, tau, k1, k2, esc);
	u = 0;
	if esc == 0
		scale([tmin, tmax, -1, 1]);
		u = 1;
	end
	(t, x) = sigExp(tmin, tmax, 0, tau);
	plot(t, (k1*t + k2).*x, 'r');
	if u == 0
		u = max((k1*t + k2).*x);
	end
	text(tmax, 0.8*u, sprintf('(%.2G*t %+.2G)*exp(-t/%.2G)', k1, k2, tau), 'r');
	line([1, 0], 0, 'c');
	line([0, 1], 0, 'c');


function seno_amort(tmin, tmax, tau, omega, k1, k2, esc)
	u = 0;
	if esc == 0
		scale([tmin, tmax, -1, 1]);
		u = 1;
	end
	(t, x) = sigExp(tmin, tmax, 0, tau);
	r = sqrt(k1*k1 + k2*k2);
	plot(t, r*x, 'c-');
	plot(t, -r*x, 'c-');
	if u == 0
		u = r*x(3);
	end
	x = x.*(k1*cos(omega*t) + k2*sin(omega*t));
	plot(t, x, 'r');
	text(tmax, 0.8*u, sprintf('exp(-t/%.2G)*[%.2G*cos(%.2G*t) %+.2G*sin(%.2G*t)]', tau, k1, omega, k2, omega), 'r');
	line([1, 0], 0, 'c');
	line([0, 1], 0, 'c');


function resp(tmin, tmax, a, R, L, C, i0, v0, esc)
	b = R/L;
	c = 1/(L*C);
	disc = b*b - 4*c;
	ii0 = -(v0 + R*i0)/L;
	if abs(disc) < 1e-4
		tau = 2/b;
		k2 = i0;
		k1 = ii0 + i0/tau;
		t_exp(tmin, tmax, tau, k1, k2, esc);
	elseif disc > 0
		disc = sqrt(disc);
		tau1 = (b == -disc)? 1e10 : 2/(b + disc);
		tau2 = (b == disc)? 1e10 : 2/(b - disc);
		k2 = (ii0 + i0/tau1)/(1/tau1-1/tau2);
		k1 = i0 - k2;
		suma_exp(tmin, tmax, tau1, tau2, k1, k2, esc);
	else
		alfa = b/2;
		omega = sqrt(-disc)/2;
		k1 = i0;
		k2 = (ii0 + alfa*i0)/omega;
		seno_amort(tmin, tmax, 1/alfa, omega, k1, k2, esc);
	end
	


function raices(R, L, C, esc)
	if esc == 0
		scale([-4, 2, -3, 3]);
	end
	scale('equal');
	b = R/L;
	c = 1/(L*C);
	line([1, 0], 0, 'c');
	line([0, 1], 0, 'c');
	plotroots([1, b, c]);
	text(2, 2.8, sprintf('Raices de s^2 %+.2G*s %+.2G', b, c), 'r');



function (t, vl, vc, il) = solve(tmin, tmax, a, R, L, C, i0, v0)
	b = R/L;
	cc = 1/(L*C);
	disc = b*b - 4*cc;
	ii0 = -(v0 + R*i0)/L;
	if abs(disc) < 1e-4
		tau = 2/b;
		k2 = i0;
		k1 = ii0 + i0/tau;
		kb = v0;
		ka = i0/C + v0/tau
		(t, x) = sigExp(tmin, tmax, 0, tau);
		il = (k1*t + k2).*x;
		vc = (ka*t + kb).*x;
		vl = L*(-k1/tau*t + (k1 - k2/tau)).*x;
	elseif disc > 0
		disc = sqrt(disc);
		tau1 = (b == -disc)? 1e10 : 2/(b + disc);
		tau2 = (b == disc)? 1e10 : 2/(b - disc);
		k2 = (ii0 + i0/tau1)/(1/tau1-1/tau2);
		k1 = i0 - k2;
		kb = (C*v0/tau1 + i0)/(C*(1/tau1 - 1/tau2));
		ka = v0 - kb;
		(t, x1) = sigExp(tmin, tmax, 0, tau1);
		(t, x2) = sigExp(tmin, tmax, 0, tau2);
		il = k1*x1 + k2*x2;
		vc = ka*x1 + kb*x2;
		vl = -L*(k1/tau1*x1 + k2/tau2*x2);
	else
		alfa = b/2;
		omega = sqrt(-disc)/2;
		k1 = i0;
		k2 = (ii0 + alfa*i0)/omega;
		ka = v0;
		kb = (i0/C + alfa*v0)/omega;
		(t, x) = sigExp(tmin, tmax, 0, 1/alfa);
		il = x.*(k1*cos(omega*t) + k2*sin(omega*t));
		vc = x.*(ka*cos(omega*t) + kb*sin(omega*t));
		vl = L*x.*((omega*k2-alfa*k1)*cos(omega*t) - (alfa*k2+omega*k1)*sin(omega*t));
	end



function potencia(t, tmin, tmax, a, esc, vl, vc, il)
	pc = vc.*il;
	pl = vl.*il;
	plot_w_p(esc, [tmin, tmax, -1, 1], t, pl, pc, a)



function energia(t, tmin, tmax, a, esc, vc, il, L, C)
	wl = (0.5*L).*il.*il;
	wc = (0.5*C).*vc.*vc;
	plot_w_p(esc, [tmin, tmax, 0, 1], t, wl, wc, a)



function plot_w_p(esc, escala, t, al, ac, a)
	if esc == 0
		scale(escala);
	end
	plot(t, al, 'r');
	plot(t, ac, 'b');
	sc = scale;
	if sc === []
		u = max(max(al), max(ac));
	else
		u = sc(4);
	end
	if a == 1
		t1 = 'Condensador';
		t2 = 'Inductor';
	else
		t1 = 'Inductor';
		t2 = 'Condensador';
	end
	tmax = escala(2);
	text(0.8*tmax, 0.8*u, t1, 'r');
	text(0.8*tmax, 0.7*u, t2, 'r');
	plot([0.85*tmax, tmax], [0.8*u, 0.8*u], 'b');
	plot([0.85*tmax, tmax], [0.7*u, 0.7*u], 'r');
	line([1, 0], 0, 'c');
	line([0, 1], 0, 'c');



function exponencial(tmin, tmax, tau)
	(t, x) = sigExp(tmin, tmax, 0, tau);
	plot(t, x);



function select(a, R, L, C, i0, v0, esc, kvar)
	krange = [0, 5; 0, 5; 0, 2; -2, 2; -2, 2];
	if kvar == 1
		u = (krange(:,2) - krange(:,1))/2;
		krange = [[R, L, C, i0, v0]' - u, [R, L, C, i0, v0]'+u];
	end
	if a == 1
		settabs('v(0) = 00.00   \t');
		slider(sprintf('R = %.2f\nL = %.2f\nC = %.2f\ni(0) = %.2f\nv(0) = %.2f', R, L, C, i0, v0), [R; L; C; i0; v0], krange, '', '', 1);
		settabs('\t\b Escala variable\t');
		button('\tEscala variable                                       Corriente en un circuito RLC serie\n\tVariables ilimitadas', [esc; kvar], 'checkmark', '', 2);
	else
		settabs('v(0) = 00.00   \t');
		slider(sprintf('G = %.2f\nC = %.2f\nL = %.2f\nv(0) = %.2f\ni(0) = %.2f', R, L, C, i0, v0), [R; L; C; i0; v0], krange, '', '', 1);
		settabs('\t\b Escala variable                                                                   \t');
		button('\tEscala variable                                       Tensi—n en un circuito RLC paralelo\n\tVariables ilimitadas', [esc; kvar], 'checkmark', '', 2);
	end



function (R, L, C, i0, v0, esc, kvar) = dragSlider(R, L, C, i0, v0, esc, kvar, id, nb, x1)
	if isempty(id)
		cancel;
	end
	if id == 1
		switch nb
		case 1
			R = round(x1*100)/100;
		case 2
			L = round(x1*100)/100;
		case 3
			C = round(x1*100)/100;
		case 4
			i0 = round(x1*100)/100;
		case 5
			v0 = round(x1*100)/100;
		end
	else
		if nb == 1
			esc = x1;
		else
			kvar = x1;
		end
	end


function (t, x) = sigExp(tmin, tmax, t0, tau)
	n = 400;

	if t0 < tmin
		t = tmin + (tmax - tmin) * (0:n) / n;
		x = exp((t0 - t) / tau);
	elseif t0 > tmax
		t = [tmin, tmax];
		x = [0, 0];
	else
		t = t0+(tmax-t0)*(0:n)/n;
		x = exp((t0-(t0+(tmax-t0)*(0:n)/n))/tau);
	end

@}