% FPEx Program to plot solutions to FPE with k/x drift term
% Part 1: Collection of solute onto electrode surface 'a'
% Part 2: Relaxation of solute from electrode surface 'a'
% D. Bakewell, 12:00 hrs Su August 22 1999
clear all;
desktop;
k=input('Please enter k>=0 DEP constant k = ')
%b=input('Please input upper BC b = ')
%a=input('Please input upper BC a = ')
KbT=1
kKbT=k/KbT
xhia=1 % viscosity constant
nu=[];
al=[];
nu=1/2+kKbT/2
al=(1-kKbT)/2
a=1
XBC=input('mult factr fr upper b=10 BC: ')
b=10*double(XBC)
SSaCo=2*al*a^(-kKbT)/(b^(2*al) - a^(2*al)); % Stdy Stat @ a fr Collctn: p(a,inf)
ICaRe=2*al*a^(-kKbT)/(b^(2*al) - a^(2*al)); % Ini Cond @ a for Relaxation: pr(a,0)
aStr=num2str(a); bStr=num2str(b); kStr=num2str(k); alStr=num2str(al);
SSaCoStr=num2str(SSaCo); ICaReStr=num2str(ICaRe);
fpgm=menu('Program(s) to run','DEP collectn & (relxtn?)','relaxation frm collctn only')
switch fpgm % Either Solute Collection (then) OR Relaxation Prgm
case 1 % ----#############---- DEP Collection Program % ----#############---- 
fcdif=2;
fcalp=2;
if abs(kKbT) < abs(eps)
   fcdif=menu('FPE Drift=0','Exact Comprisn wth Diffusn','No comparisn')
elseif kKbT==2
   fcalp=menu('alpha=1.5, k=2','Verify wth explicit Exprssn','No comparisn')
end % end if
% -------- Beta Eigenvalue determination for DET = 0 --------
reigloop=2;
reigcountr=0;
while reigloop==2; % --- REdetermining eIGenvalues LOOP ---
eigscale=1;   
reigcountr=reigcountr+1;
if reigcountr>=2;
   eigscal=menu('ENLARGE from origin X','1','2','3','4','5','6','7','8','9','10',...
      '11','12','13','14','15','16','17','18','19','20');
   eigscale=double(eigscal)
end % end if reigcountr
frescale=1;
while frescale==1;
   Beta=[];
   Det12AB=[];
	for j=1:1:200;
		be=double(j)/(10*eigscale);
   	Beta(j)=be;
		Det12AB(j)=fpexdet(be,nu,al,a,b);
	end % end computation loop
