%RSCJ equations, b_c(Rs) is used, and V/IcRs is normalization
%
global betac Ib betaL Rsg
betac= input('What is beta-c?  ');
betaL=  input('What is beta-L?  ');
Ib = input('What is the current?  ');
Rn= input('What is the Resistance Rn/Rs?  ');
Rsg= input('What is the Resistance Rsg/Rs?  ');
init = input('initial starting point   ');
t0 =0; tf=200;
x0 = [init Ib]';
[t,x] = ode45('whaneq',t0,tf,x0);
figure;
plot(t,x(:,2));
title(['voltage vs time']);
pause;
figure;
plot(t,sin(x(:,1)));
title(['phase vs time']);
pause;
%figure;
%plot(t,(x(:,1)));
%
%to have the phase portrait
%pause;
figure;
plot(rem(x(:,1),2*pi),x(:,2),'.');
%axis([-4 4 -1 1]);
title(['Junction with beta_c = ', num2str(betac), 'and I= ', num2str(Ib)]);
xlabel('theta');
ylabel('theta-dot');
