% ORUH.m program plots solutions to Ornstein Uhlenbech equation/process
% formerly called DEPDIF developed by D. Bakewell Tu/18/5/1999
% latest upgrade 12:50 Fr/27/10/2000
clear all; desktop;
line1 = 'Dave`s "ORUH" - prgrm to plot Ornstein-Uhlenbeck particle trnsport  ';
line2 = 'fr investgtion DEP accum. onto  INTerdigitated electrode surfaces.  ';
line3 = '& ch5 PhD thesis                                                    '; 
line4 = 'Please take note of the MENU BOX on the UPPER LEFT of the screen.   ';
display1 = [line1;line2;line3;line4]; 
figure; axis off; text(0,0.5,display1)
fchoice = menu('enter choice of Eqn','O-U IC=DELTA fn','O-U IC=RECT fn'); close
switch fchoice
case 1 % --- Ornstein - Uhlenbeck FPE with Initial Condition = Delta(x) ---
   k=10;    D=1;    Xo=1;    tmax=10;    tes=5; % time exponential starting point
   xmax=41; xscale=20; figure;
   for tc=1:1:tmax;
      t=10^(tc-tes); Xav=Xo*exp(-k*t); Xvar=D/(2*k)*(1-exp(-2*k*t))
		pvec=[];
      for xc=1:1:xmax; x=(xc-1)/xscale;
         p=exp(-((x-Xav)^2)/(2*Xvar))/sqrt(2*pi*Xvar);
         pvec(xc)=p; xvec(xc)=x;
      end
   	plot(xvec,pvec,'r:'); hold on; pause(1);
   end
   pvec=[]
   for xc=1:1:xmax; % --- plot at equilibrium ---
         x=(xc-1)/xscale; Xav=0; Xvar=D/(2.*k);
         p=1/sqrt(2*pi*Xvar)*exp(-((x-Xav)^2)/(2*Xvar));
         pvec(xc)=p; xvec(xc)=x;
   end
   plot(xvec,pvec,'b');
   dummy=menu('            ','Press to continue');
case 2 % --- Ornstein - Uhlenbeck FPE with Initial Condition = RECT(x-L/2)/L ---
   fplot=menu('Plot(s) of solution','component parts','entire solution');
   k=20; D=1; tmax=5; tscale=50; tes=5; % time exp. starting value
   xmax=305; Xo=100; xscale=200; L=1; figure;
   for tc=1:1:tmax;
      t=10^(tc-tes)
      %t=tc/tscale;
      A=D*(1-exp(-2*k*t))/(4*k);  B=k*t; C=L*exp(-k*t); pvec=[]; psum=0;
      for xc=1:1:xmax;
         x=(xc-Xo-1)/xscale;
         %Xav=Xo*exp(-k*t);
         %Xvar=D/(2.*k)*(1-exp(-2.*k*t));
         xd=x/(2*sqrt(A)); xdC=(x-C)/(2*sqrt(A));
         y1vec(xc)= sqrt(pi/A)*exp(-(xd^2));  % *exp(B)/L;
         y2vec(xc)=-sqrt(pi/A)*exp(-(xdC^2)); %*exp(B)/L;
         %yvec(xc)=exp(B)/(2*L)*(exp(-(xd^2)) - exp(-((xd-C)^2)));
         p=exp(B)/(2*L)*(double(erf(xd))-double(erf(xdC)));
         pvec(xc)=p; xvec(xc)=x; pmat(xc,tc)=double(p); psurf(xc,tc)=p; 
      end
      switch fplot
      case 1
      	plot(xvec,y1vec,'r:'); hold on;
      	dummy=menu('            ','Press to continue');
      	plot(xvec,y2vec,'b:');
			%plot(xvec,yvec,'k:');
      	dummy=menu('            ','Press to continue');
      case 2
      plot(xvec,pvec,'r:'); hold on;
      %dummy=menu('            ','Press to continue');
      end % end switch fplot
   end % end for tc
   dummy=menu('            ','Press to continue'); close; surf(psurf)
   hold on; dummy=menu('            ','Press to continue'); close
   for xc=1:1:xmax; % --- plot at equilibrium ---
         x=(xc-1)/xscale; Xav=0; Xvar=D/(2.*k);
      	p=1/sqrt(2*pi*Xvar)*exp(-((x-Xav)^2)/(2*Xvar));
         pevec(xc)=p; xvec(xc)=x;
   end
   plot(xvec,pevec,'b:'); hold on;
   dummy=menu('            ','Press to continue'); close;
   fratio=menu('ratio test?','No','Yes')
   switch fratio
   case 1
   	for xc=1:1:xmax;
     	 rvec(xc)=pvec(xc)/pevec(xc);
   	end
   	figure; plot(xvec,rvec,'g:')
      dummy=menu('            ','Press to continue')
   case 2
   end
   hold off; close; close
   fCt=menu('Plot time evolution?','Yes','No')
   sxmax=num2str(xmax)
   message=['Xmax is currently at: ' sxmax ]
   Cxmx=input('Please enter INTeger Xmax value for pixel summation: ')
   figure;
	switch fCt
		case 1
         for tc=1:1:tmax;
            sump=0
      		sump=(pmat(1,tc) + pmat(Cxmx))/2 % trapezoidal integration
				for xc=2:1:Cxmx-1;
            	sump=pmat(xc,tc) + sump;
         	end
            Ct(tc)=sump; time(tc)=10^(tc-tes)
         end % end pixel sum
         semilogx(time,Ct,'m:')
         dummy=menu('            ','Press to continue')
         close;
      case 2
      end
   end % end switch fchoice