%RSCJ equations to do iv with jj.m
% see jj.m for units
% Can  choose up or down ramp or both up and down together
global betac Ib
betac= input('What is beta-c?  ');
v=[]; 
I=[];
init = [0 0];
x0 = init';
% Must choose one of the other for up or down ramp
% for Ib = 0:.1:2;  % for up ramp
% for Ib = 2:-.1:0;  % for down ramp
   for Ib=[0:.1:2, 2:-.1:0]; % for up and down together  
t0 =0; tf=200;
[t,x] = ode45('jj',[t0,tf],x0);
y=x(ceil(.1*length(x)):length(x),1);
t2=t(ceil(.1*length(x)):length(x));
v=[v, trapz(t2,y)/(max(t2)-min(t2))];
I=[I, Ib];
%attempt to have the last value of the previous calculation 
% as the starting point for the next iteration to have hysteresis
x0=[x(size(x,1),:)];
end;
figure;
%v/sqrt(betac) puts V in units of IcRn
plot(v/sqrt(betac),I,'*'); % plot with point
%plot(v/sqrt(betac),I);  % plot with line
%axis([0 max(I) 0 max(I)]);
title(['Junction with beta_c = ', num2str(betac) ]);
xlabel('voltage in units of IcRn');
ylabel('current in units of Ic');