% Prova d'esame unicz del 2011-01-10
function q

clear,  format short g,  format compact

[ST,I] = dbstack;
%fprintf(2,['$$$$ ' num2str(100*rand) ' $$$$\n'])
fprintf(2,'$$$$ %0.5g $$$$ %s\n', 100*rand,ST.file)

%%
dispah(' Quesito 1 ')
clear

E = 40;    J  = 10;
R1 = 80;   R2 = 80;   R6 = 80;   R7 = 80;
R3 = 160;  R4 = 160;  R5 = 160;

I7cc = E/(R7+paral(R6,R1));
I1cc = I7cc*R6/(R6+R1);
I3cc = J*R4/(R4+R3);
Jcc  = I1cc+I3cc;
disp(table(I7cc,I1cc,I3cc,Jcc))

R0a  = R1+paral(R6,R7);
R0b  = R3+R4;
R0   = paral(R0a,R0b);

R02  = paral(R0,R2);
VAB  = R02*Jcc;
P2   = VAB^2/R2;
disp(table(R0a,R0b,R0,R02,VAB,P2))

%%
disp(' ')
dispah(' Quesito 2 ')

clear

EM = 40;  omega = 100;
R = 10;   L = 100e-3;  C = 500e-6;

dispah(' t < 0 ')
V0 = 0,  I0 = 0

dispah(' t > 0 ')
sigma = 1/2/R/C           ;  omegar = 1/sqrt(L*C);
delta = sigma^2 - omegar^2;  omegad = sqrt(-delta);
tau = 1/sigma;
disp(table(sigma,tau,omegar,delta,omegad))

% regime
% scegliamo i fasori nella rappresentazione con il seno
E   = EM;
ZL  = 1i*omega*L;
ZC  = 1/(1i*omega*C);
den = 1 - (omegar/omega)^2 - 1i*2*sigma/omega; % vedi dispensa
VRr = E/den;
VRr = (E/ZC)/(1/R+1/ZL+1/ZC); % Milllmann, ovviamente stesso risultato
VRr_mod   = abs(VRr);
alfar     = angle(VRr);
alfar_deg = alfar*180/pi;
disp(table(VRr,VRr_mod,alfar,alfar_deg))

% condizioni iniziali
dvdt0p  = 0;
dedt0p  = omega*EM;
dvRdt0p = dvdt0p+dedt0p;

B   = -imag(VRr);
A   = (dvRdt0p + sigma*B - omega*real(VRr)) / omegad;
K   = abs(A+1i*B);
Phi = angle(A+1i*B);
disp(table(dvRdt0p,B,A,K,Phi))

%% plot vR(t) e pR(t)
T = 2*pi/omega;  Td = 2*pi/omegad;  Tc = max([tau,T,Td]);
disp(' tempi caratteristici per asse tempi grafici ')
Tmax = 3*Tc; % 3 tempi caratteristici invece di 4 o 5 perche' siamo
             % interessati principalmente alla potenza, che ha tempi
             % caratteristici dimezzati
disp(table(tau,T,Td,Tc))
N   = 1000;
t   = linspace(0,Tmax,N);                 % vettore con N componenti
vRt = K*exp(-sigma*t).*sin(omegad*t+Phi); % nota il .*
vRr = VRr_mod*sin(omega*t+alfar);
vR  = vRt+vRr;
pR  = vR.^2/R;                            % nota il .^
subplot(3,1,1)
plot(t,vRt)
grid,  title('risposta transitoria'),  xlabel('t'),  ylabel('v_{Rt}(t)')
subplot(3,1,2)
plot(t,vR)
grid,  title('soluzione completa'),  xlabel('t'),  ylabel('v_R(t)')
subplot(3,1,3)
plot(t,pR)
grid,  title('potenza istantanea'),  xlabel('t'),  ylabel('p_R(t)')


function c=paral(a,b)
c=a.*b./(a+b);

function dispah(stringa)
disp(sprintf('<a href="">%s</a>', stringa))