plot(Beta,Det12AB','r');
hold on;
frescle=input('Press 1 to cntinu, or rescle valu: ')
	if frescle~=1
      eigscale=frescle;
      hold off;
      close;
	elseif frescle==1;
   	frescale=2;
	end % end if
end % end while frescale
nei=input('Enter no. of eigenvalues wanted: ')
errorlim=10^-5
tol=[];
trace=[];
eigBeta=[];
ic=0;
x=0; x1=0; x2=0;
message=['Click mouse PAST first peak ONCE for INITIAL stepping interval']
[x2,y]=ginput(1);
intval=abs(x2);
intvalscale=1;
eigv=0;
em=1;
eigloop = 1;
%if reigcountr>=2;
intvalscale=menu('enter intval scale','1','1/2','1/3','1/4',...
   '1/5','1/6','1/7','1/8','1/9','1/10');
%end % end if reigcountr
x=intval
while eigloop==1;
   eigv=0;
   eigv=fzero('fpexdet',x,tol,trace,nu,al,a,b);
   eigvd=abs(fpexdet(eigv,nu,al,a,b));
   if eigvd <= errorlim;
      uplim=eigv+intval/1000;
      lolim=eigv-intval/1000;
      eigfind=find(eigBeta>lolim); % find eigenvals already near current valu
      if any(eigfind)==0;
         ic=ic+1;
         eigBeta(ic)=eigv;
         veczero(ic)=0;
      elseif any(eigfind)==1;
         message=['Eignvalu same as last => ic=']
      end
   else
      message=['WARNING error limit exceeded: abs(eigv)= ',eigvd]
      dummy=menu('WARNING error limit exceeded','Press to continu');
   end
   eigerror(ic)=eigvd;
   %plot(eigv,0,'b+','Erasemode','none')
   if ic>=nei;
      eigloop=menu('Press to','find more eignvlus','display eignvls & cntinu');
   end %end if
   if em>=200
      eigloop=menu('Press to','find more eignvlus','display eignvls & cntinu');
   end
   em=em+1; % emergency time-out
   x=x+intval/double(intvalscale);
end % end while eigloop
if all(diff(eigBeta))>eps
	message=['eigenvalues are monotoniclly increasing']
else
   message=['WARNING eigenvalues are NOT monotoniclly increasing']
   dummy=menu('WARNING egnvls NOT MONOTNIC','Press to continu');
   diff(eigBeta)
end
plot(eigBeta,veczero,'b+','Erasemode','none')
eigBeta
merror=max(eigerror)
hold on;
xlabel('Beta value');
ylabel('Bessel function Determinant');
title('Bessel Determinant vs Beta');
kbaStr=['DEP vel. is k/x, k = ' kStr ', alpha= ' alStr ...
      ' a = ' aStr  ', b = ' bStr];
gtext(kbaStr);
print -dbitmap
reigloop=menu('p(x,t) dbitmapped','O/P WORD, PASTE, CONT','REPEAT process')
close
end % --- end while reigloop (REdetermining eIGenvalues LOOP) ---
jmax=ic
Cr=[];
Cracomp=[];
for j=1:1:jmax;
	xb=eigBeta(j);
	aB=a*xb;
	bB=b*xb;
	Cr(j)=double(BesselJ(1-nu,aB)/BesselJ(nu-1,aB)); 
   Crdiff(j)=abs(Cr(j)- (BesselJ(1-nu,bB)/BesselJ(nu-1,bB))); 
   % -- Check Eigenvalues satisfy boundary conditions --
	ga=(a^al)*(BesselJ(-nu,aB)+Cr(j)*BesselJ(nu,aB));
   dga1=al*(ga/a)+(a^al)*(xb/2)*(BesselJ(-1-nu,aB)-BesselJ(1-nu,aB));
   dga2=(a^al)*Cr(j)*(xb/2)*(BesselJ(-1+nu,aB)-BesselJ(1+nu,aB)); 
   fluxa(j)=dga1+dga2+((1-2*al)/a)*ga;
   gb=(b^al)*(BesselJ(-nu,bB)+Cr(j)*BesselJ(nu,bB));
   dgb1=al*(gb/b)+(b^al)*(xb/2)*(BesselJ(-1-nu,bB)-BesselJ(1-nu,bB));
   dgb2=(b^al)*Cr(j)*(xb/2)*(BesselJ(-1+nu,bB)-BesselJ(1+nu,bB)); 
   fluxb(j)=dgb1+dgb2+((1-2*al)/b)*gb;
end
format short e
Cr
CrMaxAbsdiff=max(Crdiff)
dummy=menu('            ','Press to continu');
maxAbsfluxa=max(abs(fluxa))
maxAbsfluxb=max(abs(fluxb))
close
% --------- C1 determination ---------
fInCo=menu('Initl Cndtn & mthd calcn:','Delta fn (b+a)/2','Unifrm - Quad8',...
   'Unifrm - Hypg1F2')
xdel=(b+a)/2;
mmax=jmax;
mmaxStr=num2str(mmax);
tol=[];
trace=[];
C1=[];
Int=[];
%switch fInCo;
%case 1 % Delta fn x=5
%fcompm=input('Please enter valu of mth eignvalu fr diagnostic:m = ')
fcompm=1
for m=1:1:mmax;
	aBet=a*eigBeta(m); 
   bBet=b*eigBeta(m);
   switch fInCo; % INitial COndition
      case 1 % IC=Delta fn x=xdel=(b+a)/2
         besarg=xdel*eigBeta(m);
         numtor=double(xdel^(1-al))*(BesselJ(-nu,besarg)+Cr(m)*BesselJ(nu,besarg));
      case 2 % IC = Uniform calctd usng Quad8 verified wth Mathematica Tu/17/8/1999
         Beta=eigBeta(m);
         BesNegNu=quad8('FpexColNegNu',a,b,tol,trace,nu,Beta);
         BesPosNu=quad8('FpexColPosNu',a,b,tol,trace,nu,Beta);
         numtor=double(BesNegNu+Cr(m)*BesPosNu)/(b-a);
      case 3 % Uniform IC = Uniform calctd usng Hypergeometric 1F2 fnctns
         Beta=eigBeta(m);
         hybe=-((bBet)^2)/4;
         hyal=-((aBet)^2)/4;
         BesNegNu=double((2^nu)*(Beta^-nu)/Gamma(1-nu))*(b*Hypg1F2(1/2,3/2,1-nu,hybe)-...
            a*Hypg1F2(1/2,3/2,1-nu,hyal)); %Hypg1F2(a1,b1,b2,z)
         BesPosNu=double((2^-nu)*(Beta^nu)/(Gamma(1+nu)*(2*nu+1)))*((b^(2*nu+1))*...
            Hypg1F2(nu+1/2,nu+3/2,1+nu,hybe)-(a^(2*nu+1))*Hypg1F2(nu+1/2,nu+3/2,1+nu,hyal));
         numtor=double(BesNegNu+Cr(m)*BesPosNu)/(b-a);
   	end % end switch
   Eqn1b=(b^2)/2*(DJBetaX(-nu,bBet)^2+(1-(nu/bBet)^2)*BesselJ(-nu,bBet)^2);
	Eqn23b=Cr(m)*(b^2)*(DJBetaX(nu,bBet)*DJBetaX(-nu,bBet)...
      +(1-(nu/bBet)^2)*BesselJ(nu,bBet)*BesselJ(-nu,bBet));
   Eqn4b=(Cr(m)^2)*(b^2)/2*(DJBetaX(nu,bBet)^2+(1-(nu/bBet)^2)*BesselJ(nu,bBet)^2);
   Eqn1a=(a^2)/2*(DJBetaX(-nu,aBet)^2+(1-(nu/aBet)^2)*BesselJ(-nu,aBet)^2);
   Eqn23a=Cr(m)*(a^2)*(DJBetaX(nu,aBet)*DJBetaX(-nu,aBet)...
            +(1-(nu/aBet)^2)*BesselJ(nu,aBet)*BesselJ(-nu,aBet));
   Eqn4a=(Cr(m)^2)*(a^2)/2*(DJBetaX(nu,aBet)^2+(1-(nu/aBet)^2)*BesselJ(nu,aBet)^2);
   Norm(m)=(Eqn1b+Eqn23b+Eqn4b)-(Eqn1a+Eqn23a+Eqn4a);
   % Norm is not used in computation of C1 - for purposes understndng only
   Int1b=double((b^2)/2)*(DJBetaX(-nu,bBet)^2+(1-(nu/bBet)^2)*BesselJ(-nu,bBet)^2);
   Int23b=double(Cr(m)*(b^2)*(BesselJ(-nu-1,bBet)*BesselJ(nu-1,bBet)...
            + BesselJ(-nu,bBet)*BesselJ(nu,bBet)));
   InX23b=double(Cr(m)*(b^2)*nu/bBet*(BesselJ(-nu,bBet)*DJBetaX(nu,bBet)...
            -BesselJ(nu,bBet)*DJBetaX(-nu,bBet)));
   Int4b=double(Cr(m)^2)*((b^2)/2)*(DJBetaX(nu,bBet)^2+(1-(nu/bBet)^2)*BesselJ(nu,bBet)^2);
   Intb=double(Int1b+Int23b+Int4b);
   Int1a=double((a^2)/2)*(DJBetaX(-nu,aBet)^2+(1-(nu/aBet)^2)*BesselJ(-nu,aBet)^2);
   Int23a=double(Cr(m)*(a^2)*(BesselJ(-nu-1,aBet)*BesselJ(nu-1,aBet)...
            + BesselJ(-nu,aBet)*BesselJ(nu,aBet)));
   InX23a=double(Cr(m)*(a^2)*nu/aBet*(BesselJ(-nu,aBet)*DJBetaX(nu,aBet)...
            -BesselJ(nu,aBet)*DJBetaX(-nu,aBet)));
   Int4a=double((Cr(m)^2)*((a^2)/2)*(DJBetaX(nu,aBet)^2+(1-(nu/aBet)^2)*BesselJ(nu,aBet)^2));
   Inta=double(Int1a+Int23a+Int4a);
   Int(m)=double(Intb-Inta);
   C1(m)=double(numtor/Int(m));
   InXDif(m)=InX23b-InX23a;
   WronskB(m)=bBet*(BesselJ(-nu,bBet)*DJBetaX(nu,bBet)-BesselJ(nu,bBet)*DJBetaX(-nu,bBet));
   WronskA(m)=aBet*(BesselJ(-nu,aBet)*DJBetaX(nu,aBet)-BesselJ(nu,aBet)*DJBetaX(-nu,aBet));
      if m==fcompm
         %Int23bdiffusion=2*Cr(m)*(-(cos(bBet)^2)/(pi*eigBeta(m)^2))
         InX23b
         InX23btest=Cr(m)*(b^2)*nu/bBet*BesselJ(nu-1,bBet)*(BesselJ(-nu,bBet)+...
               Cr(m)*BesselJ(nu,bBet))
         InX23a
         InX23atest=Cr(m)*(a^2)*nu/aBet*BesselJ(nu-1,aBet)*(BesselJ(-nu,aBet)+...
               Cr(m)*BesselJ(nu,aBet))
         WrnskB=bBet*(BesselJ(-nu,bBet)*DJBetaX(nu,bBet)-BesselJ(nu,bBet)*DJBetaX(-nu,bBet))
         WrnskA=aBet*(BesselJ(-nu,aBet)*DJBetaX(nu,aBet)-BesselJ(nu,aBet)*DJBetaX(-nu,aBet))
         InX23WrB=WrnskB*Cr(m)*nu/(eigBeta(m)^2)
         InX23WrA=WrnskA*Cr(m)*nu/(eigBeta(m)^2)
         dummy=menu('            ','Press to continu');
		end % end if diagnostic
   end % end C1 evaluation loop 
%Norm
WronskB
WronskA
NormIntDifMax=max(abs(Int-Norm))
InXDifmax=max(InXDif)
dummy=menu('            ','Press to continu');
format short e
Int
C1
% --------- Summing Fourier-Bessel series ---------
%xcchoice=input('Enter Integer value of xc for sampling = ')
xcchoice=1;
%linpveccol=input('enter plotting line color & type e.g. r--: ','s') % linecolor/type string
tmax=input('Integer no. of time points = ')
timescale=input('timescale factor (e.g. 0.1 to 10) = ');
befos=0;
bg=0;
xcmax=20;
for tc=1:1:tmax
   %tim=(tc-1)*(XBC^2);
   tim=(10^((tc-1)/10)-1)*timescale;
   pvec=[];
   tvec(tc)=tim;
   for xc=1:1:xcmax;
      x=a+(b-a)*(xc-1)/(xcmax-1);
      befos=0;
      bg=0;
      mterm=0;
      fmterm=0;
      for m=1:1:mmax;
         bg=eigBeta(m)*x;
         rhosq=double((eigBeta(m)^2)*KbT/xhia);
         extim=double(exp(-rhosq*tim));
         mterm=double(C1(m)*(x^al)*(double(BesselJ(-nu,bg))+Cr(m)*double(BesselJ(nu,bg)))*extim); 
         befos=double(befos)+double(mterm);
         if xc==xcchoice & tc==1 % sample Bessel-Fourier series @ tc=0
            termC1(m)=C1(m)*(x^al)*BesselJ(-nu,bg);
            termC2(m)=C1(m)*Cr(m)*(x^al)*BesselJ(nu,bg);
            termC3(m)=termC1(m)+termC2(m);
         end % end of check terms
         end % end m Fourier-Bessel sum
      if k==KbT
         Cs=1/log(b/a)
      else
         al2=double(2*al);
         Cs=double(al2/(b^al2-a^al2));
      end
      befos0=double(Cs*x^(-kKbT));
      xvec(xc)=x;
      pvec(xc)=double(befos0)+double(befos);
   end % end position loop
   plot(xvec,pvec,'b') % linpveccol)
   hold on;
   if tc <= 3 % enable viewer to see first 3 profiles slowly
      pause(1);
   end % end slow-view   
   colvec(tc)=pvec(1);
