% vibrationdata_circular_center_clamped_rim_free.m  ver 2.2 (2025‑08‑03)
% Author: Tom Irvine
% -------------------------------------------------------------------------
% Centre‑clamped / rim‑free solid circular plate (NASA SP‑160 roots).
% Computes natural frequencies, participation factors, and plots selected
% low‑order modes, using **polarplot3d** if available (no seam wrap, colour
% matches both sides); otherwise falls back to SURF.
% -------------------------------------------------------------------------
% Syntax
%   [fn, part, nextFig] = vibrationdata_circular_center_clamped_rim_free( ...
%                      Omega, a, h, rho, E, nu, firstFig, zFactor )
%
% Inputs
%   Omega    – spin rate,  rad/s   (0 → stationary disk)
%   a        – radius,      m
%   h        – thickness,   m
%   rho      – density,     kg/m^3
%   E        – Young’s modulus, Pa
%   nu       – Poisson’s ratio
%   firstFig – first figure number to use
%   zFactor  – vertical exaggeration multiplier (default 8×)
%
% Outputs
%   fn       – column vector of natural frequencies (Hz)
%   part     – participation Γ for each mode
%   nextFig  – the next available figure number
% -------------------------------------------------------------------------
function [fn, part, fig_num] = ...
          vibrationdata_circular_center_clamped_rim_free(Omega, a, h, ...
                                                         rho,  E,  nu, ...
                                                         fig_num, zFactor)
if nargin < 9 || isempty(zFactor);  zFactor = 8;  end  % default exaggeration

%% plate bending stiffness -------------------------------------------------
D = E * h^3 / (12 * (1 - nu^2));

%% NASA SP‑160 eigen‑value tables (m = 0:3, n = 0:3) ----------------------
K1 = [ 0     3.30   9.90  19.80 ;
       1.00  5.95  14.20  25.75 ;
       2.35  8.95  18.85  32.05 ;
       4.05 12.30  23.85  38.70 ];
K2 = [  14.10   437.30  3683    14330 ;
          0      421.20  3336    14380 ;
         28.97  1212.00  7164    23410 ;
        155.30  2839.00 11700    36274 ];

%% integration grid -------------------------------------------------------
Nr  = 400;  Nth = 240;
r   = linspace(0, a, Nr);
th  = linspace(0, 2*pi, Nth);
[R, TH] = meshgrid(r, th);
dr  = r(2)  - r(1);
dth = th(2) - th(1);
[Xcart, Ycart] = pol2cart(TH, R);   % for fallback SURF

%% header -----------------------------------------------------------------
fprintf('\n   m   n   f_NASA(Hz)     Γ        M_eff (%% tot)\n');
fprintf('  --------------------------------------------------------\n');

RM = zeros(16,2); kk = 1;

%% ========================= MAIN LOOPS ===================================
for m = 0:3
    for n = 0:3
        K1mn = K1(m+1,n+1);  K2mn = K2(m+1,n+1);
        if isnan(K1mn);  continue; end

        %% natural frequency -------------------------------------------
        fNASA = sqrt(K1mn*Omega^2 + K2mn*D/(rho*h*a^4)) / (2*pi);

        %% raw mode shape ---------------------------------------------
        if (m==1 && n==0 && K2mn==0)               % rigid coning
            Wraw = (R/a) .* cos(TH);
        elseif K2mn==0                             % membrane‑like
            Wraw = (R/a).^n .* cos(m*TH);
        else                                       % bending family
            beta = K2mn^(1/4);
            kr   = beta * R / a;
            if m==0                                % axisymmetric
                J0=besselj(0,kr);  I0=besseli(0,kr);
                J1b=besselj(1,beta); I1b=besseli(1,beta);
                Wraw = J0 - (J1b/I1b)*I0;
            else                                   % m ≥ 1
                Jm=besselj(m,kr);  Im=besseli(m,kr);
                Jmb=besselj(m,beta); Imb=besseli(m,beta);
                Wraw = (Jm - (Jmb/Imb)*Im) .* cos(m*TH);
            end
        end

        %% mass‑normalise & participation ------------------------------
        int_phi2 = rho*h * sum(sum(Wraw.^2 .* R)) * dr * dth;
        Wn       = Wraw / sqrt(int_phi2);
        G        = rho*h * sum(sum(Wn .* R)) * dr * dth;
        fprintf('  %2d  %2d   %9.2f   %7.4f   %8.2f\n', m, n, fNASA, G, (G^2)*100);

        RM(kk,:) = [fNASA, G]; kk = kk+1;

        %% plot subset (m,n ≤ 1) --------------------------------------
        if m<=1 && n<=1
            if abs(Wraw(1,1)) > 1e-12*max(abs(Wraw(:)))
                Wraw = Wraw - Wraw(1,1);          % pin hub exactly at 0
            end

            %% colour data (±target) ----------------------------------
            target = 0.05*a; if m==0&&n==0, target=0.15*a; end
            if m==0&&n>0, target=0.08*a; end
            Wc = target * Wraw / max(abs(Wraw(:)));    % colour matrix

            %% exaggerated Z‑data -------------------------------------
            Wz = zFactor * Wraw;                       % height matrix

            %% try polarplot3d first ---------------------------------
            if exist('polarplot3d','file') == 2
                Zp = Wz.';                 % polarplot3d expects θ rows
                figure(fig_num); fig_num = fig_num + 1;
                polarplot3d(Zp, ...
                             'PlotType','surfn', ...
                             'RadialRange',[0 a], ...
                             'MeshScale',0.5);                  % uniform mesh
                colormap(turbo(256));
                cMax = max(abs(Zp(:))); caxis([-cMax cMax]);
                hold on;
                contour3(Xcart, Ycart, Wz, [0 0], 'k', 'LineWidth', 4);
                hold off;
            else
                %% fallback SURF ------------------------------------
                figure(fig_num); fig_num = fig_num + 1;
                surf(Xcart, Ycart, Wz, Wc, 'EdgeColor','none', 'FaceColor','interp');
                caxis([-target target]); colormap(turbo(256));
                hold on;
                contour3(Xcart, Ycart, Wz, [0 0], 'k', 'LineWidth',4);
                hold off;
            end

            axis equal;  box on;
            xlabel('X'); ylabel('Y'); zlabel('Z');
            view(45,30);

            so=strtrim(sprintf('%8.3g Hz',Omega/(2*pi)));
            title(sprintf('f = %.2f Hz   m=%d, n=%d   Spin Rate=%s', ...
                          fNASA, m, n, so));
        end
    end
end  % loops

%% trim & sort -----------------------------------------------------------
RM(kk:end,:) = [];
RM           = sortrows(RM,1);
fn   = RM(:,1);
part = RM(:,2);

end  % =====================  function  ================================
