[V,P,t,C]=uppsalator(VOLT,VINI,NUM,Cini,NOPLOT)

function [V,P,t,C]=uppsalator(VOLT,VINI,NUM,Cini,NOPLOT)
%   An integro-differential pde solver (Dirichlet, 1-D), 
%   constant boundary values, uniform mesh.
%
%UPPSALATOR is an electroosmotical oscilLATOR developed in UPPSALA 
%by T. Teorell. Velocity V of liquid flow through a porous membrane, 
%and pressure drop P oscillate in time at constant supercritical 
%current xor voltage values. Up to now oscillations were observed 
%only at fixed currrent. In earlier theories, oscillations at fixed 
%voltage were impossible. This program demonstrates the oscillations
%at supercritical fixed voltage (here VOLT>7.69)
%
% The program solves the following boundary value problem:
% find electrolyte concentration C(T,X) as function of time T>0
% and coordinate X, 0<=X<=1. C(T,X) obeys the equation of convective
% diffusion 
%   pi^2*dCdT=d2CdX2-V*dCdX
%   dCdT and dCdX are time and spatial derivatives of C respectively, 
%   d2CdX2 is the second spatial derivative, pi^2 inherited from  [1]
%
%   Boundary conditions:
%   C(T,0)=CLEF; 0<=CLEF<1  ( not very close to 1, here CLEF=0.1)
%   C(T,1)= 1;
%
%   Initial condition is not restrictive, if we are looking
%   for steady oscillations. Here by default 
%   C(0,X)=linspace(CLEF,1,N), N = mesh size
%
%   Velocity V contains the hydraulic (P) and the electroosmotic 
%   (VOLT*F(T)) components:
%   V= -P + VOLT*F(T)   (*)
%   VOLT = constant voltage > 0;
%   P is the pressure drop across the membrane, its time 
%   derivative dPdT is proportional to velocity V:
%    dPdT=lambda*V; 
%   lambda > 0 empirical constant (here 0.2)
%
%   F(T) = electroosmotical factor, a functional of C defined as 
%   ratio of two integrals over X:
%   F(T)=integral(0,1,1/C(T,X)^1.5 dX)/integral(0,1,1/C(T,X) dX)
%
%   P(0) is found from (*) via V(0) and F(0).
%
%   Theories before [1] neglected time dependency of F(T) in eq.(*), 
%   therefore oscillations at constant voltage were "impossible".
%
%   [1] Pastushenko V.Ph. 1997. Teorells membrane system oscillates at
%   constant voltage. Bioelechtrochem. Bioenerget. 43: 143-150.
%   [2] Pastushenko V.Ph. Uppsalator's acceleration. Electrophoresis, in
%   press.
%
%Call:
%                         uppsalator;
%               [V,P,T,C]=uppsalator(VOLT);
%               [V,P,T,C]=uppsalator(VOLT,VINI);
%               [V,P,T,C]=uppsalator(VOLT,VINI,N);
%               [V,P,T,C]=uppsalator(VOLT,VINI,N,Cini);
%               [V,P,T,C]=uppsalator(VOLT,VINI,N,Cini,NOPLOT);
%Input:
%		VOLT=constant voltage, def. 7
%       VINI = initial (usually nonzero) velocity, default 1.
%		N = mesh size, odd positive integer (default 23):
%       X-coordinate is represented as a mesh  X=linspace(0,1,N)
%       so that the number of unknown concentration values is N-2
%
%       Cini = C(0,X): initial values of concentration, default or empty
%       Cini=linspace(CLEF,1,N) (here CLEF=0.1). However, for slightly 
%       supercritical voltages, where the convergency to limit cycle 
%       is slow, Cini can be selected as the last set of C values 
%       obtained in the previous run of the UPPSALATOR
%
%       NOPLOT: graph output suppressor for repetitive runs:
%       if set, no plotting
%
%Output:
%   variables:
%       V = velocity
%       P = pressure
%		T = time, here linspace(0,50,1500)
%       C = concentration, a matrix length(T) by N
%
%   graphs:
%       phase trajectory V,P;
%       time dependency of P,V,F;
%       C(T,X) - linspace(CLEF,1,N): deviation of concentration
%       from the steady-state profile
%
%   DEMO:
%    uppsalator;  %subcritical voltage: damped oscillations
%    uppsalator(8);    %supercritical voltage: limit cycle
%    uppsalator(20); %far supercritical: relaxational oscillations