end % end time loop
xlabel('Prpndclr lngth frm elctrde (um)')
ylabel('p(x,t) or [Solute] collected onto elctrde')
title('Collection p(x,t) vs lngth frm elctrde')
kbaStr=['DEP vel. is k/x, k = ' kStr ', alpha= ' alStr ...
      ' a = ' aStr  ', b = ' bStr ]
parStr=['Exact p(a,inf) = ' SSaCoStr ', no. eignvlus = ' mmaxStr];
gtext({kbaStr,parStr});
print -dbitmap
dummy=menu('p(x,t) dbitmapped, OPEN WORD & PASTE','Press to continu');
% --- Fourier sum fr Diffision case - cmpr whn FPE drift k=0 or al=1/2 ---
switch fcdif
case 1
lincompcol=input('enter plotting line color & type e.g. r--: ','s') % linecolor/type string
for tc=1:1:tmax
   tim=(tc-1)*(XBC^2);
   pvec=[];
   for xc=1:1:xcmax;
      x=a+(b-a)*(xc-1)/(xcmax-1);
      mterm=0;
      fos=0; % Fourier sum
      for m=1:1:mmax;
         rhosq=((pi*m/(b-a))^2)*KbT/xhia;
         extim=exp(-rhosq*tim);
         switch fInCo
         	case 1 % Delta fn x=(a+b)/2 case
            mterm=2/(b-a)*cos(m*pi*(a-xdel)/(b-a))*cos(m*pi*(a-x)/(b-a))*extim; 
         	case 2
         end % end switch fInCo
         fos=fos+mterm;
         %if xc==xcchoice & tc==1 % sample Bessel-Fourier series @ tc=0
         %end % end of check terms
      end % end m Fourier sum
      switch fInCo;
         case 1 % Delta fn x=5 case
         	fos0=1/(b-a);
         case 2
         end % end switch fInCo
  		pcompvec(xc)=fos0+fos;
     end % end position loop
   plot(xvec,pcompvec,lincompcol)
   hold on;
   %pause(1)
