% code to the equation w_t + WW_x=0 using
% method of lines and Method of characteristics

function gelation
global N  dx C0 k t

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;

C0 = 2/3;
L = 1; % length of the interval
k = 3 ; % number of binding sites on a monomer 
N = 101;  % number of grid points
dx = L/(N-1);
x =[0:N-1]*dx;
% first find the numerical solution with MOL

w = Wofz(x);  % initial data

tstep = 0.1;
t_end = 1.5;
tspan = [0:tstep:t_end];
    
[T,S] = ode23s(@deRHS,tspan,w);
figure(1)
plot(x,S','linewidth',2)
xlabel('z','fontsize',20)
ylabel('W(z,t)','fontsize',20)
 
if(cmykflg==1)
  print('../../figs_c/chapt_9/gel_by_mol','-deps','-cmyk')
end 
   
% Now find the exact solution using method of characteristics
for j = 1:length(T)
    t = T(j);
z = Wofz(x)*t+x;

% find W(1):
x_end = z0(t);  % find the value of x for which z = 1
 W1 =   Wofz(x_end);
 W1s(j) = W1;
 ts(j) = t;

 if (j==1)
figure(2)
plot(z,w,'--','linewidth',2)
hold on
 else
     figure(2)
plot(z,w,'linewidth',2)
set(gca,'Xtick',0:0.25:1.5)
xlabel('z','fontsize',20)
ylabel('W(z,t)','fontsize',20)
 end
 if(cmykflg==1)
  print('../../figs_c/chapt_9/gel_by_char','-deps','-cmyk')
end 
 z1 = find(z<=1);
 
 
figure(5)
plot(z(z1),w(z1),'linewidth',2)

hold on
plot(1,W1,'k*','linewidth',3)
xlabel('z','fontsize',20)
ylabel('R_g','fontsize',20)
set(gca,'Xtick',0:0.25:1.5)
end
%plot([1,1],[0,0.5],'--','linewidth',2)
axis([0 1.5 0 0.5]) 

if(cmykflg==1)
  print('../../figs_c/chapt_9/gel_dyn_1','-deps','-cmyk')
end 

 
% plot more features of the solution
N = 21;
dx = L/(N-1);

x =[0:N-1]*dx;
w = Wofz(x);
wp = C0*k*(1-(k-1)*x.^(k-2));
z = x-w./wp;
t = -1./wp;

for j = 1:length(x);
figure(3)
plot([x(j),w(j)*t_end+x(j)],[0,t_end],'linewidth',1)
xlabel('z','fontsize',20)
ylabel('t','fontsize',20)
hold on
end
% use the resultant
t = [0.01:.01: 1.5];
z = (9*C0^2*t.^2 + 6*C0*t + 1)./(12*t*C0)

plot(z,t,'--','linewidth',3)
axis( [0 1.2 0 1.5])
set(gca,'YTick',0:0.25:1.5);
 hold off
 
 if(cmykflg==1)
  print('../../figs_c/chapt_9/char_for_gel','-deps','-cmyk')
 end 


% now calculate monomer in gel vs monomer in polymer
% using the characteristic solution
for j = 1:50
 
    t = (j-1)*0.05;
    tj(j)  = t;
   
x_end = z0(t);    
%  use this to find the integral of W

tspan = [0,x_end];
    
[Z,WI] = ode23s(@deRHSw,tspan,0);
Wint(j) = WI(end);
x0(j) = x_end;
W1(j) = Wofz(x_end);
    end
  R = 3*C0./(3*C0*tj+1);  

figure(4)
plot( tj,W1,'g',ts,W1s,'k*','linewidth',2)
xlabel('t','fontsize',20)
ylabel('W(1,t)','fontsize',20)

if(cmykflg==1)
  print('../../figs_c/chapt_9/gel_dyn_2','-deps','-cmyk')
end 
figure(6)
plot( tj,W1,'g',tj,(2*Wint-W1)/(k-2),'r',...
    tj,R-W1,'k',tj,R,'k--', 'linewidth',2)
xlabel('t','fontsize',20)
legend('R_g','M_1','R_s','R', 'fontsize',20)
if(cmykflg==1)
  print('../../figs_c/chapt_9/gel_dyn_3','-deps','-cmyk')
end 
end

function s_prime=deRHS(t,w)
global N dx
W = w(2:N-1); 
Fw = -[0;W.*(w(3:N)+W)-w(1:N-2).*(W+w(1:N-2));0]/(4*dx);
 
s_prime = [Fw];

end 

function Wz = Wofz(x)
global C0 k t
Wz = C0*k*(x-x.^(k-1));
 
end

function w_prime=deRHSw(x,q)
global C0 k t
  
z = Wofz(x)*t+x;
w = Wofz(x);
wp = C0*k*(1-(k-1)*x.^(k-2))*t +1;
w_prime = [w*wp];

end 

function xend = z0(dummy)
global C0 k t
% find the z0 value for which z = 1 using bisection

 if (t<= 1/(C0*k*(k-2)))
        xend = 1;
    else
        %find x_end<1 by bisection

        xa = 0;
fa = -1;
xb = 0.99999;
fb = Wofz(xb)*t+xb-1;
 for n = 1:20
    xc = (xa+xb)/2;
    fc = Wofz(xc)*t+xc-1;
   
    if(fb*fc>0)
        fb=fc;
        xb=xc;
    else
        fa=fc;
        xa=xc;
    end
    xend = xc;
 end
 
 end
    
end




