%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This is the main program for show.m
%%
%% MtOption, ShowOption and SurfOption(for ShowOption=3 only)
%% should be provided at the beginning.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function show_main

global GeomParam MatlParam MiscParam AnisParam AnisFlag;
global MatlEvalStr MiscEvalStr LayerEvalStr;
global LayerNm MatNm;
global ShowOption MtOption SurfOption;
global M F H S w m modenum;
global VECTOR eigvalue;
global HsVect HsVectStr OmVect OmegaInp LEN hsom;
global LayerLabelFlag MessageTextFlag UIprintFlag; 
global SpacingFactor MeshLineFlag NodePointFlag SurfContourFlag;
global R Rout Rin BtmPos rcod zcod rnum pnum znum Layerflag;
global glueyes PadAngle maxn p2n Thick;


err=1.0e-10;
cla;
legend off;
zoom off;
rotate3d off;
if ShowOption==1 | ShowOption==9,

%=================mload======================
err=1.0e-10;
Rin=[];Rout=[];Thick=[];BtmPos=[];rnum=[];znum=[];glueyes=[];
pnum=length(GeomParam);
for i=1:pnum,
   Layermat(i)=str2num(GeomParam{i}{1});
   Layerflag(i)=2*(str2num(GeomParam{i}{3})-1)+str2num(GeomParam{i}{2});
  %Layerflag=1:stator 2:rotor 3:fricless stator 4:fricless rotor
   Rin(i)=str2num(GeomParam{i}{4});
   Rout(i)=str2num(GeomParam{i}{5});
   Thick(i)=str2num(GeomParam{i}{6});
end;
FunBtmPos;  %determine BtmPos
for i=1:length(MatlParam),
   Modulus(i)=str2num(MatlParam{i}{2});
   Conductivity(i)=str2num(MatlParam{i}{3});
   Expansion(i)=str2num(MatlParam{i}{4});
   Density(i)=str2num(MatlParam{i}{5});
   Specific(i)=str2num(MatlParam{i}{6});
   Poisson(i)=str2num(MatlParam{i}{7});
end;
for k=1:length(MatlParam),
   if AnisFlag(k)==1,
   for i=1:6,
     AnisD{k}(i,1)=str2num(AnisParam{k}{i,7});
     for j=i:6
     AnisC{k}(i,j)=str2num(AnisParam{k}{i,j});
     AnisC{k}(j,i)=AnisC{k}(i,j);
     end;     
   end;   
   end;
end;
Diffusivity=Conductivity./(Density.*Specific);
Mvflag=mod(Layerflag-1,2); %moving=1;stationary=0
bcid1=str2num(MiscParam{1})-1; %note: minus 1 due to popup button
bcid2=str2num(MiscParam{2})-1; %note: minus 1
meshsize=str2num(MiscParam{3});
HsVectStr=MiscParam{4};
OmegaInp=MiscParam{5};
PadAngle=str2num(MiscParam{6});
Fric=str2num(MiscParam{7});

if isempty(HsVectStr)==1,
   HsVect=0:4:20;
else
   HsVect=str2num(HsVectStr);
end;
OmVect=str2num(OmegaInp);
%%[sortK,ind]=sort(Conductivity);
fcoeff=PadAngle/360*Fric;  
%obtain equivalent f for brake problem

glueyes(1:pnum)=1; %initialize 
for p=2:pnum,     	
    for q=1:pnum,
	    if IsBound(p,q)==0, 
	       glueyes(p)=0;
	       break;   %break from q loop
	    end;
    end;
end;


%================zmesh start=========

Loop_Init_Omega=500; %% initial iteration w
if PadAngle<360, Loop_Init_Omega=50; end; %if brake
avgm=round(0.5*HsVect(1)+0.5*HsVect(length(HsVect)));
if isempty(OmegaInp)==1,
   ww=1000;  % to find skin region
   if PadAngle<360, ww=100;end;