colcomp(tc)=pcompvec(1);  
end % end time loop
print -dbitmap
dummy=menu('p(x,t) dbitmapped, OPEN WORD & PASTE','Press to continu');
case 2
end % end switch fcdif
% --- Trig-Fourier case k=2, alpha=1.5, for cmprsn wth Fourier-Bessel sum  ---
switch fcalp
case 1 
lincompcol=input('enter al=1.5 pltng lin colr & type e.g. r--: ','s') % linecolor/type string
for tc=1:1:tmax
   %tim=(tc-1)*(XBC^2);
   tim=(10^((tc-1)/10)-1);
   pvec=[];
   for xc=1:1:xcmax;
      x=a+(b-a)*(xc-1)/(xcmax-1);
      mterm=0;
      fos=0; % Fourier sum
      gmx0=0;
      gmx=0;
      for m=1:1:mmax;
         rhosq=((pi*m/(b-a))^2)*KbT/xhia;
         extim=exp(-rhosq*tim);
         switch fInCo;
         case 1 % Delta fn x=(a+b)/2 case cf p150 Book 1
            gmx0=cos(m*pi*(xdel-a)/(a-b)) + (b-a)*sin(m*pi*(xdel-a)/(a-b))/(m*pi*xdel);
            gmx=cos(m*pi*(x-a)/(a-b)) + (b-a)*sin(m*pi*(x-a)/(a-b))/(m*pi*x);
            mterm=2/(b-a)*gmx*gmx0*(xdel/x)*extim; 
         case 2
         	gmxUni=(sin(m*pi/2)/(m*pi/2))^2;
            gmx=cos(m*pi*(x-a)/(a-b)) + (b-a)*sin(m*pi*(x-a)/(a-b))/(m*pi*x);
            mterm=-2/x*gmxUni*gmx*extim; 
         end % end switch fInCo
         fos=fos+mterm;
         %if xc==xcchoice & tc==1 % sample Bessel-Fourier series @ tc=0
         %end % end of check terms
      end % end m Fourier sum
      switch fInCo;
         case 1 % Delta fn x=5 case
         	fos0=a*b/((b-a)*x*x);
         case 2
            fos0=a*b/((b-a)*x*x);
         end % end switch fInCo
   	pcompvec(xc)=fos0+fos;
	end % end position loop
	plot(xvec,pcompvec,lincompcol)
  	hold on
   %pause(1)
   colcomp(tc)=pcompvec(1);
