%RSCJ equations to do iv, with non linear subgap resistance
%
% v is in units of IcRn
% t is in units of t_J
% so that this is the Method 1 formulation 
% from jjnlin units
%
global betac Ib
betac= input('What is beta-c?  ');
Rsg =input('What is the ratio of Rsg/Rn?   ')
v=[]; 
I=[];
init = [eps eps];
x0 = init';
% for Ib = -eps:.1:2; % for going up
 for Ib = 2:-.1:-eps; % for going down
t0 =0; tf=200;
[t,x] = ode45('jjnlin',t0,tf,x0);
y=x(.1*size(x,1):size(x,1),1);
t2=t(.1*size(x,1):size(x,1));
v=[v, trapz(t2,y)/(max(t2)-min(t2))];
I=[I, Ib];
x0=[y(size(y,1)),x(size(x,1),2)]';
end;
figure;
% V is in units of IcRn
plot(v,I,'*');
%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');