else
   ww=OmVect(length(OmVect));
end;

switch meshsize
      case 1, ns=1; num=3; bsvect=[1.50,2.5];IterTol=0.08;
      case 2, ns=2; num=4; bsvect=[1.50,2.0];IterTol=0.06;
      case 3, ns=3; num=6; bsvect=[1.40,1.8];IterTol=0.04;
      case 4, ns=4; num=8; bsvect=[1.40,1.8];IterTol=0.02;
      case 5, ns=5; num=10; bsvect=[1.30,1.5];IterTol=0.015;
      case 6, ns=6; num=12; bsvect=[1.30,1.4];IterTol=0.01;
      case 7, ns=8; num=14;bsvect=[1.20,1.3];IterTol=0.01;   
end;

len=length(Conductivity);
[sortK,ind]=sort(Conductivity);
for p=1:pnum,
   if Layermat(p)==ind(len),
      bs=bsvect(1);
   else
      bs=bsvect(2);
   end;
   
   Bias(p)=1;
   for q=1:pnum,
       if IsBound(p,q)==0 & Layerflag(p)<3  %frictional contact
          Bias(p)=bs;
          break;
       end; 
   end; %q
   for q=1:pnum,
       if IsBound(q,p)==0 & Layerflag(q)<3  %frictional contact
          Bias(p)=Bias(p)/bs;
          break;
       end; 
   end; %q

   if avgm==0,
      kai=2*pi*sqrt(2*Diffusivity(Layermat(p))/(ww*1));  %banding mode
   else
      kai=2*pi*sqrt(2*Diffusivity(Layermat(p))/(ww*avgm)); %skin region
   end;
   if abs(Bias(p)-1.0)<0.01,
      znum(p)=round(0.5*Thick(p)/kai*(1+meshsize*0.3));
   else
      znum(p)=round(log(Thick(p)*(bs^ns-1)/kai+1)/log(bs));
   end;
   
   if Layermat(p)==ind(len), 
            znum(p)=round(znum(p)/2); 
   end %good conductor
   if znum(p)<2, znum(p)=2;
   elseif znum(p)>25, znum(p)=25;
   end
end; %p
 
%===============zmesh ends===========
%%determine r mesh
sortR=sort([Rin,Rout]);
len=length(sortR);
R=sortR(1);
dL=(sortR(len)-sortR(1))/num;
for i=1:len-1,
R1=sortR(i);
R2=sortR(i+1);
if abs(R1-R2)>err,
  numi=ceil((R2-R1)/dL-err);
  dRi=(R2-R1)/numi;
  R=[R,R1+dRi*[1:numi]];
end;
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Generate necessary model information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%node information
z0=0;n=0;e=0;n2p=[];p2n=[];
nbot=[];ntop=[];gbot=[];gtop=[];nofric=[];rcod=[];zcod=[];
for p=1:pnum,
   for i=1:length(R),
      if abs(Rin(p)-R(i))<err;numa=i;end;
      if abs(Rout(p)-R(i))<err;numb=i;end;
   end;
   rcod{p}=R(numa:numb);
   rnum(p)=length(rcod{p})-1;
   
   bs=Bias(p);   
   i=0:znum(p);
   if(abs(bs-1)<0.01)  %bias==1 needs special treatment
      detz=Thick(p)*i/znum(p);
   else  
      detz=Thick(p)*(power(bs,i)-1)/(power(bs,znum(p))-1);
   end
   zcod{p}=detz+BtmPos(p);
        
   for j=1:znum(p)+1,
      for i=1:length(rcod{p}),
         n=n+1;
         n2p(n,:)=[p,i,j];
         p2n(p,i,j)=n;
      end;
   end;
   maxn(p)=n;
   
   if p>1,     	
     for q=1:pnum,
       if IsBound(p,q)>-1  %contact or bound
        sortR=sort([Rin(p),Rout(p),Rin(q),Rout(q)]);
        R1=sortR(2);
        R2=sortR(3);
        for i=1:length(rcod{p}),
           if abs(R1-rcod{p}(i))<err;numa=i;end;
           if abs(R2-rcod{p}(i))<err;numb=i;end;
        end;
        nvect=p2n(p,numa,1):p2n(p,numb,1);
        if IsBound(p,q)==0,
           ntop=[ntop,nvect];
        else
           gtop=[gtop,nvect];
        end;
                    
        for i=1:length(rcod{q}),
           if abs(R1-rcod{q}(i))<err;numa=i;end;
           if abs(R2-rcod{q}(i))<err;numb=i;end;
        end;
        nvect=p2n(q,numa,znum(q)+1):p2n(q,numb,znum(q)+1);
        if IsBound(p,q)==0,
           nbot=[nbot,nvect];   
        else
           gbot=[gbot,nvect];
        end;
        if Layerflag(p)>2, %% fricless plate
           nofric=[nofric,nvect];
        end;
      end;% if isbound
     end;% for q
   end %p>1