end % end time loop
hold on;   
print -dbitmap
dummy=menu('p(x,t) dbitmapped, OPEN WORD & PASTE','Press to continu');
case 2
end % end switch fcalp
% dummy=menu('            ','Press to continu');
hold off;
close
% ------------- p(t) @ x=a plot & print & store data -------------
message=['At t=']
tvec(1)
message=['Normalised concentration or Probability at electrode edge = ']
colvec(1)
fICcorrct=menu('p(a,0) displayed','plot as is numricly','correct to therticl IC?');
switch fICcorrct;
case 1
case 2
   switch fInCo;
   case 1
      colvec(1)=0;
   case 2
      colvec(1)=1/(b-a)
   end % end fInCo switch
end % end fICcorrct
plot(tvec,colvec,'b--')
xlabel('Time (a.u.)')
ylabel('p(1,t) collection onto electrode @ a')
title('DEP solute collection onto electrode vs time')
kbaStr=['DEP vel. is k/x, k = ' kStr ', alpha= ' alStr ...
      ' a = ' aStr  ', b = ' bStr ]
parStr=['Exact p(a,inf) = ' SSaCoStr ', no. eignvlus = ' mmaxStr];
gtext({kbaStr,parStr});
hold on;
switch fcalp
	case 1
		plot(tvec,colcomp,lincompcol);
   	hold on;
   case 2
