%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INSTRUCTIONS FOR TEI.M:
%%
%% This is the main program. You are first prompted to input the 
%% Fourier numbers (hot spot numbers)  , following the general Matlab
%% format. For example, you can input 2,6,10,14 or 2 6 10 14 or 2:4:14. 
%% They are the same thing. It is OK that you hit the 'return' key 
%% without inputing anything.  The default value is 0:2:20, i.e., 
%% 0,2,4,6,...20.
%%   
%% Then you are prompted to specify the sliding speed (rad/s) you 
%% would study.  If a sliding speed is specified, then you will get
%% the growth rate and the migration speed of the perturbation, but no
%% critical speed available.  If you want to compute the critical speed 
%% instead (and the associated migration speed), just hit the 'return' 
%% key without inputing anything.    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function tei_loop

global GeomParam MatlParam MiscParam AnisParam AnisFlag;
global LayerNm MatNm;
global M F H S w m;
global VECTOR eigvalue;
global HsVect HsVectStr OmegaInp OmVect;
global modenum GrowRate MigSpeed;
global ShowOption MtOption SurfOption;
global killa killb nkeep nrem nods hsom;
global nu af E AnisC AnisD AnisFlag;
global r z r0 a b;
global textwin textwinHndl;
global LayerLabelFlag MessageTextFlag UIprintFlag;
global R Rout Rin BtmPos rcod zcod rnum pnum znum Layerflag 
global glueyes PadAngle maxn p2n Thick;
global LoadFileName MatlFileName;



timeflag=0;
hsom=[];
SurfOption=1;


%%%%%demo version;  removed in the full version
clear m w VECTOR eigvalue hsom textwin;
if ~isempty(findstr(upper(LoadFileName),'BRAKE'))   %contain string of 'brake'
   if strcmp(LoadFileName,'demo_brake.d')==0 % if not the file 'demo_brake' 
      LoadFileName='demo_brake.d';
	[GeomParam,MatlParam,MiscParam,AnisParam,AnisFlag]=Read_Param(LoadFileName,MatlFileName); 
      MtOption=1;ShowOption=1;figure(1);show_main;
    end;   
    load demo_data_brake m w VECTOR eigvalue hsom textwin;
    SurfOption=2;
else
  if strcmp(LoadFileName,'demo_clutch.d')==0   % if not the file 'demo_clutch' 
      LoadFileName='demo_clutch.d';
	[GeomParam,MatlParam,MiscParam,AnisParam,AnisFlag]=Read_Param(LoadFileName,MatlFileName); 
      MtOption=1;ShowOption=1;figure(1);show_main;
  end;   
  load demo_data_clutch m w VECTOR eigvalue hsom textwin;
  SurfOption=5;
end;
show_textwin;
%%%%%%end of demo

%=================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======================

%%%Main program starts here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%FEA algorithm is hidden in demo version
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=1;    
MtOption=0;
%%SurfOption=1; 




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 show_textwin
global textwin textwinHndl;

for i=1:7,strlst{i}=textwin{i};end;
for i=8:length(textwin),
mypause(0.13);
strlst{i}=textwin{i};
set(textwinHndl,'string',strlst,'listboxtop',length(strlst));
drawnow;
end;


function mypause(pausetime)

for i=1:round(100000*pausetime),
sin(1);
end;