end; %p

nd=n; %total node number
nc=length(nbot); 
fullnods=1:nd;
bc1=[1:rnum(1)+1];
bc2=[nd-rnum(pnum):nd];
if bcid1==5,   % cyclic bc
   gbot=[gbot,bc2];
   gtop=[gtop,bc1];
end;

%element information
for p=1:pnum,
   for j=0:znum(p)-1,
      for i=1:rnum(p),
         e=e+1;
         if p==1,
         elem(e,1)=j*(rnum(p)+1)+i;
         else
         elem(e,1)=maxn(p-1)+j*(rnum(p)+1)+i;
         end;
         elem(e,2)=elem(e,1)+1;
         elem(e,3)=elem(e,2)+rnum(p)+1;
         elem(e,4)=elem(e,1)+rnum(p)+1;
         pelm(e)=p;
      end;
   end;
   maxe(p)=e;
end;

%==================mload ends======================
end;
if ShowOption==1,   %generate mesh
   figure(1);
   view(0,90);
   maxR=max(Rout);
   minR=min(Rin);
   rmean=(minR+maxR)/2;
   %textR=maxR+[0.07,0.17]*(maxR-minR);
   spacing=0.05*(maxR-minR);
   DrawBtmPos(spacing);  %determine btmpos in draw mode (spacing is included)
   plot(rcod{1}([1,rnum(1)+1]),[0 0],'r','linewidth',1.5);
   hold on;
   
  for p=1:pnum,
      nr=rnum(p)+1;
      nlay=znum(p)+1;
      
      Rcod=rcod{p};
      Zcod=zcod{p}-zcod{p}(1)+BtmPos(p);
      
	if MeshLineFlag==1, 
      for i=2:nlay-1,plot(Rcod([1,nr]),Zcod(i)*[1,1]);end; %horizontal lines
      for i=2:nr-1,plot(Rcod(i)*[1,1],Zcod([1,nlay]));end; %vertical lines   
	end;       
      plot(rcod{p}([1,nr,nr,1,1]),Zcod([1,1,nlay,nlay,1])...
         ,'Color','r','linewidth',1.5);
      if NodePointFlag==1, 
         for i=1:nlay,
         plot(Rcod,Zcod(i)*ones(1,nr),'g.','markersize',14,'linewidth',2);%show nodes!!
         end; 
      end;
      if LayerLabelFlag==1,
         textoffset=(Rout(p)-Rin(p))*0.02;
         text(0.5*(Rin(p)+Rout(p))+textoffset,0.5*(Zcod(1)+Zcod(nlay)),num2str(p));
      end;
   end;
      
      axis([min(Rin)-0.1*spacing max(Rout)+0.1*spacing 0 Zcod(nlay)+0.1*spacing]);
      ShowStr=sprintf('Total : %d disks, %d layers, %d nodes, %d elements',...
         pnum-sum(glueyes)+1,pnum,nd,maxe(pnum));
      if MessageTextFlag==1,
         text(rmean,-(maxR-minR)*0.07,ShowStr,'HorizontalAlignment','center');
      end;
      axis equal;
      axis off;
      hold off;
      zoom on;