end % end switch fcalp
print -dbitmap
fstore=menu('p(x,t) dbitmapped, OPEN WORD/PASTE','store ColVec data/continu','Cntinu ONLY');
switch fstore
case 1
   j=1:1:tmax; % set-up COLlection MATrix for outputting
   colmat(j,1)=tvec(1,j)';
   colmat(j,2)=colvec(1,j)';
   ofcolmat=['FpxCo' aStr '_' kStr '_' bStr '_.dat'] % Output File colmat
   dlmwrite(ofcolmat,colmat, '\t');
case 2
end % end switch fstore
% ------- Summing Fourier-Bessel series LOOP with different eignvalus fr a comprison -------
fBesFoLoop=1; % Flag Bessel-Fourier SUm PLot Loop
while fBesFoLoop==1;
   mchoose=input('Tstng Eignvlu expsn @ x=a enter max m (MUST BE =< nei) = ');
	befos=0; bg=0; 
	for tc=1:1:tmax
   	%tim=(tc-1)*(XBC^2);
      tim=(10^((tc-1)/10)-1)*timescale;
      pvec=[];   tvec(tc)=tim;
   	%for xc=1:1:xcmax;    x=a+(b-a)*(xc-1)/(xcmax-1); 
   	befos=0;      bg=0;      mterm=0;      fmterm=0;
   	x=a; xc=1;
      for m=1:1:mchoose; 
         bg=eigBeta(m)*x; 
         rhosq=double((eigBeta(m)^2)*KbT/xhia);
         extim=double(exp(-rhosq*tim));
         mterm=double(C1(m)*(x^al)*(double(BesselJ(-nu,bg))+Cr(m)*double(BesselJ(nu,bg)))*extim); 
         befos=double(befos)+double(mterm);
      end % end m Fourier-Bessel sum
      if k==KbT
         Cs=1/log(b/a)
      else
         al2=double(2*al);
         Cs=double(al2/(b^al2-a^al2));
      end
      befos0=double(Cs*x^(-kKbT));
      xvec(xc)=x;
      pvec(xc)=double(befos0)+double(befos);
   	%end % end position x loop
   	colLpvec(tc)=pvec(1); % Coll
	end % end time loop
	lincol=input('enter plotting line color & type e.g. r--: ','s') % linecolor/type string
   plot(tvec,colLpvec,lincol)
   fBesFoLoop=menu('Press to - ','Bessel-Fourier usng max m again','print & end prgrm',...
         'Contiue on..')
      switch fBesFoLoop;
      case 2
         hold on;
			print -dbitmap
			dummy=menu('Bessel-Fourier p(t) @ x=a dbitmapped, OPEN WORD & PASTE','Press to cntinu');
		end % end switch fBesFoLoop print
end % end while fBesFoSuPlLoop
switch fcalp % --- Trialing & plotting few m terms for k=2, nu=3/2 case ---
case 1
   switch fInCo;
   case 2
      faprxloop = 1; 
      while faprxloop == 1;
      mchoose=input('nu=3/2, tstng Eignvlu expsn @ x=a enter max m (MUST BE =< nei) = '); 
      for tc=1:1:tmax
         %tim=(tc-1)*(XBC^2);
         tim=(10^((tc-1)/10)-1);
         mterm=0;
         fos=0;
   		for m=1:1:mchoose;
         	rhosq=((pi*m/(b-a))^2)*KbT/xhia;
         	extim=exp(-rhosq*tim);
         	gmaUni=(sin(m*pi/2)/(m*pi/2))^2;
         	gma=1;
            mterm=-2/a*gmaUni*gma*extim;
            fos=fos+mterm;
      	end  % end m loop
      	fos0=a*b/((b-a)*a*a);
      	colapprx(tc)=fos0+fos;
      end % end time loop
      lincol=input('enter plotting line color & type e.g. r--: ','s') % linecolor/type string
      plot(tvec,colapprx,lincol)
      faprxloop=menu('Press to - ','attempt apprx again','print & end prgrm',...
         'end this Collctn prgrm')
      switch faprxloop;
      case 2
         hold on;
			print -dbitmap
			dummy=menu('Approx p(t) @ x=a dbitmapped, OPEN WORD & PASTE','Press to continu');
		end % end switch faprxloop
   end % end while loop
   case 1 % case 1 fInco
   end % end fInco switch
   case 2 % case 2 fcalp 
