disp(' ');
disp(' rectangular_plate_base_accel_fea.m   ver 1.8   February 26, 2013');
disp('  ');
disp(' by Tom Irvine   Email: tomirvine@aol.com ');
disp(' ');
disp(' This script calculates the natural frequencies and ');
disp(' mode shapes of a flat, thin, rectangular plate. ');
disp(' ');
disp(' It also calculates the frequency response functions ');
disp(' for base excitation.      ');
disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
close all hidden;
%
clear mass;
clear stiff;
clear mass_unc;
clear stiff_unc;
clear mass_hold;
clear stiff_hold;
clear new_mass;
clear new_stiff;
%
clear fn;
clear omega;
clear ModeShapes;
clear MST;
clear modes;
%
clear node_matrix;
clear element_matrix;
clear constraint_matrix;
%
clear emm;
clear pdof;
%
clear accel_half_sine;
clear trans_acc_frf;
clear trans_vel_frf;
clear trans_disp_frf;
clear trans_rel_disp_frf;
%
clear fig_num;
clear con;
clear length;
clear ea;
clear nd;
%
clear acc;
%
clear forig;
%
iu=1;
%
fig_num=1;
%
[fig_num,number_nodes,node_matrix,element_matrix,nelem,...
       nodex,nodey,nx,ny,node1,node2,node3,node4,dof]=...
                rectangular_plate_nodes_elements(fig_num);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[D,rho,mu,beta,thick]=rectangular_plate_materials(iu,L,W);
%
[mass,stiff,total_volume,total_mass]=rectangular_plate_mass_stiff...
 (nelem,nodex,nodey,node1,node2,node3,node4,mu,beta,thick,dof,rho,D);   
%
mass_unc=mass;
stiff_unc=stiff;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
sz=size(mass);
dof=(sz(1));
num=dof;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
out1=sprintf('\n number of dofs =%d \n',num);
disp(out1);
%
[ea,nd,nem,idn]=rectangular_plate_drive_nodes(nx,ny);
%
[con_nodes,num_con_nodes,con_dof,num_con_dof,rflag]=...
                    rectangular_plate_constraint_nodes(nx,ny,nd,idn,nem);                
%
mass_unc=mass;
stiff_unc=stiff;
%
if(rflag==1)
%
    con_dof=sort(con_dof,'descend');
%
    for i=1:num_con_dof
        nn=con_dof(i);
%       
        mass(nn,:)=[];
        mass(:,nn)=[];
%
        stiff(nn,:)=[];
        stiff(:,nn)=[];
%
        for j=1:nem
            if(ea(j)>nn)
                ea(j)=ea(j)-1;
            end
        end
%
    end
%
end
%
%
num=max(size(mass)); % may have changed due to constraints
nff=num-nem;
%
[nff,ngw,Mwd,Mww,fn,omegan,ModeShapes,MST,part,T1,T2]=...
 partition_matrices_accel(num,nem,ea,mass,stiff);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Select modal damping input method ');
disp('  1=uniform damping for all modes ');
disp('  2=damping vector ');
idm=input(' ');
%
disp(' ');
if(idm==1)
   disp(' Enter damping ratio ');
   dampr=input(' ');
   damp(1:num)=dampr;
else
   disp(' Modal Damping Vector ');
   damp = input(' Enter the vector:  ');
end
%
[acc_frf,fstart,fend,freq]=...
 partition_matrices_accel_frf(num,nff,nem,ngw,Mwd,Mww,fn,damp,omegan,...
                                          ModeShapes,MST,part,T1,T2);
%
disp(' rectangular_plate_base_fea_post ');
%                                            
[trans_acc_frf,trans_vel_frf,trans_disp_frf,trans_rel_disp_frf,fig_num]=...
 rectangular_plate_base_fea_post(con_dof,num_con_dof,nodex,nodey,acc_frf,...
                                                 freq,fstart,fend,fig_num);
%
%
PF=part;
%
disp(' ');
disp(' Apply half-sine base input?  1=yes  2=no ');
isb=input(' ');
if(isb==1)
%
  [t,acc,rel_disp,abase]=...
    mdof_half_sine_core(fn,damp,PF,ModeShapes,T1,T2,ngw);  
%  
  [acc_half_sine,rel_disp_half_sine,fig_num]=...
      rectangular_plate_th_plot(t,acc,rel_disp,abase,...
                                  nodex,nodey,con_dof,num_con_dof,fig_num);
end
%
disp(' ');
disp(' Apply arbitrary time history base input?  1=yes  2=no ');
iar=input(' ');
if(iar==1)
%
  [t,acc,rel_vel,rel_disp,abase]=...
                mdof_arbit_base_core(fn,damp,PF,ModeShapes,T1,T2,ngw);
%  
  [acc_arb,rel_disp_arb,fig_num]=...
               rectangular_plate_th_plot(t,acc,rel_disp,abase,...
                                  nodex,nodey,con_dof,num_con_dof,fig_num);    
end
%
%
disp(' ');
disp(' Output arrays: ');
disp(' ');
disp('   trans_acc_frf');
disp('   trans_vel_frf');
disp('   trans_disp_frf');
disp('   trans_rel_disp_frf');
%
if(isb==1)
    disp('   acc_half_sine');
    disp('   rel_disp_half_sine');       
end
%
if(iar==1)
    disp('   acc_arb');
    disp('   rel_disp_arb');       
end