%
%
disp(' ');
disp(' circular_point_four.m  ver 1.0  October 3, 2011 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' This program calculates the natural frequencies of the ');
disp(' bending modes of a circular plate with four equally-spaced '); 
disp(' point supports. ');
%
clear f;
clear kk;
clear nn;
clear lambda;
clear lambda2;
clear ZZ;
clear ZZr;
iflag=0;
%
disp(' ');
disp(' Enter unit type: 1=English 2=metric');
iu=input(' ');
disp(' Enter plate type: ');
disp(' 1=homogeneous 2=honeycomb-sandwich');
ib=input(' ');
%
if(iu==1)
    disp(' ');
    diam=input(' Enter diameter (inch) ');
    disp(' ');
else
    disp(' ');
    disp(' Select diameter unit:');
    disp(' 1=meters  2=cm  3=mm ');
    ilu=input(' ');  
    if(ilu==1)
        disp(' ');
        diam=input(' Enter diameter (m) ');       
    end
    if(ilu==2)
        disp(' ');
        diam=input(' Enter diameter (cm) ');
        diam=diam/100;
    end  
    if(ilu==3)
        disp(' ');
        diam=input(' Enter diameter (mm) ');
        diam=diam/1000;
    end   
end
radius=diam/2;
%
if(ib==1) % homogeneous
%
    if(iu==1)
        disp(' ');
        h=input(' Enter thickness (inch) ');
    else
        disp(' ');
        h=input(' Enter thickness (mm) ');
        h=h/1000;
    end
%
    disp(' ');
    imat=input(' Enter material: \n 1=aluminum  2=steel  3=G10  4=other ');
%
    if(imat==1) % Al
        if(iu==1)
            rho=0.1;
            E=1.0e+07;
        else
            rho=2768;
            E=6.891e+010;
        end
        mu=0.3;
    end
    if(imat==2) % steel
        if(iu==1)
            rho=0.29;
            E=3.0e+007;
        else
            rho=7887;
            E=2.067e+011;
        end
        mu=0.3;
    end
    if(imat==3) % G10
        if(iu==1)
            rho=0.065;
            E=2.7e+006;
        else
            rho=1799;
            E=1.861e+010;        
        end
        mu=0.12;
    end
    if(imat ~=1 & imat ~=2 & imat ~=3)
        if(iu==1)
            disp(' ');   
            E = input(' Enter the elastic modulus (lbf/in^2) ');
            disp(' ');
            rho = input(' Enter the mass density (lbm/in^3) '); 
        else
            disp(' ');   
            E = input(' Enter the elastic modulus (N/m^2) ');
            disp(' ');
            rho = input(' Enter the mass density (kg/m^3) ');        
        end
        disp(' ');
        mu=input(' Enter Poisson ratio ');    
    end
%
    volume=pi*(radius^2)*h;
    De=E*h^3/(12*(1.-mu^2));
%
    smass=rho*(diam^2)*h*pi/4;
%
else % honeycomb-sandwich
    if(iu==1)
        disp(' ');
        disp(' Enter core thickness (inch)')
        c=input(' ');
%
        disp(' ');
        disp(' Enter individual face sheet thickness (inch)')
        t=input(' ');
    else
        disp(' ');
        disp(' Enter core thickness (mm)')
        c=input(' ');
        c=c/1000;
%
        disp(' ');
        disp(' Enter individual face sheet thickness (mm)')
        t=input(' ');
        t=t/1000;
    end
%
    disp(' ');
    disp(' Select face sheet material ');
    disp(' 1=aluminum  2=other ');
    imat = input(' ');
%
    if(imat==1)
        if(iu==1)
            E=1.0e+07;
            rho=0.1;
        else
            E=6.891e+010;
            rho=2768;
        end
        mu=0.3;
    else
        if(iu==1)
            disp(' ');
            disp(' Enter elastic modulus of face sheets (lbf/in^2)');
            disp(' ');
            E=input(' ');
            disp(' ');
            rho = input(' Enter the face sheet mass density (lbm/in^3) ');     
        else
            disp(' ');   
            disp(' Enter the elastic modulus of face sheets (N/m^2) ');
             E=input(' ');
            disp(' ');       
            rho = input(' Enter the face sheet mass density (kg/m^3) ');    
        end
        disp(' ');
        mu=input(' Enter the face sheet Poisson ratio ');  
    end
%
    if(iu==1)
        disp(' ');
        rho_core = input(' Enter the core mass density (lbm/in^3) ');
    else
        disp(' ');
        rho_core = input(' Enter the core mass density (kg/m^3) ');
    end
%
    area=pi*(radius^2);
    mass_face=area*(2*t)*rho;
    mass_core=area*c*rho_core;
    mass=mass_face+mass_core;
%
    H=c+t;
    De=E/(1-mu^2);
    De=(De*t*H^2)/2;
%
    thickness = c+2*t;
    h=thickness;
    volume = area*thickness;
    rho = mass/volume;
%
    smass=mass;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
if(iu==1)
    out1=sprintf(' Structural mass = %8.4g lbm',smass);
else
    out1=sprintf(' Structural mass = %8.4g kg',smass);
end
disp(out1);
disp(' ');
disp(' Add uniformly distributed non-structural mass ?');
disp(' 1=yes  2=no ')
imass=input(' ');
%
NSM=0;
if(imass==1)
  disp(' ');
  if(iu==1) 
    disp(' Enter added mass(lbm) ');
  else
    disp(' Enter added mass(kg) ');
  end
  NSM=input(' ');
end
%
mass=NSM+smass;
M=mass;
rho=M/volume;
%
disp(' ');
if(iu==1)
    out1=sprintf('     Structural mass = %8.4g lbm',smass);
    disp(out1);
    out1=sprintf(' Non-structural mass = %8.4g lbm',NSM);
    disp(out1);  
    out1=sprintf('          Total mass = %8.4g lbm ',M);
    disp(out1);
    out1=sprintf('     Total thickness = %8.4g in ',h);
    disp(out1);     
    out1=sprintf('            Diameter = %8.4g in ',2*radius);
    disp(out1);  
    out1=sprintf('              Volume = %8.4g in^3 ',volume);
    disp(out1);
    out1=sprintf('Overall mass density = %8.4g lbm/in^3 ',rho);
    disp(out1); 
    %
    disp(' ');
    out1=sprintf(' Plate Stiffness Factor  D = %8.4g lbf in',De);
    disp(out1); 
    out1=sprintf('                       D/E = %8.4g in^3',De/E);
    disp(out1);    
    %
    mass=mass/386;
    M=M/386;
    rho=rho/386;   % mass per volume
else
    out1=sprintf('     Structural mass = %8.4g kg',smass);
    disp(out1);   
    out1=sprintf(' Non-structural mass = %8.4g kg',NSM);
    disp(out1); 
    out1=sprintf('          Total mass = %8.4g kg ',M);
    disp(out1);
    out1=sprintf('     Total thickness = %8.4g mm ',h*1000);
    disp(out1);     
    out1=sprintf('            Diameter = %8.4g m ',2*radius);
    disp(out1);     
    out1=sprintf('              Volume = %8.4g m^3 ',volume);
    disp(out1);   
    out1=sprintf('Overall mass density = %8.4g kg/m^3 ',rho);
    disp(out1);     
    %
    disp(' ');
    out1=sprintf(' Plate Stiffness Factor  D = %8.4g N m',De);
    disp(out1);     
    out1=sprintf('                       D/E = %8.4g m^3',De/E);
    disp(out1);       
end
%
disp(' ');
%%%% n=input(' Enter the number of point supports  ');
n=4;
n2=n^2;
pi2=pi^2;
%
lower=n-1;
upper=n;
%
nd=360;
nr=50;
%
clear r;
clear theta;
clear omegan;
clear omegan_min;
%
omegan_min=1.0e+20;
ratio_error=1.0e+20;
dZdthrec=1.0e+20;
nde_min=1.0e+20;
%
nr=50;
nd=360;
%
dr=(radius/nr);
dtheta=(360/nd)*(nd+1);
%
r=linspace(0,nr*dr,nr);
r(1)=r(2)/100000;
rrn=r/radius;
%
theta=linspace(0,nd,dtheta);
theta=theta*pi/180;
%
clear length;
%
clear A;
clear B;
clear C;
clear T;
clear V;
clear omegan;
%
dZdth=1.0e+20;
d2rec=1.0e+20;
d3rec=1.0e+20;
%
%%% disp(' Enter number of trials ');
%%% nnn=input(' ');
%%% disp(' Enter number of cosines ');
%%% nc=input(' ');
%
nnn=1;
nc=4;
%
clear argr;
clear AAr;
%
argr=zeros(nc,1);
AAr=zeros(nc,1);
%
for(k=1:nnn)
    slopemax=-1.0e+20;
%   
%%    R=radius;
    R=1;
    
    clear AA;
    clear arg;
    clear    cos_sum;
    clear   dcos_sum;    
    clear  ddcos_sum;
    clear dddcos_sum;  
    
   cos_sum=zeros(length(r),1);    
  dcos_sum=zeros(length(r),1);        
 ddcos_sum=zeros(length(r),1);    
dddcos_sum=zeros(length(r),1);    


arg=[ 0.8709    0.0846    0.6709    0.6157 ]*pi;

AA= [1.9619   -1.3048   -3.4583    3.0771];
    
for(ijk=1:nc)
    Q=AA(ijk);
    eta=arg(ijk);
    for(j=1:length(r))
           alpha=arg(ijk)*rrn(j);
           cos_sum(j)=   cos_sum(j)+Q*cos(alpha);
          dcos_sum(j)=  dcos_sum(j)-Q*sin(alpha)*(eta/radius);
         ddcos_sum(j)= ddcos_sum(j)-Q*cos(alpha)*(eta/radius)^2;   
        dddcos_sum(j)=dddcos_sum(j)+Q*sin(alpha)*(eta/radius)^3;         
    end
end
shift=min(cos_sum);
cos_sum=cos_sum-min(cos_sum);
%
     clear A;
     clear B;
     clear x;
%      
      der=-0.0065;
      mec=3.6059;
      med=3.9948;
 
     A=[ R^(mec-1) R^med; (mec-1)*R^(mec-2) med*R^(med-1) ];
     B=[ 1; der ];
     x=A\B;
     bb=0;
     cc=x(1);
     dd=x(2);
%
cc=2.8808;
dd=-1.8808;
%%
     clear A; 
     A=0.8744;

%%
    phi=(0.1+0.8)*rand();
    lambda=acos(phi);
    ZZ=zeros(length(theta),length(r));
    T=0;
    V=0;
    tpi=2*pi;
    Mr_max=0;
    Vr_max=0;
    dZdth_max=0;
    radius2=radius^2;
    radius3=radius^3;
%
    for(i=1:length(theta))  % theta
%       
        cos_nth=cos(n*theta(i));
        dcos_nth=-n*sin(n*theta(i));
        ddcos_nth=-n^2*cos(n*theta(i));
        sin_nth=sin(n*theta(i));
%
        for(j=1:length(r))  % radius
%
             rr=r(j);
            rrp2=rr^2;
                 
            rrn2=rrn(j)^2;
            rrn3=rrn(j)^3;
            rrn4=rrn(j)^4;
%
              modu=(                         cc*rrn(j)^(mec-1) +                     dd*rrn(j)^med);
             dmodu=(                 (mec-1)*cc*rrn(j)^(mec-2) +                 med*dd*rrn(j)^(med-1))/radius;
            ddmodu=(         (mec-1)*(mec-2)*cc*rrn(j)^(mec-3) +         med*(med-1)*dd*rrn(j)^(med-2))/radius2;  
           dddmodu=( (mec-1)*(mec-2)*(mec-3)*cc*rrn(j)^(mec-4) + med*(med-1)*(med-2)*dd*rrn(j)^(med-3))/radius3;  
%        
              wave=0.5+0.5*cos_nth;
              term=   modu*wave;
             dterm=   dmodu*wave;
            ddterm=  ddmodu*wave;      
             dddterm=  dddmodu*wave;  
%
            Z=A*term + cos_sum(j); 
%           
            if(i==1)
                ZA(j)=A*term + cos_sum(j);
                ZL(j)=Z;
            end
            if(j==length(r))
                ZR(i)=term;
            end  
 %       
                    ZZ1(j,i)=Z;
                 dZdr=  A*dterm  +   dcos_sum(j);
                dZ2dr2=A*ddterm  +  ddcos_sum(j);    
                dZ3dr3=A*dddterm + dddcos_sum(j);
                
                dZdth=A*modu*(0.5*dcos_nth);
               dZ2dth2=A*modu*(0.5*ddcos_nth);
               
               dZ3drth2=A*dmodu*(0.5*ddcos_nth);
               
               dZ2drdth=A*dmodu*(0.5*dcos_nth);
                
% 
               if( j==(length(r)) )
                   
                   if(i==1)
                       keys=dZ2dr2;
                   end

                    Mr=dZ2dr2 +mu*( (1/r(j))*dZdr   +(1/r(j)^2)*dZ2dth2    );

                    R=r(j);
                    
                    Vr1=dZ3dr3 -(1/R^2)*dZdr +(1/R)*dZ2dr2 -(2/R^3)*dZ2dth2  +(1/R^2)*dZ3drth2;
                    Vr2=(2/R^3)*dZdth -(1/R^2)*dZ2drdth    -(1/R^3)*dZ2dth2  +(1/R^2)*dZ3drth2;
                    
                    Vr=Vr1+(1-mu)*Vr2;
                    
                    if(abs(Mr)>Mr_max && wave >0.2  )
                       Mr_max=abs(Mr);
                    end
                    if(abs(Vr)>Vr_max && wave >0.2  )
                       Vr_max=abs(Vr);
                    end
             
               end
               
               if( dZdr > slopemax)
                   slopemax= dZdr;
               end  
%
      
               if( abs(dZdth) > dZdth_max && rrn(j) < 0.5)
                   dZdth_max= abs(dZdth);
               end  
%
          AA1=dZ2dr2;
          AA2=dZdr/rr;
          AA3=dZ2dth2/rrp2;
          AAS= (AA1+AA2+AA3)^2;
%
          BB1=dZdr/rr;
          BB2=dZ2dth2/rrp2;
          BB=-2*(1-mu)*dZ2dr2*(BB1 + BB2);          
% 
          term=A*(((mec-2)*bb*rrn(j)^(mec-3)   +(med-1)*cc*rrn(j)^(med-2) )/radius2)*(0.5*dcos_nth);
%
          CC=2*(1-mu)*( term^2 );
%
          V=V+(AAS+BB+CC)*rr;
          T=T+(Z^2)*rr;
%
          if(T<1.0e+40)
          else
              disp(' T error ');
              gg=input(' ');
          end
          if(V<1.0e+40)
          else
              disp(' V error ');
              gg=input(' ');
          end          
%
        end
    end
%
    out1=sprintf('\n V=%8.3g  T=%8.3g  V/T=%8.3g  sqrt(V/T)=%8.3g  \n',V,T,V/T,sqrt(V/T));
    disp(out1);
    V=V*(De/2)*(dr*dtheta);
    T=T*rho*h*(dr*dtheta);
    omegan=sqrt(V/T);
%    
    if(omegan<1.0e+40)
    else
       disp(' omegan error ');
       gg=input(' ');
    end   
%
    ratio= max(ZZ1(:,1)) / max(ZZ1(50,:));
    fn=omegan/(2*pi);
    ZZr=ZZ1;
    ZAr=ZA;
    ZRr=ZR;
%    
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
out1=sprintf(' fn = %7.4g Hz  ',fn);
disp(out1)
%
figure(1);
polarplot3d(ZZr,'PlotType','surfn','meshscale',0.5,'PolarGrid',[20 72]);
%
title(out1);
figure(2);
plot(r,ZAr);
grid on;
figure(3);
plot(ZRr);
grid on;