end;          

if ShowOption==9,   %generate 3d model with lighting
   tnum=60;
   dt=2*pi/tnum;
   maxR=max(Rout);
   minR=min(Rin);
   rmean=(minR+maxR)/2;
   textR=maxR+[0.07,0.17]*(maxR-minR);
   if PadAngle<360-err %brake
       EnlargeFactor=1; % enlarge thickness in view
       spacing=0.1*(maxR-minR);
   else 
       EnlargeFactor=2;
       spacing=0.15*(maxR-minR);
   end;
   DrawBtmPos(spacing);
      
   for p=1:pnum,
   switch Layermat(p),
   case 1, facecolor=[0 1 0];
   case 2, facecolor=[1 0 0];
   case 3, facecolor=[1 1 0];
   case 4, facecolor=[1 1 0];
   case 5, facecolor=[0 1 1];
   otherwise, facecolor=[0 1 0];
   end;
   %if Layerflag(p)>2, facecolor=[0 0 1];end;  %% endplate or piston
   nodvec=[];
   facevec=[];
   rin=Rin(p);
   rout=Rout(p);
   th=Thick(p);
      
   th0=BtmPos(p);
   if (p<3 | p>pnum-2) & PadAngle<360-err, %if brake
   padtnum=round(PadAngle/360*tnum);
   for i=0:padtnum,
   nodvec=[nodvec;...
        rin*cos(dt*i),rin*sin(dt*i),th0;...
        rout*cos(dt*i),rout*sin(dt*i),th0;...
        rout*cos(dt*i),rout*sin(dt*i),th0+th;...
        rin*cos(dt*i),rin*sin(dt*i),th0+th;...
        ];
   end;
   for i=0:padtnum-1,
   facevec=[facevec;...
           [1,2,6,5;...
            4,3,7,8;...
            2,6,7,3;...
            1,5,8,4]+4*i...
        ];
   end; %i
   facevec=[facevec;[1,2,3,4]; [1,2,3,4]+padtnum*4];
   
   
   else  %not brake
   for i=0:tnum-1,

   nodvec=[nodvec;...
        rin*cos(dt*i),rin*sin(dt*i),th0;...
        rout*cos(dt*i),rout*sin(dt*i),th0;...
        rout*cos(dt*i),rout*sin(dt*i),th0+th;...
        rin*cos(dt*i),rin*sin(dt*i),th0+th;...
        ];

   if i<tnum-1,
   facevec=[facevec;...
           [1,2,6,5;...
            4,3,7,8;...
            2,6,7,3;...
            1,5,8,4]+4*i...
        ];
   else
   facevec=[facevec;...
           [1+4*i,2+4*i,6,5;...
            4+4*i,3+4*i,7,8;...
            2+4*i,6,7,3+4*i;...
            1+4*i,5,8,4+4*i]...
           ];
   end; %if
   end; %i
   end; %if brake
   nodvec(:,3)=nodvec(:,3)*EnlargeFactor;
    patch('faces',facevec,'vertices',nodvec,...
        'FaceVertexCData',facecolor,'FaceColor','flat');   
   end; %p
    
   axis equal;
   axis tight;
   axis off;
   hold off;
   view(86,16);
   %camlight('right');
   camlight headlight;
   rotate3d on;
   zoom off;
end;          


if ShowOption==6; % critV-m curve
      plot(hsom(:,1),hsom(:,2),'*','MarkerSize',10,'linewidth',1.5);
	    MyAxis=axis;
	    MyAxis([1,3])=0;
      axis(MyAxis*1.1);
      set(gca,'fontsize',12);
      xlabel('hot spot number','fontsize',14);
      ylabel('critical speed (rad/s)','fontsize',14);
end;
	
