%	Test program for numerical solution of radial Couette flow
%	Test if computed streamfunction will reproduce exact velocity and vorticity
	clear all
   clf reset
   beta=100;
   jmax=10;
   jm1=jmax-1;
   gamma=1/log(beta);
   gam2=gamma^2;
   del=1/(jmax-1);
   r=zeros(1,jmax);
   vexact=zeros(1,jmax);
	vnumr=zeros(1,jmax);
	vnumz=zeros(1,jmax);
   wnum=zeros(1,jmax);
   s=zeros(1,jmax);
   I=1:jmax;
   r=exp((I-1)*del/gamma);
   vexact=(r-1./r)./(beta-1/beta);
   wexact=2/(beta-1/beta);
   for i=2:jmax
      s(i)=s(i-1)-0.5*(vexact(i-1)+vexact(i))*(r(i)-r(i-1));
   end
   for i=2:jm1
      wnum(i)=-(s(i+1)-2*s(i)+s(i-1))*(gamma/(r(i)*del))^2;
      vnumr(i)=-(s(i+1)-s(i-1))/(r(i+1)-r(i-1));
      vnumz(i)=-(s(i+1)-s(i-1))*(gamma/r(i))/(2*del);
   end
   figure(1)
   plot(r,vexact,'r-',r,vnumr,'b-+',r,vnumz,'g-*')
   title(['Velocity Profiles, beta = ',num2str(beta),', JMAX = ',num2str(jmax)])
   xlabel('Radius')
   ylabel('Velocity')
   figure(2)
   plot([1 beta],[wexact wexact],'r-',r,wnum,'g-*')
   title(['Vorticity Profiles, beta = ',num2str(beta),', JMAX = ',num2str(jmax)])
   xlabel('Radius')
   ylabel('Vorticity')
   
      
	
   