end % end fcalp switch
hold off;
close
close
%frelax=menu('DEP collctn cmplt...','Press fr diffsn relaxtn','Press to end prgm')
%fpgm=frelax+1 % change state of fpgm
case 2 % switch fpgm=2 ----##---- Relaxation (from DEP Collection) Program % ----##----
% --- Diffusion relaxation Fourier sum (FPE drift k=0 or al=1/2) ---
%linrelax=input('enter plotting line color & type e.g. r--: ','s') % linecolor/type string
rint(1)=2;rint(2)=4;rint(3)=6;rint(4)=8;rint(5)=10;
testk=find(k==rint);
if testk>0 % for k = 2, 4, 6, 8, 10......
frlxK=menu('DEP k={2,4,..} - cmputatn choice','fast n-trig','Quad8/Gamma')
else
frlxK=2;
end % end testk
tvec=[];
relvec=[];
Smkab=[];
Cmkab=[];
xcmax=20;
tmax=input('Integer no. of time points = ');
mmax=input('Integer no. of eigenvalues = ');
mmaxStr=num2str(mmax);
switch frlxK
	case 1
      n=menu('plse confrm/chng k =','2','4','6','8','10')% n=k/2
		al=double(1/2-n);   
		for m=1:1:mmax; % pre-computation: setup Smkab, Cmkab vectors
   		del=pi*m/(b-a);
			Smkab(m)=FpexRlxKIntSin(b,n,del)-FpexRlxKIntSin(a,n,del);
   		Cmkab(m)=FpexRlxKIntCos(b,n,del)-FpexRlxKIntCos(a,n,del);
      end
   case 2
      fSCvec=menu('Smkab, Cmkab vectors','Numrcl Quad8','GammaInc')
      switch fSCvec;
      case 1;
         tol=[];
         trace=[];
         for m=1:1:mmax; % pre-computation: setup Smkab, Cmkab vectors
         	del=pi*m/(b-a);
            Smkab(m)=quad8('FpexRlxNumSin',a,b,tol,trace,al,del);
            Cmkab(m)=quad8('FpexRlxNumCos',a,b,tol,trace,al,del);
         end
      case 2;
      	for m=1:1:mmax; % pre-computation: setup Smkab, Cmkab vectors
         	%del=pi*m/(b-a);
         	ComIncGam2alA=gamma(2*al)*(1.-gammainc(-i*a*pi*m/(b-a),2*al));
         	ComIncGam2alB=gamma(2*al)*(1.-gammainc(-i*b*pi*m/(b-a),2*al));
         	exgadif=exp(i*al*pi)*(ComIncGam2alA-ComIncGam2alB);
         	Smkab(m)=FpexRlxGamSin(al,del,a,b);
   			Cmkab(m)=FpexRlxGamCos(al,del,a,b);
         end
      end % fSCvec switch
  	end % end frlx switch
frelaxloop=1;
while frelaxloop==1;
timescale=input('timescale factor (e.g. 0.1 to 10) = ')   
for tc=1:1:tmax
   %tim=(tc-1)*timescale;
   tim=(10^((tc-1)/10)-1)*timescale;
   pvec=[];
   tvec(tc)=tim;
   for xc=1:1:xcmax;
      x=a+(b-a)*(xc-1)/(xcmax-1);
      mterm=0;
      fos=0; % Fourier sum
      for m=1:1:mmax;
         rhosq=((pi*m/(b-a))^2)*KbT/xhia;
         extim=exp(-rhosq*tim);
         del=pi*m/(b-a);
         mterm=(cos(pi*a*m/(b-a))*Cmkab(m) + sin(pi*a*m/(b-a))*Smkab(m))...
            *cos(m*pi*(x-a)/(b-a))*extim; 
         fos=fos+mterm;
      end % end m Fourier sum
      fos0=1/(b-a);
      pvec(xc)=fos0+(2/(b-a))*(2*al/(b^(2*al)-a^(2*al)))*fos;
      xvec(xc)=x;
   end % end position loop
   plot(xvec,pvec,'r')% linrelax)
   hold on;
   %pause(1)
	relvec(tc)=pvec(1);  