if ShowOption==7, % migV-m curve
      xx=hsom(:,1);
      yy=hsom(:,2);
      zz=abs(hsom(:,3));
      len=length(zz);
      
      if yy(len)-zz(len)>zz(len),
      ShowStr='speed relative to stators';
      else
      zz=yy-zz;
      ShowStr='speed relative to rotors';
      end;
      
      for i=1:len,
         if xx(i)==0,zz(i)=0;end;
      end;   %%find banding mode mig speed
      
      plot(xx,zz,'*','MarkerSize',10,'linewidth',1.5);
      MyAxis=axis;
      MyAxis([1,3])=0;
      axis(MyAxis*1.1);
     %%text(0.55*MyAxis(2),1.25*MyAxis(4)-0.132*MyAxis(3),ShowStr,...
	  %%'HorizontalAlignment','center');
      set(gca,'fontsize',12);
      xlabel('hot spot number','fontsize',14);
      ylabel('migration speed(rad/s)','fontsize',14); 
      legend(ShowStr);  
end;

if ShowOption==7.1, %v-b curve
      mymark={'*-','o-','^-','+-','d-','s-','v-','p-','x-'};
      nn1=length(HsVect);
      nn2=length(OmVect);
      growMat=reshape(hsom(:,3),nn2,nn1)';
      eigMat=reshape(hsom(:,4),nn2,nn1)';
      
      for j=1:nn1,
      if HsVect(j)>=10,
      mylegend(j,:)=cat(2,'n=',num2str(HsVect(j)));
      else
      mylegend(j,:)=cat(2,'n= ',num2str(HsVect(j)));
      end;
      plot(OmVect,growMat(j,:),mymark{j});
      hold on;
      end;
      legend(mylegend,2);
      plot([0 OmVect(length(OmVect))],[0 0],'r--','linewidth',1.5);
      
      hold off;
      set(gca,'fontsize',12);
      xlabel('Sliding speed(rad/s)','fontsize',14);
      ylabel('Growth rate','fontsize',14); 
      %%axis([0 OmVect(nn2) min(min(growMat)) max(max(growMat))]);
end;

if ShowOption>1 & ShowOption<6    %show computation results

ycm; %context menu
if (modenum<1 | modenum>LEN);modenum=1;end;
if(SurfOption<1 | SurfOption>pnum);SurfOption=1;end;   
rl=real(eigvalue(modenum));
tvect=real(VECTOR(:,modenum));
%%tcon=(reshape(tvect,nr,nlay))';
if m==0
   mig=0;
else
   mig=imag(eigvalue(modenum))/m;
end;
if modenum==1 & rl<0
   fprintf('\n !! NO POSITIVE GROWTH RATES AVAILABLE. USE A LARGER SPEED.\n')
end;

