% TBI (traumatic brain injury block for bistable equation
%This code calculates the critical relationship for A_r vs Y (the
%  length for propagation block following TBI)
% for a function f(u) of your choice
function block
global alp Ar un F_int U1
 
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7);

cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;
 
alplist = [0.25,0.3,0.4];  %This is the list of alpha values to examine
for ja = 1:3
alp = alplist(ja);
%The three zeros of the bistable function are:
U0 = 0;
U1 = alp;
U2 = 1;

N = 400; %This determines the accuracy for the integral of f
du = (U2-U0)/N;
un=[U0:du:U2];
% calculate the integral of f
 F_int = cumsum(feval(un))*du;  %feval is your bistable function of choice

 %find the zero of F_int using bisection
 % guess a lower bound
 a = U0+1.e-5;
 % guess an upper bound
 b=U2;
 umx = bisect(a,b,@fint);  %This is the zero of F_int
  
%The two basic trajectories
Fa = fint(U1);
% the upper trajectory
u1=un;
F1 = fint(U2);
w1 = sqrt(2*(F1-fint(u1)));
 w1s=sqrt(2*(F1-Fa));
 
  % discretize the (U0,umx) interval
  %the lower trqjectory
  M = 50; %
 du0 = (umx-U0)/M;
u0 = [U0:du0:umx];
 Fw = fint(u0);
w0 = sqrt(-2*Fw);
 

MN = 200;
% discretize the (U0,U1) interval for starting values Us
dua = (U1-U0)/MN;
F1 = fint(U2);
Fa = fint(U1);
Uj = [U0+dua:dua:U1-dua];  %These are the starting values for the intermediate trajectory
for j = 1:length(Uj)
    Us = Uj(j);
    
 % find the starting value for W
Ws0 = sqrt(-2*fint(Us));

Ar =(F1-Fa-Ws0^2/2)/(fint(Us)-Fa); % the critical value of As
A(j) = Ar;

tspan=[0:1:200];
s = [Us Ws0];
% integrate the intermediate curve
[T,S] = ode23s(@deRHS,tspan,s,odeset('Events',@events)); 

Y(j) = T(end);
end

[miny,minj]= min(Y);
minA = A(minj);
 
figure(1)
plot(Y(1:minj),A(1:minj),Y(minj:end),A(minj:end),'--',miny,minA,'k*', ...
     'linewidth',2)
axis([0 100 0 1000])
ylabel('A_r','fontsize',20)
xlabel('Y','fontsize',20)
axis([0 50 0 50])
text(8,4.8,'\alpha = 0.4','fontsize',20)
text(18,12,'\alpha = 0.3','fontsize',20)
text(35,21.5,'\alpha = 0.25','fontsize',20)
 text(30,35,'blocked waves', 'fontsize',20)
hold on 
% for the minimal Y value, make a phase portrait
if(cmykflg==1)
  print('../../figs_c/chapt_8/tbi_crit_crv','-deps','-cmyk')
   end 

% Pick a starting value for Us < U1
Us = Uj(minj);  
% find the corresponding starting value for W
Ws0 = sqrt(-2*fint(Us));
Ar =(F1-Fa-Ws0^2/2)/(fint(Us)-Fa);

Ks = Ws0^2/2+Ar*fint(Us)
 Fmsk = (Ks-Ar*fint(u0) >= 0);
  Ws = sqrt(2* (Ks-Ar*fint(u0(Fmsk))));
 
tspan=[0:.1:50];  % hopefully this time interval is sufficiently long
s = [Us Ws0];
% integrate the intermediate curve
[T,S] = ode23s(@deRHS,tspan,s,odeset('Events',@events)); 
figure(1+ja)
plot(u0,w0,u1,w1 ,u0(Fmsk),Ws,'k--',umx,0,'k*', Us,Ws0,'k*',  ...
    U1,w1s,'k*',[U1,U1],[0,0.3],'--','linewidth',2)
xlabel('U','fontsize',20)
ylabel('W','fontsize',20)
text(0.12,0.035,'(U_0,W_0)','fontsize',18)

% clean this up:
if(ja==2)
 annotation('arrow',[.22 .24],[.53,.63])
   annotation('arrow',[.56 .63],[.75,.67])
   annotation('arrow',[.14 .18],[.17,.24])
  
    if(cmykflg==1)
  print('../../figs_c/chapt_8/tbi_pp','-deps','-cmyk')
   end 
end
end
hold off




function s_prime=deRHS(t,s)  % this is the differential equation for the intermediate section

global  Ar
 u = s(1);
 w = s(2);
 Fu = w/Ar;
 Fw = - feval(u);
 s_prime = [Fu;Fw];
    
% here is where to specify the function f(u);
function f=feval(u);
global alp
f = u.*(1-u).*(u-alp);

        
function fI = fint(u)
global un F_int

%Calculate the integral of f
% may be a vector
 
for j = 1:length(u)
    
ndx = min(find(u(j)<=un));
%linear interpolation
if(ndx==1)
    
    fI(j) = 0;
else
u1 = un(ndx-1);
f1 = F_int(ndx-1);
u2 = un(ndx);
f2 = F_int(ndx);
fI(j) = ((u(j)-u1)*f2-(u(j)-u2)*f1)/(u2-u1);
end
end

% event function
function [value,isterminal,direction]=events(t,y)
global U1

value = y(1)-U1; 
isterminal = 1;
direction = 0;

% reminder how to use event functions
%[T,S] = ode23s(@deRHS,tspan,s,odeset('Events',@events)); 

function root = bisect(a,b,feval)

ul = a;
fl = feval(ul);
uu = b;
fu = feval(uu);

% we make the assumption, without checking, that 
% fu*fl<0

% if not, the algoritm fails to find a root.

N = 20;  % number of iterates
% the main bisection algorithm
for j = 1:N
uc = (ul+uu)/2;
fc = feval(uc);
ftest = (fc.*fl>=0);
ul = ftest.*uc+(1-ftest).*ul;
fl = ftest.*fc + (1-ftest).*fl;

uu = (1-ftest).*uc+ ftest.*uu;
fu = (1-ftest).*fc + ftest.*fu;
end
root = uc;