%	Vassili Pastushenko	 March-May 2006
%==============================
%Set the UPPSALATOR configuration defined by CLEF and lambda=LAM:
CLEF=.1;  %specify left boundary concentration
LAM=0.2;  %membrane/capillary area, porosity
%===================================
NAR=nargin;
if NAR<1
    VOLT=7;
end

if NAR<2
    VINI=1;
end
 
if NAR<3
    NUM=23;
end
%========================

N=NUM;
pi2=pi^2;
MUL1=(N-1)/2/pi2;
MUL2=(N-1)^2/pi2;
PARS=struct('mul1',MUL1,'mul2',MUL2,'volt',VOLT,'lam',LAM,'num',N);
 
C0=linspace(CLEF,1,N)'; %Steady-state concentration profile
F0=fundyn(C0);%simsum(ic0.^1.5)/simsum(ic0);

if NAR<4
    Cini=C0;
elseif isempty(Cini)
    Cini=C0;     
end

Fini=fundyn(Cini);%simpsum(1./Cini.^1.5)/simpsum(1./Cini);
P0=VOLT*F0; %Pressure at focus
PINI=-VINI+VOLT*Fini;
Yini=[Cini;PINI];
 
OPT=odeset('RelTol',1e-5,'AbsTol',1e-7);
TMAX=50;
[t,Y]=ode15s(@teorders,linspace(0,TMAX,1500),Yini,OPT,PARS);

P=Y(:,N+1);
C=Y(:,1:N);
CI=1./C.';
F=simpsum(CI.^(1.5))./simpsum(CI);
V=-P+VOLT*F';

%Focus P
P=P-P0;
if NAR>4
    return
end

subplot(2,2,1)
%clf
plot(0,0,'k+',V,P,'linewidth',1.5,'markersize',12)
axis tight 
setfont
xlabel('V')
ylabel('P - P_0');
Pini=P(1);
title(['{\fontname{symbol}F = }',num2str(VOLT,3),';   \rm V_{ini} = ', num2str(VINI),';    P_{ini} - P_0 = ', num2str(Pini)])
subplot(2,2,2)
cv=diff(V,2);
in=find(cv(1:end-1).*cv(2:end)<0)+1;
plot(t,V,t,-P,':',t,V+P,'m-.',t(in),V(in),'ro',t,V,'b','linewidth',1.5)
setfont
axis tight
cline(t);
legend('V','-P','F\fontname{symbol}{F}','V inflects','Orientation','horizontal','Location','Northoutside');
xlabel('Time')

%h=cline(1i*P0,'ro');
%set(h,'linewidth',3)

subplot(2,1,2)
%figure(2)
%cla
CC=C(1:3:end,:);
SC=size(CC,1);

DC=CC-repmat(C0',SC,1);
%step=round(size(DC,1)/500);
%DCC=DC(1:step:end,:);

DCC=DC';%interp1(DC',linspace(1,N,500),'spline');

shimspm(DCC,[TMAX 1])
setfont
xlabel('Time');
ylabel('Coordinate');
title('C  -  C_0')
shg


%=========================    TERORDERS
function D=teorders(t,Y,PARS)
N=PARS.num;
C=Y(1:N);
P=Y(N+1);
%CINTER=interpol(C,9);
%CINV=1./CINTER;
F=fundyn(C);%simpsum(CINV.^(1.5))/simpsum(CINV);
V=-P+PARS.volt*F;
D2CDX2=PARS.mul2*diff(C,2);
DCDX=(C(3:N)-C(1:N-2))*PARS.mul1;
D=zeros(N+1,1);
D(2:N-1)=D2CDX2-V*DCDX;
D(N+1)=PARS.lam*V;
%========================================================

function	s=fundyn(y)
%y=interpol(y,9);
Y=1./y;
NN=length(Y);
Y32=Y.^1.5;
s1=(2*sum(Y(3:2:NN-2))+4*sum(Y(2:2:NN-1))+Y(1)+Y(NN));
s32=(2*sum(Y32(3:2:NN-2))+4*sum(Y32(2:2:NN-1))+Y32(1)+Y32(NN));
s=s32/s1;

%=================================      SHIMSPM
function shimspm(Z,S)
[R,C]=size(Z);
X=(0:C-1)/(C-1)*S(1);
Y=(0:R-1)/(R-1)*S(2);
surf(X,Y,Z);
view(0,90);
xlabel('Time')
ylabel('Coordinate')
colormap(bow(2));
h=colorbar;
set(h,'tickdir','out','fontsize',15)
shading interp;
axis tight

%==============================
function  setfont(Fontsize)
if nargin<1,Fontsize=15;end
set(gca,'Fontsize',Fontsize)