end % end time loop
xlabel('Prpndclr lngth frm elctrde (um)')
ylabel('p(x,t) or [Solute] remaining on elctrde')
title('Relaxation p(x,t) vs lngth frm elctrde')
kbaStr=['IC value: k = ' kStr ', alpha= ' alStr ...
      ' a = ' aStr  ', b = ' bStr ]
parStr=['Exact p(a,0) = ' ICaReStr ', no. eignvlus = ' mmaxStr];
gtext({kbaStr,parStr});
print -dbitmap
dummy=menu('p(x,t) dbitmapped, OPEN WORD & PASTE','Press to continu');
close;
fICcorrct=menu('p(a,0) displayed','plot as is numricly','correct to therticl IC?');
switch fICcorrct;
	case 1
	case 2
   	relvec(1)=2*al*(a^double(2*al-1))/(b^(2*al)-a^(2*al));
end % end fICcorrct
plot(tvec,relvec,'r'); %linrelax);
xlabel('Time (a.u.)')
ylabel('Pr(1,t) Solute Relaxation or diffusion from electrode @ a')
title('Relaxation of collected solute vs time (au)')
kbaStr=['IC value: k = ' kStr ', alpha= ' alStr ...
      ' a = ' aStr  ', b = ' bStr ]
parStr=['Exact p(a,0) = ' ICaReStr ', no. eignvlus = ' mmaxStr];
gtext({kbaStr,parStr});
print -dbitmap
fstore=menu('p(x,t) dbitmapped, OPEN WORD/PASTE','store colvec & cntinu','Cntinu ONLY');
switch fstore
case 1
   j=1:1:tmax; % set-up COLlection MATrix for outputting
   relmat(j,1)=tvec(1,j)';
   relmat(j,2)=relvec(1,j)';
   ofrelmat=['FpxRe' aStr '_' kStr '_' bStr '_.dat'] % Output File colmat
   dlmwrite(ofrelmat,relmat, '\t');
case 2
end % end switch fstore
frelaxloop=menu('Press to - ','different timescale','test for diff eignvlus',...
   'end this relaxation prgrm')
end % end frelaxloop
% ---- Summing Fourier series LOOP with different eignvalus fr a comprison ----
fFoLoop=frelaxloop-1; % Flag Bessel-Fourier SUm PLot Loop
while fFoLoop==1;
mchoose=input('Tstng Eignvlu expsn @ x=a enter max m (MUST BE =< nei) = ');
%timescale=input('timescale factor (e.g. 0.1 to 10) = ')   
for tc=1:1:tmax
   tim=(10^((tc-1)/10)-1)*timescale;
   pvec=[];  tvec(tc)=tim;
   %for xc=1:1:xcmax;  x=a+(b-a)*(xc-1)/(xcmax-1);
   x=a; xc=1; mterm=0; fos=0; % Fourier sum
      for m=1:1:mchoose;
         rhosq=((pi*m/(b-a))^2)*KbT/xhia;
         extim=exp(-rhosq*tim);
         del=pi*m/(b-a);
         mterm=(cos(pi*a*m/(b-a))*Cmkab(m) + sin(pi*a*m/(b-a))*Smkab(m))...
            *cos(m*pi*(x-a)/(b-a))*extim; 
         fos=fos+mterm;
      end % end m Fourier sum
      fos0=1/(b-a);
      pvec(xc)=fos0+(2/(b-a))*(2*al/(b^(2*al)-a^(2*al)))*fos;
      xvec(xc)=x;
   %end % end position loop
   %plot(xvec,pvec,'r')% linrelax);   hold on;    %pause(1)
	relLpvec(tc)=pvec(1);  
end % end time loop
hold on;
lincol=input('enter plotting line color & type e.g. r--: ','s') % linecolor/type string
plot(tvec,relLpvec,lincol)
fFoLoop=menu('Press to - ','Fourier usng max m again','print & end prgrm',...
         'End prgrm')
switch fFoLoop;
	case 2
   	hold on;
		print -dbitmap
		dummy=menu('Fourier p(t) @ x=a dbitmapped, OPEN WORD & PASTE','Press to cntinu');
	end % end switch fFoLoop print
end % end while fFoLoop
end % end switch fpgm
close;
%    query isn't it gamma(A,X)? or "A must be real"