switch ShowOption  
case 2,   %contour plot on cross section
   maxR=max(Rout);
   minR=min(Rin);
   rmean=(minR+maxR)/2;
   maxtvect=max(tvect);
   mintvect=min(tvect);
   %textR=maxR+[0.1,0.2]*(maxR-minR);
   spacing=0.05*(maxR-minR);
   DrawBtmPos(spacing);
     
   for p=1:pnum,
      nr=rnum(p)+1;
      nlay=znum(p)+1;
	Zcod=zcod{p}-zcod{p}(1)+BtmPos(p);      
      tcon=(reshape(tvect(maxn(p)-nr*nlay+1:maxn(p)),nr,nlay))';
      contourf(rcod{p},Zcod,tcon,10);
      hold on;
      caxis([mintvect,maxtvect]);
      shading flat;
      plot(rcod{p}([1,nr,nr,1,1]),Zcod([1,1,nlay,nlay,1]),'k'); 
      %plot block lines
      if p<pnum & glueyes(p+1)==1
         LineColor='g';
      elseif p==pnum,
         LineColor='k';
      else 
         LineColor='r';
      end;
      if LayerLabelFlag==1,
         textoffset=(Rout(p)-Rin(p))*0.02;
         text(0.5*(Rin(p)+Rout(p))+textoffset,0.5*(Zcod(1)+Zcod(nlay)),num2str(p));
      end;
   end;
   axis([min(Rin)-0.1*spacing max(Rout)+0.1*spacing 0 Zcod(nlay)+0.1*spacing]);
   hold off;
   axis equal;
   axis off;
   
   ShowStr=sprintf('N=%d, mode=%d, v=%8.2e(rad/s), b=%9.2e, c=%8.3e',...
      m,modenum,w,rl,abs(mig));
   if MessageTextFlag==1,
      text([rmean],[-(maxR-minR)*0.15],ShowStr,'HorizontalAlignment','center');  
   end;
      
  case 3,   %%% contour on surface
     maxabstvect=max(abs(tvect));
     cnum=120;
     p=SurfOption;
     rin=Rin(p);
     rout=Rout(p);
     rmean=(rin+rout)/2;
  	  for i=0:rnum(p),
		for j=0:cnum,
		x(i+1,j+1)=(rin+i*(rout-rin)/rnum(p))*cos(j*2*pi/cnum);
		y(i+1,j+1)=(rin+i*(rout-rin)/rnum(p))*sin(j*2*pi/cnum);
		tcon(i+1,j+1)=real(exp(m*2*pi/cnum*j*1.0j)*VECTOR(maxn(p)-rnum(p)+i,modenum));
		end
	  end

    if SurfContourFlag==1,
      surf(x,y,tcon);
      view(-74,64);
    else
      contourf(x,y,tcon,10);
      caxis([-maxabstvect,maxabstvect]);
	   shading flat;
	   THETA=0:0.002:2*pi;
	   hold on;
	   plot(rin*cos(THETA),rin*sin(THETA),'k',...   %plot circles
  		 rout*cos(THETA),rout*sin(THETA),'k','linewidth',1.);
      hold off;
      ShowStr=sprintf('N=%d, mode=%d, surf=%d, v=%8.2e(rad/s), b=%9.2e, c=%8.3e',...
       m,modenum,SurfOption,w,rl,abs(mig));
     if MessageTextFlag==1,
        text([0],[-rout-rmean*0.2],ShowStr,'HorizontalAlignment','center');   
     end;
	   axis equal;
     axis off;
    end; %surforcontour
     
  case 4,
     pvect=1:pnum;
     if PadAngle<360-err, %if brake
        midr=rcod{1}(round(rnum(1)/2));
        pvect(pnum-2)=[];pvect(pnum-3)=[];
     else
        midr=R(round(length(R)/2));
     end;
     for p=pvect,
        tcon=[];
        for i=1:rnum(p)+1,
           if abs(rcod{p}(i)-midr)<err,midnum=i;end;
        end;
        for i=1:znum(p)+1,tcon(i)=tvect(p2n(p,midnum,i));end;
        plot(zcod{p},tcon,'-',zcod{p},tcon,'m.','markersize',18,'linewidth',1.5);
        hold on;
     end;
      MyAxis=axis;
      for p=pvect(2:length(pvect)),      
      if glueyes(p)==1
          plot(zcod{p}(1)*[1,1],MyAxis(3:4),'g--');
      else 
          plot(zcod{p}(1)*[1,1],MyAxis(3:4),'r--');
      end;
      end; %interface lines
      hold off;
      axis([0,sum(Thick(pvect)),MyAxis(3),MyAxis(4)]);
      set(gca,'fontsize',12);
      xlabel('thickness direction (m)','fontsize',14);
      ylabel('temperature perturbation','fontsize',14);
      ShowStr=sprintf('N=%d, mode=%d, v=%8.2e(rad/s), b=%9.2e, c=%8.3e',...
         m,modenum,w,rl,abs(mig));
      if MessageTextFlag==1,
         text(sum(Thick(pvect))/2,.95*MyAxis(4)+0.05*MyAxis(3),ShowStr,...
            'HorizontalAlignment','center');
      end;
      %title(ShowStr,'FontSize',10);
   case 5,
      pvect=1:pnum;
      if PadAngle<360-err, %if brake
        midr=rcod{1}(round(rnum(1)/2));
        pvect(pnum-2)=[];pvect(pnum-3)=[];
      else
        midr=R(round(length(R)/2));
      end;

      tcon=[];
      for p=pvect,
        for i=1:rnum(p)+1,
           if abs(rcod{p}(i)-midr)<err,midnum=i;end;
        end;
        for i=1:znum(p),
           tcon=[tcon,tvect(p2n(p,midnum,i))];
        end;
     end;
     tcon=[tcon,tvect(p2n(p,midnum,znum(p)+1))];  %last one
     plot(1:length(tcon),tcon,'-',1:length(tcon),tcon,'m.','markersize',18,'linewidth',1.5);
     MyAxis=axis;
     hold on;
     for i=2:length(pvect),
      p=pvect(i);      
      if glueyes(p)==1
          plot(sum(znum(pvect(1:i-1)))*[1,1]+1,MyAxis(3:4),'g--');
      else 
          plot(sum(znum(pvect(1:i-1)))*[1,1]+1,MyAxis(3:4),'r--');
      end;
      end; %interface lines
         hold off;
      axis([1 length(tcon) MyAxis(3) MyAxis(4)]);
      set(gca,'fontsize',12);
      xlabel('node (equally spaced) thru thickness','fontsize',14);
      ylabel('temperature perturbation','fontsize',14);
      ShowStr=sprintf('N=%d, mode=%d, v=%8.2e(rad/s), b=%9.2e, c=%8.3e',...
         m,modenum,w,rl,abs(mig));
      if MessageTextFlag==1,
	   text(0.5+length(tcon)/2,0.95*MyAxis(4)+0.05*MyAxis(3),ShowStr,...
         'HorizontalAlignment','center');
      end;
      %titleHndl=title(ShowStr,'fontsize',10);      
