%RSCJ equations for a single junction and plots
%use with jj.m 
% see jj.m for units of voltage and time 
%
%
global betac Ib
betac= input('What is beta-c?  ');
Ib = input('What is the current?  ');
init = input('initial starting point   ');
t0 =0; tf=200;
x0 = init';
% ode23 gave too much noise
%
[t,x] = ode45('jj',[t0,tf],x0);
figure;
% To plot the phase-dot (voltage) verus time
subplot(311); plot(t,x(:,1));
title(['voltage vs time']);
% To plot the phase verus time
subplot(312);plot(t,sin(x(:,2)));
title(['phase vs time']);
%to have the phase portrait
%pause;
subplot(313);plot(rem(x(:,2),2*pi),x(:,1),'.');
%axis([0 7 -.5  2.5]);
title(['Junction with beta_c = ', num2str(betac), 'and I= ', num2str(Ib)]);
xlabel('theta');
ylabel('theta-dot');