%% PROBLEM 6_1 ALL MATLAB
%% Calls P6_1, which calls P6_1fun
%% Carol Lucas

gm=.05;cm=1;gc=.1;I1=0;I2=0;

A1=[-(gm+gc)/cm gc/cm];
A2=[gc/cm  -(gm+gc)/cm];
A=[A1;A2]
F=[I1/cm;I2/cm]
init=[1;0]
ss=-inv(A)*F
ssv=[ss;[0;0]]

lam=eig(A)
dt0=A*init+F

M=[1 1 0 0;0 0 1 1; lam.' 0 0; 0 0 lam.'];  %% Want non conjugate transpose
b=[init; dt0];
as=inv(M)*(b-ssv)

t=[0:.01:60];
x=as(1)*exp(lam(1)*t)+as(2)*exp(lam(2)*t)+ss(1,1);
y=as(3)*exp(lam(1)*t)+as(4)*exp(lam(2)*t)+ss(2,1);



[t,xs]=ode45('P6_1fun',t,init);
v1=xs(:,1);
v2=xs(:,2);


figure(1); clf;plot(t,x,t,y,t,v1,t, v2);xlabel('time (secs)');
ylabel('V1 and V2');title('Original Simulation');
mess='Hit return to continue.'
pause;
ttest=floor(length(t)/2);

%% Test system time constant

V1slope=(log(v1(end))-log(v1(ttest)))/(t(end)-t(ttest))
V2slope=(log(v2(end))-log(v2(ttest)))/(t(end)-t(ttest))

Vs=(v1+v2)/2;      Vd=(v1-v2);

figure(2);clf;subplot(211);plot(t,Vs,t,Vd);title('Vs and Vd Plots');
ylabel('Vs and Vd regular');
subplot(212);semilogy(t,Vs,t,Vd);xlabel('time (secs)');
ylabel('Vs and Vd semilog');title('Vs and Vd');
mess='Hit return to continue.'
pause;

Vsslope=(log(Vs(end))-log(Vs(1)))/(t(end)-t(1))
Vdslope=(log(Vd(end))-log(Vd(1)))/(t(end)-t(1))

AVs=[-gc/(2*cm)], fVS=[(I1+I2)/(2*cm)]
AVd=[-(2*gc+gm)/cm], fVd=[(I1-I2)/cm];

tauVs=1/AVs, tauVd=1/AVd