function fn = Timoshenko_fixed_free_beam(fn, k, E, G, area, I, L, rho)
% Fast Bruch–Mitchell solver for a Timoshenko cantilever (fixed–free).
% Inputs:
%   fn    Southwell–Dunkerley estimate for the 1st mode (Hz)
%   k     shear-correction factor (kappa)
%   E,G,area,I,L,rho  beam properties (SI)
%
% Output:
%   fn    refined Bruch–Mitchell frequency (Hz)

% ---- precompute constants
A   = area;
r2  = I/(A*L^2);
s2  = E*I/(k*A*G*L^2);
b2c = rho*A*L^4/(E*I);

% no tip mass/eccentricity in this case
Mbar = 0; sigma2 = 0;
MT_initial   = 0; %#ok<NASGU>
term_initial = 0; %#ok<NASGU>

% ---- coarse scan to find a sign-change bracket
fstart = max(1e-6, 0.60*fn);
df     = 0.5;                 % coarse step in Hz
fmax   = max(1.5*fn, fstart + 2000);   % generous initial ceiling

fgrid = fstart:df:fmax;
detg  = arrayfun(@det_real, fgrid);

ix = find(isfinite(detg(1:end-1)) & isfinite(detg(2:end)) & ...
          detg(1:end-1).*detg(2:end) <= 0, 1, 'first');

% expand if needed
expand = 2;
while isempty(ix) && expand <= 20
    fgrid2 = (fmax+df):df:(expand*fn);
    if isempty(fgrid2), break; end
    detg2  = arrayfun(@det_real, fgrid2);
    fgrid  = [fgrid, fgrid2]; %#ok<AGROW>
    detg   = [detg , detg2 ]; %#ok<AGROW>
    ix = find(isfinite(detg(1:end-1)) & isfinite(detg(2:end)) & ...
              detg(1:end-1).*detg(2:end) <= 0, 1, 'first');
    fmax = expand*fn;
    expand = expand + 1;
end

if isempty(ix)
    % no valid bracket found: return the best coarse minimum as fallback
    [~,j] = min(abs(detg));
    fn = fgrid(j);
    return
end

fL = fgrid(ix);
fR = fgrid(ix+1);

% ---- refine with fzero on the bracket
opts = optimset('TolX',1e-9);
fn   = fzero(@det_real, [fL fR], opts);

% ---------- nested determinant (real part) ----------
    function val = det_real(fHz)
        omega = 2*pi*fHz;
        b2    = b2c*(omega^2);
        b     = sqrt(b2);

        sqterm = sqrt((r2 - s2)^2 + (4/b^2));
        alpha2 = 0.5*( -(r2 + s2) + sqterm );
        beta2  = 0.5*(  (r2 + s2) + sqterm );

        % guard invalid/complex regions to keep scan finite & real
        if ~(alpha2 > 1e-20) || ~(beta2 > 0)
            val = 1e6;
            return
        end

        alpha = sqrt(alpha2);
        beta  = sqrt(beta2);

        ba = b*alpha;  bb = b*beta;

        coshba = cosh(ba);  sinhba = sinh(ba);
        cosbb  = cos(bb);   sinbb  = sin(bb);

        MC  = Mbar*b2;
        MCC = 0.5*Mbar*sigma2*b2;

        alphaC = (alpha2 + s2)/alpha;
        betaC  = (beta2  - s2)/beta;

        R1  = (b/alpha)*sinhba + MC*coshba;
        R2  = (b/alpha)*coshba + MC*sinhba;
        R3  = (b/beta)*sinbb   + MC*cosbb;
        R4  = -(b/beta)*cosbb  + MC*sinbb;

        R1P =  alphaC*( ba*coshba - MCC*sinhba);
        R2P =  alphaC*( ba*sinhba - MCC*coshba);
        R3P = -betaC *( bb*cosbb  - MCC*sinbb );
        R4P =  betaC *(-bb*sinbb  - MCC*cosbb );

        Mmat = [ 1      0      1      0
                 0   alphaC   0    betaC
                 R1     R2     R3     R4
                 R1P    R2P    R3P    R4P ];

        num = det(Mmat);
        val = real(num);
        if ~isfinite(val), val = 1e6; end
    end
end