end; %end switch
end %if ShowOp>1 & <6

function isbound=IsBound(toplayer,botlayer)
%%%isbound= 0:make contact   1:bound     -1:separate
global Layerflag BtmPos Thick Rin Rout;

distance=abs(BtmPos(botlayer)+Thick(botlayer)-BtmPos(toplayer));
if distance<2e-6 & ~(Rin(toplayer)>Rout(botlayer) | Rin(botlayer)>Rout(toplayer)),
% make contact and not separated radially
  if Layerflag(toplayer)>2 
     isbound=0;
  elseif mod(Layerflag(botlayer)-Layerflag(toplayer),2)==1
     isbound=0;
  else
     isbound=1;
  end;
else   % the two layers are separate
     isbound=-1;
end;
  
function FunBtmPos
global Rin Rout Thick BtmPos pnum;

BtmPos=[];
BtmPos(1)=0;
for i=2:pnum,
[ht, hind]=sort(BtmPos(1:i-1)+Thick(1:i-1));
for j=i-1:-1:1,
  z=hind(j);
  if ~(Rin(i)>Rout(z) | Rin(z)>Rout(i))
    BtmPos(i)=BtmPos(z)+Thick(z);
    break;
  end;
end; %j
end; %i

function DrawBtmPos(spacing)
global Rin Rout Thick BtmPos pnum glueyes SpacingFactor;

BtmPos=[];
BtmPos(1)=0;
for i=2:pnum,
[ht, hind]=sort(BtmPos(1:i-1)+Thick(1:i-1));
for j=i-1:-1:1,
  z=hind(j);
  if ~(Rin(i)>Rout(z) | Rin(z)>Rout(i))
    BtmPos(i)=BtmPos(z)+Thick(z)+(1-glueyes(i))*spacing*SpacingFactor;
    break;
  end;
end; %j
end; %i
